A random matrix model with localization and ergodic transitions

Motivated by the problem of many-body (MB) Localization [1] and the applicability of the Boltzmann's statistics in interacting disordered media [2], there was recently a revival of interest to the Anderson localization (AL) problem on hierarchical lattices such as the Bethe lattice (BL) or the random regular graph (RRG). Due to hierarchical structure of the Fock space connected by the two-body interaction, statistics of random wave functions in such models is an important playground for MB localization. In particular, the non-ergodic extended phase on disordered hierarchical lattices could model a breakdown of conventional Boltzmann statistics in interacting MB systems and an emergence of a phase of a 'bad metal' [3] or unconventional fluid phases [4] in systems of interacting particles.

However, even for the one-particle AL existence of such a phase in a finite interval of disorder strengths is a highly non-trivial issue.

According to earlier studies [5, 6] there is only one transition in such models at a disorder strength $W={W}_{\mathrm{AT}}$ which is the AL transition (AT) that separates the localized and ergodic extended states. However, recent numerical studies [7] of level statistics on RRG seem to indicate on the second transition at $W={W}_{\mathrm{ET}}\lt {W}_{\mathrm{AT}}$ which is identified as the transition between the ergodic and non-ergodic extended states (ET). Subsequent studies [8, 9] raise doubts about the existence of the second transition on RRG. Numerical results of [8, 9] indicate on the non-ergodic states on RRG in a wide range of disorder strengths down to very low disorder $W=5\ll {W}_{\mathrm{AT}}\approx 17.5$, while in [9] it is demonstrated how an apparent non-ergodic behavior for the intermediate matrix sizes N in Levy random matrix (RM) ensemble evolves into the ergodic one at larger N's. Complexity of RRG and the controversy associated with existence of the ergodic transition at $W={W}_{\mathrm{ET}}$ necessitate a search for a simpler model in which such a transition may occur.

Inspired by the success of Wigner–Dyson RM theory [10] which predictions are relevant in such seemingly different fields of physics as nuclear physics and nano- and mesoscipic physics, our goal is to search for a RM model that would be able to give a simple and universal description of all the three phases: good metal, MB insulator and 'bad metal', which are relevant in the problem of MB localization. An important heuristic argument to construct such a model is that RRG with disordered on-site energies ${\varepsilon }_{i}$ is essentially a two-step disorder ensemble. The disorder of the first level is the structural disorder due to the random structure of RRG where each of N sites of the graph is connected with the fixed number $K+1$ of other sites in a random manner. An ensemble of tight-binding models on such graphs with deterministic on-site energies ${\varepsilon }_{i}$ and hopping integrals is believed to be equivalent to the Gaussian RM ensemble [11]. The disorder of the second level is produced by randomization of ${\varepsilon }_{i}$ fluctuating independently around zero with the distribution function $p(\varepsilon )$. For numerical calculations this distribution is often taken in the form $p(\varepsilon )=\theta \left(W/2-| \varepsilon | \right)/W$, with $\theta (x)$ being the Heaviside step function.

One can expect that the following RM ensemble (the Rosenzweig–Porter (RP) ensemble [12]) is a close relative of RRG with on-site energy ${\varepsilon }_{i}$ disorder. It is an ensemble of N × N random Hermitian matrices which entries H_nm with $n\gt m$ are independent random Gaussian numbers, real for orthogonal RP model ($\beta =1$) and complex for the unitary RP model ($\beta =2$), fluctuating around zero with the variance $\langle | {H}_{{nm}}{| }^{2}\rangle =(\beta /2)\sigma$. The diagonal elements have the same properties with the variance $\langle {H}_{{nn}}^{2}\rangle =1$. The case $\sigma =1$ corresponds to the Gaussian orthogonal (GOE) or Gaussian unitary (GUE) ensembles and represents the structural disorder in RRG. The additional ${\varepsilon }_{i}$-disorder in RRG corresponds to $\sigma \lt 1$. One should also take into account that in order to significantly deviate from the GOE or GUE behavior, the ratio $\langle | {H}_{{nm}}{| }^{2}\rangle /\langle {H}_{{nn}}^{2}\rangle$ must be proportional to some negative power of the matrix size N, as the number of the off-diagonal terms is $\sim N$ times larger. Thus we consider the model:

$$\begin{eqnarray}&&\langle {H}_{{nn}}^{2}\rangle =1,\ \ \ \ \langle | {H}_{n\ne m}{| }^{2}\rangle =(\beta /2)\sigma ={\lambda }^{2}/{N}^{\gamma },\end{eqnarray} \tag{ 1 }$$

where $\lambda \sim O({N}^{0})$ and γ is the main control parameter of the problem.

