Anderson localization is a famous wave phenomenon. Since Anderson’s seminal work [1], it was known that the low energy quantum states of a Schrödinger operator with a disordered lattice potential are strongly localized within a small subregion, provided that the degree of randomness in the lattice is sufficiently large.

Here we consider high-dimensional Anderson localization for the quantum states of the following stationary Schrödinger equation with a Dirichlet, Neumann, or Robin boundary condition (Anderson’s original work [1] considers a tight-binding model which can be viewed as the discretization of the continuous model studied here):where H = −△ + KV is the Hamiltonian with V = V(x) being a disordered potential, λ > 0 is an energy level (eigenvalue), u = u(x) is the associated quantum state (eigenmode), and K > 0 is the strength of disorder. The domain ω = (0, 1)^d is the d-dimensional unit hypercube with each side being divided uniformly into N intervals. In this way, the hypercube ω is divided into N^d smaller hypercubes of the same size. In each small hypercube ω_k, the potential V is a constant with its value being sampled from a given probability distribution, which is often chosen as the Bernoulli distributionor the uniform distributionThe values of V in these small hypercubes are independent of each other. The boundary condition of the model is rather general with g ≥ 0 being a given constant, h = h(x) ≥ 0 being a given function on ∂ω, and n = n(x) being the unit outward pointing normal vector field on ∂ω. If g = 0 and h = 1, then it reduces to the Dirichlet boundary condition; if g = 1 and h = 0, then it reduces to the Neumann boundary condition; if g = 1 and h ≠ 0, then it is the Robin boundary condition.

Note that most previous papers impose the Dirichlet boundary condition in the study of Anderson localization [16]. This is equivalent to trapping a particle by imposing an infinite potential outside the relevant domain. Here we also consider other boundary conditions due to the following reason. Recall that for the time-dependent Schrödinger equation i∂_tφ = −Hφ, the probability density of the particle is defined as $\rho =| \phi {| }^{2}=\phi \bar{\phi }$ and the probability current j = j(x) of the particle is defined aswhere $\bar{\phi }$ is the complex conjugate of φ. The probability density ρ and the probability current j are linked by the continuity equation ∂_tρ + ∇ · j = 0. Note that the Dirichlet, Neumann, and Robin boundary conditions can all guarantee thatwhich means that there is no net probability current across the boundary. Hence imposing the Neumann boundary condition can confine the system within the domain without an infinite potential well, as opposed to the Dirichlet boundary condition. This is the case of perfect reflection on the boundary (see section 5.2 of [20] for a detailed explanation). The Neumann boundary condition turns out to be very useful in the R-matrix theory of scattering [21, 22] and the theory of quantum graphs [23]. The Robin boundary condition is also considered here since it builds a bridge between the Dirichlet and Neumann boundary conditions.

In fact, the above model can be generalized to a more complicated model as follows:whereis an arbitrary second-order elliptic operator, V = V(x) is an arbitrary random potential, and ${\rm{\Omega }}\subset {{\mathbb{R}}}^{d}$ is an arbitrary bounded domain. Note that the operator L may be non-symmetric. Mathematically, the operator L is called symmetric if there exists a smooth function ρ = ρ(x) such thatwhere ${C}_{c}^{\infty }({\rm{\Omega }})$ is the space of all smooth functions on ω with compact supports andWe emphasize that here the symmetry of the operator L should be distinguished from the symmetry of the matrix $A={({a}^{{ij}})}_{d\times d}$ of diffusion coefficients. The latter is always symmetric, but the former may not. It is well-known that the operator L is symmetric if and only if the vector field A⁻¹(2b − ∇ · A) is conservative (has a potential function), i.e. there exists a smooth function U = U(x) such thatwhere A = (a^ij) is the diffusion matrix, b = (bⁱ) is the drift, and ∇ · A is a vector field whose ith component is given by ( ∇ ·A)ⁱ = ∂_ja^ij [24, 25]. Clearly, the Laplace operator △ is symmetric. The eigenvalues of a symmetric operator must be all real numbers, while the eigenvalues of a non-symmetric operator may be complex numbers. In the following, we shall define the localization landscape and its valley lines for the model (3).

A recent theory [16] has shown that under the Dirichlet boundary condition (g = 0 and h = 1), the spatial location of the localized quantum states of the eigenvalue problem (1) can be predicted by the solution of an associated Dirichlet problemwhere the solution w = w(x) is called the localization landscape. In fact, the theory in [16] was developed for symmetric elliptic operators and the Dirichlet boundary condition using the technique of the Green function. The landscape theory was further developed in [26] using a probabilistic approach and generalized in [19] to the Neumann boundary condition. Here we extend the concept of localization landscape to general non-symmetric elliptic operators and more complicated boundary conditions using a different probabilistic approach.

