A multiscale reduced basis method for Schr\"{o}dinger equation with multiscale and random potentials

TL;DR

This paper introduces a multiscale reduced basis method for the semiclassical Schrödinger equation with multiscale and random potentials, effectively addressing high-frequency oscillations and enabling efficient numerical simulations in complex quantum systems.

Contribution

The paper develops a novel multiscale reduced basis approach combining optimization, proper orthogonal decomposition, and quasi-Monte Carlo methods for efficient Schrödinger equation solutions.

Findings

Method reduces spatial grid size proportional to the semiclassical parameter

Number of random samples inversely proportional to the semiclassical parameter

Numerical results confirm efficiency and accuracy of the approach

Abstract

The semiclassical Schr\"{o}dinger equation with multiscale and random potentials often appears when studying electron dynamics in heterogeneous quantum systems. As time evolves, the wavefunction develops high-frequency oscillations in both the physical space and the random space, which poses severe challenges for numerical methods. In this paper, we propose a multiscale reduced basis method, where we construct multiscale reduced basis functions using an optimization method and the proper orthogonal decomposition method in the physical space and employ the quasi-Monte Carlo method in the random space. Our method is verified to be efficient: the spatial gridsize is only proportional to the semiclassical parameter and the number of samples in the random space is inversely proportional to the same parameter. Several theoretical aspects of the proposed method, including how to determine the number of samples in the construction of multiscale reduced basis and convergence analysis, are studied with numerical justification. In addition, we investigate the Anderson localization phenomena for Schr\"{o}dinger equation with correlated random potentials in both 1D and 2D.

Figures

Tables

$H$	${\text{Error}}_{L^2}$	Order	${\text{Error}}_{H^1}$	Order
$2\pi /32$	0.09862312		0.32096262
$2\pi /64$	0.00129644	6.25	0.01449534	4.47
$2\pi /128$	0.00002892	5.49	0.00076150	4.25
$2\pi /256$	0.00000950	1.61	0.00014161	2.42

Equations

$$\{\begin{array}{rl}i\epsilon {\partial }_t{\psi }^{\epsilon }& =-\frac{{\epsilon }^2}{2}\Delta {\psi }^{\epsilon }+v^{\epsilon }(\text{x},\omega ){\psi }^{\epsilon },\text{x}\in D,\omega \in \Omega ,t\in \mathbb{R},\\ {\psi }^{\epsilon }& \in H_{\text{P}}^1(D),\omega \in \Omega ,t\in \mathbb{R},\\ {\psi }^{\epsilon }{|}_{t=0}& ={\psi }_{\text{in}}(\text{x}),\text{x}\in 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^H(\text{x}):i=1,...,N_H\}\subset H_{\text{P}}^1(D).$$

$${\varphi }_i$$

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

$$D_0$$

$$D_{\ell }$$

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

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

$${\varphi }_i^{\text{loc}}$$

$${\int }_{D_{l^{\ast }}}\varphi {\phi }_j^H\mathrm{d}\text{x}={\delta }_{i,j},\forall 1\le j\le N_H,$$

$$\varphi (\text{x})=0,\text{x}\in D\backslash D_{l^{\ast }},$$

$$\displaystyle C(\textbf{x},\textbf{y})=\sigma^{2}\exp\big{(}-\sum_{i=1}^{d}\frac{|x_{i}-y_{i}|^{2}}{2l^{2}_{i}}\big{)},$$

$$P(\omega \in \Omega |v^{\epsilon }(\text{x},\omega )\in [v_{\min },v_{\max }],\forall \text{x}\in D)=1.$$

$$v^{\epsilon }(\text{x},\omega )={\bar{v}}^{\epsilon }(\text{x})+\sum _{i=1}^{\infty }\sqrt{{\lambda }_i}{\xi }_i(\omega )v_i(\text{x}),$$

$$v_m^{\epsilon }(\text{x},\omega )={\bar{v}}^{\epsilon }(\text{x})+\sum _{i=1}^m\sqrt{{\lambda }_i}{\xi }_i(\omega )v_i(\text{x}),$$

$$\displaystyle\frac{\sum_{q=1}^{Q}\Big{|}\Big{|}\tilde{\phi}_{k}^{\textrm{loc}}(\textbf{x},\omega_{q})-\sum_{l=1}^{m_{k}}\big{(}\tilde{\phi}_{k}^{\textrm{loc}}(\textbf{x},\omega_{q}),\zeta^{k}_{l}(\textbf{x})\big{)}_{X}\zeta^{k}_{l}(\textbf{x})\Big{|}\Big{|}_{X}^{2}}{\sum_{q=1}^{Q}\Big{|}\Big{|}\tilde{\phi}_{k}^{\textrm{loc}}(\textbf{x},\omega_{q})\Big{|}\Big{|}_{X}^{2}}=\frac{\sum_{s=m_{k}+1}^{Q}\lambda_{s}}{\sum_{s=1}^{Q}\lambda_{s}},$$

$${\varphi }_k^{\text{loc}}(\text{x},\omega )\approx {\zeta }_0^k(\text{x})+\sum _{l=1}^{m_k}{c}_l(\omega ){\zeta }_l^k(\text{x}).$$