One can estimate the strength of disorder required for the AL transition as corresponding to the typical fluctuation of diagonal matrix element equal to the typical off-diagonal matrix element times the coordination number K. For the coordination number $K\sim N$ (each site is connected with any other one) this results in $\sqrt{\sigma }N\sim 1$, or $\sigma \sim 1/{N}^{2}$. However, this estimation does not take into account a random, sign-alternating character of the off-diagonal matrix elements. It is likely that for sign-alternating hopping there is another relevant coordination number $\sim \sqrt{N}$ with the critical scaling $\sigma \sim 1/N$. As we show below it corresponds to the ergodic transition. For technical reasons the most significant progress in the analytical studies of the model was achieved [13, 14, 16] for the 'unitary' RP (URP) ensemble. The conclusion was that at $\gamma =2$ the spectral form-factor (two-level correlation function) is neither of the Wigner–Dyson nor of the Poisson form [13, 15, 16] which is typical for the AL transition point. In contrast, at $\gamma =1$ the level statistics was found to be GUE [14]. The papers [13, 14, 16] have a status of classic keynote papers in the field.

The Dyson ideas of the Brownian motion of energy levels first applied to the RP ensemble in [13] were developed in the series of works [17, 18]. It was shown that the possible transitions in the level statistics are associated with the fixed points of parameter ${\rm{\Lambda }}=\sigma (N)/{[\delta (N)]}^{2}\propto {N}^{-\gamma }/{[\delta (N)]}^{2}$, where $\delta (N)$ is the mean level spacing. Then assuming $\delta (N)\propto 1/N$ established in [16] for $\gamma \gt 1$ one obtains the transition point at $\gamma =2$. If, however, the Wigner–Dyson semicircle level density is assumed with $\delta (N)\propto 1/\sqrt{N}$, then the transition would occur at $\gamma =1$ [17]. Unfortunately, $\delta (N)\propto 1/\sqrt{N}$ only at $\gamma =0$. For $0\lt \gamma \leqslant 1$ the following result is valid (see e.g. equation (170) in [19]) for the mean density of states $\rho (\varepsilon )=\sqrt{2S-{\varepsilon }^{2}}/(\pi S)$, where $S={\displaystyle \sum }_{n}\langle | {H}_{{nm}}{| }^{2}\rangle \propto {N}^{1-\gamma }$. Thus we obtain $\delta (N)=1/(N\rho (0))\propto {N}^{-(1+\gamma )/2}$ for $0\lt \gamma \leqslant 1$ and $\delta (N)\propto 1/N$ for $\gamma \gt 1$, resulting in ${\rm{\Lambda }}\propto N$ for $0\lt \gamma \leqslant 1$ and ${\rm{\Lambda }}\propto {N}^{2-\gamma }$ for $\gamma \gt 1$. We conclude that the only fixed point of Λ is possible at $\gamma =2$, and no transition at $\gamma =1$ can be obtained from the results of [17, 18].

In this paper by a more sophisticated analysis of the two-level correlations and the eigenfunctions statistics we show that the above extension of the RP model indeed contains not one but two transitions. One of them at $\gamma =2$ corresponds to the transition from the extended to the localized states. However, the extended states emerging at $\gamma \lt 2$ are not ergodic: their support set contains infinitely many ${N}^{{D}_{1}}$ sites in the $N\to \infty$ limit, which, however, is a zero fraction of all sites, since ${D}_{1}\lt 1$. Such non-ergodic extended states on RRG are recently discussed in [8]. With further decrease of γ the second transition at $\gamma =1$ happens which is a transition from the non-ergodic extended states to the ergodic extended states with ${D}_{1}=1$ similar to the eigenstates of the GOE.

We prove this statement in three steps. As the first step we use the perturbative arguments to compute the statistics of wave function amplitude $| \psi ({r}_{{\rm{o}}}){| }^{2}$ in a certain observation point r_o. We obtain a drastic change of the character of this distribution at $\gamma =1$ and $\gamma =2$ which is summarized in figure 1. This result is fully confirmed by a numerical diagonalization of the Hamiltonian (see figures 2–4). It is also confirmed by the numerical analysis of the moments of random wave functions which determine their Shannon entropy and the support set dimension D₁ (see figure 5). Then we compute numerically the overlap of amplitudes for two different wave functions with the energy difference ω and find the scaling with N of the Thouless energy ${E}_{\mathrm{Th}}\sim {N}^{-z}$ which exponent z changes abruptly at $\gamma =1$ and $\gamma =2$ (see figure 6). Finally, we perform a rigorous calculation of the spectral form-factor which also shows the transition at $\gamma =1$ and $\gamma =2$ (see figure 7). In the last section we compare the corresponding results for our model and for the RRG and demonstrate their similarity. It allows us to unify both models in a special universality class of random hierarchical models which differs from the one realized in localization transition points of two- and three-dimensional Anderson models. Further details concerning this model can be found in [20].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As the off-diagonal matrix elements in equation (1) are small, one can employ the perturbation theory for computing the distribution function of the amplitudes $x=N| \psi ({r}_{{\rm{o}}}){| }^{2}$. The first order perturbation theory gives:

$$\begin{eqnarray}&&{| {\psi }_{n}({r}_{m})| }^{2}=| {H}_{{nm}}{| }^{2}/{({H}_{{nn}}-{H}_{{mm}})}^{2},\ \ \ \ \ (n\ne m),\end{eqnarray} \tag{ 2 }$$

where the maximum of ${\psi }_{n}(r)$ is supposed to be at $r={r}_{n}$.

The perturbative series converge absolutely if the typical off-diagonal matrix element $| {({H}_{{nm}})}_{\mathrm{typ}}| \sim \lambda {N}^{-\gamma /2}$ times the coordination number N is much smaller than the typical difference of the diagonal matrix elements $| {({H}_{{nn}}-{H}_{{mm}})}_{\mathrm{typ}}| \sim W\gg \delta (N)$. Thus it converges absolutely for $\gamma \gt 2$ irrespectively of the statistics of diagonal matrix elements. For $\gamma \leqslant 2$ the convergence of the series occurs only because of the random and independently fluctuating signs of H_nm and $({H}_{{nn}}-{H}_{{mm}})$. Although it is hard to prove such a convergence rigorously, a plausible argument in its favor is that the effective coordination number of oscillatory contributions is $\sqrt{N}$ rather than N. The corresponding criterion of convergence is $\lambda \;{N}^{-\gamma /2}\;\sqrt{N}\ll W$ which is satisfied at $\gamma \gt 1$.

Consider the regular part of the characteristic function $Q(\xi )=\langle {{\rm{e}}}^{{\rm{i}}\xi \;N| \psi ({r}_{{\rm{o}}}){| }^{2}}\rangle$. For the Gaussian distribution of matrix elements equation (1) we obtain $Q(\xi )={\mathcal{Q}}(\xi N\sigma )$, where:

$$\begin{eqnarray}&&{\mathcal{Q}}(\zeta )={{\rm{e}}}^{-{\rm{i}}\zeta /2}\;{\rm{Erfc}}\left(\sqrt{-{\rm{i}}\zeta /2}\right)\approx 1-\sqrt{-2{\rm{i}}\zeta /\pi }.\end{eqnarray} \tag{ 3 }$$

The eigenfunction amplitude distribution function $P(x)={\displaystyle \int }_{-\infty }^{\infty }{{\rm{e}}}^{-{\rm{i}}\zeta \;\frac{x}{N\sigma }}\;{\mathcal{Q}}(\zeta )\displaystyle \frac{{\rm{d}}\zeta }{2\pi N\sigma }$ at $x\gg N\sigma \sim O({N}^{1-\gamma })$ is dominated at $\gamma \gt 1$ by small $\zeta \ll 1$. That is why it is only the expansion of equation (3) at small ζ what matters for P(x) at $\gamma \gt 1$. Thus we obtain for the regular part of the function P(x):

$$\begin{eqnarray}&&{\mathcal{P}}(x)={(\sqrt{2}\pi )}^{-1}\;{(N\sigma )}^{1/2}/{x}^{3/2}.\end{eqnarray} \tag{ 4 }$$

There are two normalization conditions for P(x): the normalization of probability equation (5) and the normalization of the wave function equation (6):

$$\begin{eqnarray}&&{\displaystyle \int }_{0}^{\infty }P(x){\rm{d}}x=1,\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{\displaystyle \int }_{0}^{\infty }x\;P(x){\rm{d}}x=1.\end{eqnarray} \tag{ 6 }$$

Equation (5) imposes a cut-off ${x}_{\mathrm{min}}\sim {N}^{-(\gamma -1)}$ to equation (4) at small x, while equation (6) determines the upper cut-off ${x}_{\mathrm{max}}\sim {N}^{\gamma -1}$. A caution, however, should be taken: by normalization ${\displaystyle \sum }_{i}| \psi ({r}_{i}){| }^{2}=1$ the amplitude $| \psi ({r}_{i}){| }^{2}\leqslant 1$ on any lattice site cannot exceed 1, and therefore $x\leqslant N$. One can see that the above estimation for x_max is valid only for $\gamma \lt 2$ when ${N}^{\gamma -1}\ll N$. For $\gamma \gt 2$ a correct ${x}_{\mathrm{max}}=N$. In order to compensate for the deficiency of normalization in equation (6) one has to assume a singular part of $P(x)={\mathcal{P}}(x)+A\;\delta (x-N)$. One can see that for $\gamma \gt 2$ equation (6) is dominated by the singular term, and $A={N}^{-1}$. This corresponds to the strongly localized wave functions. The mechanism of emergence of the singular term at the AL transition at $\gamma =2$ is somewhat similar to the Bose-condensation, where the singular term also appears because of the deficiency of normalization of the Bose–Einstein distribution.

One can express the distribution function equation (4) through the spectrum of fractal dimensions [8, 21]:

$$\begin{eqnarray}&&f(\alpha )=\underset{N\to \infty }{\mathrm{lim}}f(\alpha ,N)=\underset{N\to \infty }{\mathrm{lim}}\mathrm{ln}[{xN}{\mathcal{P}}(x)]/\mathrm{ln}N,\end{eqnarray} \tag{ 7 }$$

where $\alpha =1-\mathrm{ln}x/\mathrm{ln}N$ or $| \psi ({r}_{{\rm{o}}}){| }^{2}={N}^{-\alpha }$. Using equations (4)–(7) one obtains:

$$\begin{eqnarray}&&f(\alpha )=\alpha /2+1-\gamma /2,\ \ ({\alpha }_{\mathrm{min}}\lt \alpha \lt {\alpha }_{\mathrm{max}}).\end{eqnarray} \tag{ 8 }$$

The upper cutoff ${\alpha }_{\mathrm{max}}=\gamma$ corresponds to the lower cutoff x_min. The lower cutoff ${\alpha }_{\mathrm{min}}$ depends on γ. In the localized region $\gamma \gt 2$, figure 1(a), ${\alpha }_{\mathrm{min}}=0$. At the AL transition point $\gamma =2$ the function $f(\alpha )$ has the same triangular shape as at $W={W}_{\mathrm{AT}}$ on RRG, figure 1(b). In the region of the extended non-ergodic states $1\lt \gamma \lt 2$, figure 1(c), ${\alpha }_{\mathrm{min}}=2-\gamma \gt 0$. It is remarkable that in the entire region $1\leqslant \gamma \leqslant 2$ the symmetry [21, 22]$f(1+x)=f(1-x)+x$ holds. Finally, at $\gamma =1$ the two limits ${\alpha }_{\mathrm{min}}$ and ${\alpha }_{\mathrm{max}}$ collapse in one point $\alpha =1$ which marks the transition point to the ergodic state, figure 1(d).

Note that $f(\alpha )$ for $\gamma \gt 2$ (see figure 1(a)) has a singular peak at $\alpha =0$ which corresponds to the singular term ${N}^{-1}\delta (x-N)$ in P(x). This singular $f(\alpha )$ is not a limit of any convex function. However, one may easily see that all the moments $N\langle | \psi {| }^{2q}\rangle \sim {N}^{-{\tau }_{q}}$ have the same exponents ${\tau }_{q}$ as for the 'convex' $f(\alpha )$ shown by the dashed line in figure 1(a): $f(\alpha )=\alpha /\gamma$ for $0\lt \alpha \lt \gamma$. Such a triangular $f(\alpha )$ with the slope smaller than 1/2 also holds in the localized phase on RRG [8].

The numerical calculation of $f(\alpha )$ which involve the rectification and extrapolation procedures described in [8], fully confirms the above results. In figure 2 we present the results of this calculation for N = 2⁸–2¹⁴ and the extrapolated $f(\alpha )$ (shown by a solid red line) for $\gamma =3$ which perfectly coincides with the prediction of our perturbative analysis above. The similar coincidence was obtained for $\gamma =1.5$, while for $\gamma =0.75$ the distribution function P(x) is practically indistinguishable from the Porter–Thomas distribution of the GOE.

By calculating the Legendre transform ${\tau }_{q}$ [21] of $f(\alpha )$ (shown in figure 1) one finds that in the intermediate phase $1\lt \gamma \lt 2$ all fractal dimensions ${D}_{q}={\tau }_{q}/(q-1)$ for $q\gt 1/2$ are the same and equal to:

$$\begin{eqnarray}&&{D}_{q}=D=2-\gamma ,\ \ \ \ (q\gt 1/2).\end{eqnarray} \tag{ 9 }$$

Thus the support set of a typical wave function is a fractal containing ${N}^{{D}_{1}}={N}^{2-\gamma }$ sites. As ${N}^{2-\gamma }\to \infty$ in the limit $N\to \infty$ it is an extended state. However the support set contains a fraction of all sites $F={N}^{1-\gamma }$ tending to zero in this limit. Thus it is a non-ergodic state.

In the localized phase $\gamma \gt 2$ (including the critical point $\gamma =2$) we obtained:

$$\begin{eqnarray}\tau (q)=\left\{\begin{array}{cc}\gamma q-1, & q\lt 1/\gamma \\ 0, & q\gt 1/\gamma \\ \end{array}.\right.\end{eqnarray} \tag{ 10 }$$

One can see that the fractal dimensions ${D}_{q}=0$ only at $q\gt 1/\gamma$, while they are non-zero and negative for $0\lt q\lt 1/\gamma$. This is not exactly the behavior of the typical Anderson insulator where all fractal dimensions with $q\gt 0$ are equal to zero. The behavior similar to equation (10) were found in certain two-dimensional random Dirac models [25–28]. Such a quasi-localized phase is referred to as the frozen phase and the corresponding transition is known as the freezing transition. In such a phase a typical wave function amplitude has several sharp peaks separated by valleys where $| \psi {| }^{2}$ is not exponentially but only power-law small in N ($| \psi {| }_{\mathrm{typ}}^{2}\gt {N}^{-\gamma }$ in our case). The same behavior is also found for the RRG [8] with $W\gt {W}_{\mathrm{AT}}$.

In order to check the existence of the intermediate phase numerically we computed the average $\langle x\mathrm{ln}x\rangle$ which is directly related with the Shannon entropy and the dimension D₁ of the support set of fractal wave functions [23]. The results are shown in the inset of figure 5 where N span from 256 up to 32768. The corresponding values of D₁ extracted from the linear in $\mathrm{ln}N$ fit are shown in figure 5 which are consistent with the transitions at $\gamma =2$ and $\gamma =1$. The deviation from the expected ${D}_{1}=2-\gamma$ shown by a dashed line in figure 5 is a finite-size effect. Indeed, the correlation volume N_c close to the localization transition at $\gamma =2$ is exponentially large ${N}_{{\rm{c}}}\propto {e}^{c/| 2-\gamma | }$. This follows from equation (17) of section 5 where the Poisson limit is reached only for ${N}^{\gamma -2}\gg O(1)$ or $\mathrm{ln}N\gg \mathrm{ln}{N}_{{\rm{c}}}\sim 1/(\gamma -2)$. The similar exponential dependence ${N}_{{\rm{c}}}\sim {e}^{c/\sqrt{| W-{W}_{{\rm{c}}}| }}$ of the correlation volume was obtained on the BL [24]. For system sizes $N\ll {N}_{{\rm{c}}}$ one should see the properties of the critical point $\gamma =2$ where ${D}_{1}=0$. Thus in the vicinity of the transition point the support set dimension extracted from the finite-size simulations should show a tendency towards smaller values as in figure 5. However, for $\gamma \lt 1.5$ at our systems sizes the support set dimension D₁ approaches the values expected in $N\to \infty$ limit (dashed line in figure 5). This fact implies that for $\gamma \lt 1.5$ we reached $N\gg {N}_{{\rm{c}}}$ and thus it may serve as numerical evidence of convergence and existence of non-ergodic extended phase in the thermodynamic limit.

We also introduced the $1/N$ corrections to the fit with its magnitude ${C}_{1/N}$ being a measure of the global curvature of the $\langle x\mathrm{ln}x\rangle$ versus $\mathrm{ln}N$ dependence. Remarkably, ${C}_{1/N}$ changes sign at both the transition points $\gamma =1$ and $\gamma =2$ (though the positive ${C}_{1/N}$ is very small for $\gamma \gt 2$). We also checked that it changes sign at the localization transition point of the 3D Anderson model (not shown). We believe that the changing of sign of ${C}_{1/N}$ is a convenient way to identify the points of both localization and ergodic transitions.

Next we compute numerically the overlap of different wave functions $K(\omega )=N{\displaystyle \sum }_{r}\langle | {\psi }_{\varepsilon }(r){| }^{2}| {\psi }_{\varepsilon +\omega }(r){| }^{2}\rangle$. Note that for the ergodic wave functions of GOE $K(\omega )=1$ is independent of the energy difference ω, as in this case the overlap is always 100%. Our results presented in figure 6 show that for $\gamma \gt 1$ the overlap $K(\omega )$ has a plateau at $\omega \lt {E}_{\mathrm{Th}}$ which is followed by a fast decrease $K(\omega )\propto 1/{\omega }^{2}$ for $\omega \gt {E}_{\mathrm{Th}}$. The Thouless energy E_Th [29] that separates the GOE-like behavior (plateau) from the system specific behavior ($K(\omega )\propto {\omega }^{-2}$), depends on N as a power-law ${N}^{-z}$. However, the scaling exponent z is different in all the three phases (see figure 6). In the localized phase figure 6(b) we obtained $z=\gamma /2$ which corresponds to rare resonances when $\omega \lt | {H}_{n\ne m}| \sim {N}^{-\gamma /2}$. In the extended non-ergodic phase figure 6(a) we found ${E}_{\mathrm{Th}}\sim {N}^{1-\gamma }\sim \delta {N}^{D}$, where $\delta =1/(N\;p(0))$ is the mean level spacing. This corresponds to all N^D sites in the support set being in resonance with each other. The corresponding N^D states produced by linear combinations of basis states localized on resonant sites form a mini-band of levels of the width ${E}_{\mathrm{Th}}\sim \delta \;{N}^{D}$. Clearly, the states inside such a mini-band should have the GOE-like correlations. On the contrary, the states separated by the energy distance $\omega \gt {E}_{\mathrm{Th}}$ should belong to different support sets which poorly overlap with each other. At the ergodic transition at $\gamma =1$ and in the entire extended ergodic state at $\gamma \lt 1$ we obtain ${E}_{\mathrm{Th}}\sim O({N}^{0})$, and the plateau extends to entire spectral band-width. The emergence of such a plateau that survives the limit $N\to \infty$ is a signature of the ergodic state [30].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Surprisingly, the overlap function $K(\omega )\sim {N}^{1-\gamma }/{\omega }^{2}\to 0$ as $N\to \infty$ at any fixed ω and $\gamma \gt 1$. This phenomenon of 'repulsion of wave functions' [30] is a peculiar feature of our model. The non-ergodic fractal states at the localization transition points of the two and three- dimensional Anderson models of the Dyson symmetry classes, as well as those of the power-law banded random matrices [21, 30, 31] show a different behavior. In these models, the Thouless energy for fractal states is proportional to the mean level spacing ${E}_{\mathrm{Th}}\sim \delta$ and the behavior for $\omega \gt {E}_{\mathrm{Th}}(N)$ is described by the conventional Chalker's scaling [32, 34, 35]:

$$\begin{eqnarray}&&K(\omega )\sim {\omega }^{-1+{D}_{2}}.\end{eqnarray} \tag{ 11 }$$

Only at a very large energy separations ω of the order of the total spectral band-width, the 'repulsion of wave functions' was observed [30].

A remarkable feature of the present model is that the Thouless energy in the region of extended non-ergodic states is much larger than the mean level spacing:

$$\begin{eqnarray}&&{E}_{\mathrm{Th}}\sim \delta \;{N}^{{D}_{2}}\propto {N}^{-z}.\end{eqnarray} \tag{ 12 }$$

One can interpret this relationship as a non-trivial dynamical scaling exponent

$$\begin{eqnarray}&&z=1-{D}_{2}\lt 1.\end{eqnarray} \tag{ 13 }$$

For non-interacting systems in two or three-dimensions in the point of Anderson transition z = 1 for all Dyson symmetry classes. A non-trivial z is known only in two-dimensional systems described by the Dirac equation with random vector-potential which belong to chiral symmetry classes [25–28] where the freezing transition is observed.

In terms of the dynamical exponent z the leading power-law term in the Chalker's scaling for $\omega \gg {E}_{\mathrm{Th}}$ can be rewritten as [28]:

$$\begin{eqnarray}&&K(\omega )\propto {\omega }^{-\mu },\ \ \ \ \ \mu =(1-{D}_{2})/z.\end{eqnarray} \tag{ 14 }$$

In our model we have:

$$\begin{eqnarray}&&K(\omega )\sim \displaystyle \frac{{E}_{\mathrm{Th}}}{{\omega }^{2}},\ \ \ \ (\omega \gg {E}_{\mathrm{Th}}).\end{eqnarray} \tag{ 15 }$$

One can consider equation (15) as a particular case of expansion in ${E}_{\mathrm{Th}}/\omega \ll 1$ with the leading term exponent $\mu =1$:

$$\begin{eqnarray}&&K(\omega )=\displaystyle \frac{1}{\omega }\left[{c}_{0}+{c}_{1}\;\displaystyle \frac{{E}_{\mathrm{Th}}}{\omega }+...\;\right],\end{eqnarray} \tag{ 16 }$$

in which the coefficient c₀ is zero. We will see below that equations (12) and (16) hold for the RRG too. However in this case the coefficient ${c}_{0}\sim 1$ is non-zero. Thus one can speak of the special universality class of models with a non-trivial $z\ne 1$ and $\mu =1$ which the present model belongs to together with the RRG model.

Finally, we present the results of a rigorous calculation of the spectral form-factor $C(t,{t}^{\prime })={\displaystyle \sum }_{n\ne m}{{\rm{e}}}^{{\rm{i}}{{tE}}_{m}+{\rm{i}}{t}^{\prime }{E}_{n}}$, with a set of eigenvalues $\{{E}_{n}\}$ of H. To this end we generalize to the case $\gamma \gt 1$ the results of [16] where $C(t,{t}^{\prime })$ was derived for URP model equation (1) with $\gamma =2$ using the Itzykson–Zuber formula of integration over unitary group. The final result (see details of the derivation in [20]) for $C(t,{t}^{\prime })=2\pi \delta (t+{t}^{\prime })\left[S\left(\frac{{E}_{\mathrm{Th}}}{2}(t-{t}^{\prime })\right)-1\right]$ is given by equation (17)

$$\begin{eqnarray}S(u) & = & 1+{{\rm{e}}}^{-2\pi {{\rm{\Lambda }}}^{2}u}{{\rm{e}}}^{-{{\rm{\Lambda }}}^{2}{u}^{2}{N}^{\gamma -2}}\left[\displaystyle \frac{2{I}_{1}(\kappa {u}^{3/2})}{\kappa {u}^{3/2}}\right.\\ & & -\left.\displaystyle \frac{1}{4\pi }\kappa {u}^{5/2}\;{N}^{\gamma -2}\;{\displaystyle \int }_{0}^{\infty }\displaystyle \frac{x\;{\rm{d}}x}{\sqrt{x+1}}\;{I}_{1}(\kappa {u}^{3/2}\sqrt{x+1}){{\rm{e}}}^{-x\;{u}^{2}{{\rm{\Lambda }}}^{2}{N}^{\gamma -2}}\right],\end{eqnarray} \tag{ 17 }$$

with the modified Bessel function ${I}_{1}(x)$, ${E}_{\mathrm{Th}}=\delta \;{N}^{2-\gamma }$, $\kappa =\sqrt{8\pi {N}^{\gamma -2}}{{\rm{\Lambda }}}^{2}$ and ${\rm{\Lambda }}=\lambda p(0)$. The unfolded spectral form-factor $R({\mathfrak{t}})$ with ${\mathfrak{t}}=\frac{1}{2}(t-{t}^{\prime })\delta$ is given by $R({\mathfrak{t}})=S({\mathfrak{t}}{N}^{2-\gamma })$. It follows from equation (17) that for $\gamma \gt 2$, $R({\mathfrak{t}})\to 1$ in the $N\to \infty$ limit, which corresponds to completely uncorrelated energy levels and the exact Poisson statistics. Another important feature of equation (17) is that $R(0)=1$ for all $\gamma \gt 1$.

In figure 7 we plot the unfolded spectral form-factor $R({\mathfrak{t}})$ for the two phases: (i) the critical phase of the AL transition at $\gamma =2$ and (ii) the intermediate phase at $1\lt \gamma \lt 2$. One can see that while at $\gamma =2$ the function $R({\mathfrak{t}})$ has a non-trivial $N\to \infty$ limit, for $1\lt \gamma \lt 2$ the limit coincides with that of the GUE, except for the point ${\mathfrak{t}}=0$ where there is a jump in $R({\mathfrak{t}})$. This jump is a hallmark of the intermediate phase. To demonstrate this more clearly we blow up the region of small ${\mathfrak{t}}$ by re-scaling the variable ${\mathfrak{t}}\Rightarrow u$ $=\displaystyle \frac{1}{2}(t-{t}^{\prime })\delta \;{N}^{2-\gamma }=\displaystyle \frac{1}{2}(t-{t}^{\prime }){E}_{\mathrm{Th}}$. Note that ${E}_{\mathrm{Th}}=\delta {N}^{2-\gamma }$ appears again as the characteristic scale where the level repulsion is taken over by the level attraction in figure 7. In this new variable $R({{uN}}^{\gamma -2})=S(u)$ has a non-GUE $N\to \infty$ limit:

$$\begin{eqnarray}&&S(u)=\mathrm{exp}(-2\pi {{\rm{\Lambda }}}^{2}\;u),\ \ \ \ 1\lt \gamma \lt 2,\end{eqnarray} \tag{ 18 }$$

which is shown in the inset to figure 7. The true GUE form factor is just identically zero in this limit. The existence of the new scale E_Th and a non-GUE $N\to \infty$ limit equation (18) in the variable $u=\frac{1}{2}(t-{t}^{\prime }){E}_{\mathrm{Th}}$ have been overlooked in [14].

Equation (18) holds for $u\gt {E}_{\mathrm{Th}}\sim {N}^{1-\gamma }$, and S(u) is saturated at $S(0)\approx {{\rm{e}}}^{-2\pi {{\rm{\Lambda }}}^{2}\;{N}^{1-\gamma }}$ for $u\lt {N}^{1-\gamma }$, which corresponds to $| t-{t}^{\prime }|$ smaller than the inverse total spectral band-width. For $\gamma =1$ we have ${N}^{1-\gamma }=1$ . Thus the value $S(0)\sim {{\rm{e}}}^{-2\pi \;{{\rm{\Lambda }}}^{2}}$ at $\gamma =1$ is smaller than 1. So, in addition to the specific critical behavior of S(u) at $\gamma =2$ (shown by the red curve in figure 7) one obtains yet another critical behavior of S(u) at the ergodic transition $\gamma =1$ (shown by the dashed yellow line in figure 7) which is stable in the $N\to \infty$ limit and is characterized by $S(0)\lt 1$. For $\gamma \lt 1$, ${N}^{1-\gamma }$ is increasing with N, making $S(0)=\mathrm{max}S(u)\to 0$ as $N\to \infty$. This is how the GUE limit $S(u)\equiv 0$ is reached.

Note that the fact that ${E}_{\mathrm{Th}}\gg \delta$ affects the level number variance ${\rm{var}}(n)$ (n and ${\rm{var}}(n)$ are the average number of levels and the level number variance in a certain spectral window, respectively) in which a new scale ${E}_{\mathrm{Th}}/\delta \sim {N}^{2-\gamma }$ appears for n at $1\lt \gamma \lt 2$:

$$\begin{eqnarray}{\rm{var}}(n)=\left\{\begin{array}{ll}\sim \mathrm{ln}n, & 1\ll n\ll {N}^{2-\gamma }\\ n, & {N}^{2-\gamma }\ll n\ll N.\end{array}\right.\end{eqnarray} \tag{ 19 }$$

The level compressibility χ [33, 34] is ill-defined in our model because of the jump in the spectral form factor $R({\mathfrak{t}}){| }_{N\to \infty }$ at ${\mathfrak{t}}=0$. Formally it can take any value from $\chi =0$ to $\chi =1$ depending on the parameter $n/{N}^{2-\gamma }$. This is in contrast to other RM models with multifractal eigenstates, e.g. PLBRM and Moshe–Neuberger–Shapiro models (see [34] and references therein) where the level compressibility is well defined and takes a definite value $0\lt \chi \lt 1$ in the $N\to \infty$ limit.

In the introduction we mentioned a heuristic relation between the Anderson model on a hierarchical RRG and the RP model which has no apparent hierarchical structure. It is instructive now to compare the main results of this paper with the corresponding results for RRG.

First of all we recall (see figure 1 and the corresponding explanations in the text) that all the moments $| \psi {| }^{2q}\rangle$ in the localized and the AT critical phases of our model have exactly the same $q$-dependence as in the corresponding phases of the RRG. The $N$- and $\gamma$- dependence of the moments is also very similar (see figure 5 with figure 8) to the corresponding N and W-dependences for the random wave functions obtained by the exact diagonalization of the Anderson model on RRG with the branching number $K=2\;\mathrm{and}$ $N=2000-16000$. However, there is an important difference: we found only one point of changing the sign of ${C}_{1/N}$ on RRG which corresponds to the known point of the AL transition at $W\approx$ 17.5.

Zoom In Zoom Out Reset image size

In figure 9 we demonstrate that in the case of RRG the scaling of the Thouless energy with the system size follows the same equation (12) as for the present RP model and is different from the standard Chalker's scaling. The falling part of $K(\omega )\;\mathrm{at}\;\omega \gg {E}_{\mathrm{Th}}$ can be described by the unified expansion equation (16) both for RRG and our model, albeit the coefficient ${c}_{0}$ is zero for the present model and is non-zero for RRG. It is important that for RRG the data for properly re-scaled $F(y)={N}^{-z}K({{yN}}^{-z})$ at different $N$ collapse on one same scaling curve $F(y)\equiv {F}_{W}(y)$ that, however, depends on $W$, as well as the exponent $z=1-{D}_{2}$ (see inset of figure 9). This is very different from the usual scaling $K(\omega )={N}^{a}\;F(\omega {N}^{z},{N}_{{\rm{c}}}/N)$ in the vicinity of a single critical point $W={W}_{{\rm{c}}}$ where the exponents $a\;\mathrm{and}\;z$ are determined by the property of this critical point and not by the distance $| W-{W}_{{\rm{c}}}|$ from this point which determines only the correlation volume ${N}_{{\rm{c}}}(| W-{W}_{{\rm{c}}}| )$. In our opinion, this implies that there is a line of critical points at $W\lt {W}_{{\rm{c}}}$ which determine the behavior of the system at least in a parametrically large interval of sizes ${N}_{{\rm{c}}2}\gg N\gg {N}_{{\rm{c}}}$, with the second characteristic size scale ${N}_{{\rm{c}}2}\gg {N}_{{\rm{c}}}$.

Zoom In Zoom Out Reset image size

Our conclusion is that the localized and AT critical states are very similar for our model and RRG. The extended states show non-ergodicity in a broad interval of $\gamma \;\mathrm{and}\;W$ and are characterized by the Thouless energy which in both models is much larger than the mean level spacing. However, the existence or non-existence of the ergodic transition is more subtle and depends on tiny features of the model. It exists in our model and most probably does not exist on RRG with the branching number K = 2. Nonetheless our study largely confirms expectation on the similarity between the RRG and RP models. This allows us to speak on the special class of models with the explicit (RRG) or hidden (RP) hierarchical structure.

We are grateful to B L Altshuler, E Bogomolny, G Biroli, M S Foster, L B Ioffe, J P Pekola and O Yevtushenko for stimulating discussions. We would like to specially thank Pragya Shukla for letting us know about her forthcoming work [36], where the second transition at $\gamma =1$ was mentioned. IMK acknowledges the support of the European Union Seventh Framework Programme: INFERNOS (FP7/2007–2013) under grant agreement no. 308850 and of Academy of Finland (Project Nos. 284594, 272218). VEK acknowledges the hospitality of Aalto University and at LPTMS of University of Paris Sud at Orsay and support under the CNRS grant ANR-11-IDEX-0003-02 Labex PALM, project MultiScreenGlass.