The key to our approach is to find the probabilistic representation for the solution of the eigenvalue problem (3). For the Dirichlet boundary condition, the solution can be represented by the expectation of a functional of a Brownian motion [26]. However, for Neumann and Robin boundary conditions, the solution can no longer be represented by a Brownian motion but should be represented by a reflecting diffusion. To see this, recall that for a bounded domain ${\rm{\Omega }}\subset {{\mathbb{R}}}^{d}$ with a unit outward pointing normal vector field n = n(x) on ∂ω, the operator L given in (4) is the infinitesimal generator of a reflecting diffusion $X={({X}_{t})}_{t\geqslant 0}$ with drift b = (bⁱ) and diffusion matrix a = (a^ij). This reflecting diffusion is the solution to the Skorokhod stochastic differential equationwhere $B={({B}_{t})}_{t\geqslant 0}$ is a d-dimensional standard Brownian motion and F_t is a continuous nondecreasing process that increases only when X_t ∈ ∂ω [27]. In particular, if L = △ is the Laplace operator, then the solution of (6) is a reflecting Brownian motion. It can be proved that when g > 0, the reflecting diffusion X_t can never exit the domain ω; once X_t touches ∂ω, it will be reflected into ω again due to the effect of the random force F_t [27]. With the aid of the reflecting diffusion, it can be shown that the solution of the problem (3) has the following probabilistic representation (see appendix A for the proof):where ${{\mathbb{E}}}_{x}$ denotes the conditional expectation given that X₀ = x. For problem (3), we define its localization landscape w = w(x) to be the solution of the following boundary value problem:Note that the boundary value problem (8) is obtained by setting the right-hand side of the eigenvalue problem (3) to be 1. Similarly, the landscape w also has a probabilistic representation using reflecting diffusion, which is given by (see appendix A for the proof)The probabilistic representations (7) and (9) are closely related. If we normalize the eigenmode u such that ∥u∥_∞ = 1, then we obtain the Filoche–Mayboroda inequalityThis inequality shows that when the energy level λ is not large, the eigenmode u must be small at those points where the landscape w is small. Hence the low energy eigenmodes must be localized to those subregions where the landscape is large. To see this, we illustrate the graphs of the landscape w and the first four eigenmodes u/λ for a one-dimensional problem when K is large (figure 1(a)). It can be seen that the inequality (10) indeed provides an accurate upper bound for the eigenmodes. Furthermore, the spatial locations of the confining subregions of these eigenmodes are perfectly predicted by the peak positions of the landscape. Recent numerical and theoretical studies [19, 28] have shown that the order of the eigenmodes is closely related to the height of the peaks of the landscape. This phenomenon was observed in our simulations in figure 1(a), where the first four eigenmodes are localized around the highest four peaks of the landscape.

In the one-dimensional case, the local minima of the landscape w divide the domain ω into many subintervals (figure 1(a)). Such subintervals give possible spatial locations where the eigenmodes can be localized. In the two-dimensional case, the confining subregions are determined by the valley lines of the landscape [16]. In fact, the valley lines are defined as the lines of steepest descent, starting from the saddle points of the landscape and going to its local minima. In this way, the valley lines separate the domain ω into many subregions. Along the valley lines, the values of w are expected to be small. From the inequality (10), the values of u are also small along the valley lines and thus the eigenmodes must be localized within the subregions enclosed by these lines. Following [16, 17], we use the watershed algorithm proposed in [29] to compute the valley lines numerically. In higher dimensions, valley lines will become hypersurfaces and the results are similar.

Figures 1(b), (c) illustrate the landscape, eigenmodes, and valley lines for a two-dimensional problem under the Neumann boundary condition. It is clear that the valley lines of the landscape separate the domain into many subregions and each eigenmode is exactly located in one of these subregions. Under the Dirichlet boundary condition, the eigenmodes must be located in the interior of the domain since they vanish on the boundary [16]. However, under the Neumann boundary condition, some eigenmodes may be localized on the boundary (see the first, second, and fourth ones in figure 1(c)). The boundary values of these eigenmodes are large and their local maxima may even appear on the boundary. This phenomenon will never appear for the Dirichlet boundary condition.

Here the eigenvalue problem (3) and the boundary value problem (8) are solved numerically by using the spectral element method instead of the classical finite difference method (FDM) or finite element method (FEM). In the spectral element method, the solution of the associated partial differential equation (PDE) is approximated by piecewise high-order polynomials and the Legendre polynomials are used as a basis to approximate the solution in each small hypercube ω_k. The degree of the Legendre polynomials is chosen to be 10 in the one-dimensional case and 6 in the two-dimensional case. In fact, classical FDM or FEM only have second-order convergence with respect to the mesh size, while the spectral element method can achieve exponential convergence with respect to the degree of the polynomials. Therefore, the spectral element method can obtain higher numerical accuracy with less computational costs than the other two methods.

