Localized computation of eigenstates of random Schr\"odinger operators

TL;DR

This paper introduces a new numerical scheme for approximating localized low-energy eigenstates of random Schrödinger operators, leveraging a preconditioned inverse iteration with multigrid methods for efficient local computation.

Contribution

The paper presents a novel localized numerical method for eigenstates of random Schrödinger operators, combining inverse iteration with multigrid solvers for improved efficiency.

Findings

Method effectively approximates localized eigenstates.

Numerical experiments demonstrate good performance in 2D and 3D.

Applicable to nonlinear random Schrödinger operators.

Abstract

This paper concerns the numerical approximation of low-energy eigenstates of the linear random Schr\"odinger operator. Under oscillatory high-amplitude potentials with a sufficient degree of disorder it is known that these eigenstates localize in the sense of an exponential decay of their moduli. We propose a reliable numerical scheme which provides localized approximations of such localized states. The method is based on a preconditioned inverse iteration including an optimal multigrid solver which spreads information only locally. The practical performance of the approach is illustrated in various numerical experiments in two and three space dimensions and also for a non-linear random Schr\"odinger operator.

Equations

$$-\Delta u+Vu=\lambda u$$

$$-△u+Vu+\delta |u{|}^{2}u=\lambda u.$$

$$a(u,v):={\int }_D\nabla u(x)\cdot \nabla v(x)+V(x)u(x)v(x)\text{d}x=\lambda (u,v)$$

$$\Vert v{\Vert }_V^2:=(Vv,v)={\int }_{D}V(x)|v(x){|}^2\text{d}x.$$

$$|||v{|||}^2:=a(v,v)=\Vert \nabla v{\Vert }^2+\Vert v{\Vert }_V^2.$$

$$\lambda (v):=\frac{a(v,v)}{\Vert v{\Vert }^2}=\frac{|||v{|||}^2}{\Vert v{\Vert }^2}.$$

$$|||u{|||}_{D\setminus B_{p\epsilon L}^{\infty }(\mathrm{supp}f)}\le c{\gamma }^p|||u|||$$

$$a(\psi ,v)=(1,v)$$

$$A_h{u}_h={\lambda }_h{M}_h{u}_h.$$

$$(A_h{)}_{ij}:=a({\phi }_{z_i}^h,{\phi }_{z_j}^h),(M_h{)}_{ij}:=({\phi }_{z_i}^h,{\phi }_{z_j}^h).$$

$$\Vert D^{2}u\Vert \le \Vert \Delta u\Vert =\Vert \lambda u-Vu\Vert \le \lambda +\beta ,$$

$$\frac{{\lambda }_{h,1}-{\lambda }_1}{{\lambda }_1}\lesssim \frac{h^2}{{\epsilon }^2}+\frac{h^4}{{\epsilon }^4},$$

$${\psi }_h=A_h^{-1}{M}_{h}1,$$

$$|u_1(x)|\le {\lambda }_1|\psi (x)|.$$

$$f_{p}(k)=\log\big{(}\|(A_{h}^{-1}M_{h})^{p}e_{k}\|_{2}\big{)}.$$

$$f_{R,p}(k)=\log\big{(}\|e_{k}^{T}(A_{h}^{-1}M_{h})^{p}R\|_{2}\big{)}.$$

$${\ell }_{\epsilon }=7,{\ell }_h=9,{\ell }_H=6,\alpha =1,\beta =5\cdot {\epsilon }^{-2}.$$

$${\ell }_{\epsilon }=5,{\ell }_h=6,{\ell }_H=4,\alpha =1,\beta =5\cdot {\epsilon }^{-2}.$$

$$-\Delta u+Vu+\delta |u{|}^{2}u=\lambda u$$

$$A_h{u}_h+\delta N_h(u_h)u_h=\lambda M_h{u}_h,$$

$$u_h^{n+1}=(A_h+\delta N_h(u_h^n){)}^{-1}{M}_h{u}_h^n$$

Full text

Keywords. random potential, iterative eigensolver, multigrid solver

AMS subject classifications. 65N25, 65N55, 47B80

1. Introduction

This paper concerns the numerical approximation of essentially localized eigenstates of the linear Schrödinger eigenvalue problem

[TABLE]

on a bounded domain $D\subseteq\mathbb{R}^{d}$ with homogeneous Dirichlet boundary conditions. The non-negative variable coefficient $V$ represents an external potential reflecting a high degree of disorder. This apparently simple problem is relevant in the context of quantum-physical processes related to ultracold bosonic or photonic gases, known as Bose-Einstein condensates [Bos24, Ein24, DGPS99, PS03]. A Bose-Einstein condensate (BEC) is an extreme state of matter formed by a dilute gas of bosons at ultra-cold temperatures, very close to absolute zero. In a BEC, individual particles (i.e., their wave packages) overlap, lose their identity, and form one single super atom. BECs allow to study macroscopic quantum phenomena such as superfluity (i.e., the frictionless flow of a fluid) on an observable scale. When BECs are trapped in a highly oscillatory high amplitude potential $V$ that exhibits a sufficiently large degree of disorder, the low-energy stationary states essentially localize in the sense of an exponential decay of their moduli. This localization is a universal wave phenomenon referred to as Anderson localization [And58]. For BECs, it has been observed experimentally in [FFI08].

