Anderson transition on the Bethe lattice: an approach with real energies

TL;DR

This paper investigates the Anderson localization transition on the Bethe lattice using real energy propagators, introducing a new stability criterion, providing precise transition estimates, and connecting analytic approximations with asymptotic analysis.

Contribution

It presents a novel stability-based criterion for localization transition, accurate transition point estimates, and links forward approximation with asymptotic results.

Findings

New stability criterion for localization transition

Precise numerical estimate of critical disorder

Connection between forward approximation and asymptotic analysis

Abstract

We study the Anderson model on the Bethe lattice by working directly with propagators at real energies $E$. We introduce a novel criterion for the localization-delocalization transition based on the stability of the population of the propagators, and show that it is consistent with the one obtained through the study of the imaginary part of the self-energy. We present an accurate numerical estimate of the transition point, as well as a concise proof of the asymptotic formula for the critical disorder on lattices of large connectivity, as given in [P.W. Anderson 1958]. We discuss how the forward approximation used in analytic treatments of localization problems fits into this scenario and how one can interpolate between it and the correct asymptotic analysis.

Equations

$$H=\sum _i{ϵ}_i|i⟩⟨i|+t\sum _{⟨i,j⟩}\left(|i⟩⟨j|+|j⟩⟨i|\right),$$

$${\mathcal{G}}_{ii}(E)=⟨i|\frac{1}{E-H}|i⟩,$$

$${\mathcal{G}}_{ii}(E)=\frac{1}{E-{ϵ}_i-{\sum }_{k=1}^{K+1}{G}_{k,k}(E)},$$

$${\mathcal{G}}_{ii}(E)=\frac{1}{E-{ϵ}_i-{\sum }_{k=1}^{K+1}{\sum }_{j=1}^{K+1}{G}_{k,j}(E)},$$

$$G_{k,j}=\sum _{\alpha }\frac{{\psi }_k^{\alpha }{\psi }_j^{\alpha }}{E-E^{\alpha }},$$

$$G_i(E)=\frac{t^2}{E-{ϵ}_i-{\sum }_{k=1}^K{G}_k(E)},$$

$$P(\hat{R},\hat{I})\text{mbox}andP(G).$$

$$\rho (E)=\int d\hat{R}d\hat{I}P(\hat{R},\hat{I})\hat{I}\text{mbox}and\underset{G\to \infty }{\lim }{G}^{2}P(G)=\rho (E).$$

$$P(G)=\int d\hat{R}d\hat{I}P(\hat{R},\hat{I})\frac{\hat{I}}{(\hat{R}-G{)}^2+{\hat{I}}^2}.$$

$$G_i=\frac{-1}{{\sum }_{k=1}^K{G}_k},$$

$$P(G)=\frac{(\pi \sqrt{K}{)}^{-1}}{K^{-1}+G^2}.$$

$$\mathcal{P}(\mathcal{G})=\frac{b/\pi }{b^2+{\mathcal{G}}^2},$$

$$P(G)$$

$$|i⟩=(H-E)\sum _j{\mathcal{G}}_{ij}^{RRG}|j⟩,$$