Next, we focus on the limit behavior of the eigenmodes when the strength K of disorder is very large. For simplicity, we assume that the potential V restricted to each small hypercube ω_k is sampled from the Bernoulli distribution given in (2), which can only take the value of 0 or 1. To proceed, letdenote the set of points at which the potential vanishes. We first consider the behavior of the eigenmodes outside D. Let x ∈ ω be the initial position of the reflecting diffusion X_t which solves the Skorokhod equation (6). When x ∉ D, we have V(X_s) = 1 when s is small, which implies thatIt thus follows from (7) and (9) and the dominated convergence theorem thatThis shows that the landscape and eigenmodes must vanish outside D in the limit of K → ∞. In other words, the eigenmodes can only be localized in the region where the potential attains its minimum. This coincides with the intuition that the quantum states tend to be localized in the region with low potential energy.

We next focus on the behavior of the eigenmodes inside D. Clearly, D can be decomposed as the disjoint union of several connected components (subregions), i.e. D = D₁ ∪ D₂ ∪ ⋯ ∪ D_M. Since the potential V vanishes in D, in the limit of K → ∞, the eigenmode u must be the solution to the following local eigenvalue problem in each subregion D_k:Therefore, when K is very large, the spectrum of the Hamiltonian H = −L + KV is composed of the local eigenvalues of the operator −L in each subregion D_k. If an eigenvalue of H coincides with one of the local eigenvalues of −L in D_k, then the corresponding eigenmode will be localized in D_k; conversely, if an eigenvalue of H coincides with neither one of the local eigenvalues of −L in D_k, then the corresponding eigenmode will not be localized in D_k. Hence the limits in (11) enable us to decompose the eigenvalues problem (3) with random potential in the whole domain ω into some local eigenvalue problems (12) with zero potential in the subregions D_k. When K is not very large, the landscape and eigenmodes are small but not zero on ∂D_k\∂ω and thus the domain ω can be divided into the many weakly coupled subregions [16].

To gain a deeper insight, we depict the potential, landscape, and eigenmodes for a one-dimensional problem under different values of K (figure 2(a)). When K is large, the landscape and eigenmodes are only localized in the region where the potential vanishes. However, this is not the case when K is relatively small. Figure 2(b) illustrates the potential and valley lines for a two-dimensional problem under different values of K. When K is small, the confining subregions enclosed by the valley lines may include many connected components of D. However, for sufficiently large K, the connected components of D are exactly separated by valley lines.

We have seen that under the Neumann or Robin boundary condition, the low energy quantum states are very likely to be localized on the boundary (see the second eigenmode in figure 1(a) and the first, second, and fourth eigenmodes in figure 1(c)). Such a phenomenon will never take place for the Dirichlet boundary condition. When an eigenmode is localized on the boundary, its boundary value must be very large and the local maxima of the eigenmode may even appear on the boundary. In what follows, each eigenmode is always normalized such that ∥u(x)∥_∞ = 1. To better characterize the boundary effect, we introduce the following definition. In the one-dimensional case, we say that an eigenmode u is localized on the boundary if it satisfiesIn the two-dimensional case, similarly, we say that an eigenmode u is localized on the boundary if it satisfiesand we say that it is localized in the corner if it satisfiesSince the potential V is stochastic, the eigenmodes may or may not be localized on the boundary. For simplicity, we only consider the first eigenmode. The probability for the first eigenmode to be localized on the boundary is denoted by P_b and the probability for the first eigenmode to be localized in the corner is denoted by P_c. These two probabilities are collectively referred to as boundary probabilities. Clearly, the boundary probabilities are both 0 for the Dirichlet boundary condition. However, for the Neumann or Robin boundary condition, these probabilities are strictly positive with their values depending on the parameters h, K, and p, where the function h is chosen to be a constant for simplicity. Figures 3(a)–(c) illustrate the boundary probabilities as h, K, and p vary under the Robin boundary condition. Here the boundary probabilities P_b and P_c are computed numerically using the Monte Carlo method. Specifically, we generate the random potential 1000 times, solve the corresponding eigenvalue problems, and then count the frequency that the first eigenmode appears on the boundary or in the corner. The spectral element method enables us to solve the eigenvalue problem thousands of times in an acceptable time. It can be seen that the boundary probabilities decrease with the three parameters in both the one- and two-dimensional cases. When h is small, the first eigenmode is localized on the boundary with a relatively high probability. In particular, the probability P_b is exceptionally large in the two-dimensional case. When h ≫ 1, the Robin boundary condition reduces to the Dirichlet one and thus all boundary probabilities tend to zero (figure 3(a)). This explains why the boundary probabilities decrease as h increases.

We next focus on the dependence of the boundary probabilities on K. When K = h = 0, the problem (1) reduces to the Laplacian eigenvalue problem with the Neumann boundary condition and thus the first eigenmode is a constant. According to the definition, a constant eigenmode is localized on the boundary (since its maximal boundary value is larger than 0.5). As K increases, the eigenmode may be localized in the interior of the domain and thus the boundary probabilities decrease. When K ≫ 1, compared to the one-dimensional case, we find that the eigenmodes for a two-dimensional problem are more likely to be localized on the boundary (figure 3(b)). The effect of p on the boundary probabilities is much more complicated and will be discussed later.