On a mathematical level, the formation of stationary quantum states can be modeled by the slightly more involved Gross-Pitaevskii eigenvalue problem (GPEVP). In non-dimensional form, the GPEVP seeks $L^{2}$-normalised eigenfunctions $u\in H_{0}^{1}(D)$ and corresponding eigenvalues (so-called chemical potentials) $\lambda\in\mathbb{R}$ such that

[TABLE]

The parameter $\delta\geq 0$ resembles the strength of repulsive particle interactions depending on physical properties of the particles that form the BEC. The linear case  (1.1) corresponds to the regime of vanishing particle interaction ($\delta=0$). Since the interactions are weak for BECs, the linear case provides very good approximations of actual BEC states [RS13]. This is why we will mostly consider the numerical solution of the linear eigenvalue problem (1.1). However, the numerical techniques can be generalized to the nonlinear setting which will be demonstrated through numerical experiments at the end of this paper.

The numerical approximation of localized Schrödinger eigenstates has recently caused a large interest in the fields of computational physics and scientific computing. Fundamentally novel methodologies have been developed to cope with the intrinsic difficulty coming from the multiscale nature of the potential [ADJ*+*16, Ste17, ADF*+*19, XZO18]. They are based on the link of the groundstate $u_{1}$, normalized by $\|u_{1}\|_{L^{\infty}(D)}=1$, and the so-called landscape function $\psi\in H^{1}_{0}(D)$ defined through of the homogeneous elliptic equation $-\Delta\psi+V\psi=1$, cf. [FM12]. The landscape function $\psi$ gives rise to surprisingly sharp eigenvalue bounds and its local maxima indicate regions where localization may occur. While the new techniques based on landscape functions are empirically successful in predicting the eigenvalues they lack any control on the accuracy of the approximation. In particular the approximation of the states is rather sketchy.

The present paper aims at a more sophisticated method leading not only to the regions of localization but also to actual approximations of the lowermost eigenstates. The starting point is the recent rigorous a priori prediction of exponentially localized low-energy states caused by the interplay of disorder (randomness) and high amplitude (contrast) of the potential trap $V$ [AHP18]. The results of [AHP18] employ the convergence analysis of a preconditioned inverse iteration in the spirit of iterative numerical homogenization [KY16] and thereby provides a role model for efficient simulation. The main idea is to exploit the localization property, meaning that these eigenfunctions can be approximated in a sophisticated manner by functions supported in the union of only a few small sub-domains. The main contribution of this paper is the finding of an appropriate starting subspace and the construction of a preconditioned eigensolver which only relies on local operations. In this way, we are able to approximate the eigenstates of lowest energy roughly at the same costs as the computation of the landscape function $\psi$. In addition, the algorithm is applicable for the nonlinear case.

The first step is a finite element discretization of the eigenvalue problem (1.1) which we discuss in Section 2. For this, we consider a uniform refinement of the $\varepsilon$-mesh on which the potential $V$ is defined. This automatically provides a mesh hierarchy for which we define a multigrid based preconditioner. Using only a small number of smoothing steps and no direct solves, one multigrid step preserves locality in the sense that the support of the resulting function is only slightly larger than the initial function. This then leads to a localization preserving eigensolver introduced in Section 3. In other words, the simple yet efficient trick is to consider a multigrid preconditioner on a hierarchy of meshes starting from the $\varepsilon$-level. This sufficed to be optimal in the sense of an $\mathcal{O}(1)$ condition number of the preconditioned operator. If no direct solver is used on the coarsest level, this preconditioner prevents the otherwise global communication of a standard multigrid involving coarser levels. The overall algorithm consists of two parts. First, we apply in parallel the localization preserving iteration scheme to finite element basis functions on a coarse grid. Based on the energies we select the most promising candidates indicating the regions of localization. Second, we consider a Ritz-Rayleigh iteration on the remaining local functions, i.e., we project the eigenvalue problem on a very small subspace. This procedure is applicable in any dimension and thus, allows to compute eigenfunctions for $d=3$ where standard eigensolvers break. This and further numerical experiments are subject of Section 4.

2. Schrödinger Eigenvalue Problem

This section is devoted to the linear Schrödinger eigenvalue problem (1.1) and the introduction of needed finite element spaces. Further, we discuss localization effects of the lowermost eigenfunctions if the potential contains disorder.

2.1. Model problem

In this paper, we consider the $d$-dimensional linear eigenvalue problem of Schrödinger type with a potential being highly oscillatory and of large amplitude. In particular, we assume that the potential includes some kind of disorder such that the first few eigenfunctions of lowest energy localize.

The variational formulation corresponding to the eigenvalue problem (1.1) reads as follows: Given a non-negative potential $0\leq V\in L^{\infty}(D)$, find non-trivial eigenpairs $(u,\lambda)\in\mathcal{V}\times\mathbb{R}$ with search and test space $\mathcal{V}:=H^{1}_{0}(D)$ such that

