The effective potential of an M-matrix

The fundamental premises of quantum physics guarantee that a potential V induces exponential decay of the eigenfunctions of the Schrödinger operator −Δ + V (on either a continuous domain $Rd$ or a discrete lattice $Zd$⁠) as long as V is larger than the eigenvalue E outside of some compact region. This heuristic principle has been established with mathematical rigor by Agmon¹ and has served as a foundation to many beautiful results in semiclassical analysis and other fields (see, e.g., Refs. 2–5 for a glimpse of some of them). Roughly speaking, the modern interpretation of this principle is that the eigenfunctions decay exponentially away from the “wells” {x : V(x) ≤ E}.

In 2012, two of the authors of the present paper introduced the concept of the localization landscape. They observed in Ref. 6 that the solution u to the equation (−Δ + V)u = 1 appears to have an almost magical power to “correctly” predict the regions of localization for disordered potentials V and to describe a precise picture of their exponential decay. For instance, if V takes the values 0 and 1 randomly on a two-dimensional lattice $Z2$ (a classical setting of the Anderson–Bernoulli localization), the eigenfunctions at the bottom of the spectrum are exponentially localized, that is, exponentially decaying away from some small region, but this would not be detected by the Agmon theory because the region {V ≤ E} could be completely percolating and there is no “room” for the Agmon-type decay, especially if the probability of V = 0 is larger than the probability of V = 1. Indeed, the phenomenon of Anderson localization is governed by completely different principles, relying on the interferential rather than confining impact of V. On the other hand, looking at the landscape in this example, we observe that the region ${1u≤E}$ exhibits isolated wells and that the eigenmodes decay exponentially away from these wells. It turns out that, indeed, the reciprocal of the landscape, $1u$⁠, plays the role of an effective potential, and in Ref. 7, Arnold, David, Jerison, and the first two authors proved that the eigenfunctions of −Δ + V decay exponentially in the regions where ${1u>E}$ with the rate controlled by the so-called Agmon distance associated with the landscape, a geodesic distance in the manifold determined by $(1u−E)+$⁠. The numerical experiments in Ref. 8 and physical considerations in Ref. 9 show an astonishing precision of the emerging estimates, although mathematically speaking in order to use these results for factual disordered potentials, one has to face, yet again, a highly non-trivial question of resonances—see the discussion in Ref. 7. At this point, we have only successfully treated Anderson potentials via the localization landscape in the context of a slightly different question about the integrated density of states.¹⁰ 

However, the scope of the landscape theory is not restricted to the setting of disordered potentials. In fact, all results connecting the eigenfunctions to the landscape are purely deterministic, and one of the key benefits of this approach is the absence of a priori assumptions on the potential V, which already in Ref. 7 allowed us to rigorously treat any operator −div A∇ + V with an elliptic matrix of bounded measurable coefficients A and any non-negative bounded potential V, a level of the generality not accessible within the classical Agmon theory. These ideas and results have been extended to quantum graphs,¹¹ to the tight-binding model,¹² and, perhaps most notably, to many-body localization in Ref. 13.

This paper shows that the applicability of the landscape theory, in fact, extends well beyond the scope of the Schrödinger operator or, for that matter, even the scope of Partial Differential Equations (PDEs), at least in the bottom of the spectrum where the region ${1u≤E}$ exhibits isolated potential wells. Indeed, let us now consider a general real symmetric positive definite N × N matrix $A=(aij)i,j∈[N]$⁠, which one can view as a self-adjoint operator on the Hilbert space ℓ²([N]) on the domain [N] ≔ {1, …, N}. In certain situations, one expects A to exhibit “localization” in the following two related aspects, which we describe informally as follows:

• (Eigenvector localization) Each eigenvector¹⁴ $ϕ=(ϕk)k∈[N]$ of A is localized to some index i of [N] so that |ϕ_k| decays when |k − i| exceeds some localization length L ≪ N.

• (Poisson statistics) The local statistics of eigenvalues λ₁, …, λ_N of A asymptotically converge to a Poisson point process in a suitably rescaled limit as N → ∞.