Based on the simulations in figures 3(a)–(c), we have made two crucial observations: (i) when h is small, the eigenmodes are localized on the boundary with relatively high probability; (ii) when K is large, compared to the one-dimensional case, the eigenmodes for a two-dimensional problem are more likely to be localized on the boundary. Here we provide a detailed explanation for these two observations when h = 0 (Neumann boundary condition) and K ≫ 1 (strong disorder).

In fact, the boundary effect can be explained intuitively using the energy functional theory. To see this, recall that for the Dirichlet and Neumann boundary conditions, the eigenmodes of (1) are actually the critical points of the energy functional [20]where the first term is the kinetic energy and the second term is the potential energy. When K ≫ 1 is large, the potential energy is dominant whenever V ≠ 0 and thus the low energy quantum states must be localized in the region $D={\bigcup }_{k=1}^{M}{D}_{k}$ where the potential V vanishes. Hence we only need to focus on the kinetic energy on each subregion D_k, i.e.Hence the energy functional perspective also guides us to focus on the local eigenvalue problem in each subregion D_k. Intuitively, for the Dirichlet boundary condition, the function u must tend to zero on the boundary, which leads to a large ∣∇u∣ and thus a large kinetic energy J_k(u) when D_k appears on the boundary. For the Neumann boundary condition, however, the function u does not have to be zero on the boundary, which makes the kinetic energy J_k(u) much smaller. This explains why the boundary effect is significant for the Neumann boundary condition.

We next examine the boundary effect in more detail. For simplicity, we first focus on the one-dimensional case. For large K, we have shown that the eigenmode of the problem (1) restricted to each subregion D_k is approximately the solution to the local problem (12), which is a Laplacian eigenvalue problem. The first eigenmode must be localized in the subregion D_k where the local eigenvalue is the smallest. Suppose that D_k = [x_k, x_k + L], where L is the length of D_k. When ${D}_{k}\cap \partial {\rm{\Omega }}=\varnothing $, i.e. x_k > 0 and x_k + L < 1, the local problem (12) can be rewritten aswhich has the Dirichlet boundary condition. The first eigenvalue and eigenmode of the system are given by λ₁ = (π/L)² and ${u}_{1}(x)=\sin (\pi (x-{x}_{k})/L)$, respectively. Clearly, the longer the subregion D_k, the smaller the local eigenvalue. If the original system has the Dirichlet boundary condition, then the first eigenmode must be localized in the longest subregion when K ≫ 1.

When D_k ∩ ∂ω ≠ φ, i.e. x_k = 0 or x_k + L = 1, the local problem (12) has a mixed boundary condition. For simplicity, we only consider the case of x_k = 0. In this case, the local problem can be rewritten aswhere the left-hand side of the interval has the Neumann boundary condition and the right-hand side has the Dirichlet one. The first eigenvalue and eigenmode of the system are given by λ₁ = (π/2L)² and ${u}_{1}(x)=\cos (\pi x/2L)$, respectively. We then make a crucial observation that the first eigenpair of this system can be regarded as that of another system in the extended subregion ${D}_{k}^{(e)}=[-L,L]$:which has the Dirichlet boundary condition (see figure 3(d) for an illustration of the extended subregion and the first eigenmodes in the original and extended subregions). In other words, when D_k is around the boundary, the first eigenmode of the local system (19) with mixed boundary condition can be viewed as that of an equivalent system (20) with Dirichlet boundary condition in the extended subregion ${D}_{k}^{(e)}$ which is twice as long as the original subregion D_k. Clearly, for a subregion of a fixed length, the local problem (19) has a smaller eigenvalue if the subregion appears on the boundary since the length of the subregion should be ‘doubled’. From the energy functional perspective, compared to the Dirichlet boundary condition, the kinetic energy J_k(u) for the Neumann boundary condition is smaller by a factor of 2 when D_k appears on the boundary (note that the factor may not be 2 in higher dimensions). This explains why compared to the interior of the domain, the eigenmodes are more likely to be localized on the boundary when h is small and K is large.

The above considerations can also be generalized to higher dimensions. For each subregion D_k, we can extend it to another subregion ${D}_{k}^{(e)}$ via mirror symmetry along the boundary (figure 3(e)). If D_k does not appear on the boundary, then the extended subregion ${D}_{k}^{(e)}$ is the same as D_k; if D_k appears on the boundary, then ${D}_{k}^{(e)}$ is strictly larger than D_k. Suppose that D_k appears on the boundary. Then the local problem (12) has the Neumann boundary condition on ∂D_k ∩ ∂ω and the Dirichlet boundary condition on ∂D_k\∂ω. In appendix B, we have proved that the first eigenmode of the local system with mixed boundary conditions is exactly the same as that of an equivalent system with Dirichlet boundary conditions in the extended subregion ${D}_{k}^{(e)}$. In particular, in the two-dimensional case, if D_k does not lie in the corner, then the extended subregion is twice as large as the original subregion. However, if D_k lies in the corner since there are two edges enclosing D_k, the extended subregion should be reflected along both edges and thus is four times as large as the original subregion. Recall that a Laplacian eigenvalue problem tends to have smaller eigenvalues on a larger domain with strong symmetry. Clearly, if D_k appears on the boundary, especially in the corner, then it has a larger extended subregion which also has strong symmetry (figure 3(e)). Therefore, subregions around the boundary tend to have smaller local eigenvalues and thus the eigenmodes are more likely to be localized on the boundary.