[TABLE]

for all test functions $v\in\mathcal{V}$. Here, $(\cdot,\cdot)$ denotes the $L^{2}$-inner product on $D$ and eigenfunctions are assumed to be normalized in the $L^{2}$-norm. Further, we assume that the eigenvalue are ordered, i.e., $0<\lambda_{1}<\lambda_{2}\leq\dots$.

As far as the theory is concerned, we focus on a representative class of potentials which are piecewise constant with respect to a mesh $\mathcal{T}^{\varepsilon}$ consisting of cubes with side length $\varepsilon\ll 1$. For the sake of simplicity we assume $\varepsilon=2^{-\ell_{\varepsilon}}$ where $\ell_{\varepsilon}$ denotes the level of the $\varepsilon$-scale. On each cube in $\mathcal{T}^{\varepsilon}$, the potential takes (randomly) a value in the interval $[\alpha,\beta]$ with moderate $0\leq\alpha\approx 1$ and large $\beta$ in the sense of $\beta\gtrsim\varepsilon^{-2}$. We emphasize that the resulting potentials are highly oscillatory due to the underlying mesh on $\varepsilon$-scale. Examples of such random potentials are shown in Figure 2.1. In the field of matter waves, one often considers disorder potentials created optically by using speckle patterns [BMRF*+*99, LFM*+*05, FFI08]. This means that a laser beam is transmitted trough a diffusive plate forming randomized and high-contrast patterns. Within this paper, we imitate such speckles by piecewise constant potentials with respect to a quadrilateral mesh consisting of cubes with side length $\varepsilon$. On each of these cubes the potential takes random values following certain statistical assumptions, cf. [Goo75, DKW08].

Within this paper, $\|\cdot\|:=\sqrt{(\cdot,\cdot)}$ denotes the standard $L^{2}$-norm on $D$, whereas the $V$-weighted $L^{2}$-norm is given by

[TABLE]

Corresponding to the Schrödinger operator, we define the energy norm as

[TABLE]

The energy of a function is characterized by the Rayleigh quotient, namely

[TABLE]

In case $v$ is an eigenfunction of the Schrödinger eigenvalue problem (2.1), the energy $\lambda(v)$ is equal to the corresponding eigenvalue. The eigenfunction of minimal energy $\lambda_{1}:=\min_{v\in\mathcal{V}\setminus\{0\}}\lambda(v)$ is called the groundstate.

2.2. Exponential decay of the Green’s functions

For certain classes of potentials it is known that the corresponding Green’s function of the Schrödinger operator decays exponentially. For constant potentials $V\equiv\beta$ with sufficiently large $\beta$ this was shown in [Glo11, Lem. 3.2].

Piecewise constant potentials with two values $\alpha$ and $\beta$ were considered in [AHP18]. For this, an operator preconditioner was constructed depending on the geometric structure of the potential, i.e., on the interaction of $\alpha$-valleys and $\beta$-peaks. It was designed in such a way that the application of this operator preconditioner only enlarges the support by a small amount of $\varepsilon$-layers within $\mathcal{T}^{\varepsilon}$. With this, one can show that the weak solution $u\in\mathcal{V}$ of the variational problem $a(u,\cdot)=(f,\cdot)$ decays exponentially fast around the support of $f$ for potentials under certain statistical assumptions (including, e.g., 1D-tensorized potentials as in Figure 2.1). More precisely, this means that

[TABLE]

for constants $c>0$ and $0<\gamma<1$ independent of $\varepsilon$, $L\approx\log(1/\varepsilon)$, and $B^{\infty}_{r}(z)$ denoting the ball of radius $r$ around $z$. We emphasize that this result only relies on $\varepsilon$ being small (oscillatory) and $\beta$ being large (high amplitude). Thus, disorder does not play a role for the exponential decay of the Green’s function.

2.3. Localization of eigenfunctions

For the localization of eigenfunctions, the potential needs to include a certain degree of disorder. For a periodic potential the lowermost part of the spectrum of the Schrödinger operator is clustered. In the one-dimensional case, which is illustrated in Figure 2.2, one can prove that there exist about $\varepsilon^{-d}/2$ (the number of potential valleys) eigenvalues in the energy range of $\varepsilon^{-2}$. Disorder changes the picture dramatically and leads to significant spectral gaps already within the first few eigenvalues. In one space dimension there is a one-to-one correspondence between the largest $\alpha$-valleys of the potential and the smallest eigenvalues. This is illustrated by dotted vertical lines which correspond to the largest $\alpha$-valleys and coincide with the spectral gaps for the random potential. This then leads to the exponential decay of the lowermost eigenfunctions, cf. Figure 2.2.