$$P(G)\simeq\frac{1}{1+4K/W^{2}}\left\{\begin{array}[]{@{\kern 2.5pt}lL}\hfill(WG^{2})^{-1}&if $|G|\geq 2/W$\\ K/W&if $|G|<2/W$,\end{array}\right.$$

$${\chi }_p=\frac{\partial G_j}{\partial {ϵ}_i}=\left(\prod _{k=1}^L\frac{\partial G_{k+1}}{\partial G_k}\right)\frac{\partial G_i}{\partial {ϵ}_i}\propto \prod _{k=0}^L{G}_k^2,$$

$${\chi }_L=\sum _{p=1}^{K^L}{\chi }_p\propto \sum _{p=1}^{K^L}\prod _{k\in p}{G}_k^2.$$

$$\Lambda =\underset{L\to \infty }{\lim }\frac{1}{L}\ln {\chi }_L.$$

$$⟨{\chi }^s⟩={\int }_0^{\infty }d\chi Q(\chi ){\chi }^s=⟨\prod _{k=1}^L{G}_k^{2s}⟩=C_L\lambda (s{)}^L,$$

$$Q(\chi )=C_L{\int }_B\frac{ds}{2\pi i}{\chi }^{-s-1}{e}^{L\mu (s)},$$

$$Q(\chi)=C_{L}\int_{B}\frac{ds}{2\pi i}e^{-L[(s+1)\theta-\mu(s)]}\Big{|}_{L\theta=\log\chi}.$$

$$Q(\theta )=C_L{\int }_B\frac{ds}{2\pi i}{e}^{-L(s\theta -\mu (s))}.$$

$${\mu }^{\prime }(s^{\ast })=\theta .$$

$$Q(\theta )\simeq e^{-L(s^{\ast }\theta -\mu (s^{\ast }))}.$$

$${\theta }_{\mathrm{typ}}=\underset{s\to 0}{\lim }\frac{\partial \mu (s)}{\partial s},$$

$$\ln K+\mu (s^{\ast })=0,$$

$$K\lambda (s^{\ast }=1/2)=1,$$

$$Q(\chi )=e^{-L(\ln K+\frac{3}{2}\theta +\frac{{\theta }^2}{2{\mu }^{\prime \prime }(1/2)})}=\frac{K^{-L}}{{\chi }^{3/2}}{e}^{-\frac{(\ln \chi {)}^2}{2L{\mu }^{\prime \prime }(1/2)}}$$

$$G_{k+1}=-\frac{1}{ϵ+\zeta +G_k},$$

$$P(G_{k+1}|G_k)={\mathbb{E}}_{ϵ,\zeta }\left[\delta \left(G_{k+1}+\frac{1}{ϵ+\zeta +G_k}\right)\right]$$

Full text

• May 2019

1 Introduction

Anderson Localization [1, 2] is one of the deepest and most universal quantum mechanical implications for the dynamics of disordered systems. Under certain circumstances (small spatial dimension, large disorder) quantum mechanical effects can be so strong as to completely hinder transport in models of non-interacting particles hopping on lattices (the so-called Anderson models). When this happens, one speaks of an Anderson localized phase. Although the localized phase has come under almost complete mathematical control [3, 4], the complementary (“delocalized”) phase, and in particular the phase transition between them has proved elusive to most analytic control, despite the constantly improving abundance and quality of numerical results [5]. Some exceptions to this is provided by the $\epsilon$ expansion of Efetov [6], for the Anderson model on a $d$-dimensional square lattice $\mathbb{Z}^{d}$ around $d=2$. However, these results fail to predict the behavior of the critical point for $d\gtrsim 3$ (see for example [7]).

Following the similarity with conventional thermodynamic transitions, one is then led to consider loop-less lattices as proxies for the $d\to\infty$ limit. The simplest example of loop-less geometry is probably the Bethe lattice and several works, starting with [8] and including [9], have dealt with the Anderson model on the Bethe lattice. The Bethe lattice, as the infinite limit of a regular random graph [10], might look like a peculiar geometry to study a nearest-neighbor hopping problem, but it is quite common in spin glasses models for two apparently disconnected reasons: on one hand, due to absence of short loops in the lattice [11, 12], one can make analytic progress in writing a mean field theory in terms of recursion equations, see Ref. [8]; on the other hand this lattice structure naturally arises in random optimization problems [13] and therefore the results have great relevance for applications [14]. For Anderson Localization (AL) a further reason that is behind the recent renaissance of interest in the problem lies in the connection, first emphasized in Ref. [15], between AL on the Bethe lattice and localization effects in interacting many-body systems, a topic now known as Many-Body Localization (MBL)[16, 17]. MBL has been at the center of a prolific research activity [18, 19, 20, 21, 22], which, in addition to the absence of transport (a trait in common with AL), has uncovered some peculiarities like the behavior of entanglement measures [23, 24, 25], emergent integrability [26, 27, 28, 29], and protection of symmetries at high temperature [30]. In a recent line of research, some of the MBL phenomenology has been found in highly excited states of systems without disorder [31, 32, 33] and even lattice gauge theories [34].

Pushing ahead on the initial motivation for connecting the Anderson model on the Bethe lattice with MBL (where regions of sub-diffusing transport are seen to exist within the ergodic phase [35, 36, 37, 38, 39, 40]), a conjecture mostly based on numerical evidence has been put forward [41] (see also the previous works [15, 42]), and consequently made into a phenomenology [43, 44]. According to these papers, the eigenstates of the Anderson model on finite $N$-vertices Regular Random Graphs (RRG) are non-ergodic, in that they possess anomalous dimensions, in a way similar to what is found in models of random matrices [45, 46]. More recently, a corresponding putative dynamical transition in the return probability [47] has been observed. This picture has been challenged by other numerical [48, 49] and analytical studies [50, 51], that support the scenario according to which ergodicity is restored for system-sizes exceeding a certain critical volume. This issue can be roughly summarized into the question: “are all the properties of the eigenfunctions and of the dynamics encoded in the recursion equations introduced in Ref. [8]?” While the question of ergodicity or lack thereof concerns finite lattices, and in particular regular random graphs with $N$ vertices having plenty of loops of length $O(\ln N)$, we feel that a more detailed numerical and analytical knowledge of the recursion relations in Ref. [8] is welcome, at least to dispel claims that some of the features observed in the RRG numerics are directly implied by the recursion equations.

In this work we provide a new perspective on the recursion relations and derive implications on the above and related issues. Our strategy will be to work with propagators computed directly at real values of the energy, at variance with the more conventional approaches in which the energies are treated as complex variables. This choice makes the numerics more convenient and more readily comparable to the exact diagonalization results, while the analytical calculations are not more complicated than in the conventional approach.

The organization of the article is as follows: in Section 2 we discuss how the Anderson transition can be described in terms of the susceptibilities of the populations of propagators, without introducing an imaginary part for the energy. In Section 3 we show how the integral equation introduced in [8] to determine the stability of the Anderson phase can be recovered from the probability distribution of the susceptibilities. We discuss some features of its solution for two different distributions of the disorder: the uniform distribution (conventionally considered in the numerics) and the Cauchy distribution, for which analytic results can be derived. Using the formalism based on the integral equation, we derive analytically the weak-disorder behaviour of the fractal dimension $D_{1}$ defined as in Ref. [44]. In Section 4 we revisit the large-connectivity limit of the problem, bringing up a new point of view on the integral equation. We present a concise proof of the known asymptotic formula [1, 8, 52] for the critical disorder $W_{c}$ at large connectivity, which relies on the identification of the eigenvector corresponding to the largest eigenvalue of the integral kernel. We then discuss why the asymptotic formula is not correctly recovered within the forward or “upper limit” approximation of Ref. [1], and introduce a generalized integral equation that interpolates between the forward approximation [53] and to the full Anderson problem at large connectivity.

2 Real-energy propagators: the criterion for localization

2.1 On the inverted thermodynamic limit

We consider the Anderson model on a Random Regular Graph with connectivity $Z=K+1$ ($K$ being the branching number). In the infinite volume limit, any finite section of this graph becomes tree-like. The Anderson model has the Hamiltonian:

[TABLE]

where $i$ labels the vertices of the graph, the second sum is over nearest neighboring vertices, and $\epsilon_{i}$ are identically distributed independent random variables. We mostly focus on two distributions of disorder $\rho(\epsilon)$: a uniform distribution with $\epsilon_{i}\in[-W/2,W/2]$, and the Cauchy distribution $\rho(\epsilon_{i})=(W/\pi)/(W^{2}+\epsilon_{i}^{2}).$ The latter is particularly interesting as it allows to derive analytic results, mainly due to its stability with respect to the iteration of the recursion relations. We notice that the ubiquity of the Cauchy distribution for random Schrödinger operators has been discussed in the mathematical physics literature in the past [4] and cannot possibly be overemphasized.

In the following, we set $t=1$ for simplicity, and we focus on the case $E=0$, $E$ being the energy, unless stated otherwise. In the numerics we always consider $K=2$. The diagonal part of the resolvent, i.e., the Fourier transformed survival amplitude:

[TABLE]

can be computed from the approximated equations [8]

[TABLE]

where the sum is over all sites $k$ that are nearest neighbors of $i$, and the $G_{k,k}(E)$’s are cavity Green functions: they correspond to the Hamiltonian restricted to the sub-tree rooted in site $k$, once the link connecting sites $i$ and $k$ is removed. These equations are expected to become exact in the infinite volume limit. The exact relation, valid for any size and any graph, is of the form:

[TABLE]

which reduces to the previous equation if the off-diagonal terms $G_{k,j}$ with $k\neq j$ are negligible for large $N$.

Why should the off-diagonal terms be negligible? We first note that the points $k$ and $j$, in absence of $i$, are not directly connected: they are not connected at all on a tree geometry while they are at distance of the order of $\ln(N)/\ln(K)$ on a RRG. We also notice that

[TABLE]

where the $\psi^{\alpha}_{k}$ and the $E^{\alpha}$ are the eigenvectors and the eigenvalues of the cavity Hamiltonian. The typical values of the quantities $\psi^{\alpha}_{k}\psi^{\alpha}_{j}$ are of order $1/N$ in the case of ergodic states, and are exponentially small in the distance between $j$ and $k$ in the case of localized states. When $E$ has an imaginary part the off-diagonal terms are always small. The same happens for real $E$ and finite $N$, where problems may arise when the value of $E$ belongs to the support of the spectrum, as the denominators vanish. In the case of localized states the off-diagonal terms are small; also in the case of extended non-ergodic states the off-diagonal terms can be neglected, while they have to be considered in the case of extended ergodic states. Under the hypothesis that the diagonal terms can be neglected, by following the usual cavity construction one can easily show that the cavity Green functions at energy $E$ satisfy the recursion

[TABLE]

where the sum is over the $K$ cavity neighbors of site $i$ and we introduced the shorthand notation $G_{k}=G_{k,k}$ for the diagonal terms. These equations are defined here for a finite tree geometry; the usual route to analyze them consists in adding a small imaginary part to $E\to E+i\eta$, which is sent to zero after the thermodynamic limit is taken. In the following we shall follow an alternative route and consider $\eta\to 0$. The inversion of the limits $\eta\to 0$ and $N\to\infty$, dubbed “inverted thermodynamic limit” in Ref. [44], has been considered in [44, 54, 55] to deal with the eigenfunction statistics. It is clear that in the inverted thermodynamic limit the $G$’s are always real; in the standard thermodynamic limit, in the extended case the equation admits one solution where the $G$’s are real, that is however unstable with respect to the addition of imaginary parts.

The use of the simplified cavity equation without the off-diagonal terms may look unfounded, as far as it is easy to convince ourselves that the off-diagonal terms cannot be neglected in the extended ergodic phase. Indeed, the equations cannot be used in this phase unless one proves that the extended states are non ergodic. Of course we can (and will) prove that these equations can be used to show that the extended states are not ergodic. Before proceeding, we have to face a serious consistency problem. In the standard thermodynamic limit in the extended non-ergodic case we have a probability distribution of a complex $G=\hat{R}+i\hat{I}$ that satisfies the same standard cavity equation in the inverted thermodynamic limit with a probability distribution of the real $G$. At the end of the day we are left with two probability distributions:

[TABLE]

Both probability distributions describe the physics of the problem. The consistency problem becomes more acute if we notice that the spectral density $\rho(E)$ can be extracted by both equations, as:

[TABLE]

This consistency problem has a simple solution: the two probability functions are not independent. Indeed, $P(G)$ can be written as a simple function of $P(\hat{R},\hat{I})$

[TABLE]

This relation was already given in [56], and it is quite possible that it has its root in the hyperbolic geometry of the $O(n|n)$ group in the $n\to 0$ limit. The proof is quite direct. It makes use the basic properties of the Cauchy distribution: it is quite simple to check that if $P(\hat{R},\hat{I})$ satisfies the complex cavity equation, the function $P(G)$ defined in Eq. (9) satisfies the real cavity equation. Notice that, as mentioned before, the inconsistency exists only in the extended non-ergodic case: in the localized case the distribution $P(\hat{R},\hat{I})$ is concentrated on $\hat{I}=0$, and in the extended ergodic case there is no simple equation for $P(G)$.

Having checked that the results of the inverted thermodynamic limit are consistent with those of the standard thermodynamic limit, we will devote the rest of the paper to the study of this limit.

2.2 Vanishing disorder: the emergence of the Cauchy distribution

Consider first the limiting case of vanishing disorder $W=0$. Setting $E=0$, the equation for $G_{i}$ reduces to

[TABLE]

which has no unique real solution (assuming $G_{i}=G$). This is natural, since the resolvent is undefined for values of $E$ belonging to the spectrum of the corresponding operator (in this case, the adjacency matrix of the graph). However, the above equations admit a solution in distribution in the form of a Cauchy-Lorentz distribution of width $1/\sqrt{K}$

[TABLE]

This, in turn, induces a Cauchy-Lorentz distribution for the $\mathcal{G}_{ii}$

[TABLE]

where $b=\sqrt{K}/(K+1)$. This agrees with numerical results, as shown in Figure 1, where the distribution of $\mathcal{P}$ in (12) is compared with numerical solutions of Eq. (2) on an ensemble of RRGs. The fact that (for random operators) the diagonal matrix elements of the resolvent evaluated on the real axis are Cauchy variables is a property that holds in surprising generality [57, 58], irrespectively of the statistics of the matrix eigenvalues and of the possible strong correlations between them, as exemplified by Coulomb gas models [59, 60, 58]. Moreover, the fact that the numerics on RRG is so accurately described by a cavity equation means that the correlations along closed paths of length $\ln_{K}N$ which loop around the RRG are irrelevant for this quantity, and that the $G_{k}$’s can be considered as independent, identically distributed variables. This is due to the presence of a gap in the spectrum of the integral equation associated with the susceptibility, as also shown below.

2.3 Finite disorder: population dynamics and exact numerics on RRG

In presence of disorder, the stable probability distribution that describes the random variables $G_{i}$ satisfies the cavity equation

[TABLE]

We determine the distribution of the cavity propagators $G_{i}$ and of $\mathcal{G}_{ii}$ by using a (by now standard) population dynamics algorithm: we iterate Eq. (6) over a population of size $n$ (up to $n=10^{7}$), initially drawn from a uniform distribution. We perform $I=100$ sweeps of the population to achieve convergence, replacing each element of the population (the new elements are thus available to be drawn before the sweep is complete). We then sample the population $s=10$ times at iteration $i=100,115\dots 250$. This gives the same result as iterating a population of size $N=ns$ directly.

To compare with exact results on the RRG geometry, we also determine $\mathcal{G}_{ii}$ by numerically solving the linear equation

[TABLE]

using for the kinetic term in $H$ the adjacency matrix of the regular random graph (with hopping $t=1$), which takes loops into consideration. The above is the equation for a vector $v=\sum_{j}\mathcal{G}^{RRG}_{ij}\left|j\right\rangle$ in the canonical form $v_{0}=Mv,$ where $v_{0}$ and $M$ are known. Since linear solvers are more efficient than exact-diagonalization algorithms, a lot of statistics can be accumulated in a short time.

The distribution of $\mathcal{G}^{RRG}_{ii}$ and the one of $\mathcal{G}_{ii}$ obtained from the analysis of the population dynamics agree within the precision of the numerical computations. For $W>0$, the distribution $P(G)$ for boxed disorder shows a depletion around $G=0$ (see Fig. 1, as well as the comparison with the exact numerical RRG result in Fig. 2). A consequence of it is that the maximum of $P(G)$ moves from $G=0$ to $G>0$ discontinuously at a certain $W\simeq 2.2$, see Fig. 2. For large $W$, as we argue in the following, $P(G)\simeq G^{-2}\rho(1/G)$, giving $P(G)\propto G^{-2}$ for large-enough $G$. More precisely, a good large-disorder approximation of P(G) is given by:

[TABLE]

where the small-$G$ behaviour is obtained using the tails of the distribution of $G_{i}$ together with $G\sim(\sum_{k=1}^{K}G_{i})^{-1}$. This yields the location of the maximum at $G=2/W$ for box disorder distribution. This is a good estimate at sufficiently high disorder (say $W>6$ for $K=2$), as shown in the inset of Fig. 2. The tails of the distribution behave as in Eq. (12), i.e. $a/\pi(a^{2}+G^{2})$, for any disorder (with $a\to K^{-1/2}$ as $W\to 0$). The coefficient $a$ can be obtained from a direct fit of the tails of $P(G)$; an alternative, more convenient way to extract it numerically is through the distribution $p(l)$ of the quantity $l=(1+W/2G)^{-1}$, as $a=\lim_{l\to 0}p(l)/W$. The coefficient $a$ can then be used (with proper normalization) to extrapolate the distribution $P(G)$ on the full real axis.

By looking at real populations one cannot directly see if the system is in a localized or delocalized phase. Rather, the populations must be used as input for a more refined analysis, which is done in the next sections.

2.4 A criterion for the transition based on susceptibilities

When $E=E_{r}+i\eta$, the localized (delocalized) phase is characterized by the (in)stability of the phase in which $\Im G$ is concentrated on its limiting value, related to the local density of states. For real $E$, one can formulate the problem in terms of the stability of a given population, or at the decay of the susceptibility. In the localized region, the change of the value of an energy level at position $i$ should not result in a modified dynamics of a site $j$, when the distance between $i$ and $j$, $d(i,j)=L$, goes to infinity (in the sense that $d\gg\xi$, the localization length). Therefore one is led to consider the decay with $L$ of the telescopic identity associated with the shortest unique path $p$ from $i$ to $j$

[TABLE]

where $k$’s are all the sites on the (single) path between $j$ and $i$, where $k=0$ corresponds to the site $i$ and $k=L+1$ to the end site $j$.

If one considers a RRG, the number of paths going from one site to the other is one for $L\lesssim\ln N/\ln K$ but quickly becomes $K^{L}$ when $L\gtrsim\ln N/\ln K$ (see [10, 61]). The susceptibility is then obtained by summing (16) over $K^{L}$ paths,

[TABLE]

The Lyapunov exponent governing the large-$L$ behavior of the susceptibility is defined by

[TABLE]

While the susceptibility $\chi_{L}$ is a wildly fluctuating random number, the Lyapunov exponent $\Lambda$ is instead self-averaging. The system is localized whenever its typical value is negative: the values of parameters at which $\Lambda=0$ correspond to the transition point.

To get an expression for the Lyapunov exponent $\Lambda$ let us first consider the distribution $Q(\chi)$ of a single path contribution $\chi_{p}$ (we drop for convenience the index $p$). This can be extracted from its $s$-th moment, interpreted as a Mellin transform of $Q(\chi)$:

[TABLE]

where we have absorbed the exponential dependence in $\lambda(s)$, and thus assumed that $C_{L}$ grows less than exponentially in $L$. We introduce a generalized exponent $\mu(s)$ through the equality $\lambda(s)=e^{\mu(s)}$ and invert the Mellin transform, getting

[TABLE]

where the integral is done over the Bromwich path, parallel to the imaginary axis. Using the notation $\theta=(\ln\chi)/L$ we have

[TABLE]

Equivalently, the distribution $Q(\theta)$ of the rescaled variable $\theta$ reads

[TABLE]

For large $L$, the integrals can be computed by the saddle point method, by passing through the real saddle point $s^{*}$, in a direction parallel to the imaginary axis, where $s^{*}=s^{*}(\theta)$ is an implicit function of $\theta$ defined by the condition:

[TABLE]

Plugging this back in (22) we have

[TABLE]

The value $\theta=0$ discriminates between an exponentially growing ($\theta>0$) and an exponentially decaying ($\theta<0$) susceptibility along a single path. The typical value of the path susceptibility is governed by $\theta_{\rm typ}$ satisfying

[TABLE]

implying the vanishing of the large-deviation function for $\theta$, see Eq. (22). The transition between the localized and delocalized phases is actually governed by atypically large fluctuations of the single-path susceptibilities; the typical decay (25) has however been considered in the literature111Notice that these references consider the typical value of the squared susceptibility.: it appears in a sufficient criterion for delocalization in Ref. [62], and it has been conjectured to be related to the transition between ergodic and non-ergodic phases in the delocalized phase in Ref. [44]. We discuss this more extensively in the next section.

A change of variables from $\theta$ to $\chi$ shows that the distribution of the product $\prod_{k\in p}G_{k}^{2}$ has long tails, implying that the sum over paths in (17) is dominated by the maximal among the summands. The typical value of the total susceptibility thus solves the equation $K^{L}Q(\chi_{L}^{\rm typ})=1$. The transition corresponds to this typical value becoming equal to one (corresponding to $\Lambda=0$), and thus it is obtained setting $\theta=0$ in (24). In this case, one is left with the equation $\mu^{\prime}(s^{*})=0$ to solve. As first noted in Ref. [8] (and in a more general setting in Ref. [63]), the symmetry of the equation relating $G_{k}$ to $G_{k-1}$ implies that $\mu^{\prime}(s)=0$ iff $s=1/2$. Keeping this in mind, the criterion for the transition reduces to:

[TABLE]

or

[TABLE]

where $\mu(s)=\ln\lambda(s)$, which is exactly the resonant criterion derived in Ref. [4]. This is not surprising, as the total susceptibility appears naturally when linearizing the self-consistent equations for the self energy around the solution having zero imaginary part [4], which is stable in the localized phase. Close to the critical point, and for $\chi=O(1)$, one can even write the full probability distribution of $\chi$ by noticing that, close to $\theta=0$, by writing $s^{*}=1/2+\epsilon$ we have $\epsilon=\theta/\mu^{\prime\prime}(1/2)$ and so

[TABLE]

is a log-normal distribution.

There are two ways of tackling the calculation of $\lambda(s)$: one can directly compute the distance between two populations initially differing by a small quantity, or one can link $\lambda(s)$ to the eigenvalue of an integral kernel, mapping explicitly to the formalism of Ref. [8].

3 The integral equation for $\lambda(s)$ and the Anderson transition

The product of $G_{i}^{2s}$ along a path is a product of correlated random variables. In the forward or upper limit approximation, these correlations are neglected and the propagators are treated as independent. To go beyond this approximation, we notice that two consecutive propagators along a path are related by the equation

[TABLE]

where $\epsilon$ is the onsite energy at site $k$, and $\zeta=\sum_{j=1}^{K-1}G_{j}$ with $G_{j}$ i.i.d. random variables with distribution $P(G)$. The conditional probability

[TABLE]

defines an integral operator $\hat{K}$ with the non-symmetric kernel $K(y,x)\equiv P(y|x).$ The operator is Markovian, since

[TABLE]

It thus defines a continuous state Markov chain. Note that the chain is non-reversible, since no detailed balance condition holds; however, the following symmetry [8, 63] holds:

[TABLE]

The stable distribution $P(G)$ solving (13) is a right eigenvector of $\hat{K}$ with eigenvalue $1$; the corresponding left eigenvector is the function with constant value $1$, as seen from Eq. (31).

The Kernel $K$ can be used to compute $\lambda(s)$ by writing

[TABLE]

Introducing the vector space notation:

[TABLE]

we have

[TABLE]

Following the usual arguments on transfer matrices, we see that $\lambda(s)$ is the largest eigenvalue of the integral Kernel

[TABLE]

i.e., it is the largest solution of the equation

[TABLE]

This integral equation is exactly the same one that appears in Ref. [8], where it was derived by investigating the tails of the probability density of the imaginary part of the self energies (equivalently, the behavior of its characteristic function close to the origin).

The integral operator can also be exploited to determine the correlation functions between subsequent cavity propagators along a path in the lattice. Since $K$ is a contraction, being the average of contractions, then $\lambda=1$ is its largest eigenvalue. Let $\lambda_{1}$ be the second largest eigenvalue, with right and left eigenvectors $\phi_{1},\psi_{1}$ respectively and assume there is a finite gap between the two:

[TABLE]

The connected part of the correlation function

[TABLE]

decays exponentially with $r$, as follows from Eq.(37), with leading behavior

[TABLE]

that is dictated by the second largest eigenvalue of the kernel $K$, irrespectively of the value of $s$, assuming of course that the quantity on the right-hand side of (39) is finite (which sets the bounds $0\leq s<1/2$). One can define a correlation length $\xi=-1/\ln\lambda_{1}$, and show (analytically for the Cauchy-Lorentz disorder and numerically for box disorder) that it is always $\lambda_{1}<1$, which means that $\xi$ is always finite. In this respect the Anderson model is not different from other statistical models on the Bethe lattice in that the correlation length does not diverge.

3.1 Box disorder distribution

For box disorder in $[-W/2,W/2]$, we compute the kernel by first determining the probability distribution function $P(G)$ of the cavity propagators with population dynamics, as outlined in Sec. 2.3. We use a cubic spline interpolator, which is efficient to integrate numerically, plus the Cauchy-like tails $\sim G^{-2}$ for $P(G)$. We then compute (35) on a discretized $(x,y)$ grid:

[TABLE]

where $\mathcal{P}_{\zeta}(\cdot)$ denotes the distribution of $\zeta=\sum_{j=1}^{K-1}G_{j}$. The discretization is parametrized with $(x,y)=(\tan\theta_{1},\tan\theta_{2})$ where $\theta\in[-\pi/2,\pi/2]$ with an equal spacing in both directions $\Delta\theta=10^{-3}$.222We have changed $\Delta\theta$ over the range $[0.5\cdot 10^{-3},2\cdot 10^{-2}]$ and for population size $N=3\cdot 10^{8}$. The dependence on $\Delta\theta$ is quite small (like $\Delta\theta^{3}$) and the estimated error we associate to the discretization is of one part in $10^{4}$.

We compute the largest eigenvalue $\lambda^{(N)}$ (for a population of given size $N$) of the discretized matrix $K_{s}$ by means of exact diagonalization, with the Arnoldi iterative algorithm. For $s=0$ and $1$, we know that $\lambda=1$; we use this as a criterion for the goodness of the discretization of $K$ (other errors, coming from the numerical integration of $P(G)$ and the exact diagonalization, are otherwise well controlled), resulting in an error less than $0.1\%$ for the chosen value of $\Delta\theta=10^{-3}$ and population size $N=ns=10^{9}$. The dependence of the resulting $\lambda^{(N)}$ on the population size $N$ (which led in a previous version of this manuscript to overestimate $W_{c}$ by about $1.5\%$) is significantly reduced if one recognizes that asymptotically $P(G)\propto a/G^{2}$ and that this asymptote is quickly reached for $G=O(1)$. The coefficient $a$ can then be extracted from a fit to typical values of $G$ with high accuracy. Values of $sN=10^{9}$ used here are then more than sufficient to reach the precision quoted in this paragraph.333If one neglects this information about the functional form of the tail of $P(G)$, and fits considering not typical but extreme values of $G$, a dependence of $\lambda^{(N)}=c_{0}+c/\log(N)^{2}$ is observed for $N<N^{*}$ which however saturates to $\lambda^{\infty}$ for $N\gg N^{*}$. This dependence can be traced back to the existence of a critical volume which scales like $e^{C/|W_{c}-W|^{1/2}}$ as $W\to W_{c}$.

In Fig. 3 we plot the resulting largest eigenvalue $\lambda(s)$ as a function of the exponent $s$. Since the minimum at $s=1/2$ is the value of interest for the transition, we also plot $\lambda(s=1/2)$ as a function of disorder. Here, the crossing of $\lambda=1/2$ happens at $W_{c}=18.11\pm 0.02$, which is the localization transition point. This is in agreement with the most recent numerical results in [64].

We now consider the behavior of the typical value $\theta_{\rm typ}=\partial\lambda(s)/\partial s|_{s=0}$ at small values of disorder. For box disorder $\partial\theta_{\rm typ}/\partial W=0$ at $W=0$, and we conjecture that this happens for all disorder distributions which fall to zero faster than $1/\epsilon^{2}$ at large $\epsilon$. We can show this by computing $\partial P(G)/\partial W$ at small $W$. The integral equation to be solved is

[TABLE]

This is a linear integral equation (a tangent map) for $v(y)=\frac{\partial P}{\partial W}(y)$. For $W=0$ we can use the result (11) and compute the integrals explicitly. The second integral can be performed easily (as $\rho(\epsilon,W=0)=\delta(\epsilon)$) while in the first one it is convenient to keep $\rho(\epsilon)=W^{-1}\theta(W/2-|\epsilon|)$ in the integral, and then expand for small $W$.

After some algebra one finds, for $K=2$ and $W\ll 1$:

[TABLE]

We have

[TABLE]

Without solving the equation (42) one can see that $v(x)=O(W)$ so we have $\theta_{\rm typ}(W)\propto W^{2}$ for small $W$ (up to constants that do not depend on $W$). Notice that this is only true for box disorder: for Cauchy disorder the corresponding first term in the right-hand side of (42) is $O(1)$ instead of $O(W)$ and this gives $\theta_{\rm typ}(W)\propto W$, see the following section.

3.2 Cauchy on-site randomness

In the case of Cauchy on-site randomness,

[TABLE]

the integral operator with Kernel (30) (corresponding to $s=0$) can be explicitly diagonalized. In this case the equilibrium distribution of the cavity fields at $E=0$ is itself a Cauchy distribution

[TABLE]

with width

[TABLE]

Plugging this into the definition of the Kernel one finds

[TABLE]

where the two constants $\alpha,\overline{\alpha}$ are related by:

[TABLE]

The Kernel, being not symmetric, has distinct left and right eigenvectors. As we show in the Appendix, the eigenvalues are:

[TABLE]

for $n=0,1,\cdots$, where $\lambda_{0}\equiv\lambda=1$. The connected component of the correlation function between the square modulus of the cavity Green functions along a path decays over a correlation length that is related to the spectral gap and reads

[TABLE]

In the clean case $W=0$, we recover $\alpha=1/\sqrt{K}$.

For $n=1,2,\cdots$ the right and left eigenfunctions ($\Phi^{(n)}(x)$ and $\Psi^{(n)}(y)$, respectively) can be compactly written as follows:

[TABLE]

and

[TABLE]

as we show in the Appendix. In (51), the normalization constant $C_{n}$ is chosen in such a way that $\int_{-\infty}^{\infty}dx\Phi^{(n)}(x)\Psi^{(n)}(x)=1$.

The typical value of the path susceptibility is readily obtained as:

[TABLE]

which agrees with the expression in [4] relating the typical value in presence and in absence of disorder in the case of a Cauchy distribution of the disorder. Notice that since $\alpha\leq 1/\sqrt{K}$ (the equality being valid only at zero disorder), at $E=0$, $\theta_{\rm typ}>\ln K$ in the deep delocalized phase. In particular, for Cauchy disorder $\theta_{\rm typ}$ grows linearly with $W$,

[TABLE]

We will discuss in the following section how, if $\theta_{\rm typ}$ is connected to the fractal dimension $D_{1}$ of the wave function of a finite tree-like lattice as suggested in Ref. [44], this would imply that the wave function is never ergodic.

3.3 On the fractal dimensions of non-ergodic extended states

In several recent works [42, 41, 44] it has been suggested that the delocalized phase of the Anderson model on the RRG with $N$ sites has two distinct phases: one in which the participation ratios of the eigenfunctions is $\propto N$ and one in which it is $\propto N^{\alpha}$ for some $\alpha<1$. This was questioned in several other papers [48, 49, 50] (however a genuine multi-fractal region was shown to exist for the finite Cayley tree with $N$ vertices [54, 65, 66], and to control the RRG properties below a given critical volume [67]); to the best of our understanding the question is whether this phase exists on the RRG or it is washed away when going from the Cayley tree to the RRG.

This question can not be addressed within the framework discussed in this work. However, we can point out the connection between some arguments presented in Ref. [44] and the quantities investigated here, in particular the largest eigenvalues $\lambda(s)$ of the integral equation (36). We will do it assuming the definition of the fractal dimension $D_{1}$ given in Eq. (76) of Ref. [44]:

[TABLE]

with

[TABLE]

where $\mu(s)=\ln\lambda(s)$ has already appeared in Eq. (20). The extremization procedure with respect to $s$ follows from the interpretation of $\Lambda$ as a free energy of a directed polymer model [42], which undergoes a 1-step RSB transition. The cavity propagators $G_{i}$ in this equation are not independent, but follow the recursion relation (6), and one can find $\lambda(s)$ by solving the integral equation for $K_{s}$.

For $W=0$, the minimum is attained at $s^{*}=1$; since $\mu(1)=\mu(0)=0$ (for any $W$) this gives $D_{1}=1$, which is interpreted as a condition for a fully ergodic state. In [44] it is shown that, when using a slightly modified version of the large disorder forward scattering approximation444The conventional forward scattering approximation consists in neglecting the term $\sum G_{k}(E)$ (which are self-energy corrections) in the denominator of the recursive equation for the cavity Green functions, thus decoupling the consecutive propagators along the path and leading to $G_{i}(E)\sim(E-\epsilon_{i})^{-1}$. In [44] this approximation is performed, but the distribution of the on-site randomness $\epsilon_{i}$ is modified introducing a cut-off around the origin, in order to partially account for the regularization effect produced by the self energies (or, equivalently, of for the anti-correlations between consecutive propagators along the path), in the same spirit of the lower limit approximation in [1], the same remains true for sufficiently small $W\leq W_{erg}\simeq 5.74$ (for $K=2$), meaning that the minimum is still achieved at $s^{*}=1$ for all these values of disorder (and thus $D_{1}=1$).

As $W$ is increased the minimum moves away to $s=s^{*}<1$, and we have

[TABLE]

which gives

[TABLE]

Notice that, due to the symmetry $\mu(s)=\mu(1-s)$ we have, for $s^{*}\simeq 1$ and $W\gtrsim W_{erg}$,

[TABLE]

Now, one can argue that $W_{erg}>0$ only if $\theta_{\rm typ}=\ln K$ at $W=W_{erg}$. This is indeed what happens within the approximation exploited in [44]. However, this result is a feature of that approximation, as we illustrate here. We first notice that, for the Cauchy distribution, $\theta_{\rm typ}=\ln K$ only at $W=0$ (see Eq.(53)) and that for $W>0$ the strict inequality $\theta_{\rm typ}>\ln K$ holds. Second, we numerically check, by using the equality $\theta_{\rm typ}=-\langle\ln G^{2}\rangle$, that the same holds also for the box distribution, although $\partial\theta_{\rm typ}/\partial W=0$ at $W=0$ for box disorder (see Fig. 4). From this we are led to conjecture that (provided that Eq. (55) is used to define the fractal dimension), $D_{1}<1$ for all $W>0$ and for all disorder distributions. We also conjecture that for all disorder distributions which fall faster than $1/\epsilon^{2}$ at large $\epsilon$, $\theta_{\rm typ}=\ln K+cW^{2}$.

Under these hypothesis we now show that

[TABLE]

as numerically found in Ref.[44] (for box disorder we find numerically $\kappa=(5.0\pm 0.6)\cdot 10^{-5}$). In fact using $\lambda(s)=\lambda(1-s)$ one can write, close to $s=1$, and by defining $s=1-\sigma$,

[TABLE]

and considering that $\mu(0,W)=0$, $\mu^{\prime}(0,W)\equiv\frac{\partial\mu(\sigma,W)}{\partial\sigma}|_{\sigma\to 0}\simeq-\ln K-cW^{2}$, and $\mu^{\prime\prime}(0,W)>0$ (numerically we find $c=(4.1\pm 0.2)\cdot 10^{-3}$), for small $W$ one has

[TABLE]

or

[TABLE]

Notice that the constant $c$ can be computed using perturbation theory in $s$ and $W$ around the $W=0,s=0$ point. As the problem for $W=0$ is reduced to the Cauchy distribution, this can be done using the eigenvectors and eigenvalues of the kernel reported in the Appendix [68].

In conclusion, we believe that the quantity $D_{1}$, defined in [44] by using quantities computed on the tree, is always analytic, with only 3 vanishing derivatives at $W=0$. One is led to conclude that the full inclusion of the self-energy corrections in the computation of $\lambda(s)$ breaks the analogy with directed polymers with independent on-site energies, so that the RSB phase present in that directed polymer calculation is not present in the Anderson model at any finite $K$. The fact that $D_{1}$ is very flat at the origin can be considered a vestige of that phase and it is somewhat surprising that a simple computation in the approximation of [44] captures the qualitative feature of a very flat $D_{1}$ at $W=0$.

4 The large disorder, large connectivity limit: an interpolation formula

In this section we give a concise proof of the asymptotic formula $W_{c}=4K\ln K$ for the Anderson transition point $W_{c}$ (this result is specialized for box disorder, the general formula is given in Eq. (77)). This formula is known since the original work of Anderson [1] and has been rigorously proved in Ref. [52], to which we refer for a more comprehensive list of references. Additionally, we provide an interpolating formula between the exact result and the “upper bound” formula $W_{c}=2eK\ln K$ [1] obtained within the forward approximation (FA) [41, 53, 44]. Due to the central importance of the FA for analytic treatments of many-body localization [17, 28, 69, 29] we believe it is important to understand how one could go beyond it in a controlled expansion at large disorder/large connectivity.

To begin with, we notice that the distribution $P(G)$ changes continuously with $W$ for any $K$: what signals the transition is the divergence with $L$ of the typical value of the susceptibility $\chi$ summed over $K^{L}$ paths. Therefore, one can make some assumption on $P(G)$ and use it in the kernel $K_{s}$. In the limit of large connectivity, $K\gg 1$, if we are interested in the transition we have $W=O(K)$ and therefore some simplifications occur. One can self consistently neglect $\sum_{i=1}^{K}G_{i}$ with respect to $\epsilon_{i}$, since the former turns out to be of $O(K/W)$, while the latter is of $O(W)$. This is licit as long as $W\gg\sqrt{K}$. Notice that this was the scale of disorder on which $D_{1}$ decayed from its maximum value $D_{1}=1$ at $W=0$. Numerically, this is also not far from the disorder value at which the maximum of $P(G)$ moves away from $G=0$. Under these hypotheses, the solution for $P(G)$ is

[TABLE]

In the case of Cauchy distribution this functional form is exact and this approximation amounts to replacing $\alpha\to 1/W$:

[TABLE]

In the case of box distribution, Eq. (65) turns out to be a very good approximation already for $K=2,W\gtrsim 10$.

Let us now show that, for $W\gg 1$, $P(G)$ is also the only eigenvector corresponding to a non-zero eigenvalue of the kernel $K_{1/2}$. Following the same reasoning as above, we can neglect the factor $\zeta=\sum_{j=1}^{K-1}G_{j}$ in the definition of the Kernel: this corresponds to replacing $\mathcal{P}_{\zeta}$ in Eq. (40) with a delta function in zero, so that the Kernel simplifies to:

[TABLE]

We then need to solve the integral equation

[TABLE]

where in the last line we have integrated with respect to $\epsilon$.

In order to prove that $P(G)$ is an eigenvector of this integral equation, after factoring out $1/y^{2},$ we need to show that $\mathcal{I}\propto\rho(1/y)$ where

[TABLE]

and the proportionality coefficient is $\lambda(1/2)$. As a first step, it is convenient to pass to the Fourier transform, by defining

[TABLE]

to get

[TABLE]

The integral can be done in terms of Bessel functions:

[TABLE]

where for $kk^{\prime}<0$ we must take only the real part of the Bessel function $Y_{0}$ of complex argument.

If the disorder $W$ is the only scale of the distribution, appearing in the functional form

[TABLE]

then we have

[TABLE]

Assuming that $\tilde{\rho}$ decays naturally on scales of $O(1)$, for large $W$ the integrals over $k,k^{\prime}$ are cut-off at $1/W\ll 1$. We can then expand the Bessel function to find the leading divergence at small $k,k^{\prime}$, which is logarithmic

[TABLE]

The integral over $k,k^{\prime}$ is done by changing variables to $\theta=kW,\theta^{\prime}=k^{\prime}W$, keeping $Wy=O(1)$, while sending $W\to\infty$, considering that $\int\frac{d\theta}{2\pi W}e^{i\theta/Wy}\tilde{\rho}(\theta)=\rho(1/y)$ and $\int\frac{d\theta}{2\pi W}\tilde{\rho}(\theta)=\rho(0)$:

[TABLE]

So, to leading order in $\log(W)$, the eigenvalue is $\lambda(1/2)=4\log(W_{c})\rho(0)$ and the transition point is located at

[TABLE]

For the box distribution this gives the familiar expression

[TABLE]

which is sometimes written in the form (equivalent to leading $O(\log(W))$)

[TABLE]

We can now pause and discuss the physical meaning of the fact that $\phi_{1/2}(x)=P(x)$ for large disorder. At the transition, the eigenvector $\phi_{1/2}$ is the conditioned probability of the propagator $G$, given that one of the $K^{L}$ paths gives a resonance 555An alternative way to interpret this comes from realizing that the eigenvector $\phi_{1/2}(x)$ is proportional to $P(x)$, and the proportionality factor is the expectation value of the paths susceptibility $\chi_{p}$, conditioned to the fact that the first cavity propagator along the path has amplitude $x$. A similar observation for the adjoint kernel is given in [52]. The limit $\phi_{1/2}(x)\to P(x)$ translates into the fact that at the transition the expectation is dominated by instances where the susceptibility is $O(1)$ due to resonances, and for these events correlations between distant propagators do not matter (thus the conditioning is immaterial and its dependence on $x$ drops).. The fact that $\phi_{1/2}=P$, which is the same distribution obtained by the forward approximation, means that one does not have to go to rare events to get a resonance at large $W$. Why then the FA does not give the correct limit for the transition point?

The reason is that the kernel corresponding to the FA, and that of $K_{1/2}$ share the eigenvector corresponding to the largest eigenvalue (in the large $W,K$ limit), but the eigenvalue itself is different (as well as the other eigenvectors and eigenvalues).

Let us clarify this point by introducing the modified kernel

[TABLE]

where $\eta\in[0,1]$. We have $K^{\eta\to 0}=K^{\textrm{FA}}$, corresponding to the forward approximation, since we are not considering the self-energy correction in the denominator. By contrast, $K^{\eta\to 1}=K$ is the full problem we have considered in this article. For any $\eta\neq 1$ the symmetry $\lambda(s)=\lambda(1-s)$ is not respected. So the transition point in the saddle point approximation

[TABLE]

is not achieved at $s^{*}=1/2$, but rather at a given point depending on $W$, see Fig. 5. The simplest case, $\eta=0$ yields

[TABLE]

which has again the eigenvector $\phi^{0}_{s}(x)=\frac{1}{x^{2}}\rho(1/x)$. The eigenvalue is

[TABLE]

Using the form $\rho(\epsilon)=\frac{1}{W}\rho(\epsilon/W)$ we have

[TABLE]

We need to minimize this expression with respect to $s$, and this yields the same result of the forward scattering approximation $\lambda(s^{*})=2e\log(W/2)/W$, which for large $K$ gives the familiar

[TABLE]

In the general case $\eta>0$, we should in principle solve the integral equation. However we can make a simplifying assumption: since $\frac{1}{x^{2}}\rho(1/x)$ is the eigenvector both for $\eta=0$ and $\eta=1$, one can assume that it is the eigenvector for any $\eta\in[0,1]$. In this way the eigenvalue is obtained by:

[TABLE]

For box disorder this reads

[TABLE]

Notice that the symmetry $\lambda(s)=\lambda(1-s)$ holds only for $\eta=1$.

Notice that, while for $\eta=0$ this yields an overestimate of $W_{c}$, for $\eta=1$ it yields an underestimate, which however becomes exact when $W\to\infty$. One can ask then what is the value of $\eta$ corresponding to the numerically obtained value for the smallest $K=2$. From our numerical investigations $W_{c}(K=2)=18.11$, which corresponds to $\eta=0.836$, quite close to $\eta=1$.

As a final result, notice that, for box distribution and $\eta=1$,

[TABLE]

where

[TABLE]

From this, we can recover the full $P(\chi)$ for any $W$. We quote here just the result for $\chi=O(1)$ and $W=W_{c}$, where we find the distribution of the single path susceptibility $\chi$,

[TABLE]

which is a log-normal distribution, in agreement with the general discussion in Sec. 2.4. The typical values of this log-normal distribution are of order $\ln\chi\sim\sqrt{L}/\ln(W/2)$.

In the general case

[TABLE]

and therefore the integral representation is

[TABLE]

where $\tilde{\theta}=\ln(W/2)\ln\chi/L$. A real saddle point exists only for $-2<\tilde{\theta}<2$, so the main range of variability666For $\tilde{\theta}$ outside this range there is a pair of complex conjugate saddle points and the resulting integral is much smaller. So we can speak of “main range” for $\chi$. of the susceptibility is $\chi\in[e^{-\frac{2L}{\ln(W/2)}},e^{+\frac{2L}{\ln(W/2)}}].$

5 Summary and perspectives

In this article we have revisited the Anderson transition on the Bethe lattice in terms of real-energy propagators. We have first determined the distribution of these propagators through population dynamics and compared it with the one resulting from exact numerical methods on Random regular Graphs, finding agreement. We have formulated a novel criterion for the transition in terms of the distribution of the susceptibilities, writing an integral equation for the quantity $\lambda(s)$ determining the Mellin transform of the distribution; we have shown that this criterion is consistent with the previously known one, given that this equation coincides with the one originally introduced in [8] to determine the stability of the Anderson phase.

In the case of Cauchy disorder, we have determined analytically the full set of eigenvalues and eigenvectors of the integral kernel at $s=0$, and we have shown explicitly that the gap between the largest eigenvalues of the integral kernel is always $>0$, implying that the cavity propagators have exponentially decaying correlations. For box disorder we have derived an accurate estimate of the transition point for lattices with connectivity $k=2$, and we have shown analytically that the fractal dimension $D_{1}$ defined as in Ref. [44] satisfies $D_{1}=1-\kappa W^{4}<1$, as was also obtained numerically in that paper.

Finally, we have revisited the derivation of the critical disorder in the limit of large disorder (equivalently, of large connectivity of the lattice), showing explicitly that the forward approximation fails to correctly capture the spectrum of the integral operator because it neglects relevant correlations between consecutive propagators along a path. We notice that for single particle problems, the asymptotic formula for $W_{c}$ has been shown to be in very good agreement with the numerical results, for dimensions $d\gtrsim 4$ [7]; therefore, it would be desirable to extend our analytic progress on the integral equation to the analytic treatment of MBL.

A natural continuation of this work would be the systematic derivation of a perturbative expansion for the eigenvalue $\lambda(s)$ of the integral kernel $K_{s}$, given by

[TABLE]

This expansion is convergent because the spectrum of $K_{s=0}$ is gapped. In the case of a Cauchy distribution of disorder, this can be worked out using the explicit expressions for the eigenstates and for the spectrum of the operator $K_{s=0}$ reported in this work. In particular, in the large disorder limit it should be possible to identify which contributions in the expansion have to be retained to recover the asymptotic correct behavior at $s=1/2$; this should give indications on how to properly account for the self-energy contributions in more conventional path expansions in the localized phase, in both tree-like and finite dimensional lattices. We leave these calculations for future work [68].

Acknowledgments

The authors would like to thank V. Kravtsov, L.Ioffe, A.Mirlin and K.Tikhonov for interesting discussions. FP thanks G. Lemarié for discussions and acknowledges the support of the project THERMOLOC ANR-16-CE30-0023-02 of the French National Research Agency (ANR). VR acknowledges the support of the Simons Foundation collaboration Cracking the Glass Problem (No. 454935 to G. Biroli). SP is partially supported by Istituto Nazionale di Fisica Nucleare (INFN) through the project “QUANTUM”. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement 694925).