$$\{\begin{array}{rl}& \underset{\mathbf{c}}{\min }\frac{1}{2}{\mathbf{c}}^{T}Q\mathbf{c},\\ & \text{s.t. }A\mathbf{c}=\mathbf{b},\end{array}$$

$$Q_{ij}=\frac{{\epsilon }^2}{2}(\nabla {\phi }_i^h,\nabla {\phi }_j^h)+(v^{\epsilon }(\text{x},{\omega }_q){\phi }_i^h,{\phi }_j^h).$$

$$\Vert \varphi (\cdot ,{\omega }_1)-\varphi (\cdot ,{\omega }_2){\Vert }_{L^{\infty }(D)}\le \frac{C}{{\kappa }^6}{\epsilon }^{-2}\Vert v^{\epsilon }(\cdot ,{\omega }_1)-v^{\epsilon }(\cdot ,{\omega }_2){\Vert }_{L^{\infty }(D)},$$

$$\mathbf{c}=Q^{-1}{A}^T(A{Q}^{-1}{A}^T{)}^{-1}\mathbf{b}.$$

$${\left(\delta V\right)}_{ij}=\left(\left(v(\cdot ,{\omega }_1)-v(\cdot ,{\omega }_2)\right){\phi }_i^h,{\phi }_j^h\right),$$

$$\Vert \delta V{\Vert }_{\infty }\le h^d\Vert v^{\epsilon }(\cdot ,{\omega }_1)-v^{\epsilon }(\cdot ,{\omega }_2){\Vert }_{L^{\infty }(D)}.$$

$$Q_2^{-1}=\sum _{n=0}^{\infty }{\left(Q_1^{-1}\delta V\right)}^n{Q}_1^{-1},$$

$${\mathbf{c}}_2-{\mathbf{c}}_1$$

$$=Q_1^{-1}\delta V{Q}_1^{-1}{A}^T(A{Q}_1^{-1}{A}^T{)}^{-1}\mathbf{b}$$

$$-Q_2^{-1}{A}^T(A{Q}_1^{-1}{A}^T{)}^{-1}(A{Q}_1^{-1}\delta V{Q}_1^{-1}{A}^T)(A{Q}_1^{-1}{A}^T{)}^{-1}\mathbf{b}+o(\Vert \delta V{\Vert }_{\infty }),$$

$$=Q_1^{-1}\delta V{Q}_1^{-1}{A}^T(A{Q}_1^{-1}{A}^T{)}^{-1}\mathbf{b}$$

$$-Q_1^{-1}{A}^T(A{Q}_1^{-1}{A}^T{)}^{-1}(A{Q}_1^{-1}\delta V{Q}_1^{-1}{A}^T)(A{Q}_1^{-1}{A}^T{)}^{-1}\mathbf{b}+o(\Vert \delta V{\Vert }_{\infty }).$$

$$|{\mathbf{c}}_2-{\mathbf{c}}_1{|}_{\infty }$$

$$\Vert A{\Vert }_{\infty }\le C{h}^d,|\mathbf{b}{|}_{\infty }=1,\Vert Q_1^{-1}{\Vert }_{\infty }\le C{h}^{-2},\Vert Q_1{\Vert }_{\infty }\le C\max \{{\epsilon }^2,h^2\}\le C{\epsilon }^2,$$

$$|{\mathbf{c}}_2-{\mathbf{c}}_1{|}_{\infty }$$

Full text

Keyword: random Schrödinger equation; multiscale reduced basis function; optimization method; quasi-Monte Carlo method; Anderson localization.

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

1 Introduction

The semiclassical Schrödinger equation describes electron dynamics in the semiclassical regime. Applications of such an equation can be found in Bose-Einstein condensation, graphene, semiconductors, topological insulators, etc. When propagating in a (quasi-)periodic microstructure, electrons experience a multiscale potential. As a consequence, the electron wavefunction develops high-frequency oscillations, which poses severe challenges from the numerical perspective. Brute-force methods are very costly and asymptotics-based methods have been proposed in the literature; see [jin2011mathematical] for review and references therein.

In [Anderson:58], Anderson proposed to study localized eigenstates in a tight-binding model with random potentials. This model was soon to be generalized to the random Schrödinger equation, i.e., the Schrödinger equation with a random potential. In this case, electrons are found to be localized provided that the strength of randomness is sufficiently large. The randomness can be realized in an experiment by enhancing the disorder of impurities in a material. Due to the importance of this model, Anderson was awarded the Nobel Prize in physics in 1977. In the presence of multiscale and random potentials, the electron wavefunction develops high-frequency oscillations in both the physical space and the random space, making numerical approximations even more difficult.

In this paper, we study the following Schrödinger equation with random potential 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},\omega)$ is the given random potential, $\psi^{\varepsilon}=\psi^{\varepsilon}(t,\textbf{x},\omega)$ is the electron wavefunction, and $\psi_{\textrm{in}}(\textbf{x})$ is the initial data. 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}\}$.

Equation (1) can be used to model electron transport in a disordered medium in a single-electron picture where the electron interaction is ignored. It is customary to write the semiclassical Schrödinger equation and the multiscale and random potential with a single parameter $\varepsilon$. But there is no reason that the parameter of the multiscale and random potential should be the same as the semiclassical parameter; see §5 for details on the parameterization of the multiscale and random potential $v^{\varepsilon}(\textbf{x},\omega)$.

The existence of Anderson localization is closely related to the electron wavefunction in (1). To be specific, assume $\psi^{\varepsilon}(t,\textbf{x},\omega)$ has zero mean with respect to the measure $\rho$ induced by $v^{\varepsilon}(\textbf{x},\omega)$ and denote $A(t)=\mathbb{E}[\int_{R^{d}}|\textbf{x}|^{2}|\psi^{\varepsilon}|^{2}d\textbf{x}]_{\rho}$ the second-order moment of the position density. When the strength of disorder is small, an electron undergoes a diffusion process with $A(t)=2C_{d}t$, $C_{d}>0$. In the presence of a strong disorder, however, $A(t)$ converges to a time-independent quantity, i.e., $\lim_{t\rightarrow\infty}A(t)=C$, which implies the localization of the electron and the system undergoes a metal-insulator transition [mott1990metal, Erdos:10]. When $d=1$, localization always occurs for (1) with random potential [Anderson:58]. When $d\geq 2$, the situation becomes complicated. Some analytical results show that localization occurs when the strength of disorder is large [frohlich1983absence, aizenman1993localization]. This motivates us to study Anderson localization in the presence of correlated random potentials [nosov2019correlation].

When the potential is deterministic, i.e., $v^{\varepsilon}(\textbf{x},\omega)=v^{\varepsilon}(\textbf{x})$, many numerical methods have been proposed; see [bao2002time, faou2006poisson, tanushev2007mountain, jin2008gaussian, faou2009computing, delgadillo2016gauge, chen2019multiscale] for example. When the potential is random, few works have been done; see [wu2016bloch, jin2019gaussian]. As mentioned above, the major difficulty is that the wavefunction $\psi^{\varepsilon}$ develops high-frequency oscillations in both the physical space and the random space, which requires tremendous computational resources.