In [AHP18] the localization of eigenfunctions was rigorously proven for the regime $\beta\gtrsim\varepsilon^{-2}$ and certain statistical assumptions on the potential. The key ingredient was to prove the existence of spectral gaps in the presence of disorder in combination with a preconditioned block inverse iteration. For this, particular finite element spaces and quasi-interpolation operators were considered. Further, the connection of valley sizes and eigenvalues was used. More precisely, we considered the first eigenfunctions of the Laplacian in the largest $\alpha$-valleys with homogeneous Dirichlet boundary conditions to obtain eigenvalue estimates of the Schrödinger operator. These Laplace eigenfunctions also serve as stating functions within the preconditioned block inverse iteration. The claimed decay then follows from the exponential convergence with respect to the number of iteration steps.

The exponential decay may also be characterized a posteriori in terms of the landscape function $\psi\in\mathcal{V}$ defined through

[TABLE]

for all test functions $v\in\mathcal{V}$, cf. [FM12]. The connection to the eigenvalue problem (1.1) gets visible by the possible reformulation as an eigenvalue problem with the effective confining potential $1/\psi$, which encodes the decay of the eigenfunction in some Agmon measure [ADJ*+*16]. This means that the eigenfunction reduces by a certain factor when crossing a valley of $1/\psi$. This interplay has also been analyzed in [Ste17] using averages over local Brownian motion paths.

2.4. Finite element discretization

In order to approximate the (localized) eigenstates of the Schrödinger operator, we consider a finite element discretization of (2.1). For this, we consider a uniform refinement of $\mathcal{T}^{\varepsilon}$, namely $\mathcal{T}^{h}$, consisting of cubes with side length $h=2^{-\ell_{h}}\leq\varepsilon$. The corresponding set of nodes is denoted by $\mathcal{N}^{h}$ and the set of interior nodes by $\mathcal{N}^{h}_{0}$.

Based on $\mathcal{T}^{h}$, we define $V_{h}:=Q_{1}(\mathcal{T}^{h})\cap\mathcal{V}\subseteq\mathcal{V}$ as the conforming $Q_{1}$-finite element space, cf. [BS08, Sect. 3.5]. This space has the dimension $n:=|\mathcal{N}^{h}_{0}|$ and is spanned by the standard $Q_{1}$ hat functions, which we denote by $\varphi^{h}_{z}$ for $z\in\mathcal{N}^{h}_{0}$. Recall that $\varphi^{h}_{z}$ is a piecewise polynomial of partial degree one with $\varphi^{h}_{z}(z)=1$ and $\varphi^{h}_{z}(w)=0$ for any other node $w\in\mathcal{N}^{h}\setminus\{z\}$. Clearly, hat functions can also be defined on the original mesh $\mathcal{T}^{\varepsilon}$ w.r.t. its set of interior nodes $\mathcal{N}^{\varepsilon}_{0}$. For these hat functions we write $\varphi^{\varepsilon}_{z}$ for $z\in\mathcal{N}^{\varepsilon}_{0}$.

We now apply a Galerkin ansatz to the eigenvalue problem (2.1). Thus, we restrict the trial and test space to $V_{h}$, which then leads to a discrete eigenvalue problem of the form

[TABLE]

Here, $A_{h}\in\mathbb{R}^{n,n}$ and $M_{h}\in\mathbb{R}^{n,n}$ denote the symmetric stiffness and mass matrices defined through

[TABLE]

Note that the stiffness matrix already includes the potential. Eigenpairs of the matrix eigenvalue problem (2.3) approximate eigenpairs of the PDE eigenvalue problem (2.1). Due to the min-max principle of eigenvalues based on the Rayleigh quotient, we know that $\lambda_{1}\leq\lambda_{h,1}$.

As mentioned above, we consider highly oscillatory potentials meaning that $\varepsilon$ is small. Further, the mesh size $h$ needs to sufficiently small compared to $\varepsilon$ in order to resolve the oscillations of the potential and thus, guarantees a reasonable approximation of the eigenpairs. Such a condition on the minimal resolution are justified by standard a priori error analysis. If $D$ is convex, any eigenfunction $u\in H^{1}_{0}(D)$ is $H^{2}$ regular. When normalized in $L^{2}(D)$, $u$ satisfies the bound

[TABLE]

where $\lambda$ is the eigenvalue it corresponds to. For the lowermost eigenvalue in the regime of [AHP18] where $\lambda_{1}\approx\beta\approx\varepsilon^{-2}$ this leads to an error bound

[TABLE]

see, e.g., [SF73, Lem. 6.1]. While the hidden constant in this bound depends only on the quasi-uniformity of the mesh, multiplicative constants in bounds for larger eigenvalues or any eigenfunction may deteriorate with the distance to neighboring eigenvalues in the case of clustered eigenvalues [Yse13]. In this sense the resolution condition $h\lesssim\varepsilon$ is minimal for the approximation of the lowermost eigenvalue and may be much more restrictive in other cases. As a result, the eigenvalue problem (2.3) is of large dimension, which makes the solution computationally costly.

As we are interested in the smallest eigenvalues and the corresponding localized eigenstates, we aim to construct a preconditioner based on local operations. Such a localization preserving preconditioner then allows parallel computations and is subject of the following section.

3. Localization Preserving Preconditioner