Empirically, phenomena (i) and (ii) are observed to occur in the same matrix ensembles A; intuitively, the eigenvector localization property (i) implies that A “morally behaves like” a block-diagonal matrix, with the different blocks of A supplying “independent” sets of eigenvalues, thus leading to the Poisson statistics in (ii). However, the two properties (i) and (ii) are not formally equivalent; for instance, conjugating A by a generic unitary matrix will most likely destroy property (i) without affecting property (ii).

We now focus on the question of establishing eigenvector localization (i). Can one deduce any uniform bound on the eigenvectors of a general matrix A depicting, in particular, a structure of the exponential decay similarly to the aforementioned considerations for a matrix of the Schrödinger operator −Δ + V? An immediate objection is that there is no “potential” that could play the role of V. Even aside from the fact that the proof of the Agmon decay relies on the presence of both kinetic and potential energy, as well as on many PDE arguments, it is not clear whether there is a meaningful function, analogous to V, which governs the behavior of eigenvectors of a general matrix. The main result of this paper is that, surprisingly, the landscape theory still works, at least in the class of real symmetric Z-matrices (matrices with non-positive entries off the diagonal). Furthermore, when A is a real symmetric non-singular M-matrix (a positive semi-definite Z-matrix), the reciprocal $1u$ of the solution to Au = 1 gives rise to a distance function ρ on the index set [N], which predicts the exponential decay of the eigenvectors.

We introduce an Agmon-type distance ρ on the index set [N] ≔ {1, …, N} associated with an N × N matrix A, a N × 1 “landscape” vector u, and an additional spectral parameter $E∈R$⁠.

It is easy to see that ρ is a pseudo-metric in the sense that it is symmetric and obeys the triangle inequality with ρ(i, i) = 0, although without further hypotheses²² on A, u, E, it is possible that ρ(i, j) could be zero or infinite for some i ≠ j. One can view ρ as a weighted graph metric on the graph with vertices [N] and edges given by those pairs (i, j) with a_ij ≠ 0 and with weights given by the right-hand side of (2.3). We discuss the comparison between ρ and the Euclidean metric in the beginning of Sec. V.

We recall that a Z-matrix is any N × N matrix A such that a_ij ≤ 0 when i ≠ j, and a M-matrix is a Z-matrix with all eigenvalues having non-negative real part. Our typical setup is the case when A is a real symmetric non-singular M-matrix, i.e., a positive definite matrix with non-positive off-diagonal entries, and in that case, we will choose u as the landscape function, i.e., the solution to Au = 1, with 1 denoting a vector with all values equal to 1. We say that a matrix A has connectivity at most W_c if every row and column has at most W_c non-zero non-diagonal entries. If A is a real symmetric non-singular M-matrix, all the principal minors are positive (see, e.g., Ref. 23), and hence, by Cramér’s rule, all the coefficients of the landscape u will be non-negative. In this case, a simple form of our main results is as follows:

Informally, the above inequality ensures that an eigenvector φ experiences exponential decay away from the wells of the effective potential $V̄=(1uk)k∈[N]$ cutoff by the energy level E. This is what typically happens for the Schrödinger operator −Δ + V (according to some version of the Agmon theory); see, for instance, Ref. 7 (Corollary 4.5). However, the existence of such an effective potential for an arbitrary M-matrix is perhaps surprising.

In fact, our results apply to the larger class of real symmetric Z-matrices A and more general vectors u and can handle “local” eigenvectors as well as “global” ones. We first introduce some more notations.

To avoid confusion, we shall sometimes refer to the original notion of an eigenvector as a global eigenvector; this is the special case of a local eigenvector in which M = [N].

We can now state a more general form of Theorem 2.5.

There are two terms on the left-hand side of (2.4), corresponding to two different lines in the display, and they serve different purposes. The bound for the term in the second line [which, in particular, yields (2.5)] asserts, roughly speaking, that the eigenvector φ_k experiences exponential decay in the regime where k is far from $KĒ$ in the sense that $ρ(k,KĒ)≫Wc$⁠. Note that Theorem 2.5 is the special case of (2.5) when A is a M-matrix and u = A⁻¹1.