With the aid of the concept of extended subregions, we are able to transform a local eigenvalue problem with a mixed boundary condition into an equivalent eigenvalue problem with the Dirichlet boundary condition in an extended subregion. When h is small and K is large, the first eigenmode will be localized in the extended subregion with the smallest Dirichlet eigenvalue. In particular, in the one-dimensional case, the first eigenmode will be localized in the longest extended subregion. If the longest two extended subregions are the same in length, then the first eigenmode may be localized in these two subregions simultaneously. This phenomenon will be discussed in detail later.

Note that in the one-dimensional case, there are at most two subregions around the boundary and the extended subregion is twice as long as the original one. However, in the two-dimensional case, the proportions of subregions that need to be extended are much higher (figure 3(e)). Moreover, if a subregion appears in the corner, then its area should be ‘magnified’ four times. These two facts explain why compared to the one-dimensional case, the eigenmodes for a two-dimensional problem have a much higher probability to be localized on the boundary.

In the previous discussion, we have explained why the boundary probabilities decrease as h and K increase. Here we focus on the dependence of the boundary probabilities on p. For simplicity, we only focus on the one-dimensional case for small h and large K.

In the one-dimensional case, the domain ω = (0, 1) is divided uniformly into N subintervals of the same length 1/N. In each subinterval, the potential V can only take the value of 0 or 1 with p being the probability of V = 0 and q = 1 − p being the probability of V = 1. Then the domain ω can be decomposed into many subregions with V = 0 and V = 1 alternatively (see the white and gray regions in figure 2(a)). Intuitively, when N is large, the lengths of these subregions are approximately independent and the number of subintervals included in each subregion with V = 0 (V = 1) approximately follows a geometric distribution with parameter p (q). For small h and large K, we have shown that the first eigenmode must be localized in the longest extended subregion with V = 0. Therefore, the probability that the first eigenmode is localized on the boundary is approximately equal to the probability that the longest extended subregion with V = 0 appears on the boundary, which is given by (see appendix C for the proof)This formula reveals a complex relationship between the boundary probability P_b on the parameter p in the one-dimensional case.

To test our theory, we compare the theoretical prediction given in (21) with numerical simulations where the boundary probability P_b is computed by generating the random potential 1000 times and then solving the associated eigenvalue problem (1) (figure 3(f)). Interestingly, while the analytical expression is very complicated, it is in perfect agreement with simulation results. Moreover, we find that P_b decreases with p. When p ≈ 0.5, we have P_b ≈ 0.26, which means that the first eigenmode has a one in four chance of being localized on the boundary. In the two-dimensional case, the dependence of the boundary probabilities P_b or P_c on the parameter p is similar (figure 3(c)), but it is very difficult to obtain their analytical expressions.

Thus far, the boundary effects were examined when the lattice potential V within each hypercube is randomly sampled from a Bernoulli distribution. From figure 1, it is clear that the boundary effect is also significant when the random potential V is sampled from a uniform distribution. A natural question is whether the boundary effect is also significant for other types of random potentials.

To answer this, we compare the boundary effects for four types of random potentials sampled from the Bernoulli distribution, the normal distribution, the gamma distribution, and the uniform distribution, respectively, as shown in figure 4. Here the strength of disorder is chosen to be K = 10⁴. For each distribution type, we illustrate the boundary probabilities P_b and P_c in the one- and two-dimensional cases for random potentials with the same mean μ = 1/2 and three different standard deviations, where the standard deviations σ within each hypercube are chosen to be σ = μ, $\sigma =\mu /\sqrt{3}$, and σ = μ/3, respectively (see the four columns of figure 4). It is clear that the boundary effect is significant for all types of random potentials. Interestingly, we find that the boundary effect is insensitive to the standard deviation of the random potential within each hypercube. For the Bernoulli and gamma distributions, the boundary probabilities P_b and P_c become slightly higher as the standard deviation decreases, i.e. as the distribution becomes more concentrated. However, the normal and uniform distributions give rise to the opposite effect, i.e. the boundary probabilities become slightly lower with the decrease of the standard deviation. Furthermore, we find that whether the low energy quantum states are localized on the boundary is remarkably affected by the distribution type. The random potentials sampled from the normal distribution lead to a much weaker boundary effect than those sampled from the other three types of distributions.