In [AHP18] we have presented a locally operating preconditioner based on a domain decomposition of $D$, which was a key ingredient for the proof of the exponential decay of the first eigenstates. The construction, however, was tailored to theoretical needs and checkerboard potentials. It included exact projections of local sub-domains. For actual computations, it turns out that already a simple multigrid preconditioned iteration may act as a localized version of a block inverse power method leading to sophisticated approximation results. In order to preserve locality, the coarsest level used within the multigrid cycle is the mesh $\mathcal{T}^{\varepsilon}$ on which the potential is defined. This still yields an optimal preconditioner due to the properties of the Green’s function of the Schrödinger operator.

Further, we discuss how to find a low-dimensional starting subspace of local functions with which we can initialize the eigenvalue iteration. This is a crucial step in order to obtain local approximations of the lowermost eigenfunctions.

3.1. Multigrid preconditioner

Recall that we have assumed that the finite element mesh $\mathcal{T}^{h}$ is defined by a uniform refinement of the mesh on $\varepsilon$-level on which the potential is defined. Thus, we have a mesh hierarchy given by $\mathcal{T}^{h}$, possible intermediate meshes, and $\mathcal{T}^{\varepsilon}$. For this hierarchy, we define a standard geometric multigrid method without direct solve on the coarsest level, which will serve as a preconditioner for the eigenvalue iteration. A major aspect is that the multigrid method maintains locality, i.e., for a local right-hand side $b$ the resulting approximate solution of $A_{h}^{-1}b$ is supported on a domain which is only slightly larger than the support of $b$.

Given the mesh hierarchy, we consider a so-called V-cycle with only one smoothing step on each level, cf. [Hac85, McC87, Yse93]. Starting with a vanishing starting vector, we consider the residuals of $A_{h}x=b$ on different levels of the hierarchy. On the coarsest level, i.e., on $\mathcal{T}^{\varepsilon}$, we run one single relaxation step. This means that we do not use a direct solver on $\varepsilon$-level with the stiffness matrix $A_{\varepsilon}$ but use instead a single step of the Jacobi iteration. This yields a reasonable approximation, since the condition number of $A_{\varepsilon}$ is of order $1$. The reason for this is the shift by the potential in the bilinear form $a$ in (2.1) and the fact that $\beta\gtrsim\varepsilon^{-2}$. A computational study of the condition numbers of the stiffness matrices $A_{h}$ and $A_{\varepsilon}$ is shown in Figure 3.1. We emphasize that a direct solver would ruin the locality immediately.

On each level of the multigrid cycle the support of the iterate grows slightly due to the application of the stiffness matrix on the respective level. In total, the application of this preconditioner enlarges the support by strictly less than three layers of $\varepsilon$-cubes.

It is well-known that multigrid solvers with a reduced tolerance (we only perform one V-cycle with a single relaxation step per level) can be used efficiently as a preconditioner within an external iterative solver, leading to methods such as pcg or lopcg, cf. [KN03]. To keep things simple, we concentrate on pcg. We denote the application of $j$ steps of pcg with the multigrid preconditioner by $\mathcal{P}_{j}$. Thus, $\mathcal{P}_{j}(b)$ yields an approximation of $A_{h}^{-1}b$, where the accuracy can be controlled by the number of iteration steps. We emphasize that $\mathcal{P}_{j}$ preserves locality in the following sense.

Theorem 3.1** (Localization preserving property of $\mathcal{P}_{j}$).**

Assume that $x\in\mathbb{R}^{n}$ is the coefficient vector of a local function $u_{h}\in V_{h}\subset\mathcal{V}$. Then, $\mathcal{P}_{j}(x)$ is the representative of a function in $V_{h}$ with support being at most $3j$ layers of $\varepsilon$-cubes larger than $\operatorname{supp}(u_{h})$.

Proof.

As already mentioned, the application of a V-cycle as described above affects less than three layers of $\varepsilon$-cubes. More precisely, the support enlarges at most within a ball of radius ${3\varepsilon-h}$ around $\operatorname{supp}(u_{h})$. The cg step involves another application of the stiffness matrix on the finest level and thus, adds only one layer of $h$-cubes to the support. In total this leads to a growth of $3$ layers with side length $\varepsilon$ per pcg step. ∎

Remark 3.2*.*

Note that the choice of the preconditioner is by far not unique and may be replaced, e.g., by local variants of hierarchical basis [Yse86] or BPX [BPX90] preconditioners or Gamblets [XZO18]. Here, local means that $\mathcal{T}^{\varepsilon}$ is the coarsest level in the underlying hierarchy of meshes.

3.2. Starting subspace

For an efficient eigenvalue iteration we also need a suitable starting subspace. Starting for example with a vector which has only positive entries, we will approach the groudstate but every iteration step is ’global’. Instead, we search for a number of local functions at positions where we expect localization of the first eigenfunctions. For the detection of a suitable starting subspace we discuss several approaches.

3.2.1. Landscape function

One possibility to find spots of localization is to compute the already mentioned landscape function $\psi\in\mathcal{V}$, cf. [FM12]. This function is defined as the solution of the Schrödinger source problem with right-hand side $1$ and homogeneous Dirichlet boundary conditions, cf. equation (2.2). In other words, $\psi$ is the outcome of one step of the inverse power method in the PDE setting. The corresponding discrete approximation $\psi_{h}$ satisfies

