A multiscale finite element method for the Schr\"{o}dinger equation with multiscale potentials

TL;DR

This paper introduces a multiscale finite element method for efficiently solving the Schrödinger equation with complex multiscale potentials in the semiclassical regime, overcoming limitations of asymptotic approaches.

Contribution

The paper develops a novel multiscale finite element approach that constructs basis functions via sparse Hamiltonian compression, enabling accurate solutions without specific potential form assumptions.

Findings

Method achieves convergence rates of first-order in $H^1$ norm and second-order in $L^2$ norm.

Numerical tests demonstrate robustness across various multiscale potentials in 1D and 2D.

Spatial and temporal discretizations are independent of the semiclassical parameter $\

Abstract

In recent years, an increasing attention has been paid to quantum heterostructures with tailored functionalities, such as heterojunctions and quantum matematerials, in which quantum dynamics of electrons can be described by the Schr\"{o}dinger equation with multiscale potentials. The model, however, cannot be solved by asymptoics-based approaches where an additive form of different scales in the potential term is required to construct the prescribed approximate solutions. In this paper, we propose a multiscale finite element method to solve this problem in the semiclassical regime. The localized multiscale basis are constructed using sparse compression of the Hamiltonian operator, and thus are "blind" to the specific form of the potential term. After an one-shot eigendecomposition, we solve the resulting system of ordinary differential equations explicitly for the time evolution. In our approach, the spatial mesh size is $ H=\mathcal{O}(\epsilon) $ where $\epsilon$ is the semiclassical parameter and the time stepsize $ k$ is independent of $\epsilon$. Numerical examples in one dimension with a periodic potential, a multiplicative two-scale potential, and a layered potential, and in two dimension with an additive two-scale potential and a checkboard potential are tested to demonstrate the robustness and efficiency of the proposed method. Moreover, first-order and second-order rates of convergence are observed in $H^1$ and $L^2$ norms, respectively.

Equations