From the simulations in figure 2(a), we can see that for different values of K, the eigenmodes may be localized in completely different subregions. When K ≫ 1, the first eigenmode is localized in the extended subregion with the smallest local eigenvalue. In particular, in the one-dimensional case, the first eigenmode is localized in the longest extended subregion with V = 0. When K is relatively small, however, the first eigenmode may be localized in somewhere else. This suggests that as K increases, the localization subregion for a given eigenmode may change and a bifurcation phenomenon may take place. Note that in figure 2(a), while the first eigenmode is not localized in the longest extended subregion with V = 0 when K is small, it is localized in the union of some long subregions with V = 0 that are separated by some short barriers with V = 1. The union of these subregions is even longer than the longest extended subregion. In other words, when K is large, the localization behavior of the system is very sensitive to the short potential barriers, while it is insensitive to these barriers when K is small. We emphasize that while figure 2 only shows a bifurcation for the first eigenmode, similar phenomena can be also observed for other eigenmodes.

To gain a deeper insight into the bifurcation phenomenon, we consider a one-dimensional toy model with periodic boundary conditions (figure 6(a)). In this model, the potential V takes the values of 0 and 1 alternately in the intervals with lengths L₂/2, L₁, L₂, L₃, L₄, L₃, and L₂/2 respectively, i.e.

The periodic boundary condition guarantees that the system has stronger symmetry and thus the theory will become much simpler. The potential V is composed of two ‘wells’ with lengths L₁ and 2L₃ + L₄ respectively, while there is a short barrier with length L₄ in the middle of the right well. The distance between the two wells is L₂. For convenience, the left and right wells are denoted by W₁ = [x₁, x₂] and W₂ = [x₃, x₅], respectively.

We next investigate the localization of the quantum states of the stationary Schrödinger equationHere we require that the left well W₁ is the longest subinterval with V = 0 and the right well W₂ is the union of two relatively long subintervals with V = 0 that are separated by a very short subinterval with V = 1. The total length of the two relatively long subintervals should be greater than the length of the longest subinterval so that bifurcation can occur. To ensure that the two wells do not interfere with each other, we further require the distance between the two wells to be large enough. To be more specific, the parameters L₁, L₂, L₃, and L₄ should satisfy: (i) L₁ > L₃, which grantees that W₁ is the longest subinterval; (ii) L₁ < 2L₃, which grantees that W₂ is longer than W₁; (iii) L₄ < L₃/2, which guarantees that the barrier in W₂ is short enough; (iv) L₂ > L₁ + 2L₃ + L₄, which grantees that the two wells are far enough; (v) L₁ + 2L₂ + 2L₃ + L₄ = 1, which guarantees that the length of the whole interval is 1.

For the toy model, we illustrate the first eigenmode of the system S given in (24) under different values of K (figure 6(b)). Since the right well is longer than the left one and the barrier in the right well is short, when K is small, the first eigenmode is mainly localized in the right well. Since the left well is the longest subinterval with V = 0, when K is large, the first eigenmode is mainly localized in the left well. With the increase of K, the left peak of the eigenmode becomes higher and the right peak becomes lower. At a critical threshold of K, the two peaks are equal in height and bifurcation takes place—the confining subregion of the eigenmode transitions from the right well to the left one. For the first eigenmode u, we introducewhich represents the height of the left peak relative to the total height of the two peaks. In general, the relative height F of the left peak for the first eigenmode increases with K (figure 6(c)). The critical threshold of K is defined as the value at which F = 0.5. At the critical value K = K_c, the relative height F of the left peak changes sharply from a low to a high value, corresponding to the occurrence of bifurcation.

When L₂ is large, the two wells are far apart and they have little influence on each other. Hence we can decompose the system S into two subsystems S₁ and S₂ with S₁ corresponding to the left peak of the eigenmode and with S₂ corresponding to the right peak. Specifically, we introduce two new potentials V₁ and V₂, where V₁ is obtained from V by setting the values in the right well to 1 and V₂ is obtained from V by setting the values in the left well to 1 (figure 6(a)), i.e.

Then the original system S can be decomposed into two subsystems S₁ and S₂, which are defined as

Here the subsystem S₁ has potential V₁ and the subsystem S₂ has potential V₂. Intuitively, for the first two eigenmodes of the system S, one is localized in the left well W₁ and the other is localized in the right well W₂. The eigenmode localized in W₁ (W₂) can be approximated by the first eigenmode of subsystem S₁ (S₂). Figure 6(d) illustrates the first two eigenmodes of the original system and the first eigenmodes of the two subsystems. It can be seen that the first eigenmode of each subsystem almost coincides with one of the first two eigenmodes of the original system. Moreover, we illustrate the first two eigenvalues of the original system and the first eigenvalues of the two subsystems under different values of K in figure 6(e). Clearly, the first eigenvalues of the two subsystems serve as good approximations to the first two eigenvalues of the original system.

We now provide a detailed explanation for the occurrence of bifurcation. From figure 6(e), the first eigenvalue of each subsystem increases with K. However, the eigenvalues of both subsystems increase with K at different rates (see two lines of different slopes in figure 6(e)). When K is small, the eigenvalue of subsystem S₂ is smaller than that of subsystem S₁, and thus the first eigenmode of the original system S is localized in the right well W₂. When K is large, the eigenvalue of S₂ exceeds that of S₁, and thus the first eigenmode of S is localized in the left well W₁. At the critical threshold K = K_c, the localization subregion of the first eigenmode transitions from W₂ to W₁. In analogy to multimodality, bifurcation occurs when the first eigenvalues of both subsystems are approximately equal, i.e. when the first and the second eigenvalues of the original system are approximately equal. This critical eigenvalue is denoted by λ = λ_c.