[TABLE]

where $\mathbf{1}=[1,\ 1,\ \dots,1]^{T}\in\mathbb{R}^{n}$. The landscape function shall indicate where the first eigenfunctions may localize. To see this, let $u_{1}$ denote the normalized ground state of the Schrödinger eigenvalue problem, i.e.,  $\|u_{1}\|_{L^{\infty}(D)}=1$. Then, it is shown in [FM12] that the landscape function satisfies pointwise

[TABLE]

This bound implies that in regions where $\psi$ is small compared to the smallest eigenvalue $\lambda_{1}$, the eigenfunction $u_{1}$ needs to be small as well. Considering the peaks of the landscape function then provides empirically accurate predictions where to expect localization. However, since we know that $\lambda_{1}\gtrsim\varepsilon^{-2}$, the estimate may degenerate quickly to $|u_{1}(x)|\leq 1=\|u_{1}\|_{L^{\infty}(D)}$. Thus, this approach does not allow rigorous predictions a priori but may serve as an indicator.

Although the estimate (3.1) only includes the groundstate, the landscape function has been used to locate a larger number of eigenfunctions by analyzing the local minima of $\psi$, cf. [ADF*+*19]. Further, it may be used to partition the computational domain $D$ into a network of valleys, cf. Figure 3.2. Yet another possibility is to consider multiple applications of the inverse Schrödinger operator, i.e., $(A_{h}^{-1}M_{h})^{k}\mathbf{1}$ in the discrete setting, cf. [Ste17].

3.2.2. Randomized landscape function

A similar approach to detect positions of localization was introduced in [LS18]. Therein, the inverse power method is applied to the unit vectors of dimension $n$. More precisely, for a fixed parameter $p\in\mathbb{N}$ and $e_{k}\in\mathbb{R}^{n}$ denoting the $k$-th unit vector, one computes for $k=1,\dots,n$

[TABLE]

In [LS18] it was shown that highly localized eigenfunctions correspond to metastable states of the power iteration and thus, are directly connected to local maxima of the function $f_{p}$. Note that in the given setting $f_{p}$ computes the inverse iteration to all hat functions on the $h$-scale, which is expensive due to the size of the eigenvalue problem.

A variant is the following randomized version: Given a random matrix $R\in\mathbb{R}^{n,m}$ with $m\ll n$ we compute

[TABLE]

Thus, we apply the inverse power method to random vectors and consider the means in each component. In other words, we replace the application of $A_{h}^{-1}M_{h}$ to $n$ local functions by the application to $m$ global functions [LS18]. The outcome of this approach is displayed in Figure 3.3. It shows a landscape function which seems more smoothed than the previous approach of Section 3.2.1. Further, an increase of the number of iterations does not focus purely on the groundstate as this would happen by applying the inverse power method to the vector $M_{h}\,\mathbf{1}$.

3.2.3. Largest valleys

Assume that the potential $V$ is given in form of a random checkerboard, i.e., the function values of $V$ are either $\alpha$ or $\beta$, each with probability $0.5$. Then, motivated by the theoretical findings in [AHP18], we may start with a number of hat functions distributed in the largest valleys. Here, a valley denotes a cube in which the potential is constant $\alpha$. For special classes of such checkerboard potentials, this approach guarantees to yield an appropriate starting subspace, meaning that the number of initial functions is $\mathcal{O}(1)$, that these functions are local, and that the inverse power method converges quickly, cf. [AHP18].

However, the search for valleys comes with the drawback that one first needs to analyze the structure of the given potential. Further, a generalization of the term valley is needed to consider general potentials. Nevertheless, the idea motivates the subsequent approach using uniformly distributed hat functions.

3.2.4. Uniformly distributed hat functions

The current approach combines the benefits of the previous subsections. In this manner, we obtain a suitable subspace of local functions which is highly eligible to serve as a starting point within a (preconditioned) eigenvalue iteration. At the same time, we restrict the computational costs such that it is comparable to a single step of the inverse power method.

Let $\mathcal{T}^{H}$ be a mesh of cubes with side length $H=2^{-\ell_{H}}\geq\varepsilon$ such that $\mathcal{T}^{\varepsilon}$ is a refinement of $\mathcal{T}^{H}$ and thus, $\ell_{H}\leq\ell_{\varepsilon}$. For each interior node $z\in\mathcal{N}^{H}_{0}$ there exists a $Q_{1}$ hat function $\varphi_{z}^{H}$. Similar to the approach in Section 3.2.2 we may now apply several steps of the inverse power method. Instead, we apply a single pcg iteration step as introduced in Section 3.1 to each hat function. Note that all these computations may be performed in parallel and that the application of $\mathcal{P}_{1}$ maintains locality in the sense that the support of $\mathcal{P}_{1}(\varphi_{z}^{H})$ is only slightly larger than the support of $\varphi_{z}^{H}$, cf. Theorem 3.1. Although we perform only one pcg step, we gain sufficient information in order to rank the importance of the basis function. For this we compute the energies $\lambda(\mathcal{P}_{1}(\varphi_{z}^{H}))$ and sort them in an increasing order. Afterwards, we only keep the functions of lowest energy defined through a fixed ratio $0<\eta<1$. This procedure is then repeated until the number of functions is small enough.