Our work is motivated by the multiscale finite element method (FEM) for solving elliptic problems with multiscale coefficients [hou1999convergence, efendiev2009multiscale]. 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 multiscale FEM basis functions. Recently, several relevant works on constructing localized basis functions that approximate the elliptic operator with heterogeneous coefficients have been proposed. In [maalqvist2014localization], 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, owhadi2017multigrid] reformulates the multiscale problem from the perspective of decision theory using the idea of gamblets as the modified basis. Hou et.al. [hou2017sparse] extend these works such that localized basis functions can also be constructed for higher-order strongly elliptic operators. Recently, Hou, Ma, and Zhang propose to build localized multiscale stochastic basis to solve elliptic problems with multiscale and random coefficients [hou2018model].

In this paper, we propose a multiscale reduced basis method to solve the Schrödinger equation with random potentials in the semiclassical regime. Our method consists of offline and online stages. In the offline stage, we apply an optimization approach to systematically construct localized multiscale reduced basis functions on each patch associated with each coarse gridpoint. These basis functions provide nearly optimal approximation to the random Schrödinger operator. In the online stage, we use these basis functions to approximate the physical space of the solution and the quasi-Monte Carlo (qMC) method to approximate the random space of the solution, respectively. We find the proposed method is efficient in the sense that the number of basis functions is only proportional to $\varepsilon$ and the number of samples in qMC is inversely proportional to $\varepsilon$. Under some conditions, we conduct the convergence analysis of the proposed method with numerical verifications. Moreover, we study how to determine the number of samples in qMC such that the corresponding multiscale reduced basis functions provide accurate approximation of the solution space. Finally we investigate the existence of Anderson localization for correlated random potentials.

The rest of the paper is organized as follows. For completeness, in §2, we introduce multiscale basis functions for the deterministic Schrödinger equation in semiclassical regime and discuss some properties of the basis functions. In §3, we propose a multiscale reduced basis method to solve the random Schrödinger equation. Analysis results are presented in §4 and numerical experiments, including both 1D and 2D examples, are conducted to demonstrate the convergence and efficiency of the proposed method in §5. Conclusions and discussions are drawn in §LABEL:sec:Conclusion.

2 Multiscale basis functions for deterministic Schrödinger equations

In this section, we briefly review the construction of multiscale basis functions based on an optimization approach to solve the Schrödinger equation with a deterministic potential. Some properties of the multiscale basis functions are also given.

2.1 Construction of multiscale basis functions

In the deterministic case, we consider the following problem

[TABLE]

$\psi_{\textrm{in}}(\textbf{x})$ is the initial data over $D$. Defining the Hamiltonian operator $\mathcal{H}(\cdot)\equiv-\frac{\varepsilon^{2}}{2}\Delta(\cdot)+v^{\varepsilon}(\textbf{x})(\cdot)$ and introducing 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 [hou1997multiscale, babuska2011optimal, maalqvist2014localization, owhadi2017multigrid, hou2017sparse] cannot be used for the Schrödinger equation. In fact, we shall utilize the similar idea to construct localized multiscale 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 multiscale 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. On 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, owhadi2017multigrid] and piecewise linear basis functions in [maalqvist2014localization]. In [li2017computing, hou2018model], it is found that the usage of FEM 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_{H}$ be the number of vertices. For every vertex $\textbf{x}_{i}\in\mathcal{N}$, let $\varphi^{H}_{i}(\textbf{x})$ denote the corresponding nodal basis function, i.e., $\varphi^{H}_{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 notation simplicity and the periodic boundary condition is incorporated into the above optimization problem through the solution space $H_{\textrm{P}}^{1}(D)$.

In general, one cannot solve the above optimization problem analytically. Therefore, we use numerical methods to solve it. Specifically, we partition the physical domain $D$ into a set of non-overlapping fine triangles with size $h\ll\varepsilon$. Then, we use standard FEM to discretize $\phi_{i}(\textbf{x})$, $\varphi^{H}_{j}(\textbf{x})$, $1\leq i,j\leq N_{H}$. In the discrete level, the optimization problem (4)-(5) is reduced to a constrained quadratic optimization problem; see (19) in Section 3.3, which can be efficiently solved using Lagrange multiplier methods. Finally, with these multiscale FEM basis functions $\{\phi_{i}(\textbf{x})\}_{i=1}^{N_{H}}$, we can solve the Schrödinger equation (2) using the Galerkin method.

Remark 2.1*.*

In analogy to the multistate FEM [hou1999convergence, efendiev2009multiscale], the multiscale basis functions $\{\phi_{i}(\textbf{x})\}_{i=1}^{N_{H}}$ are defined on coarse elements with mesh size $H$. However, they are represented by fine-scale FEM basis with mesh size $h$, which can be pre-computed and done in parallel.

Remark 2.2*.*

The notation $||\cdot||_{V}$ in (3) does not define a norm. However, as long as the potential $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, li2017computing] for details.

2.2 Exponential decay of the multiscale finite element basis functions