Clearly, the critical threshold of K plays a crucial role in the bifurcation phenomenon. Numerically, in order to compute the critical value, we need to solve the original system under many different values of K, which is very time-consuming. Here we provide an analytical theory of the critical threshold and the associated critical eigenvalue.

To this end, we first simplify the subsystems S₁ and S₂. For each subsystem S_i, due to the periodic boundary condition, we can shift the potential V_i along the interval ω = (0, 1) without changing the eigenvalues. In particular, we move the well W_i to the middle of the interval (figure 6(a)) and the shifted potentials for the two subsystems are defined as

The subsystem with the shifted potential ${\hat{V}}_{i}$ is denoted by ${\hat{S}}_{i}$, which can be written explicitly as

Note that the shifted potential ${\hat{V}}_{i}$ satisfies ${\hat{V}}_{i}(x)={\hat{V}}_{i}(1-x)$. Thus the eigenmodes of the subsystem ${\hat{S}}_{i}$ satisfy u(x) = u(1 − x) and $u^{\prime} (x)=-u^{\prime} (1-x)$. Taking x = 1/2, we obtain $u^{\prime} (1/2)=0$. Moreover, taking x = 0 and using the periodic boundary condition $u^{\prime} (0)=u^{\prime} (1)$, we obtain $u^{\prime} (0)=0$. Then the eigenvalues of the subsystem ${\hat{S}}_{i}$ coincide with the eigenvalues of the following equivalent subsystem in the interval (0, 1/2) with the Neumann boundary condition:

where the potential ${\tilde{V}}_{i}$ of the equivalent subsystem is defined as (see figure 6(a) for an illustration)

After the above simplification, the eigenvalues of each subsystems S_i are consistent with those of the equivalent subsystem ${\tilde{S}}_{i}$. For a fixed value of K, each eigenvalue λ of the subsystem ${\tilde{S}}_{1}$ satisfies the equation (see appendix D for the proof)

and each eigenvalue λ of the subsystem ${\tilde{S}}_{2}$ satisfies the equation

where $\alpha =\sqrt{\lambda }$, $\beta =\sqrt{K-\lambda }$, t₀ = L₁/2, t₁ = L₄/2, t₂ = L₄/2 + L₃, and t₃ = 1/2. At the critical threshold K = K_c, the first eigenvalues of the two subsystems are approximately equal. Then the critical threshold K_c and the critical eigenvalue λ_c approximately satisfy the following system of algebraic equations:This system of equations can be used to analytically predict the critical threshold K_c and the critical eigenvalue λ of bifurcation.

To test our theory, we sample the parameters L₁, L₂, L₃, and L₄ randomly from a large swathe of the parameter space. Specifically, the parameters are chosen as L₃ ∼ U[0.03, 0.055], L₁ ∼ U[1.3L₃, 1.7L₃], L₄ ∼ U[0.2L₃, 0.4L₃], and L₂ is determined so that L₁ + 2L₂ + 2L₃ + L₄ = 1, where U[a, b] denotes the uniform distribution over the interval [a, b]. We then use (28) to predict K_c and λ_c. For 1000 sample points, the mean relative prediction errors ((predicted value − real value)/real value) for K_c and λ_c are only 1.52 × 10⁻⁴ and 1.61 × 10⁻⁴, respectively. This shows that our analytical theory can indeed capture the critical values of bifurcation accurately without carrying out numerical simulations.

We have seen that the critical threshold K_c and critical eigenvalue λ_c satisfy the system of algebraic equations (28). Clearly, the value of K_c depends on the parameters L₁, L₂, L₃, and L₄. However, the relationship between them is very complicated. Here we investigate their relationship numerically and reveal the key factors of bifurcation.

Since L₁ + 2L₂ + 2L₃ + L₄ = 1, there four parameters only have three degrees of freedom. For simplicity, we definewhere P₁ represents the total length of the two wells relative to the length of the whole interval, P₂ represents the length of the left well relative to the total length of the two wells, and P₃ represents the length of the barrier relative to the length of the right well. The relationship between K_c and P₁, P₂, and P₃ is complicated. However, we can obtain an approximate relationship between them using our analytical theory. Figures 7(a)–(c) illustrate the critical threshold K_c as P₁, P₂, and P₃ vary. Here we generate different values of P₁, P₂, or P₃ randomly from a large swathe of the parameter space and keep the other two parameters invariant. Interestingly, we find that K_c depends on P₁ in a power law form asTo check this, we plot $\mathrm{log}({K}_{c})$ as a function of $\mathrm{log}({P}_{1})$ in figure 7(a) and find that there is an excellent linear relationship between them with an exceptionally high R² and a slope of −2. Similarly, we find that K_c depends on P₂ in an exponential form as (figure 7(b))and K_c depends on P₃ in a power law form as (figure 7(c))These results indicate that K_c decreases with P₁ and P₃ at a power law speed and decreases with P₂ at an exponential speed. In particular, short potential wells (small P₁), a short left well relative to the right well (small P₂), and a short barrier in the right well (small P₃) are more capable of producing a large critical value of K. The exponential decay of K_c with P₂ suggests that the relative length of the left well affects the critical threshold more severely than the other two factors.