$$\{\begin{array}{rl}i\epsilon {\partial }_t{\psi }^{\epsilon }& =-\frac{{\epsilon }^2}{2}\Delta {\psi }^{\epsilon }+v^{\epsilon }(\text{x}){\psi }^{\epsilon },\text{x}\in {\mathbb{R}}^d,t\in \mathbb{R},\\ {\psi }^{\epsilon }{|}_{t=0}& ={\psi }_{\text{in}}(\text{x}),\text{x}\in {\mathbb{R}}^d,\end{array}$$

$$\{\begin{array}{rl}i\epsilon {\partial }_t{\psi }^{\epsilon }& =-\frac{{\epsilon }^2}{2}\Delta {\psi }^{\epsilon }+v^{\epsilon }(\text{x}){\psi }^{\epsilon },\text{x}\in D,t\in \mathbb{R},\\ {\psi }^{\epsilon }& \in H_{\text{P}}^1(D),\\ {\psi }^{\epsilon }{|}_{t=0}& ={\psi }_{\text{in}}(\text{x}).\end{array}$$

$$||{\psi }^{\epsilon }|{|}_V=\frac{1}{2}(\mathcal{H}{\psi }^{\epsilon },{\psi }^{\epsilon })=\frac{1}{2}{\int }_D\frac{{\epsilon }^2}{2}|\nabla {\psi }^{\epsilon }{|}^2+v^{\epsilon }(\text{x})|{\psi }^{\epsilon }{|}^2\mathrm{d}\text{x}.$$

$$V^H=\{{\phi }_i(\text{x}):i=1,...,N_x\}\subset H_{\text{P}}^1(D),$$

$${\varphi }_i$$

$${\int }_D\varphi {\phi }_j\mathrm{d}\text{x}={\delta }_{i,j},\forall 1\le j\le N_x.$$

$$D_0$$

$$D_{\ell }$$

$$\sqrt{V_0}H/\epsilon \lesssim 1,$$

$$||\nabla {\varphi }_i(\text{x})|{|}_{D\backslash D_{\ell }}\le C{\beta }^{\ell }||\nabla {\varphi }_i(\text{x})|{|}_D,$$

Full text

Keyword: Schrödinger equation; localized basis function; operator compression; optimization method; multiscale potential.

AMS subject classifications. 65M60, 74Q10, 35J10

1 Introduction

In solid state physics, one of the most popular models to describe electron dynamics is the Schrödinger equation in the semiclassical regime

[TABLE]

where $0<\varepsilon\ll 1$ is an effective Planck constant describing the microscopic and macroscopic scale ratio, $d$ is the spatial dimension, $v^{\varepsilon}(\textbf{x})$ is the given electrostatic potential, $\psi^{\varepsilon}=\psi^{\varepsilon}(t,\textbf{x})$ is the wavefunction, and $\psi_{\textrm{in}}(\textbf{x})$ is the initial data. In the community of mathematics, there has been a long history of interest from both mathematical and numerical perspectives; see for example the review paper [JinActa:2011] and references therein.

In the simplest situation, $v^{\varepsilon}(\textbf{x})=u(\textbf{x})$, where $u(\textbf{x})$ is an (external) macroscopic potential. $\psi^{\varepsilon}(t,\textbf{x})$ propagates oscillations with a wavelength of $\mathcal{O}(\varepsilon)$, a uniform $L^{2}-$approximation of the wavefunction requires the spatial mesh size $h=o(\varepsilon)$ and the time step $k=o(\varepsilon)$ in finite element method (FEM) and finite difference method (FDM) [BaoJinMarkowich:2002, JinActa:2011]. If the spectral time-splitting method is employed, a uniform $L^{2}-$approximation of the wavefunction requires the spatial mesh size $h=\mathcal{O}(\varepsilon)$ and the time stepsize $k=o(\varepsilon)$ [BaoJinMarkowich:2002]. For a perfect crystal, in the presence of an external field, $v^{\varepsilon}(\textbf{x})=v(\frac{\textbf{x}}{\varepsilon})+u(\textbf{x})$, where $v(\frac{\textbf{x}}{\varepsilon})$ describes the electrostatic interaction of ionic cores. A number of methods have been proposed by taking advantage of the periodic structure of $v(\frac{\textbf{x}}{\varepsilon})$, such as the Bloch decomposition based time-splitting spectral method [Huangetal:2007, Huang:2008], the Gaussian beam method [Jin:08, Jin:2010, Qian:2010, Yin:2011], and the frozen Gaussian approximation method [DelgadilloLuYang:2016]. The Bloch decomposition based time-splitting spectral method requires a mesh strategy $h=\mathcal{O}(\varepsilon)$ and $k=\mathcal{O}(1)$ for the uniform $L^{2}-$approximation of the wavefunction. The Gaussian beam method and the frozen Gaussian approximation method are based on asymptotic analysis, and thus are especially efficient when $\varepsilon$ is very small.

With recent developments in nanotechnology, a variety of material devices with tailored functionalities have been fabricated, such as heterojunctions, including the ferromagnet/metal/ferromagnet structure for giant megnetoresistance [Zutic:2004], the silicon-based heterojunction for solar cells [Louwenetal:2016], and quantum metamaterials [Quach:2011]. A basic feature of these devices is the combination of dissimilar crystalline structures, which results a heterogeneous interaction from ionic cores with different lattice structures. Therefore, when travelling through a device, electrons experience a potential $v^{\varepsilon}(\textbf{x})$ which cannot be written in the abovementioned form. Consequently, all the available methods based on asymptotic analysis cannot be applied. Moreover, direct methods, such as FEM and FDM, are extremely inefficient with strong mesh size restrictions. This motivates us to design efficient numerical methods for (1) in the general situation.

Our work is motivated by the multiscale FEM for solving elliptic problems with multiscale coefficients [HouWuCai:99, EfendievHou:09]. The multiscale FEM is capable of correctly capturing the large scale components of the multiscale solution on a coarse grid without accurately resolving all the small scale features in the solution. This is accomplished by incorporating the local microstructures of the differential operator into the finite element basis functions. We remark that in the past four decades, many other efficient methods have been developed for the multiscale PDEs in the literature; see [Babuska:94, Hughes:1998, ChenHou:02, Jenny:03, EngquistE:03, Kevrekidis:2003, OwhadiZhang:07, HanZhangCMS:12] for example and references therein.

Recently, several works relevant to the compression of elliptic operator with heterogeneous and highly varying coefficients have been proposed. In [Peterseim:2014], Malqvist and Peterseim construct localized multiscale basis functions using a modified variational multiscale method. The exponentially decaying property of these modified basis has been shown both theoretically and numerically. Meanwhile, Owhadi [Owhadi:2015, Owhadi:2017] reformulates the multiscale problem from the perspective of decision theory using the idea of gamblets as the modified basis. In particular, a coarse space of measurement functions is constructed from Bayesian perspective, and the gamblet space is explicitly constructed. In addition, the gamblets are still proven to decay exponentially such that localized computation is made possible. Hou and Zhang [hou2017sparse] extend these works such that localized basis functions can also be constructed for higher-order strongly elliptic operators.

In this paper, we propose a multiscale FEM to solve the Schrödinger equation in the semiclassical regime. The localized multiscale basis are constructed using sparse compression of the Hamiltonian operator, and thus are “blind” to the specific form of the potential. After an one-shot eigendecomposition, we can solve the resulting system of ordinary differential equations explicitly for the time evolution. In our approach, $H=\mathcal{O}(\varepsilon)$ and $k$ is independent of $\varepsilon$. Numerical examples in one dimension with a periodic potential, a multiplicative two-scale potential, and a layered potential, and in two dimension with an additive two-scale potential and a checkboard potential are tested to demonstrate the robustness and efficiency of the proposed method. Moreover, first-order and second-order rates of convergence are observed in $H^{1}$ and $L^{2}$ norms, respectively.

The rest of the paper is organized as follows. In §2, we introduce a multiscale FEM for the semiclassical Schrödinger equation and discuss the properties of the proposed method. Numerous numerical results are presented in §LABEL:sec:NumericalExamples, including both one dimensional and two dimensional examples to demonstrate the efficiency of the proposed method. Conclusions and discussions are drawn in §LABEL:sec:Conclusion.

2 Multiscale finite element method for the semiclassical Schrödinger equation

In this section, we construct the multiscale finite element basis functions based on an optimization approach, and use these basis functions as the approximation space in the Galerkin method to solve the Schrödinger equation. A couple of properties of the proposed method are also given.

2.1 Construction of multiscale basis functions

Recall that the Schrödinger equation (1) is defined in $\mathbb{R}^{d}$. However, numerically we can only deal with bounded domains, thus artificial boundary condition is needed here. For the sake of brevity, we shall restrict ourselves to a bounded domain with prescribed boundary condition. In fact, artificial boundary condition can also be combined with the proposed approach which will be investigated in a subsequent work. Therefore we consider the following problem

[TABLE]

Here $D=[0,1]^{d}$ is the spatial domain and $H_{\textrm{P}}^{1}(D)=\{\psi|\psi\in H^{1}(D)\textrm{ and }\psi\textrm{ is periodic over D}\}$. $\psi_{\textrm{in}}(\textbf{x})$ is the initial data over $D$. Define the Hamiltonian operator $\mathcal{H}(\cdot)\equiv-\frac{\varepsilon^{2}}{2}\Delta(\cdot)+v^{\varepsilon}(\textbf{x})(\cdot)$ and introduce the following energy notation $||\cdot||_{V}$ for Hamiltonian operator

[TABLE]

Note that (3) does not define a norm since $v^{\varepsilon}$ usually can be negative, and thus the bilinear form associated to this notation is not coercive, which is quite different from the case of elliptic equations. However, this does not mean that available approaches [HouWu:97, BabuskaLipton:2011, Peterseim:2014, Owhadi:2017, hou2017sparse] cannot be used for the Schrödinger equation. In fact, we shall utilize the similar idea to construct localized multiscale finite element basis functions on a coarse mesh by an optimization approach using the above energy notation $||\cdot||_{V}$ for the Hamiltonian operator.

To construct such localized basis functions, we first partition the physical domain $D$ into a set of regular coarse elements with mesh size $H$. For example, we divide $D$ into a set of non-overlapping triangles $\mathcal{T}_{H}=\cup\{K\}$, such that no vertex of one triangle lies in the interior of the edge of another triangle. In each element $K$, we define a set of nodal basis $\{\varphi_{j,K},j=1,...,k\}$ with $k$ being the number of nodes of the element. From now on, we neglect the subscript $K$ for notational convenience. The functions $\varphi_{i}(\textbf{x})$ are called measurement functions, which are chosen as the characteristic functions on each coarse element in [hou2017sparse, Owhadi:2017] and piecewise linear basis functions in [Peterseim:2014]. In [LiZhangCiCP:18, hou2018model], it is found that the usage of nodal basis functions reduces the approximation error and thus the same setting is adopted in the current work.

Let $\mathcal{N}$ denote the set of vertices of $\mathcal{T}_{H}$ (removing the repeated vertices due to the periodic boundary condition) and $N_{x}$ be the number of vertices. For every vertex $\textbf{x}_{i}\in\mathcal{N}$, let $\varphi_{i}(\textbf{x})$ denote the corresponding nodal basis function, i.e., $\varphi_{i}(\textbf{x}_{j})=\delta_{ij}$. Since all the nodal basis functions $\varphi_{i}(\textbf{x})$ are continuous across the boundaries of the elements, we have

[TABLE]

Then, we can solve optimization problems to obtain the multiscale basis functions. Specifically, let $\phi_{i}(\textbf{x})$ be the minimizer of the following constrained optimization problem

[TABLE]

The superscript $\varepsilon$ is dropped for notational simplicity and the periodic boundary condition is incorporated into the above optimization problem through the solution space $H_{\textrm{P}}^{1}(D)$. With these multiscale finite element basis functions $\{\phi_{i}(\textbf{x})\}_{i=1}^{N_{x}}$, we can solve the Schrödinger equation (2) using the Galerkin method.

Remark 2.1*.*

Note that the energy notation $||\cdot||_{V}$ in (3) does not define a norm. However, as long as $v^{\varepsilon}(\textbf{x})$ is bounded from below and the fine mesh size $h$ is small enough, the discrete problem of (4) - (5) is convex and thus admits a unique solution; see [hou2017sparse, LiZhangCiCP:18] for details.

2.2 Exponential decay of the multiscale finite element basis functions

We shall show that the multiscale basis functions $\{\phi_{i}(\textbf{x})\}_{i=1}^{N_{x}}$ decay exponentially fast away from its associated vertex $x_{i}\in\mathcal{N}_{c}$ under certain conditions. This allows us to localize the basis functions to a relatively smaller domain and reduce the computational cost.

In order to obtain localized basis functions, we first define a series of nodal patches $\{D_{\ell}\}$ associated with $\textbf{x}_{i}\in\mathcal{N}$ as

[TABLE]

Assumption 2.1**.**

We assume that the potential term $v^{\varepsilon}(\textbf{x})$ is bounded, i.e., $V_{0}:=||v^{\varepsilon}(\textbf{x})||_{L^{\infty}(D)}<+\infty$ and the mesh size $H$ of $\mathcal{T}_{H}$ satisfies

[TABLE]

where $\lesssim$ means bounded from above by a constant.

Under this assumption, many typical potentials in the Schrödinger equation (2) can be treated as a perturbation to the kinetic operator. Thus, they can be computed using our method. Then, we can show that the multiscale finite element basis functions have the exponentially decaying property.

Proposition 2.2** (Exponentially decaying property).**

Under the resolution condition of the coarse mesh, i.e., (8), there exist constants $C>0$ and $0<\beta<1$ independent of $H$, such that

[TABLE]

for any $i=1,2,...,N_{x}$.

Proof of (9) will be given in [ChenMaZhang:prep]. The main idea is to combine an iterative Caccioppoli-type argument [Peterseim:2014, LiZhangCiCP:18] and some refined estimates with respect to $\varepsilon$. To demonstrate the exponentially decaying property of multiscale basis functions, we use the multiscale basis function centered at $x=1/2$ in Example LABEL:example2 for a sequence of $\varepsilon$ from $1/40$ to $1/160$ for illustration. The left figure plots $\nabla\phi/||\phi||_{L_{2}}$ with respect to the distance to $x=1/2$, which shows both the exponential decay and the $\varepsilon$ dependence with respect to the distance. The right figure plots $E_{\textrm{relative}}=\dfrac{||\nabla\phi||_{D}-||\nabla\phi||_{D_{\ell}}}{\max(||\nabla\phi||_{D}-||\nabla\phi||_{D_{\ell}})}$ with respect to the patch size $\ell$, which shows the decay rate of $E_{\textrm{relative}}$ with respect to $\ell$ is independent of $\varepsilon$, and thus the estimate in (9) is sharp. $||\phi||_{L_{2}}$ and $\max(||\nabla\phi||_{D}-||\nabla\phi||_{D_{\ell}})$ in the denominators are used such that $\nabla\phi/||\phi||_{L_{2}}$ and $E_{\textrm{relative}}$ with respect to $\varepsilon$ are in a similar range of magnitudes.