Appendix A Eigenfunctions of the Cauchy Kernel at s=0

In this Appendix we derive the expressions for the right and left eigenfunctions of the Kernel Eq. (47), given in Eqs. (51), (52) in the main text. For the purpose of the derivation, it is convenient to consider the Kernel in Fourier space, which reads:

[TABLE]

where we assume $\theta(0)=0$, $\theta(x)\delta(x)=0$, and where $J_{1}(\cdot)$ is a Bessel function of the first kind. We denote with $\tilde{\Psi}^{(n)}(k^{\prime})$ and $\tilde{\Phi}^{(n)}(k)$ the left and right eigenvectors of the Kernel, respectively, with $n=0,1,...$ and with eigenvalues $\lambda_{n}=(-1)^{n}\alpha^{2n}$. We aim at showing that the left eigenvectors (up to normalization) are equal to:

[TABLE]

where

[TABLE]

$L_{m}^{(1)}(\cdot)$ are generalized Laguerre polynomial, and

[TABLE]

For $n=0$, the eigenvalue equation is trivially satisfied by $\delta(k^{\prime})$. For $n=1,2,\cdots$, the eigenvalue equation reads:

[TABLE]

where

[TABLE]

while

[TABLE]

and

[TABLE]

To compute this integral, we exploit the following identity, holding for any real $t$:

[TABLE]

Choosing $t=\alpha|k|$, we get

[TABLE]

Using the orthogonality condition for the generalized Laguerre polynomials:

[TABLE]

we have

[TABLE]

Now,

[TABLE]

and thus

[TABLE]

This implies that

[TABLE]

where

[TABLE]

so that

[TABLE]

meaning that the eigenvalue equation is satisfied, given that $\alpha+\overline{\alpha}=\alpha^{-1}$.

Following the same procedure, one can show that the right eigenvectors read (up to normalization):

[TABLE]

with the same coefficients $c_{m}^{(n)}$ defined above.

The expressions in Eqs. (51), (52) are obtained by performing the inverse Fourier transform. For the left eigenvectors we get:

[TABLE]

We now exploit the generating function of the generalized Laguerre polynomials:

[TABLE]

so that

[TABLE]

to write

[TABLE]

and thus

[TABLE]

which is Eq. (52) in the main text, up to a normalization factor.

For the right eigenvector, we need to evaluate the following integral:

[TABLE]

Exploiting once more the generating function (114), we get

[TABLE]

and so

[TABLE]

which is Eq. (51) in the main text, up to a normalization factor.