We emphasize that the support of each function only grows gently in each step whereas the number of functions decreases by the prescribed factor $\eta$. In other words, the proposed algorithm is approximately as costly as computing the landscape function $\psi$, cf. Section 3.2.1. Here, however, we directly obtain a number of candidates where the first eigenfunctions will localize without the need of constructing a network structure. An illustration of the algorithm is given in Figure 3.4.

3.3. Approximation of eigenfunctions

In this subsection we finally introduce an algorithm to approximate the eigenfunctions of lowest energy of the linear Schrödinger eigenvalue problem (2.1) under disorder potentials. Thus, we do not only intend to find possible regions of localization but actually compute approximations of the smallest eigenvalues and their corresponding eigenstates.

Recall the definition of $\mathcal{P}_{j}$ in Section 3.1 containing $j$ pcg steps preconditioned my a multigrid cycle. Further, we fix the following parameters.

[TABLE]

The proposed algorithm consists of two parts: In the pre-iteration, we apply the selection process described in Section 3.2.4. For this, we consider hat functions on a (coarser) mesh of level $\ell_{H}$. Using the hat functions corresponding to $\mathcal{T}^{H}$ rather than $\mathcal{T}^{\varepsilon}$ reduces the number of initial functions. Further, we truncate these hat functions in order improve the locality. After $K_{\text{pre}}$ iteration steps with ratio $\eta$ we then have a much smaller number of functions, which already give a rough approximation of certain eigenfunctions. Most notably, however, they indicate with high precision where localization will take place without the need of constructing an additional mesh or finding local minima as with the landscape approach. If we are only interested in the groundstate, then one may reduce the number of functions until only a single function is left. Otherwise, we need to keep a sufficiently large amount of functions within the iteration process in order to prevent that only the groundstate is approximated. Further, since we only perform rough approximations of the inverse power method, we cannot guarantee that the functions of lowest energy correspond to the lowermost eigenstates. Therefore, we always keep at least $3\,N_{\text{eigs}}$ functions within the pre-iteration.

In the post-iteration, we combine all functions which were selected by the pre-iteration in form of a Ritz-Rayleigh approximation. This means that we first perform three preconditioned cg steps to each function and then project the mass and stiffness matrices to the span of these functions and solve a (small) eigenvalue problem of the form $\mathbf{M}\alpha=\mu\mathbf{A}\alpha$. Note that the dimension of the eigenvalue problem depends on the number of pre-iteration steps. The first eigenvector $\alpha_{1}$ contains the coefficients of the optimal linear combination of the ansatz functions and provides the approximate groundstate. The second eigenvector $\alpha_{2}$ then defines the function which serves as approximation of the second eigenstate and so on. Further, we decrease the number of ansatz functions by cutting off the candidates of highest energy. Similar to the pre-iteration, we always keep a certain number of functions, namely twice the amount of eigenfunctions we are interested in, i.e., $2\,N_{\text{eigs}}$. Note that, in contrast to the parallel steps within the pre-iteration, we consider here a combined iteration of all remaining candidates.

A summary of the procedure is given in Algorithm 1 and obtained numerical results for the two- and three-dimensional case are discussed in the subsequent section.

Remark 3.3*.*

The simple pcg steps in the pre-iteration of Algorithm 1 (line 6) may also be replaced by $K_{\text{pre}}$ steps of lopcg, cf. [Kny00, Kny01]. Since the main purpose of the pre-iteration is the detection of regions of localization, we consider here only this simple variant.

Remark 3.4*.*

In order to accelerate the iteration and improve the stability of the method, one may consider additional cut-offs. In particular, one may set components of $\alpha_{j}$ in line 18 of Algorithm 1 to zero if they are beneath a certain threshold. This also increases the level of locality.

The procedure of Algorithm 1 assumes that the first eigenfunctions are indeed localized. Thus, global eigenstates as they would appear for periodic potentials are not well-approximated by this method. However, one may adapt the selection process in the algorithm to such an extent that periodic structures are at least detected. For this, one may substitute the fixed-ratio criterion by a selection step which is based on the energies of the considered functions. In a periodic structure, where all hat functions would have a comparable energy, we would then not cut off any functions, meaning that the parameter $m$ does not decrease and that the eigenvalue problem in line 17 of the algorithm is still of large dimension.

4. Numerical Examples

We perform several numerical tests in order to explore the performance and feasibility of the proposed method. In the first experiment we consider the two-dimensional random potential, which was already used in Section 3 for the illustration of the approaches to find regions of localization. Second, we consider a three-dimensional speckle potential leading to huge eigenvalue problems which quickly overcharge standard eigenvalue solvers. Finally, we apply the method to the nonlinear Gross-Pitaevskii eigenvalue problem.