By taking advantage of the term in the first line of (2.4), we can proceed further and demonstrate an approximate diagonalization, or decoupling, of A on the collection of disjoint subregions defined by the landscape function u by following the arguments from Ref. 7. The details are too technical to be put in the Introduction, and we refer the reader to Sec. V. In short, the idea is that viewing [N] as a graph induced by A (with the vertices connected whenever a_ij ≠ 0), we can define a Voronoi-type splitting of this graph into subgraphs, Ω_ℓ, each containing an individual connected component of $KĒ$ (or sometimes merging a few components if convenient). Then, A can be essentially decoupled into smaller matrices $A|Ωℓ$ with the strength of coupling exponentially small in the ρ_A,u distance between individual “wells.” Related to this, the spectrum of A will be exponentially close to the combined spectrum of $A|Ωℓ$’s.

Note how the geometry of the metric ρ is sensitive to the spatial distribution of the matrix A and, in particular, to the connectivity properties of the graph induced by the locations of the nonzero locations of A. For instance, conjugating A by a generic orthogonal matrix will almost certainly destroy the localization of the eigenvectors φ but will also heavily scramble the metric ρ (and most likely also destroy the property of being an M- or Z-matrix). On the other hand, conjugating A by a permutation matrix will simply amount to a relabeling of the (pseudo-)metric space ([N], ρ) and not affect the conclusions of Corollary 2.5 and the decoupling results in Theorem 5.2 and Corollary 5.5 in any essential way.

We will show some results of the numerical simulations in Sec. III and then pass to the proofs, but let us say a few more words about the particular cases that would perhaps be of most interest.

Here, the connectivity is W_c = 2W. Strictly speaking, the random Gaussian band matrix models A considered in Example 1.1 do not fall under the scope of Corollary 2.5 because the matrices will not be expected to have non-positive entries away from the diagonal nor will they be expected to be positive definite. However, one can modify the model to achieve these properties (at least with high probability) by replacing the Gaussian distributions by distributions supported on the negative real axis and then shifting by a suitable positive multiple of the identity to ensure positive definiteness with high probability. These changes will likely alter the semicircle law for the bulk distribution of eigenvalues, but in the spirit of the universality phenomenon, one may still hope to see localization of eigenvectors, say, in the bulk of the spectrum, as long as the width W of the band matrix is small enough (in particular, when W ≪ N^1/2). In this case, Corollary 2.5 entails exponential decay of the eigenvectors governed by the landscape $1u$⁠, and Theorem 5.2 and Corollary 5.5 yield the corresponding diagonalization of A. Of course, the key question is the behavior of the landscape. If the set $KĒ$ of wells is localized to a short interval, then this corollary will establish localization in the spirit of (i) above; however, if $KĒ$ is instead the union of several widely separated intervals, then an eigenvector could, in principle, experience a resonance in which non-trivial portions of its ℓ² energy were distributed amongst two or more of these intervals. Whether or not this happens is governed to some extent by Theorem 5.2 and Corollary 5.5. These results indicate that the resonances have to be exponentially strong in the distance between the wells, and our numerical experiments suggest that such strong resonances are, in fact, quite rare.

When A is a matrix of the tight-binding Schrödinger operator (a standard discrete Laplacian plus a potential) in a cube in $Zd$⁠, the connectivity parameter W_c is now the number of nearest neighbors, 2d, and the size of the matrix is the sidelength of the cube to the power d. If the potential is non-negative, A is an M-matrix with the entries a_ij equal to −1 whenever i ≠ j corresponds to the nearest neighbors in the graph structure induced by $Zd$⁠, and a_ii = 2d + V_i. This particular case has been considered in Ref. 12, and our results clearly cover it. However, the tight-binding Schrödinger is only one of many examples, even when concentrating on applications in physics. We can treat any operator in the form −div A∇ + V on any graph structure, provided that the signs of the coefficients yield an M-matrix. We can also address long range hopping for a very wide class of Hamiltonians.