Under the Neumann boundary condition, the first eigenmode is localized in the longest extended subregion with V = 0 when K is large. Therefore, the event that the first eigenmode is localized on the boundary is equivalent to the event that the longest extended subregion with V = 0 appears on the boundary. This event can be classified into the following three situations.

First, we consider the situation where V(0) = V(1) = 0, whose probability is q². In this case, the lengths of the extended subregions with V = 0 are given by 2X₁, X₂, ⋯ , X_M−1, 2X_M. Then the probability that the longest extended subregion appears on the boundary can be computed explicitly as

Second, we consider the situation where V(0) = 0 and V(1) = 1, whose probability is pq. In this case, the lengths of the extended subregions with V = 0 are given by 2X₁, X₂, ⋯ , X_M. Then the probability that the longest extended subregion appears on the boundary can be computed explicitly as

Third, we consider the situation where V(0) = 1 and V(1) = 0, whose probability is pq. In this case, the lengths of the extended subregions with V = 0 are given by X₁, X₂, ⋯ , 2X_M. In analogy to the second case, the probability that the longest extended subregion appears on the boundary is given by

Finally, using the total probability formula, the probability that the first eigenmode is localized on the boundary is given by

which gives (21) in the main text.

Under the Dirichlet boundary condition, the first eigenmode is localized in the longest subregion with V = 0 when K is large. If there is a unique longest subregion with V = 0, then the first eigenmode is unimodal; if there are two or more longest subregions with V = 0, then the first eigenmode is multimodal. Therefore, the event that the first eigenmode is unimodal is equivalent to the event that there is a unique longest subregion with V = 0. This probability can be computed explicitly as

which gives (22) in the main text.

Under the Neumann boundary condition, the first eigenmode is localized in the longest extended subregion with V = 0 when K is large. Therefore, the event that the first eigenmode is unimodal is equivalent to the event that there is a unique longest extended subregion with V = 0. This event can be classified into the following four situations.

First, we consider the situation where V(0) = V(1) = 0, whose probability is q². In this case, the lengths of extended subregions with V = 0 are given by 2X₁, X₂, ⋯ , X_M−1, 2X_M. Then the probability that there is a unique longest extended subregion is given by

where [x] represents the largest integer less than or equal to x. Second, we consider the situation where V(0) = 0 and V(1) = 1, whose probability is pq. In this case, the length of extended subregions with V = 0 are given by 2X₁, X₂, ⋯ , X_M. Then the probability that there is a unique longest extended subregion is

Third, we consider the situation where V(0) = 1 and V(1) = 0, whose probability is pq. In this case, the probability that there is a unique longest extended subregion is given by p₃ = p₂. Finally, we consider the situation in which V(0) = V(1) = 1, whose probability is p². In this case, the length of the extended subregions with V = 0 are given by X₁, X₂, ⋯ , X_M. Then the probability that there is a unique longest extended subregion is given by p₄ = 1 − p_D, where p_D is the multimodal probability for the Dirichlet boundary condition. In summary, the probability that the first eigenmode is unimodal is given by

which gives (23) in the main text.

Here we compute the eigenvalues of the subsystem ${\tilde{S}}_{1}$ given in (27). For convenience, we define

In the interval [0, t₀], the solution of the subsystem ${\tilde{S}}_{1}$ is given by

and in the interval [t₀, 1/2], the solution of (27) is given by

where A, B, C, D are four constants to be determined. Under the Neumann boundary condition, due to the continuity of the eigenmodes, we have

which can be rewritten in matrix form as

Therefore, λ is an eigenvalue of the subsystem ${\tilde{S}}_{1}$ if and only if the above coefficient matrix M is nondegenerate, i.e.

which can be rewritten as

Here we compute the eigenvalues of the subsystem ${\tilde{S}}_{2}$ given in (27). For convenience, we define

In the interval [0, t₁], the solution of the subsystem ${\tilde{S}}_{2}$ is given by

in the interval [t₁, t₂], the solution of the subsystem ${\tilde{S}}_{2}$ is given by

and in the interval [t₂, t₃], the solution of the subsystem ${\tilde{S}}_{2}$ is given by

where A, B, C, D, E, F are six constants to be determined. Under the Neumann boundary condition, due to the continuity of the eigenmodes, we have

which can be rewritten in matrix form as

Therefore, λ is an eigenvalue of the subsystem ${\tilde{S}}_{1}$ if and only if the above coefficient matrix M is nondegenerate, i.e.

which can be rewritten as