All computations have been performed with Matlab on an HPC Infiniband cluster (1.7 TB RAM, 2 Tesla V100 32GB GPUs). The reference solutions, which are used for the error plots, are obtained by the Matlab solver eigs with tolerance $10^{-12}$ for the 2D examples and $10^{-10}$ in 3D, both with a maximum of $100$ iterations. The code is available as supplementary material.

4.1. Random potential in 2D

In this first numerical experiment we consider a random potential in two space dimensions with parameters

[TABLE]

We are interested in the first five eigenstates of lowest energy. The outcome of the pre-iteration with $\eta=0.2$ was already illustrated in Figure 3.4, leading to a small number of regions in which we expect the smallest eigenstates to localize. Applying the Ritz-Rayleigh method and the additional selection process, we end up in a ten-dimensional subspace of $V_{h}$, where all basis functions have local support. The five functions with lowest energy then serve as approximations of the eigenfunctions $u_{1},\dots,u_{5}$. The results are very convincing and depicted in Figure 4.1. For a comparison with the actual eigenstates we refer to Figure 3.2 (left).

As mentioned above, the pre-iteration only provides very rough approximations of the eigenstates and serves more the finding of an appropriate starting subspace for the Rayleigh-Ritz method. As a result, the assignments of the localization regions to the actual eigenstates may be inaccurate in the beginning. In the present example, four Rayleigh-Ritz steps were needed until the fifth eigenstate was correctly assigned. This effect can also be observed in the convergence plot in Figure 4.2. Therein, the relative error in $\lambda_{5}$ drops significantly from step $6$ to step $7$. Such effects are directly influenced by gaps within the spectrum of the Schrödinger operator. Here, the fifth and sixth eigenvalues are given by $2.6076\cdot 10^{4}$ and $2.6093\cdot 10^{4}$ and thus, at close quarters.

Figure 4.3 shows that the proposed methods also works if we are interested in a larger number of eigenvalues and eigenfunctions. If we only consider the pre-iteration, the relative error of the eigenvalues seems to be almost constant which is in agreement with the observations in [ADF*+*19]. In contrast to the landscape function approach we are able to control the accuracy of the method by the number of post-iterations. Further, the method is flexible enough to recognize when there are several eigenfunctions located in a single valley.

4.2. Speckle potential in 3D

The second example considers the linear Schrödinger eigenvalue problem (1.1) in three space dimensions with a speckle-like potential as it is used in actual experiments [FFI08]. As parameters we choose

[TABLE]

An illustration of the potential is given in Figure 4.4. As for the previous example we compare the outcome of Algorithm 1 with the results obtained by the eigs solver in Matlab. As a consequence, we are not able to exceed level $\ell_{h}=6$, as this marks the limit of the integrated eigenvalue solver.

The approximation of the first eigenfunction and the corresponding result obtained by eigs are compared in Figure 4.4. The displayed approximations are computed with $K_{\text{pre}}=3$ pre- and $K_{\text{post}}=5$ post-iteration steps. A comparison of computation times for various mesh levels of $\mathcal{T}^{\varepsilon}$ and $\mathcal{T}^{h}$ are summarized in the following table:

[TABLE]

4.3. Nonlinear Gross-Pitaevskii eigenvalue problem

In this final example we consider the nonlinear version of the Schrödinger equation, describing quantum-physical processes with interaction. This GPEVP has the form

[TABLE]

with $\delta\geq 0$ regulating the nonlinearity. For moderate values of $\delta$ similar localization results for the groundstates can be observed [APV18]. As a result, the here presented localization preserving preconditioner promises good approximation results. We emphasize that this nonlinear case is not easily treated by the landscape function approach. Considering again a finite element discretization as in Section 2.4, we obtain the nonlinear matrix eigenvalue problem

[TABLE]

where $N_{h}(u_{h})$ equals the finite-dimensional approximation of the nonlinearity. Inspired by the inverse iteration for the linear case, i.e., $u^{n+1}_{h}=A_{h}^{-1}M_{h}u^{n}_{h}$, one may consider the iteration

[TABLE]

with an additional normalization step. This then leads to the results shown in Figure 4.5, illustrating that the groundstate can be well approximated by local functions also in the nonlinear case. It goes without saying that the iteration scheme may be replaced by more advanced methods, e.g., based on gradient flows [BD04, HP18, AHP19].

5. Conclusion

In this paper, we have constructed a novel iterative scheme using analytical inside and adapting numerical analysis techniques used in localization proofs of the eigenstates [AHP18]. The novel method is solely based on local operations (relative to the oscillation/correlation length of the potential) and, hence, able to approximate eigenfunctions of the random Schrödinger operator up to high precision. For this, we exploit the fact that oscillatory, high-amplitude random potentials lead to a localization of the lowermost eigenstates, which allows an efficient computation.

The numerical experiments prove the applicability of the method also for the three-dimensional case as well as the corresponding nonlinear eigenvalue problem. Future research aims to further develop this approach in view of other applications with similar localization effects (such as light in a disordered medium [WBLR97]) as well as time-dependent problems.