Much more generally, in statistical physics, the probability distribution over all possible microstates (or the density matrix in the quantum setting) of a given system evolves through elementary jumps between microstates. This evolution is a Markov process whose transition matrix is a Z-matrix that is symmetric up to a multiplication by a diagonal matrix. For a micro-reversible evolution, the matrix A is symmetric and is akin to a weighted Laplacian on the high-dimensional indirect graph whose vertices are the microstates and whose edges are the possible transitions.

One essential result of statistical physics is that under the condition of irreducibility of the transition matrix, the system eventually reaches thermodynamical equilibrium. Our approach might open the way to unravel the structure of the eigenvectors of the Markov flow and thus to understand how localization of these eigenvectors can induce a many-body system to remain “frozen” for mesoscopic times out of equilibrium. This effect is referred to as many-body localization. A first successful implementation of the landscape theory in this context has been recently achieved by Balasubramanian et al.¹³ for a many-body system of spins with nearest-neighbor interaction. In this work, the authors cleverly use the ideas of Ref. 24 to transfer the problem to the Fock space and to deduce an Agmon-type decay governed by the corresponding effective potential. Once in the Fock space, their results are also a particular case of Theorem 2.7 and Theorem 5.2. From that point, however, the authors of Ref. 13 go much farther to discuss, based on physical considerations, deep implications of such an exponential decay on many-body localization, but in the present paper, we restrict ourselves to mathematics and will not enter those dangerous waters.

Finally, we would like to mention that the idea of trying the localization landscape and similar concepts in the generality of random matrices has appeared before, e.g., in Refs. 25 and 26. However, the authors relied on a different principle, extending the inequality |φ| ≤ Eu from Ref. 6 to these more general contexts, which by itself, of course, does not prove exponential decay. Reference 25 actually deals with a different proxy for the landscape and different inequalities, but we (and the authors) believe that these are related to the landscape and that, again, they do not prove exponential decay estimates. However, we would like to mention that the importance of M-matrices was already suggested in Ref. 26, and it was inspiring and reassuring to arrive at the same setting from such different points of view.

We ran numerical simulations to compute the localization landscape u, the effective potential $1u$⁠, and the eigenvectors for several realizations of random symmetric M-matrices. The diagonal coefficients are random variables that follow a centered normal law of variance 1. The off-diagonal coefficients belonging to the first W_c/2 diagonals of the upper triangle of the matrix are minus the absolute values of random variables following the same law. The remaining off-diagonal coefficients of the upper triangle are taken to be zero, and the lower triangle is completed by symmetry. This creates A₀, a Z-matrix of bandwidth W_c + 1 (and connectivity W_c). To ensure positivity, we add a multiple of the identity

with λ₀ being the smallest eigenvalue of A₀ and ɛ = 0.1. The smallest eigenvalue of the resulting matrix A is thus ɛ. The matrices A and A₀ clearly have the same eigenvectors, and their spectra differ only by a constant shift.

Below are the results of several simulations. We take N = 10³. Figures 1–4 correspond to random symmetric M-matrices constructed as above of connectivity W_c = 2, 6, 20, and 32. Each figure consists of two frames.

The top frame displays the localization landscape u superimposed with the first five eigenvectors plotted in log₁₀ scale. The exponential decay of the eigenvectors can clearly be observed on this frame for W_c = 2, 6, and 20. One can see that, as expected, it starts disappearing around W_c = 32 (W_c being in this case roughly equal to $N$⁠). It is important to observe that in all cases, the eigenvectors decay exponentially except for the wells of $1u$ (equivalently, the peaks of u) where they stay flat. This is exactly the prediction of Theorem 2.7.

The bottom frame displays the effective potential $1u$ superimposed with the first five eigenvectors plotted in linear scale. The horizontal lines indicate the energies of the corresponding eigenvectors. One can clearly see the localization of the eigenvectors inside the wells of the effective potential.