It can be proved that the multiscale basis functions $\{\phi_{i}(\textbf{x})\}_{i=1}^{N_{H}}$ decay exponentially fast away from its associated vertex $\textbf{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. 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 the potential $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 resolution assumption for the coarse mesh, 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_{H}$.

Proof of (9) will be given in [ChenMaZhang:prep]. The main idea is to combine an iterative Caccioppoli-type argument [maalqvist2014localization, li2017computing] and some refined estimates with respect to $\varepsilon$.

The exponential decay of the basis functions enables us to localize the support sets of the basis functions $\{\phi_{i}(\textbf{x})\}_{i=1}^{N_{H}}$, so that the corresponding stiffness matrix is sparse and the computational cost is reduced. In practice, we define a modified constrained optimization problem as follows

[TABLE]

where $D_{l^{*}}$ is the support set of the localized multiscale basis function $\phi_{i}^{\textrm{loc}}(\textbf{x})$ and the choice of the integer $l^{*}$ depends on the decaying speed of $\phi_{i}^{\textrm{loc}}(\textbf{x})$. In (11) and (12), we have used the fact that $\phi_{i}(\textbf{x})$ has the exponentially decaying property so that we can localize the support set of $\phi_{i}(\textbf{x})$ to a smaller domain $D_{l^{*}}$. In numerical experiments, we find that a small integer $l^{*}\sim\log(L/H)$ will give accurate results, where $L$ is the diameter of domain $D$. Moreover, the optimization problem (10)-(12) can be solved in parallel. Therefore, the exponentially decaying property significantly reduces our computational cost in constructing basis functions and computing the solution of the Schrödinger equation (2).

With the localized multiscale finite element basis functions $\{\phi_{i}^{\textrm{loc}}(\textbf{x})\}_{i=1}^{N_{H}}$, we can approximate the wavefunction by $\psi^{\varepsilon}(\textbf{x},t)=\sum_{i=1}^{N_{H}}c_{i}(t)\phi_{i}^{\textrm{loc}}(\textbf{x})$ using the Galerkin method.

3 Multiscale reduced basis functions for the random Schrödinger equation

3.1 Parametrization of the random potential

The random potential $v^{\varepsilon}(\textbf{x},\omega)$ is used to model the disorder in a given material. Specifically, we assume $v^{\varepsilon}(\textbf{x},\omega)$ is a second order random field, i.e., $v^{\varepsilon}(\textbf{x},\omega)\in L^{2}(D,\Omega)$, with mean $\mathbb{E}\left[v^{\varepsilon}(\textbf{x},\omega)\right]=\bar{v}^{\varepsilon}(\textbf{x})$ and covariance kernel $C(\textbf{x},\textbf{y})$. For example, we can choose the covariance kernel as

[TABLE]

where $\sigma$ is a constant and $l_{i}$’s are the correlation lengths in each dimension. We also assume that the random potential $v^{\varepsilon}(\textbf{x},\omega)$ is almost surely bounded, namely there exist $v_{\mathrm{max}}$ and $v_{\mathrm{min}}$, such that

[TABLE]

Circulant embedding method [circulantEmbed1997] and Karhunen-Loève (KL) expansion method [Karhunen:47, Loeve:78] are commonly used to generate samples of $v^{\varepsilon}(\textbf{x},\omega)$, and the latter will be used in the current work. The KL expansion of $v^{\varepsilon}(\textbf{x},\omega)$ reads as

[TABLE]

where $\xi_{i}(\omega)$’s are mean-zero and uncorrelated random variables, i.e., $\mathbb{E}\left[\xi_{i}\right]=0$, $\mathbb{E}\left[\xi_{i}\xi_{j}\right]=\delta_{ij}$, and $\{\lambda_{i},v_{i}(\textbf{x})\}_{i=1}^{\infty}$ are the eigenpairs of the covariance kernel $C(\textbf{x},\textbf{y})$. Generally, $\lambda_{i}$’s are sorted in a descending order and their decay rates depend on the regularity of the covariance kernel. It has been proven that an algebraic decay rate, i.e. $\lambda_{i}=\mathcal{O}(i^{-\gamma})$, is achieved asymptotically if the covariance kernel is of finite Sobolev regularity, and an exponential decay rate is achieved, i.e., $\lambda_{i}=\mathcal{O}(e^{-\gamma i})$ for some $\gamma>0$, if the covariance kernel is piecewisely analytic [schwab2006karhunen].

In practice, we truncate the KL expansion (15) into its first $m$ terms and obtain a parametrization of the random potential as

[TABLE]

which will be used in both analysis and numerics in the remaining part of the paper.

Remark 3.1*.*

In general, the decay rate of $\lambda_{i}$ depends on the correlation lengths $l_{i}$, $i=1,...,d$ of the random field $v^{\varepsilon}(\textbf{x},\omega)$. Small correlation length results in slow decay of the eigenvalues. When the correlation lengths approach zero, the random field $v^{\varepsilon}(\textbf{x},\omega)$ becomes a spatially white noise, which is the case used in the original physics paper [Anderson:58].

3.2 Construction of the multiscale reduced basis functions

For the random Schrödinger equation (1), it is prohibitively expensive to construct multiscale basis functions for each realization of the random potential using (10) - (12). To address this issue, we use a model reduction method to build a small number of reduced basis functions that enable us to obtain multiscale basis functions in a cheaper way without loss of approximation accuracy.

For every $\textbf{x}_{k}\in\mathcal{N}$, we first compute a set of samples of multiscale basis functions associated to the vertex $\textbf{x}_{k}$. Specifically, let $\{v^{\varepsilon}(\textbf{x},\omega_{q})\}_{q=1}^{Q}$ be samples of the random potential that are obtained using Monte Carlo (MC) method or qMC method, where $Q$ is the number of samples. Denote $\zeta_{0}^{k}(\textbf{x})=\frac{1}{Q}\sum_{q=1}^{Q}\phi_{k}^{\textrm{loc}}(\textbf{x},\omega_{q})$ the sample mean of the basis functions, and $\tilde{\phi}_{k}^{\textrm{loc}}(\textbf{x},\omega_{q})=\phi_{k}^{\textrm{loc}}(\textbf{x},\omega_{q})-\zeta_{0}^{k}(\textbf{x})$ is the fluctuation of the $k-$th basis function.

We apply the proper orthogonal decomposition (POD) method [berkooz1993proper, sirovich1987turbulence] to $V=\{\tilde{\phi}_{k}^{\textrm{loc}}(\textbf{x},\omega_{q})\}_{q=1}^{Q}$ and build a set of basis functions $\{\zeta^{k}_{1}(\textbf{x}),\zeta^{k}_{2}(\textbf{x}),...,\zeta^{k}_{m_{k}}(\textbf{x})\}$ with $m_{k}\ll Q$ that optimally approximates $V$. Quantitatively, we have the following approximating property.

Proposition 3.1**.**

*Let $\lambda_{1}\geq\lambda_{2}\geq...\geq\lambda_{m_{k}}\geq\lambda_{m_{k+1}}\geq...>0$ be positive eigenvalues of the covariance kernel associated with the snapshot of the fluctuations $V$ and the corresponding eigenfunctions are $\zeta_{1}^{k}(\textbf{x})$, …, $\zeta_{m_{k}}^{k}(\textbf{x})$,…. Then, the reduced basis functions $\{\zeta_{l}^{k}(\textbf{x})\}_{l=1}^{m_{k}}$ have the following approximation property *

[TABLE]

where $X=L^{2}(D)$ or $X=H^{1}(D)$ and the number $m_{k}$ is determined according to the ratio $\rho=\frac{\sum_{s=1}^{m_{k}}\lambda_{s}}{\sum_{s=1}^{Q}\lambda_{s}}$.

In practice, we choose the first $m_{k}$ dominant reduced basis functions such that $\rho$ is close enough to $1$ to achieve a desired accuracy, say $\rho=99\%$. More details of the POD method can be found in [berkooz1993proper, sirovich1987turbulence]. Notice that reduced basis functions $\zeta^{k}_{0}(\textbf{x})$ and $\zeta^{k}_{l}(\textbf{x})$, $l=1,...,m_{k}$ approximately capture the mean profile and the fluctuation of multiscale basis functions associated with $\textbf{x}_{k}$, respectively. Thus, it is expected that for each realization of the random potential the associated multiscale basis functions can be approximated by the reduced basis functions, i.e.,

[TABLE]

Remark 3.2*.*

To construct the multiscale reduced basis functions, we partition the coarse grids $D^{k}$ into fine-scale quadrilateral elements with meshsize $h\ll\varepsilon$, which requires additional computational cost in the offline stage. However, the precomputed reduced basis functions can be used repeatedly to solve (1) for each realization of the random potential and different initial data, which results in considerable savings.

3.3 Estimation of the number of learning samples

We shall study the continuous dependence of multiscale basis functions on the random potential, which provide a guidence on how to determine the number of samples in the construction of multiscale basis functions. For notational simplification, we carry out the analysis for multiscale basis functions without localization.

Let $\varphi^{h}_{s}(\textbf{x})$, $s=1,...,N_{h}$ denote the finite element basis functions defined on fine mesh with size $h$ and $N_{h}$ is the number of fine-scale finite element basis functions. When we numerically solve (4)-(5), we represent the multiscale basis function as $\phi_{i}(\textbf{x})=\sum_{s=1}^{N_{h}}c_{s}\varphi^{h}_{s}(\textbf{x})$ and obtain the following quadratic programming problem with equality constraints

[TABLE]

where $\boldsymbol{c}=[c_{1},...,c_{N_{h}}]^{T}$ is the coefficients and $Q$ is a symmetric positive definite matrix on the fine triangularization $\mathcal{T}_{h}$ with the $(i,j)$ component

[TABLE]

In (19), $A$ is an $N_{h}$-by-$N_{H}$ matrix with $A_{ij}=(\varphi^{h}_{i},\varphi^{H}_{j})$ and $\boldsymbol{b}$ an $N_{h}$-by-$1$ vector with only the $i-$th entry being $1$ and others being [math].

The following result states the continuous dependence of multiscale basis functions on the random potential.

Theorem 3.2**.**

Assume the random potential $v^{\varepsilon}(\textbf{x},\omega)$ is almost surely bounded, i.e. (14) is satisfied and mesh size of the fine-scale triangles is small such that: (1) $h/\varepsilon=\kappa$ is small; and (2) $h^{d}\lVert v^{\varepsilon}(\cdot,\omega_{1})-v^{\varepsilon}(\cdot,\omega_{2})\rVert_{L^{\infty}(D)}<1$. Then for two realizations $\omega_{1}$ and $\omega_{2}$ of the random potential $v^{\varepsilon}(\textbf{x},\omega)$, the corresponding multiscale basis functions satisfy

[TABLE]

where the constant $C$ is independent of $h,\varepsilon$, and $\lVert v^{\varepsilon}(\cdot,\omega_{1})-v^{\varepsilon}(\cdot,\omega_{2})\rVert_{L^{\infty}(D)}$.

Proof.

Under the assumptions that $v^{\varepsilon}(\textbf{x},\omega)$ is almost surely bounded and $h/\varepsilon=\kappa$ is small, we know that $Q$ is a positive definite matrix. Moreover, we know that $A$ has full rank, i.e., $\mathrm{rank}(A)=N_{H}$. Therefore, the quadratic optimization problem (19) has a unique minimizer, satisfying the Karush-Kuhn-Tucker condition. Specifically, the unique minimizer of (19) can be explicitly written as

[TABLE]

For two realizations $\omega_{1}$ and $\omega_{2}$, we define $\delta V=Q_{1}-Q_{2}$. Then

[TABLE]

and thus

[TABLE]

We choose $h$ to be small enough such that $\lVert\delta V\rVert_{\infty}\leq 1$, and have

[TABLE]

and thus

[TABLE]

Therefore,

[TABLE]

By their definitions, we have

[TABLE]

and thus

[TABLE]

We complete the proof since $h/\varepsilon=\kappa$ and $\lVert\phi(\cdot,\omega_{2})-\phi(\cdot,\omega_{1})\rVert_{L^{\infty}(D)}\leq\lvert\boldsymbol{c}_{2}-\boldsymbol{c}_{1}\rvert_{\infty}$. ∎

Equipped with Theorem 3.2, we can estimate the number of samples in the construction of multiscale reduced basis functions. Suppose the random potential is of the form (16). For any $\delta>0$, we choose an integer $Q_{\delta}$ and a set of random samples $\{v^{\varepsilon}(\textbf{x},\omega_{q})\}_{q=1}^{Q_{\delta}}$ such that

[TABLE]

where the expectation is taken over the random variables in $v^{\varepsilon}_{m}(\textbf{x},\omega)$ of the form (16). We can give a way to choose the random samples $\{v^{\varepsilon}(\textbf{x},\omega_{q})\}_{q=1}^{Q_{\delta}}$ since the distribution of the random variables $\xi_{i}(\omega)$, $i=1,...,m$ is known.

For every $\textbf{x}_{k}\in\mathcal{N}$, let $\{\phi_{k}(\textbf{x},\omega_{q})\}_{q=1}^{Q_{\delta}}$ be the samples of multiscale basis functions associated with $\textbf{x}_{k}$. Then, we have

[TABLE]

Given parameters $\varepsilon$ and $h$, we choose $\delta$ and $Q_{\delta}$ so that the right-hand side of (26) is small. Then the space of multiscale basis functions can be well approximated by the samples of multiscale basis functions $\{\phi_{k}(\textbf{x},\omega_{q})\}_{q=1}^{Q_{\delta}}$ with controllable accuracy and the POD method is further applied to construct multiscale reduced basis functions.

3.4 Derivation of our method based on the multiscale reduced basis functions

In this section, we present our method for solving the random Schrödinger equation: in the physical space, we use the multiscale reduced basis functions obtained in §3.2; in the random space, we use the qMC method.

The implementation of the qMC method is fairly easy. For instance, given a set of qMC samples, expectation of the solution is approximated by

[TABLE]

where $n$ is the number of qMC samples. Details of the generation of qMC samples and its convergence analysis will be discussed in §4.

Now, we focus on how to approximate the wavefunction in the physical space for each qMC sample $\omega_{s}$. For each node point $\textbf{x}_{k}\in\mathcal{N}$, we have constructed a set of multiscale reduced basis functions $\{\zeta^{k}_{i}\}_{i=0}^{m_{k}}$ and represent the wavefunction by

[TABLE]

where $m_{k}$ is the number of multiscale reduced basis functions associated with the node $\textbf{x}_{k}$. In the Galerkin formulation, we have the following weak form

[TABLE]

where $\mathcal{H}(\textbf{x},\omega_{s})$ is a deterministic operator. To numerically solve (29), we introduce some notations. Let $S$, $M$, and $V(\omega_{s})$ be matrices with dimension $\sum_{k=1}^{N_{H}}(m_{k}+1)\times\sum_{k=1}^{N_{H}}(m_{k}+1)$. Their entries are given by

[TABLE]

Then, we can reduce the weak formulation (29) into the following ODE system

[TABLE]

where the column vector $\textbf{c}(t,\omega_{s})=(c^{1}_{0}(t,\omega_{s}),...,c^{1}_{m_{k}}(t,\omega_{s}),...,c^{N_{H}}_{0}(t,\omega_{s}),...,c^{N_{H}}_{m_{k}}(t,\omega_{s}))^{T}$ consisting of all expansion coefficients of the solution $\psi^{\varepsilon}(t,\textbf{x},\omega_{s})$ onto multiscale reduced basis functions. We can further rewrite (30) as

[TABLE]

where $B(\omega_{s})=M^{-1}A(\omega_{s})$ and $A(\omega_{s})=\frac{\varepsilon^{2}}{2}S+V(\omega_{s})$. In the end, we can solve the above ODE system using existing ODE solvers.

Before ending this section, we shall explain why we choose the qMC method to approximate the random space of the electron wavefunction. Since the parameterization of a random potential may have high dimension, i.e., $m$ is large in (15), non-intrusive methods, such as sparse grid method [bungartz2004sparse] and stochastic collocation method [nobile2008sparse], become prohibitively expensive to solve PDEs with random coefficients. Polynomial chaos expansion (PCE) methods [ghanem2003stochastic, xiu2003modeling] are also frequently used in the literature to solve PDEs with random coefficients. This type of methods is useful if the solution is sufficiently smooth in the random space with small dimensionality. The performance of MC method does not depend on the dimension of the random space. However, its convergence rate is merely $O(\frac{1}{\sqrt{n}})$. The convergence rate of the qMC method is better both theoretically and numerically; see (45) in Theorem 4.5. Therefore, we choose the qMC method and its implementation is almost the same as the MC method.

4 Convergence analysis

We shall analyze the approximation error of the proposed method, where the emphasis is put on computing functionals of the wavefunction.

4.1 Regularity of the wavefunction with respect to the random variables

Since the potential $v^{\varepsilon}(\textbf{x},\omega)$ in (1) is parametrized by $m$ random variables $\xi_{i}(\omega)$, $i=1,...,m$ in (16), i.e., $v_{m}^{\varepsilon}(\textbf{x},\omega)=v^{\varepsilon}(\textbf{x},\xi_{1}(\omega),...,\xi_{m}(\omega))$. The wavefunction $\psi^{\varepsilon}_{m}(t,\textbf{x},\omega)$ satisfies

[TABLE]

The Doob-Dynkin’s lemma implies the wavefunction $\psi^{\varepsilon}_{m}(t,\textbf{x},\omega)$ in (32) can also be represented by a functional of these random variables, i.e., $\psi^{\varepsilon}(t,\textbf{x},\omega)=\psi^{\varepsilon}(t,\textbf{x},\xi_{1}(\omega),...,\xi_{m}(\omega))$.

First of all, we analyze the error introduced by the parameterization of the random potential. We have the following estimate result.

Lemma 4.1**.**

The difference between wavefunctions to (32) and (1) satisfies

[TABLE]

Proof.

The difference $\delta\psi=\psi_{m}^{\varepsilon}-\psi^{\varepsilon}$ satisfies

[TABLE]

By a direct calculation, we have

[TABLE]

where $\rho(\omega)$ is the probability measure induced by the randomness in the potential (16) and thus

[TABLE]

Therefore, we obtain

[TABLE]

which completes the proof. ∎

To analyze the qMC method, it is crucial to bound the mixed first derivatives of $\psi_{m}^{\varepsilon}$ with respect to $\xi_{i}(\omega)$. Denote $\boldsymbol{\xi}(\omega)=(\xi_{1}(\omega),\cdots,\xi_{m}(\omega))^{T}$ for convenience. Let $\boldsymbol{\nu}=(\nu_{1},\cdots,\nu_{m})$ denote a multi-index of non-negative integers, with $\lvert\boldsymbol{\nu}\rvert=\sum_{j=1}^{m}\nu_{j}$ and $\lvert\boldsymbol{\nu}\rvert_{\infty}=\max_{1\leq j\leq m}\nu_{j}$ . The value of $\nu_{j}$ determines the number of derivatives to be taken with respect to $\xi_{j}$, and $\partial^{\boldsymbol{\nu}}\psi_{m}^{\varepsilon}$ denotes the mixed derivative of $\psi_{m}^{\varepsilon}$ with respect to all variables specified by the multi-index $\boldsymbol{\nu}$.

Lemma 4.2**.**

For any $\omega\in\Omega$, any time $T$, and for any multi-index $\boldsymbol{\nu}$ with $\lvert\boldsymbol{\nu}\rvert<\infty$, the partial derivative of $\psi_{m}^{\varepsilon}(t,\textbf{x},\omega)$ satisfies the following a-priori estimate

[TABLE]

Proof.

When $\lvert\boldsymbol{\nu}\rvert=1$, we take the derivative of (32) with respect to $\xi_{j}(\omega)$. Let $\partial_{j}\psi_{m}=\partial_{\xi_{j}}\psi_{m}^{\varepsilon}$ and $\partial_{j}v_{m}=\partial_{\xi_{j}}v_{m}^{\varepsilon}$, we have

[TABLE]

Thereafter, we have the following estimate by a direction calculation

[TABLE]

and

[TABLE]

When $\lvert\boldsymbol{\nu}\rvert\geq 2$, we have

[TABLE]

According to the definition of the random potential (15), we have $\partial^{\boldsymbol{\nu}-\boldsymbol{\mu}}v_{m}^{\varepsilon}=0$ if $\lvert\boldsymbol{\nu}-\boldsymbol{\mu}\rvert\geq 2$. Thus the above equation can be simplified as

[TABLE]

Similarly, we obtain

[TABLE]

and

[TABLE]

Now we are ready to prove the theorem by mathematical induction. From (35), we know that (34) holds for $\lvert\boldsymbol{\nu}\rvert=1$. Assume that (34) holds for $\boldsymbol{\mu}$ with $\lvert\boldsymbol{\nu}-\boldsymbol{\mu}\rvert=1$. Substituting this into (36) yields the desired estimate for the $\boldsymbol{\nu}$ case. ∎

Remark 4.1*.*

The above derivation is similar to that in [jin2019gaussian], where an estimate in $L^{2}(D,\Omega)$ norm is obtained. Here, for each random realization $\omega$, we have the esitmate (34) in $L^{2}(D)$ norm, which will be used to prove the convergence in qMC.

4.2 Main result of the error analysis

In the framework of uncertainty quantification, we are interested in computing some statistical quantities of the electron wavefunction. As such, we shall present the error analysis of our method in computing functionals of $\psi_{m}^{\varepsilon}$.

Let $\mathcal{G}(\cdot)$ be a continuous linear functional on $L^{2}(D)$, then there exists a constant $C_{\mathcal{G}}$ such that

[TABLE]

for all $u\in L^{2}(D)$. Consider the following integral

[TABLE]

with $F(\boldsymbol{\xi})=\mathcal{G}(\psi_{m}^{\varepsilon}(\cdot,\boldsymbol{\xi}))$. We approximate the integral over the unit cube by randomly shifted lattice rules

[TABLE]

where $\boldsymbol{z}\in\mathbb{N}^{m}$ is the (deterministic) generating vector and $\boldsymbol{\Delta}\in[0,1]^{m}$ is the random shift which is uniformly distributed over $[0,1]^{m}$. Notice that $m$ is the dimension of the random vector $\boldsymbol{\xi}$ in the random potential and $n$ is the number of the sample point in implementing the qMC method. The interested reader is referred to [dick2013high] for more details of the randomly shifted lattice rules in the qMC method.

Lemma 4.3**.**

Let $F$ be the integrand in (37). Given $m,n\in\mathbb{N}$ with $n\leq 10^{30}$, weights $\boldsymbol{\gamma}=(\gamma_{\boldsymbol{\mathfrak{u}}})_{\boldsymbol{\mathfrak{u}}\subset\mathbb{N}}$, a randomly shifted lattice rule with $n$ points in $m$ dimensions can be constructed by a component-by-component algorithm such that, for all $\lambda\in(1/2,1]$,

[TABLE]

with

[TABLE]

Proof.

The proof of this result is essentially an application of the Koksma-Hlawka inequality, which is the same as the proofs of Theorem 15, Theorem 16, and Theorem 17 in [graham2015quasi], or Theorem 5.10 in [dick2013high] with the following modification of estimates:

[TABLE]

with $\zeta(x)=\sum_{j=1}^{\infty}j^{-x}$ the Riemann zeta function, and $C^{*}=\lVert\mathcal{G}\rVert_{L^{2}(D)}$. ∎

To analyze the error of our method, we need to make some assumptions on the regularity of the eigenfunctions and the decay rate of the eigenvalues in the KL expansion (16) of the random potential.

Assumption 4.4**.**

1. (a)

\let

0

There exist $C>0$ and $\Theta>1$ such that $\lambda_{j}\leq Cj^{-\Theta}$ for $j\geq 1$; 2. (b)

\let

0

The Karhunen-Loéve eigenfunctions $v_{j}(\textbf{x})$ are continuous and there exist $C>0$ and $\eta\in[0,\frac{\Theta-1}{2\Theta})$ such that $\lVert v_{j}\rVert_{C^{0}(\bar{D})}\leq C\lambda_{j}^{-\eta}$ for $j\geq 1$; 3. (c)

\let

0

The sequence defined by $\sqrt{\lambda_{j}}\lVert v_{j}\rVert_{C^{0}(\bar{D})},\;j\geq 1$ satisfies $\sum_{\begin{subarray}{c}j\geq 1\end{subarray}}\left(\sqrt{\lambda_{j}}\lVert v_{j}\rVert_{C^{0}(\bar{D})}\right)^{p}<\infty$ for some $p\in(0,1]$, and $\sum_{\begin{subarray}{c}j\geq 1\end{subarray}}\sqrt{\lambda_{j}}\lVert v_{j}\rVert_{C^{0}(\bar{D})}<\frac{\varepsilon}{T}\sqrt{\varrho(\lambda)}$ for $\lambda\in(1/2,1]$.

Recall that $\psi^{\varepsilon}$ and $\psi^{\varepsilon}_{m}$ are solutions to (1) and (32), respectively. Denote $\psi_{H,m}^{\varepsilon}$ the solution obtained by our method using the multiscale reduced basis functions in the physical space and the qMC method in the random space. Under the assumptions for the random potential, we have the error estimate.

Theorem 4.5**.**

Consider the approximation of $\mathbb{E}[\mathcal{G}(\psi^{\varepsilon})]$ via qMC multiscale finite element methods, denoted by $Q_{m,n}(\cdot;\mathcal{G}(\psi_{H,m}^{\varepsilon}))$, where we assume $\psi^{\varepsilon}\in L^{2}(\Omega;H^{2}(D))$. A randomly shifted lattice rule $Q_{m,n}$ is applied to $\mathcal{G}(\psi_{m}^{\varepsilon})$. Then, we can bound the root-mean-square error with respect to the uniformly distributed shift $\boldsymbol{\Delta}\in[0,1]^{m}$ by

[TABLE]

for $0<\chi\leq(1/2-\eta)\Theta-1/2$, and with $r=1/p-1/2$ for $p\in(2/3,1]$ and $r=1-\delta$ for $p\leq 2/3$, with $\delta$ arbitrarily small. Here the constant $C$ is independent of $\varepsilon$, $m$, and $n$ but depends on $T$.

Proof.

The linearity of operator $\mathcal{G}$ implies

[TABLE]

Under the assumption $\psi^{\varepsilon}\in L^{2}(\Omega;H^{2}(D))$, we have, see for example [ChenMaZhang:prep],

[TABLE]

Under the assumptions (b) and (c) in Assumption 4.4, we have, based on Lemma 4.1,

[TABLE]

for all $0<\chi\leq(1/2-\eta)\Theta-1/2$. Detailed derivation is essentially the same as the proof of Theorem 8 in [graham2015quasi].

Finally, when applying the qMC method to (42), we need to analyze the error in the qMC method. We adopt the standard framework, i.e., the Koksma-Hlawka inequality. Under Assumption 4.4, we have, based on Lemma 4.2 and Lemma 4.3,

[TABLE]

where $r=1/p-1/2$ for $p\in(2/3,1]$ and $r=1-\delta$ for $p\leq 2/3$, with $\delta$ arbitrarily small. Detailed derivation is essentially the same as the proof of Theorem 20 in [graham2015quasi]. A combination of above estimates completes the proof. ∎

Remark 4.2*.*

The term $\frac{m^{-\chi}}{\varepsilon}$ in the error estimate (41) can be viewed as a modeling error. When the $m$-term KL truncation potential $v^{\varepsilon}_{m}(\textbf{x},\omega)$ in (16) provides an accurate approximation to the potential $v^{\varepsilon}(\textbf{x},\omega)$, the term $\frac{m^{-\chi}}{\varepsilon}$ becomes small or negligible.

Remark 4.3*.*

In §5, we will show the proposed method works well for a large class of random potentials, even when the eigenvalues in the KL expansion have a relatively slow decay rate. Therefore, Assumption 4.4 is a rather technical assumption for the convergence analysis of the proposed method.

Remark 4.4*.*

In the error analysis for the qMC method, we assume $\boldsymbol{\xi}=(\xi_{1},...,\xi_{m})\in[0,1]^{m}$ for notational convenience; see (37), where $\xi_{i}$ are i.i.d. uniform random variables. In the KL expansion (16) representation for $v_{m}^{\varepsilon}(\textbf{x},\omega)$, we choose $\xi_{i}\in[-\sqrt{3},\sqrt{3}]$, $i=1,...,m$ so that the conditions $\mathbb{E}\left[\xi_{i}\right]=0$, $\mathbb{E}\left[\xi_{i}\xi_{j}\right]=\delta_{ij}$ are satisfied. The same convergence result can be obtained with little modification of the current proof.

5 Numerical examples

In this section, we conduct numerical experiments to test the accuracy and the efficiency of our method. Specifically, we will present convergence tests with respect to the physical meshsize, the number of multiscale reduced basis functions, and the number of qMC samples. In addition, we will investigate the existence of Anderson localization in both 1D and 2D. For convenience, we first introduce $L^{2}$ norm and $H^{1}$ norm as

[TABLE]

In what follows, we compare the relative error between expectations of the numerical solution $\psi^{\varepsilon}_{\textrm{num}}$ and the reference solution $\psi^{\varepsilon}_{\textrm{ref}}$ in both $L^{2}$ norm and $H^{1}$ norm

[TABLE]

Here $\mathds{E}[\psi^{\varepsilon}_{\textrm{num}}]=\int_{\Omega}\psi^{\varepsilon}_{\textrm{num}}(t,\textbf{x},\omega)\mathrm{d}\rho(\omega)$, $\mathds{E}[\psi^{\varepsilon}_{\textrm{ref}}]=\int_{\Omega}\psi^{\varepsilon}_{\textrm{ref}}(t,\textbf{x},\omega)\mathrm{d}\rho(\omega)$, $\Omega$ is the random space, and $\rho(\omega)$ is the probability measure induced by the randomness in (16). The reference solution refers to the numerical wavefunction using a very fine mesh and a large amount of qMC samples. In numerical experiments, we use the MATLAB’s Statistics Toolbox to generate the Sobol sequence to implement the qMC method. When we use the POD method to construct multiscale reduced basis functions, we observed similar decay behaviors of the associated eigenvalues at each coarse grid point. Therefore, we choose the same reduced basis number $m_{k}$ for all the coarse grid points.

5.1 Convergence in the physical space

Consider the 1D Schrödinger equation over $D=[-\pi,\pi]$

[TABLE]

where the periodic condition is imposed, the initial data $\psi_{\textrm{in}}(x)=(\frac{10}{\pi})^{1/4}e^{-20(x-0)^{2}}$, and the random potential $v^{\varepsilon}(x,\omega)$ is defined as

[TABLE]

In the random potential (47), $\sigma$ is used to control the strength of the random potential, and $\xi_{j}(\omega)$’s are mean-zero and independent random variables uniformly distributed in $[-\sqrt{3},\sqrt{3}]$. Moreover, we choose $\varepsilon=\frac{1}{16}$, $\sigma=1$ and $E=[\frac{1}{9},\frac{1}{13},\frac{1}{11}]$, i.e., the characteristic length scale of randomness is different from the semiclassical parameter.

Convergence with respect to the coarse mesh size $H$. In our numerical test, we set the final computational time $T=1$. For the reference solution, we choose the fine mesh to be $h=\frac{2\pi}{2048}$ and the qMC sample number to be $n=16000$. For our method, we choose the POD modes $m_{k}=3$, the sampling number in the offline training stage to be $200$ and the number of qMC samples in the online stage to be $2560$.

In Table LABEL:table_physical, we compute the relative errors of the expectation of the wavefunction in both $L^{2}$ norm and $H^{1}$ norm for a series of coarse meshes with meshsize ranging from $H=\frac{2\pi}{32}$ to $H=\frac{2\pi}{256}$. Nice convergence in the physical space is observed.