Figure 5 provides numerical evidence for finer effects encoded in Theorems 2.7 and 5.2. The two Theorems combined prove that exponential decay away from the wells of the effective potential governed the Agmon distance associated with 1/u, at least in the absence of resonances. In Fig. 5, we display, for several values of connectivity W_c and several eigenvectors, the values −ln |ψ_i| against the distance ρ_A,u,E(i, i_max), taking as the origin the point i_max where |ψ| is imal and using the corresponding eigenvalue as the threshold E. The linear correspondence down to e⁻⁴⁰ is quite remarkable and shows that $e−cρA,u,E(i,imax)$ is not only an upper bound but actually an approximation of the eigenfunction and that the resonances are indeed unlikely. On the other hand, the constant c does not appear to be equal to $1/Wc$⁠, which means that in this respect, our analysis is probably not optimal, at least in the class of random matrices. Indeed, we believe that the application of the deterministic Schur test in the proof does not yield the best possible constant for random coefficients, but since we emphasize the universal deterministic results, this step cannot be further improved.

Finally, Fig. 6 shows that Hypothesis 5.1 is actually fulfilled in some realistic situations. The top frame displays the example already presented in Fig. 2. Superimposed to the eigenvectors, the set K_E introduced in Definition 2.1 is also drawn (gray rectangles) for E = 0.7 (horizontal dashed red line). The middle frame displays the plot of the Agmon distance to K_E. Thresholding this plot at S = 2 (green horizontal line) allows us to draw the S-neighborhood of K_E (the orange rectangles). The bottom frame shows a possible partition of the entire domain into five subdomains (Ω₁, …, Ω₅), each subdomain containing at least one well of the effective potential 1/u. The distances ρ(∂⁻Ω_ℓ, K_ℓ) defined in Hypothesis 5.1 are here, respectively, 48.1584, 2.8093, 4.2169, 3.6784, and 9.3756. They all are larger than S, thus satisfying Hypothesis 5.1.

To be more precise, let us turn to the exact statements.

In this section, we prove Theorem 2.7. We will use a double commutator method. Let [A, B] ≔ AB − BA denote the usual commutator of N × N matrices and ⟨,⟩ denote the usual inner product on ℓ²([N]). We observe the general identity

whenever $A=(aij)i,j∈[N]$is a matrix, D = diag(d₁₁, …, d_nn) is a diagonal matrix, and $u=(ui)i∈[N]$is a vector. In particular, we have

whenever A is a Z-matrix and the entries of u have constant sign. It will be this negative definiteness property that is key to our arguments. One can compare (4.1) and (4.2) to the Schrödinger operator identity

for any (sufficiently well-behaved) functions $V,g,u:Rd→R$⁠.

To exploit (4.1), we will use the following identity:

We can now conclude as follows:

The strategy is then to apply this corollary with a sufficiently slowly varying function G so that one can hope to mostly control the right-hand side of (4.3) by the left-hand side.

Let the notation and hypotheses be as in Theorem 2.7. We abbreviate $ρ=ρA,u,Ē$ and $K=KĒ$⁠. To illustrate the decoupling phenomenon, we place the following hypothesis on the potential well set K:

We remark that axioms (i)–(iii) imply that the full boundary ∂Ω_ℓ, defined as the union of the inner boundary ∂⁻Ω_ℓ and the outer boundary ∂⁺Ω_ℓ consisting of those j ∉ Ω_ℓ such that a_kj ≠ 0 for some k ∈ Ω_ℓ, stays at a distance at S from K since every element of an outer boundary ∂⁺Ω_ℓ either lies in the inner boundary of another Ω_ℓ′ or else lies outside of all of the Ω_ℓ′. We also remark that axiom (iii) is a strengthening of axiom (ii), since if there was an element k in $Ωℓc$ at a distance less than S from K_ℓ, then by taking a geodesic path from K_ℓ to k, one would eventually encounter a counterexample to (iii), but we choose to retain the explicit mention of axiom (ii) to facilitate the discussion below.

Informally, to obey Hypothesis 5.1, one should first partition K into “connected components” K_ℓ, concatenating two such components together if their separation ρ is too small so that the separation $S̄:=infℓ≠ℓ′ρ(Kℓ,Kℓ′)$ is large and then perform a Voronoi-type partition in which Ω_ℓ consists of those k ∈ [N] that lie closer to K_ℓ in the ρ metric than any other K_ℓ′. The axioms (i) and (ii) would then be satisfied for any $S<S̄/2$⁠, thanks to the triangle inequality, and when $S̄$ is large, one would expect axiom (iii) to also be obeyed if we reduce S slightly. It seems plausible that one could weaken axiom (iii) and still obtain decoupling results comparable to those presented here, but in this paper, we retain this (relatively strong) axiom in order to illustrate the main ideas.

We have already demonstrated in Sec. III non-vacuousness of Hypothesis 5.1, at least in typical numerical examples. Let us say a few more words in this direction. Recall the simulations in Sec. III. Much as there, let us assume for the moment that we are working with a band matrix and W is the band width, that is, a_ij = 0 whenever |i − j| > W. One can deduce a rather trivial lower bound for the Agmon distance associated with v as in Definition 2.1. If v_i ≥ v_min for all i in an interval $I=[i1,iq]∪N$ of length q, then the Agmon distance between two components of the complement of I is

Here, the lower bracket as usual stands for the floor function. The above inequality follows directly from the observation that the number of non-trivial components such that $viℓ≠0$ and $viℓ+1≠0$ and $aiℓiℓ+1≠0$ of the path from i₁ − 1 to i_q + 1 is at least $iq−i1+1−WW$⁠. Going back to our definitions and fixing some $Ē>E$ and the respective partition of $KĒ=∪ℓKℓ$ into disjoint components, we denote by d the minimal “Euclidean” distance between the components, i.e.,

It is in our interest to make this distance (or rather the corresponding Agmon distance) substantial, so we might combine several disjoint components into one K_ℓ. With this at hand, we choose Ω_ℓ to be imal possible neighborhoods of K_ℓ, which are still disjoint. Since the inner boundary ∂⁻Ω_ℓ consists of i ∈ Ω_ℓ such that j ∉ Ω_ℓ and a_ij ≠ 0, that is, has “width” at most W, one can see that with the aforementioned choices, the “Euclidean” distance between K_ℓ and ∂⁻Ω_ℓ

or, to be more precise, $d+12−W$⁠. By design, the complement of K_ℓ in Ω_ℓ consists of points such that $vi>Ē$⁠, that is, there is v_min > 0 such that v > v_min in Ω_ℓ\(∂⁻Ω_ℓ ∪ K_ℓ). Hence, at the very least, for this v_min > 0, we have

Clearly, we could take a smaller subinterval of Ω_ℓ\(∂⁻Ω_ℓ ∪ K_ℓ) and make v_min > 0 larger, not to mention that this is a trivial lower estimate that does not take into account high values of v. In any case, this demonstrates that Hypothesis 5.1 is non-vacuous.

Let ψ^j denote the complete system of orthonormal eigenvectors of A on [N] with eigenvalues λ_j. Let Ψ_(a,b) denote the orthogonal projection in ℓ²([N]) onto the span of eigenvectors ψ^j with the eigenvalue λ_j ∈ (a, b). For a fixed ℓ, let φ^ℓ,j denote a complete orthonormal system of the local eigenvectors of A on Ω_ℓ with eigenvalues μ_ℓ,j and let Φ_(a,b) be the orthogonal projection onto the span of the eigenvectors φ^ℓ,j with the eigenvalue μ_ℓ,j ∈ (a, b) over all ℓ and j.

The goal of this section is to prove that under the assumption of Hypothesis 5.1, A can be almost decoupled according to $⋃$_ℓΩ_ℓ, with the coupling exponentially small in S. More precisely, we have the following result, which is an analog of Ref. 7 (Theorem 5.1) in the M-matrix setting:

We are thankful to Guy David for many discussions in the beginning of this project and Sasha Sodin for providing the literature and context with regard to the edge state localization.

M.F. was supported by the Simons Foundation (Grant No. 601944). S.M. was partly supported by the NSF RAISE-TAQS (Grant No. DMS-1839077) and the Simons Foundation (Grant No. 563916). T.T. was supported by the NSF (Grant No. DMS-1764034) and a Simons Investigator Award.

The data that support the findings of this study are available from the corresponding author upon reasonable request.