Robustness of delocalization to the inclusion of soft constraints in long-range random models

TL;DR

This paper investigates how soft constraints affect localization-delocalization transitions in long-range random matrix models, revealing unexpected robustness of delocalized phases through numerical and analytical methods.

Contribution

It introduces a new analytical approach to understand the robustness of delocalized phases under soft constraints in long-range random models.

Findings

Delocalized phases are robust to correlations in long-range hopping.

Analytical methods confirm numerical observations of robustness.

Eigenstate calculation method is developed for correlated random matrices.

Abstract

Motivated by the constrained many-body dynamics, the stability of the localization-delocalization properties to the inclusion of the soft constraints is addressed in random matrix models. These constraints are modeled by correlations in long-ranged hopping with Pearson correlation coefficient different from zero or unity. Counterintuitive robustness of delocalized phases, both ergodic and (multi)fractal, in these models is numerically observed and confirmed by the analytical calculations. First, matrix inversion trick is used to uncover the origin of such robustness. Next, to characterize delocalized phases a method of eigenstate calculation, sensitive to correlations in long-ranged hopping terms, is developed for random matrix models and approved by numerical calculations and previous analytical ansatz. The effect of the robustness of states in the bulk of the spectrum the inclusion of to soft constraints is generally discussed for single-particle and many-body systems.

Figures

Equations

$$\hat{H}=\hat{\epsilon }+\hat{j},\hat{\epsilon }$$

$$\hat{j}$$

$$j_{mn}/j_{m^{\prime }{n}^{\prime }}=f(m,n)/f(m^{\prime },n^{\prime })>0.$$

$$0<⟨j_{mn}{j}_{m^{\prime }{n}^{\prime }}⟩<\sqrt{⟨j_{mn}^2⟩⟨j_{m^{\prime }{n}^{\prime }}^2⟩}.$$

$$⟨{\epsilon }_n⟩=0,⟨{\epsilon }_n^2⟩=W^2.$$

$$j_{mn}=j_{nm}^{\ast }=j^0(|m-n|)+h_{mn}\tilde{j}(|m-n|),$$

$$⟨j_{mn}⟩=j^0(|m-n|),{\sigma }_j=\tilde{j}(|m-n|),$$

$${\tilde{j}}_{PLRBM}=\frac{1}{|m-n{|}^{\tilde{a}}}\text{ and }{\tilde{j}}_{RP}=N^{-\tilde{\gamma }/2},$$

$$j_{\mathrm{BM}}^0=\frac{1}{|m-n{|}^{a_0}}\text{, and }{j}_{\mathrm{YS}}^0=N^{-{\gamma }_0/2}.$$

$$f(\alpha )=1-\alpha +\underset{N\to \infty }{\lim }\ln [P(|{\psi }_E(n){|}^2=N^{-\alpha })]/\ln N,$$

$$|{\psi }_E(n){|}_{typ}^2\equiv \exp \left[⟨\ln |{\psi }_E(n){|}^2⟩\right]\sim |n-n_0{|}^{-a_{\mathrm{eff}}}.$$

$$f(\alpha)=\left\{\begin{array}[]{lc}1+(\alpha-\tilde{\gamma})/2,&\alpha_{\min}<\alpha<\tilde{\gamma}\\ -\infty,&\alpha<\alpha_{\min}\text{ or }\alpha>\tilde{\gamma}\\ \end{array}\right.\ ,$$

$${\gamma }_{\mathrm{YS}}=\max ({\gamma }_0,2)\ge 2.$$

$$|{\psi }_E(n)|\propto |{\psi }_E(n_0)|/|n-n_0{|}^a$$

$$a_{\mathrm{BM}}=\max (a,2-a)\ge 1.$$

$$j_{mn}=h_{mn}\tilde{j}(|m-n|)+h_{mn}^0{\tilde{j}}^0(|m-n|),$$

$$j_{mn}=h_{mn}\max \left[\tilde{j}(|m-n|),{\tilde{j}}^0(|m-n|)\right]$$

$${\gamma }_{\mathrm{eff}}$$

$$a_{\mathrm{eff}}$$

$$\hat{j}^{0}=\sum_{\left\langle m,n\right\rangle}j^{0}(|m-n|)\left|m\right\rangle\left\langle n\right|\equiv\sum_{p}j^{0}_{p}\left|p\right\rangle\left\langle p\right|\$$

$$\hat{M}=(1+{\hat{j}}^0/E_0{)}^{-1}=\sum _p\frac{|p⟩⟨p|}{1+j_p^0/E_0}.$$

$$\left[(1+{\hat{j}}^0/E_0{)}^{-1}(E-\hat{\epsilon }+E_0)-E_0\right]|{\psi }_E⟩=0.$$

$$\left[(1+{\hat{j}}^0/E_0{)}^{-1}(E-\hat{\epsilon }-\hat{\tilde{j}}+E_0)-E_0\right]|{\psi }_E⟩=0,$$

$$\left[\hat{\epsilon }+\hat{\tilde{j}}+{\hat{\tilde{j}}}^0+\hat{r}\right]|{\psi }_E⟩=E|{\psi }_E⟩.$$

$${\tilde{j}}_{mn}^0=-\frac{E-{\epsilon }_n+E_0}{N+E_0{N}^{{\gamma }_0/2}}\sim N^{-{\gamma }_{\mathrm{YS}}/2},$$

$$r_{mn}=\frac{{\sum }_l{\tilde{j}}_{ln}}{N+E_0{N}^{{\gamma }_0/2}}\sim \frac{N^{(1-\tilde{\gamma })/2}}{N+E_0{N}^{{\gamma }_0/2}}\ll {\tilde{j}}_{mn}.$$

$$(E_n-{\epsilon }_k){\psi }_{E_n}(k)=j^0\sum _l{\psi }_{E_n}(l)+\tilde{j}\sum _l{h}_{kl}{\psi }_{E_n}(l),$$

$${\psi }_{E_n}(k)$$

$$J_{kn}$$

$$a_n$$

$$P_{kn}$$

$${⟨{\left|{\psi }_{E_n}(k)\right|}^2⟩}_{h_{kn}}\sim \frac{⟨|{\psi }_{E_n}(n)J_{kn}{|}^2⟩}{({\epsilon }_n-{\epsilon }_k+\Delta E_n{)}^2+{\Gamma }_n^2}.$$

Full text

I Introduction

Absence of thermalization in interacting many-body quantum systems has attracted significant interest and boosted numerous studies of different possibilities to violate eigenstate thermalization hypothesis (ETH) both in static and driven systems. The first and most developed way to do this is to randomize system parameters by including disorder. This phenomenon is called many body localization (MBL) Basko et al. (2006); Gornyi et al. (2005). Like in single-particle case of Anderson localization Anderson (1958) disorder induces destructive interference and provokes emergent local integrals of motion Serbyn et al. (2013); Huse et al. (2014) blocking the excitation transport.

An alternative way to break ETH in absence of disorder is to add some hard constraints to the many-body system, that crucially reduce the Hilbert space by separating the Hamiltonian into the disjoint sub-block structure, see Fig. 1(a). These hard constraints can be realized either by infinitely strong interactions Turner et al. (2018a, b); Choi et al. (2019); Khemani et al. (2019), additional integrals of motion Sala et al. (2019); Pai and Pretko (2019); Pai et al. (2019), or gauge invariance Brenes et al. (2018).

As a result, such hard constraints produce special low-entanglement states (such as many-body scars) in the bulk of the spectrum Turner et al. (2018a, b), giving significant contribution to the typical infinite-temperature states and revealing themselves via infinitely long-lived oscillations in quenched observables Bernien et al. (2017). However, in real life none of barriers is infinite. The effect of soft constraints on the thermalization in such systems is non-trivial and under hot debate nowadays as the finite-energy barriers between disjoint Hilbert space sub-blocks might be prevailed at high temperature, Fig. 1(c,e).

For this reason it is of fundamental importance to study a simple model in which hard and soft constraints can be easily realized and corresponding localization properties can be precisely investigated. This would provide efficient criteria to characterize the effects of soft constraints in generic cases.

The straightforward analog of hard constraints in single-particle systems is given, e.g., by fully-correlated long-ranged hopping in random matrix models Burin and Maksimov (1989); Owusu et al. (2009); Modak et al. (2016); Deng et al. (2018); Nosov et al. (2019). Indeed, in these models the complete correlations of all hopping terms, Fig. 1(b), impose the hard constraints like in the case of many-body scars and localize the states in the bulk of the spectrum. However, due to the single-particle nature of these systems, all states in the spectral bulk become localized TI- . Soft constraints in this case can be easily realized by considering partial correlations of long-ranged hopping with non-integer Pearson correlation coefficient. In this work we consider such single-particle disordered models with partially correlated long-ranged hopping, Fig. 1(d, f), investigate both numerically and analytically their eigenstate statistics and phase diagrams, and reveal unexpected robustness of delocalized phases to soft constraints, Fig. 1(f).

In most random matrices the delocalized side of the localization transition (corresponding to ETH in many-body systems) is represented by Gaussian Wigner-Dyson ensembles Mehta (2004). One of most well-known examples of a $d$-dimensional random matrix model confirming this statement and demonstrating the Anderson localization transition (ALT) at any $d$ (including $d\leq 2$) is the power-law random banded matrix ensemble (PLRBM) Mirlin et al. (1996),

[TABLE]

written in $d$-dimensional basis of $N$ lattice sites $\left|n\right\rangle$. This model is characterized by the independent Gaussian distributed hopping terms $j_{mn}$, with the standard deviation $\langle j_{mn}^{2}\rangle^{1/2}\propto|m-n|^{-a}$ power-law decaying at large distances $|m-n|\gg 1$ and its ALT is governed by the decay exponent $a$ responsible for the ratio of the on-site disorder $\varepsilon_{n}=H_{nn}$ On- and hopping terms. The system shows ergodic (localized) wave-function statistics for $a<d$ ($a>d$), while the ALT point $a=d$ is characterized by so-called nonergodic extended (multifractal) wave functions typical for the ALT phase diagram at the criticality MF- ; PLR .

However recently there have been found several models showing the whole nonergodic extended phases, see, e.g., Kravtsov et al. (2015); Nosov et al. (2019). The milestone random matrix example in this row is the Rosenzweig-Porter ensemble (RP) Rosenzweig and Porter (1960). This nominally $1$d model has infinitely long-ranged independent Gaussian distributed hopping elements with the $N$-dependent variance $\langle j_{mn}^{2}\rangle\propto N^{-\gamma}$ and apart from the ALT transition at $\gamma=2$ Pandey (1995); Brézin and Hikami (1996); Guhr (1996); Altland et al. (1997); Kunz and Shapiro (1998); Shukla (2000, 2005), it exhibits an ergodic transition (ET) at $\gamma=1$ from the ergodic phase ($\gamma<1$) to a whole phase of nonergodic extended states ($1<\gamma<2$) characterized by a non-trivial fractal support set De Luca et al. (2014); Kravtsov et al. (2018) of wave functions Kravtsov et al. (2015) with the fractal Hausdorff dimension $D=2-\gamma$, $0<D<1$. This behavior has been further confirmed by several analytical and numerical papers Facoetti et al. (2016); Truong and Ossipov (2016); Amini (2017); von Soosten and Warzel (2018, 2019); Monthus (2017); Bogomolny and Sieber (2018a); De Tomasi et al. (2019); RP- .

The question of constraints (hopping correlations) imposed in both above mentioned models has been considered recently in Nosov et al. (2019). Indeed, the new paradigm of the ALT suggested there states that hopping correlations $\left\langle j_{mn}j_{m^{\prime}n^{\prime}}\right\rangle-\left\langle j_{mn}\right\rangle\left\langle j_{m^{\prime}n^{\prime}}\right\rangle\neq 0$ shrink in general an ergodic phase towards smaller disorder strengths extending both localized and multifractal phases. In the case when all hopping integrals are fully-correlated (with unit Pearson’s coefficient) the localization at any disorder strength is restored Burin and Maksimov (1989); Deng et al. (2018); Celardo et al. (2016) similar to the case of the short-range Anderson model in $d=1,2$ Anderson (1958); Abrahams et al. (1979). An example of such random matrix models with fully-correlated hopping elements $j_{m\neq n}=C|m-n|^{-a}$, decaying with the distance $|m-n|$ as a power-law like in PLRBM, has been suggested in a seminal paper by Burin and Maksimov (BM) Burin and Maksimov (1989). The infinitely long-ranged limit of this model (analogous to RP) with complete correlations between hopping terms $j_{mn}=CN^{-\gamma/2}$ has been shown to be exactly integrable by Yuzbashyan and Shastry (YS) Owusu et al. (2009); Modak et al. (2016). Both these models demonstrate localization for all eigenstates, except measure zero, for all values of parameters $a$ and $\gamma$. Note that the statistics of the site-independent scalar $C\sim 1$ does not play any role here.

A representative type of correlations considered in Ref. Nosov et al. (2019) is the hard constraint (correlations with Pearson’s coefficient $+1$) of certain pairs, $(m,n)$ and $(m^{\prime},n^{\prime})$, of hopping terms,

[TABLE]

Here $f(m,n)>0$ is the deterministic function of indices $m$, $n$, and possibly of the system size $N$. For uncorrelated models (like PLRBM and RP) the pairs are only $(m^{\prime},n^{\prime})=(n,m)$, while in fully-correlated examples (BM and YS) all pairs of $(m^{\prime},n^{\prime})$ and $(m,n)$ are involved. In the intermediate translation-invariant case $|m^{\prime}-n^{\prime}|=|m-n|$ TI- .

In this paper we address a complimentary aspect of soft constraints, namely partial hopping correlations of all pairs of hopping terms

[TABLE]

For this we consider both PLRBM and RP models with finite hopping average values that interpolate between original uncorrelated PLRBM and RP ensembles and their fully-correlated counterparts, BM and YS models 111Here we take the scalar $C$ in fully-correlated models to be constant, as its statistics does not play any role, and consider cross-correlations $<j_{mn}j_{m^{\prime}n^{\prime}}>$ instead of cross-covariances $<j_{mn}j_{m^{\prime}n^{\prime}}>-<j_{mn}><j_{m^{\prime}n^{\prime}}>\neq 0$ without loss of generality..

The unexpected stability of the delocalized phases to soft constraints in both cases is demonstrated. The delocalization is shown to survive even for relatively narrow distributions with mean values $j^{0}$ much larger than the width $\tilde{j}$, see Fig. 1(f). The positions of the ALT and possible ET are shown to be governed solely by the distribution width $\tilde{j}$, but not by relative fluctuations $\tilde{j}/j^{0}$, see Fig. 2 and $P_{\rm eff}$ in Fig. 1(f). This brings us to the conclusion that soft constraints added to the initially delocalized phase do not break delocalization of any state in the bulk of the spectrum, Fig. 1(e), even if the hard constraint does, Fig. 1(a).

In order to uncover the origin of this counterintuitive result we first use the matrix inversion trick suggested in Nosov et al. (2019) to rewrite the eigenproblem in the coordinate basis in an alternative way. Furthermore we develop the self-consistent method of eigenvector calculation based on the averaging over off-diagonal matrix elements, allowing one to access wave-function statistics and, in particular, confirming the phenomenological ansatz known in the literature for RP ensemble Monthus (2017); Bogomolny and Sieber (2018a); De Tomasi et al. (2019) (see also Truong and Ossipov (2016)). Unlike the standard renormalization group analysis Levitov (1990); Mirlin et al. (1996) or the Wigner-Weisskopf approximation Monthus (2017) used in the literature before this self-consistent method is sensitive to the hopping correlations. In the current problem the full ALT diagram of previously mentioned models is calculated with help of these methods.

The rest of the paper is organized as follows. In Sec. II we formulate the random matrix models in the focus. Sec. III shows how the naive guess for the behavior of these models fails, and provide numerical results along with localization-delocalization phase diagrams. Sec. IV describes the matrix inversion trick which explains the behavior of these models and allows us to uncover the origin of unexpected stability of delocalized phases. In Sec. V we demonstrate the self-consistent method of eigenfunction calculation on the example of RP ensemble with finite mean hopping values. In Conclusion we sum up or results and give an outlook.

II Models

Throughout the text we focus on the generalized Anderson’s single-particle model with long-ranged hopping terms, represented by the Hamiltonian (1). The uncorrelated diagonal disorder is given by independent identically distributed random on-site energies $\varepsilon_{n}$ with zero mean and fixed variance

[TABLE]

The summation in (1b) is taken over pairs of sites $m,\leavevmode\nobreak\ n$ coupled by hopping integrals

[TABLE]

characterized by the distant-dependent mean and standard deviation values, respectively

[TABLE]

where $\sigma_{j}^{2}=\left\langle j_{mn}^{2}\right\rangle-\left\langle j_{mn}\right\rangle^{2}$. For simplicity we restrict our consideration to $d=1$, unless stated otherwise. Here and further we denote by $h_{mn}$ i. i. d. random variables with zero mean $\left\langle h_{mn}\right\rangle=0$ and unit variance $\left\langle|h_{mn}|^{2}\right\rangle=1$.

The PLRBM and RP ensembles correspond to $j^{0}=0$ and

[TABLE]

respectively, while for BM and YS models in contrast $\tilde{j}=0$,

[TABLE]

As infinitely long-ranged models (like RP and YS) do not have the notion of distance, the main tool used to characterize the properties of their delocalized and localized phases is the standard multifractal (MF) analysis. This analysis is based on the spectrum of fractal dimensions Evers and Mirlin (2008)

[TABLE]

defined via the distribution of wave-function amplitudes $P(|\psi_{E}(n)|^{2})$ with the $N$-scaling of eigenfunctions $|\psi_{E}(n)|^{2}\propto N^{-\alpha}$, see Fig. 3.

Other long-ranged models (like PLRBM and BM) provide additional tools, e.g., the spatial decay of eigenfunctions with the distance $|n-n_{0}|$ from its maximal value at $n=n_{0}$ Deng et al. (2018); Nosov et al. (2019) given by the typical wave-function decay, see Fig. 4,

[TABLE]

Here $\langle\ldots\rangle$ denotes the average over disorder and over eigenstates in the middle of the spectrum. The energy level statistics (see, e.g., Atas et al. (2013)) as basis-invariant characteristics gives the definite information about the fully-ergodic (Wigner-Dyson) phase and the phase localized in a certain basis with the Poisson level statistics Nosov et al. (2019); r-s .

For RP-model the spectrum of fractal dimensions $f(\alpha)$ is shown to be linear in $\alpha=-\ln|\psi_{E}(n)|^{2}/\ln N$ for $\tilde{\gamma}\geq 1$, with the slope $1/2$ Kravtsov et al. (2015)

[TABLE]

and an additional point $f(0)=0$ for $\tilde{\gamma}>2$. Here $\alpha_{\min}=\max(0,\;2-\tilde{\gamma})$. The $f(\alpha)$ in the ergodic phase, $\tilde{\gamma}<1$, coincides with the one at $\tilde{\gamma}=1$ and is represented by the only point $f(1)=1$ Kravtsov et al. (2015).

The Yuzbashyan-Shastry (YS) ensemble (or as sometimes called the Type-1 model) Owusu et al. (2009); Modak et al. (2016); Ossipov (2013); Celardo et al. (2016); Deng et al. (2018) characterized by deterministic infinitely long-range hopping terms $j_{mn}=j^{0}g_{m}g_{n}$, with the constants $g_{m}\propto O(1)$ of order of one, is exactly integrable Owusu et al. (2009); Modak et al. (2016) and known to have all localized states for $j^{0}<1/N$ and all, except one, localized states for $j^{0}>1/N$ Modak et al. (2016); Ossipov (2013); Celardo et al. (2016); Nosov et al. (2019). The generalization of this ensemble to $N$-dependent hopping elements (8) shows a single-site localized phase for $\gamma_{0}>2$ and a critical behavior at $\gamma_{0}<2$ with the spectrum of fractal dimensions given by (11) with $\gamma_{0}$ replaced by the following expression Nosov et al. (2019) (see Appendix A for analytical derivations)

[TABLE]

PLRBM undergoes the ALT at $a=1$, showing ergodic behavior at $a<1$ and power-law localization at $a>1$ MF- ; PLR , with the decay rate equal to the parameter $a$

[TABLE]

at the large distance $|n-n_{0}|\gg 1$ from the maximal point $n_{0}$, $\max_{n}|\psi_{E}(n)|=|\psi_{E}(n_{0})|$.

In BM-model Burin and Maksimov (1989); Rodriguez et al. (2000, 2003); Balagurov et al. (2004); de Moura et al. (2005); Celardo et al. (2016); Deng et al. (2018); Nosov et al. (2019), the fully-correlated counterpart of PLRBM, which is determined by the Hamiltonian (1) with the hopping elements (8), all, except measure zero of the states, are power-law localized (13) in the entire region of the parameter $a$. However, the power-law decay rate $a_{\rm BM}$ is not equal to $a$, but instead is always larger than one (see Deng et al. (2018); Nosov et al. (2019) for details)

[TABLE]

III Intuitive guess and numerical results

What would be the phase diagram for general models (1 - 6) with both finite mean $j^{0}$ and fluctuating $\tilde{j}$ hopping terms? For the first glance, it is natural to expect that the behavior of uncorrelated models (7) should be dominant as soon as $\tilde{j}\gg j^{0}$ ($\tilde{\gamma}<\gamma_{0}$ or $\tilde{a}<a_{0}$) as the distribution of hopping elements is relatively wide and nearly centered at zero, see Fig. 1(d), and vice versa the models with deterministic hopping (8) should dominate at $\tilde{j}\ll j^{0}$ ($\tilde{\gamma}>\gamma_{0}$ or $\tilde{a}>a_{0}$) when the distribution is relatively narrow and its width can be neglected, Fig. 1(b,f).

However, this is not the case. Indeed, from numerical calculations one can see that these models undergo the ALT (and the ET for RP case) at the same points as their uncorrelated counterparts: $\tilde{\gamma}=2$ ($\tilde{\gamma}=1$) and $\tilde{a}=1$ irrespective to the amplitude $j^{0}$ as if all mean values are zero, $j^{0}=0$, Fig. 2. Moreover, the wave-function statistics of such models in all phases, Figs. 3 and 4, coincides with the one of a simple mixture of two uncorrelated long-ranged hopping models of the type (7), with hopping terms

[TABLE]

where the function $\tilde{j}^{0}(|m-n|)$ is described by (8) with the parameters $\gamma_{0}$ and $a_{0}$ of fully-correlated models (YS and BM) replaced by their effective values (12) and (14), respectively. Thus, in the leading approximation one can replace the mixture of hopping elements by their maximum

[TABLE]

and map the model with partial correlations to the one with zero mean (equivalent to RP and PLRBM) and parameters $\gamma$ and $a$ in Eqs. (11, 13) replaced by effective ones

[TABLE]

This is the main result of this paper shown here numerically, Figs. 3 and 4, and confirmed further analytically.

The qualitative explanation of this unexpected stability of delocalized states originates from the fact that the models with deterministic long-range hopping terms (like YS and BM) demonstrate only localized eigenfunction statistics in the bulk of the spectrum, but never delocalized. As a result, in the mixture these fully-correlated models can compete with their uncorrelated counterparts only in the localized phase, $\tilde{a}>1$ or $\tilde{\gamma}>2$, Fig. 2, affecting the wave-function spatial profile in the locator expansion Anderson (1958); Levitov (1989, 1990) (see also Nosov et al. (2019) for more detailed discussion). In terms of soft constraints this means that as soon as typical states in the spectral bulk are considered the infinite temperature corresponding to them prevail over the finite barrier of soft constraint between previously disjoint sub-blocks of the Hamiltonian and brings the system to the phase where it was before imposing constraints.

To understand the origin of the above mentioned behavior of partially-correlated models (3), summarized in Eqs. (15 –17a), in the next section we describe a matrix inversion trick Nosov et al. (2019) providing an alternative representation of the eigenproblem and apply it to the mixture of YS and RP model as an example.

IV Matrix inversion trick

Here we restrict our consideration of the matrix inversion trick to the case of the mixture of RP and YS models (for the mixture of PLRBM and BM models please see Appendix B). For the first time this method has been suggested by us in Ref. Nosov et al. (2019) to analytically prove the duality of the eigenfunction power-law decay in the $1$d BM-model  (14) numerically discovered in Ref. Deng et al. (2018) and to generalize both Anderson localization Anderson (1958); Levitov (1989, 1990) and Mott ergodicity Mott (1966) principles for the models with correlated hopping. However recently there have been found several many-body Luitz and Bar Lev (2019); Maksymov et al. (2019) and higher-dimensional, $d>1$, single-particle models Cantin et al. (2018); Deng et al. (2019), applied to which this method easily uncovers their phase diagrams and the wave-function structure by the extension the locator expansion validity range.

Let’s first consider the pure deterministic (BM or YS) model (8). The matrix inversion trick is based on the spectral properties of the hopping matrix

[TABLE]

diagonalized in a certain basis $\left|p\right\rangle$ (momentum basis for BM model). If the spectrum $j^{0}_{p}$ of this matrix diverges from either side in the thermodynamic limit (e.g., $j^{0}_{\max}=\max_{p}j^{0}_{p}\rightarrow\infty$ for $N\to\infty$ and $j^{0}_{\min}=\min_{p}j^{0}_{p}$ is finite) or even from both sides, but has a finite spectral gap, see Fig. 5(a), one can diminish the effect of these divergent terms to the hopping elements, inverting the matrix $(1+\hat{j}^{0}/E_{0})\equiv\hat{M}^{-1}$

[TABLE]

This matrix inversion sends the diverging top-spectrum (or edge-spectrum) terms close to $j^{0}_{\max}$ to the denominator of the sum while the condition $E_{0}+j^{0}_{\min}>0$ avoids the divergence of the contributions from the bottom (or states close to the gap) of the spectrum, see Fig. 5(b). The optimization of $E_{0}\sim N^{\beta}$ over the parameter $\beta$ Nosov et al. (2019) gives the smallest effective hopping terms at $\beta=0$ in the whole parameter range (please see Appendix B for details).

After the matrix inversion trick the problem takes the form

[TABLE]

The diagonal part $M_{0}$ of the matrix $(1+\hat{j}^{0}/E_{0})^{-1}_{m,n}=M_{m-n}$ forms effective on-site disorder $M_{0}\varepsilon_{n}$ and eigenvalue $M_{0}(E+E_{0})-E_{0}$ of the problem, while the hopping terms are formed by $M_{m-n\neq 0}(E-\varepsilon_{n}+E_{0})$.

The main idea behind this matrix inversion trick uses the fact that the eigenstates with large hopping energies $|j^{0}_{p}|\gg 1$ are barely affected by the disorder $\hat{\varepsilon}$. Thus, they nearly coincide with those hopping matrix eigenstates $\left|p\right\rangle$ that give the main contribution to (18). All other eigenstates corresponding to small hopping energies $|j^{0}_{p}|\sim O(1)$ are orthogonal to these large-energy states at the spectral edge and thus almost orthogonal to the main contribution to the hopping matrix given by them. As a result the states in the bulk of the spectrum “see” the hopping terms $M_{m-n\neq 0}(E-\varepsilon_{n}+E_{0})$ which are significantly reduced compared to the initial ones $j^{0}_{m-n}$.

For the case of YS model the rank-$1$ matrix $\hat{j}^{0}$ has the only non-zero eigenvalue $j^{0}_{p=0}=N^{1-\gamma_{0}/2}$ corresponding to the zero-momentum state $\left\langle n|0\right\rangle_{p}=N^{-1/2}$ with arbitrary rest basis states $\left|g_{k\neq 0}\right\rangle$ orthonormal to $\left|0\right\rangle_{p}$ 222For the general YS model the latter should be modified by $<n|g_{0}>=g_{n}/G$, with $G^{2}=\sum_{n}|g_{n}|^{2}$..

The same matrix inversion trick (20) can be applied as well to the model (6) with partial correlations (3)

[TABLE]

where we inverse only the deterministic hopping part with the semi-infinite spectrum.

As a result Eq. (21) takes the predicted form of Eq. (15)

[TABLE]

with the fluctuating elements of the matrix $\hat{\tilde{j}}^{0}$

[TABLE]

which do not break the locator expansion and the residual term $\hat{r}$ small compared to $\hat{\tilde{j}}$

[TABLE]

In Eq. (23) $\gamma_{\rm YS}$ is given by (12). This derivation confirms our numerical observation (15) and concludes this section. The consideration of the mixture of PLRBM and BM models is addressed in Appendix B.

To sum up, in this section we have shown that the effective locator expansion result Nosov et al. (2019) analogous to (13) is applicable in the case of the RP-YS mixture to calculate the wave function in the whole localized phase coinciding with the one of the uncorrelated model ($\tilde{\gamma}>2$) at any $\gamma_{0}$, Fig. 2.

To calculate the wave-function statistics in all phases, including delocalized ones, in the next section we develop a self-consistent method sensitive to hopping correlations.

V Self-consistent eigenstate calculation

In this section we consider the self-consistent method of the wave-function calculation, generalizing the perturbation theory. For simplicity we restrict our consideration to the mixture of RP and YS ensembles. For more general analysis please see Appendix C. Separating hopping terms $j_{mn}=j^{0}+h_{mn}\tilde{j}$ (5) in the Hamiltonian (1) into deterministic $j^{0}=j^{0}_{YS}=N^{-\gamma_{0}/2}$ (8) and fluctuating $\tilde{j}=\tilde{j}_{RP}=N^{-\tilde{\gamma}/2}$ (7) parts in the eigenproblem

[TABLE]

one can formally write its solution

[TABLE]

in terms of the sums

[TABLE]

and $E_{n}-\varepsilon_{n}=J_{nn}$ given by (26b), that should be calculated self-consistently. Here and further we choose the index $n$ in the energy $E_{n}$ in such a way that in absence of off-diagonal elements $j^{0}=\tilde{j}=0$ the wave function is localized at $k=n$, $\psi_{E_{n}}(k)=\delta_{k,n}$, with the energy $E_{n}=\varepsilon_{n}$.

Averaging (26) over the hopping elements $h_{mn}$ with fixed bare energies $\varepsilon_{n}$ gives the following expression for the wave-function intensity

[TABLE]

Here we assume that the sums $P_{mn}$ and $a_{n}$ are self-averaging and thus they are uncorrelated from each other and from $h_{kl}$. For RP model itself this approximation of self-averaging over hopping terms has been used in several papers Facoetti et al. (2016); Truong and Ossipov (2016) and confirmed there by other methods. The fluctuating energy shift $J_{nn}=E_{n}-\varepsilon_{n}$ after this averaging leads both to the energy shift $\Delta E_{n}$ and the level broadening $\Gamma_{n}$ in Eq. (28). Both $\Delta E_{n}$ and $\Gamma_{n}$ are of the same order as $E_{n}-\varepsilon_{n}$ and in principle contain contributions from all cumulants of $P_{nn}$ and $a_{n}$, however for our analysis it is enough to consider only mean values and the first-order perturbation theory term $\tilde{j}h_{nn}$ with standard deviation $\tilde{j}$ taken into account (for further details see Appendix C)

[TABLE]

Here $\Gamma_{n}$ plays a role of the level broadening. This level broadening determines the size of the miniband of almost fully-correlated eigenfunctions like in RP-model Kravtsov et al. (2015); De Tomasi et al. (2019). The factor $\left\langle|\psi_{E_{n}}(n)J_{kn}|^{2}\right\rangle$ in the numerator of Eq. (28) guarantees the wave-function normalization and not important for the wave-function statistics.

Focusing on the $N$-scaling of $\Gamma_{n}\sim\Delta E_{n}$ one can show that the localized state realizes at $\Gamma_{n}$ smaller than the mean level spacing $\delta\simeq W/N$ of the model without hopping, the ergodic state corresponds to $\Gamma_{n}$ large compared to the bare band of the system $W\sim O(1)$, while the fractal phase appears at intermediate values:

[TABLE]

In order to estimate the scaling of the level broadening and identify the corresponding phases one can make use of the self-consistent equations for $P_{kn}$ and $a_{n}$ which could be obtained by substituting expressions (26) for $\psi_{E_{n}}(k)$ and $E_{n}$ to (27). The resulting equations read as

[TABLE]

The next essential step is to average these equations over $h_{nm}$ relying on the above-mentioned self-averaging properties of these sums. As a result, Eqs. (31) take the form

[TABLE]

with the sum $S_{1}$ given by

[TABLE]

The latter can be calculated, e.g., in two limiting cases of (i) the completely rigid spectrum $\varepsilon_{n}-\varepsilon_{m}=(n-m)\delta$, with $\delta$ being a bare mean level spacing without hopping, and (ii) the Poisson statistics of $\varepsilon_{n}$. Up to prefactors unimportant for the multifractal analysis $S_{1}$ takes the form

[TABLE]

Considering cases (30) and substituting (34) into (32) and (29) one easily obtains the result up to prefactors $O(1)$

[TABLE]

with $\gamma_{\rm YS}=\max(\gamma_{0},2)$ given by (12).

Comparison of the resulting $\Gamma_{n}$ with the one of the RP model shows the same $N$-scaling with the parameter $\tilde{\gamma}$ replaced by $\gamma_{\rm eff}$ from Eq. (17a). One can also check it by the direct calculation of the multifractal spectrum. Indeed, the resulting approximate wave function (28) scales as

[TABLE]

where the scaling of $\Delta E_{n}\sim\Gamma_{n}$ is given by (35) and we parameterized $(\varepsilon_{n}-\varepsilon_{k})\sim N^{-p/2}$ with the parameter $p$. For a certain distribution of the diagonal terms $\varepsilon_{n}$ with the width $W\sim O(1)$ the scaling of the probability reads as

[TABLE]

for positive $p>0$ as $P(\varepsilon_{n}\to 0)\sim O(1)$. For negative $p<0$ the probability is at least exponentially small in $\varepsilon_{n}$ and thus, for multifractal analysis one should neglect it focusing on $p\geq 0$. As a result the spectrum of fractal dimensions

[TABLE]

can be found from the expression

[TABLE]

where the maximization in r.h.s. is taken with respect to the condition $N^{-\alpha}\sim(N^{-\tilde{\gamma}}+N^{-\gamma_{\rm YS}})/({N^{-p}+\Gamma_{n}^{2}})$. The maximal probability is given by the condition $0<p<-2\ln\Gamma_{n}/\ln N$ leading to $p=\gamma_{\rm eff}-\alpha$ and to the result (11) with $\gamma_{\rm eff}$ given by (17a). Moreover the boundaries $0<p<-2\ln\Gamma_{n}/\ln N$ provide the correct bounds for $\alpha$

[TABLE]

This analysis confirms numerical results shown in Fig. 3 and concludes this section. Note that the method developed in this section is powerful and accurate for the multifractal analysis as one can take into account cumulants of any order of the sums (31) fluctuating with the hopping terms and sensitively distinguish models with slightly different hopping correlations.

VI Conclusion and discussions

To sum up, in this work we address the effect of soft constraints on the phase diagram of random matrix models with long-ranged correlated hopping. We demonstrate unexpected robustness of the delocalized phases to partial hopping correlations imposed by soft constraints and determine wave-function statistics and corresponding phase diagrams of milestone disordered long-range models, power-law banded random matrix and Rosenzweig-Porter ensembles. This main result (17) is confirmed by both numerical calculations and two analytical approaches. The matrix-inversion trick developed in Nosov et al. (2019) uncovers the effective Hamiltonian (22) and confirms the wave-function behavior in the localized phase. The self-consistent method allows to calculate wave-function statistics in delocalized phases as well and confirm the main result of the paper.

A parallel drown between constrained random-matrix models and many-body systems, brings us to the conclusion that in general soft constraints added to the initially delocalized phase do not break delocalization of any typical (infinite temperature) state, even if the hard constraint does. Indeed, the infinite temperature corresponding to the typical states in the spectral bulk prevail over the finite barrier of soft constraint between previously-disjoint sub-blocks of the Hamiltonian and brings the system to the phase where it was before imposing constraints. However, the relations of slow-dynamics phenomena Lezama et al. (2019) to hard and soft constraints both in many-body systems Brenes et al. (2018); Turner et al. (2018a, b); Choi et al. (2019); Khemani et al. (2019); Sala et al. (2019); Pai and Pretko (2019); Pai et al. (2019); Torres-Herrera and Santos (2015); Távora et al. (2017); Torres-Herrera et al. (2018); Santos and Torres-Herrera (2018) and in closely related single-particle disordered models De Tomasi et al. (2019); Bera et al. (2018) is still under debates and consideration. The question of the dynamics and relaxation of highly non-local operators Khaymovich et al. (2019) is also in the focus of the research in the community.

Another intriguing question is how the interplay between partial correlations in hopping and interaction amplitudes could affect localization properties in many-body systems with either or both hopping and interaction terms being long-ranged in the coordinate space. Specifically, the limiting fully-correlated case of the interacting version of Burin-Maksimov model was recently analyzed in Yao et al. (2014); Singh et al. (2017); Luitz and Bar Lev (2019); Botzung et al. (2018, 2018); Roy and Logan (2019); Nag and Garg (2019), and the opposite situation without any constraints was considered in details in Burin and Kagan (1994); Burin (2015a, b); Gornyi et al. (2017); Tikhonov and Mirlin (2018); De Tomasi (2019). However, the intermediate regime represented by both finite means and dispersion in distribution functions of matrix elements needs deep consideration.

Finally, in order to more precisely investigate the role of partial correlations in long-range random matrix models, it would be insightful to construct the corresponding effective field-theoretic description with imposed soft constraints. For instance, this approach may naturally incorporate ergodicity-breaking phenomena in gauge-invariant lattice models Brenes et al. (2018). Although the desirable theory is more technically involved than usual super-symmetric non-linear sigma model due to finite means of hopping and effective non-locality, several attempts were made to describe systems with correlations in off-diagonal terms Ossipov (2013); Ossipov and Kravtsov (2006). Moreover, so-called “virial expansion” for almost diagonal random matrices developed in Yevtushenko and Kravtsov (2003); Yevtushenko and Ossipov (2007) seems to be a suitable candidate for an appropriate representation of constrained models with disjoint sub-blocks in Hilbert space.

Acknowledgements.

We are grateful to V. E. Kravtsov for stimulating discussions. P. N. appreciates warm hospitality of the Max-Planck Institute for the Physics of Complex Systems, Dresden, Germany, extended to him during his visits when this work was done. P. N. acknowledges funding by the RFBR, Grant No. 17-52-50013, and the Foundation for the Advancement to Theoretical Physics and Mathematics BASIS Grant No. 17-11-107. I. M. K. acknowledges the support of German Research Foundation (DFG) Grant No. KH 425/1-1 and the Russian Foundation for Basic Research Grant No. 17-52-12044.

Appendix A Yuzbashyan-Shastry model

A.1 Model

Consider an ensemble of random $N\times N$ Hermitian matrices

[TABLE]

with a random real vector $\varepsilon=(\varepsilon_{1},\ldots,\varepsilon_{N})$ and deterministic complex vector $g=(g_{1},\ldots,g_{N})$ with $N$-independent elements. $\varepsilon_{i}$ are statistically independent entries with a zero mean and a unit variance

[TABLE]

while $j^{0}$ scales with the matrix size $N$ as $j^{0}\propto N^{-\gamma_{0}/2}$. Note that the translation-invariant case of $g=(1,\ldots,1)$ corresponds to the finite mean $j^{0}$ of the hopping elements (plus additional unimportant shift of energy).

As Type-1 Hamiltonian Owusu et al. (2009); Modak et al. (2016) $H_{0}$ provides an exact eigenproblem solution Owusu et al. (2009)

[TABLE]

with the (possibly $N$-dependent) normalization constant $C_{n}=j^{0}\sum_{m}g_{m}\psi_{E_{n}}(m)$,

[TABLE]

and the secular equation for the spectrum $E=E_{n}$

[TABLE]

giving all eigenvalues (except the highest one $E_{N}>\varepsilon_{N}$) lying between adjacent bare levels

[TABLE]

Here we assumed $\varepsilon_{m}$ to be ordered in ascending order $\varepsilon_{m}<\varepsilon_{m+1}$, $m=\overline{1,N-1}$. Note that this model includes the limiting case $a=0$ of long-range deterministic hoppings considered in Deng et al. (2018) both for positive ($g_{m}=1$) and staggering ($g_{m}=(-1)^{m}$) hopping elements.

A.2 Spectral statistics

To find the spectrum $E_{n}$ one should consider the secular equation in more details. Let’s assume that all $g_{k}$ are of the same order $g_{k}\simeq g\sim N^{0}$ and consider the variation of $E_{n}$ from its bare value $\varepsilon_{n}$ as

[TABLE]

Separating positive and negative summands of the sum (45) for $E=E_{n}$ we obtain

[TABLE]

Now we have to estimate a typical value of the sum taking into account the inequality (47). To do so we consider two limiting cases.

(i) In the limit of a completely rigid spectrum $\varepsilon_{n}=n\delta$, with the mean level spacing $\delta\sim 1/[\rho(0)N]$ and the density of states at the Fermi level $\rho(0)$, the sum (48) can be taken explicitly

[TABLE]

giving

[TABLE]

where

[TABLE]

Analogously the normalization constant governed by (44) takes the form

[TABLE]

In the latter equality we neglect prefactors, focus only on the $N$-scaling and take into account that $\Gamma_{n}\lesssim\delta$.

(ii) In the opposite limit of uncorrelated eigenstates one can calculate (48) as follows

[TABLE]

where we take into account (47) in the lower limits of integration. Here the density of states is

[TABLE]

Again considering only $N$-scaling of $\Gamma_{n}$ in (53) one can obtain

[TABLE]

Analogously the normalization constant governed by (44) takes the form of (52)

[TABLE]

As the scaling (50, 55) of the energy deviation $\Gamma_{n}$ (47) and of the normalization constant $C_{n}$ (52, 56) are the same in both limiting cases, we conclude that these parameters weakly depends on the statistics of bare levels $\varepsilon_{n}$.

A.3 Eigenstate statistics

Using the results (43) and (52) one can calculate the spectrum of fractal dimensions $f(\alpha)$ for the wave function intensity

[TABLE]

Indeed, as the probability of $\varepsilon_{n}-\varepsilon_{k}\sim N^{-p/2}$ for $p>0$ is

[TABLE]

one can easily find

[TABLE]

and

[TABLE]

for $p=\gamma_{\rm YS}-\alpha>0$, giving the spectrum of fractal dimensions of the form of Kravtsov et al. (2015)

[TABLE]

with $\gamma_{\rm YS}$ given by (51).

As a result, unlike the Rosenzweig-Porter model, the YS model (41) shows only localized and critical wave functions of all eigenstates, except the only top energy state at $\gamma_{0}<2$ Modak et al. (2016). At $\gamma_{0}>2$ the wave-function statistics (61) coincides with the one of the RP, with $\gamma_{\rm YS}=\gamma_{0}$, while at all $\gamma_{0}<2$ instead of the delocalized phases YS model shows the critical localization with $\gamma_{\rm YS}=2$.

Appendix B Matrix inversion trick

In this Appendix we first give general estimate of the optimal parameter $\beta$, $E_{0}\sim N^{\beta}$ and then apply the matrix inversion trick Nosov et al. (2019) to the PLRBM model with partial correlations (3)

B.1 Optimization over $E_{0}$ in the matrix inversion trick

Let’s focus for simplicity on the case of $E\sim W\sim N^{0}$. The effective hopping matrix $J_{m-n}=M_{m-n}(E+E_{0}-\varepsilon_{n})/M_{0}$ in Eq. (20)

[TABLE]

can be estimated as follows

[TABLE]

Here we neglected the term $\varepsilon_{n}$ for simplicity as $E\sim\varepsilon_{n}\sim N^{0}$ and divided Eq. (20) by

[TABLE]

in order to have diagonal disorder in the standard form of $\varepsilon_{n}$. Here $j_{p,\min}^{0}$ is the typical hopping energy level and $E_{0}$ is taken to be in the gap of the hopping spectrum or beyond it (see Fig. 5).

Each term in the sum (63) can be minimized over $E_{0}$ giving

[TABLE]

However, one should take into account that $E_{0}$ has to be beyond the spectrum $|E_{0}|\gtrsim|j^{0}_{p,\min}|$, leading to the final result

[TABLE]

In the case of BM or YS deterministic hopping $\hat{j}^{0}$, the typical energy level $j^{0}_{p,\min}\sim N^{0}$, thus, the optimal $E_{0}\sim N^{0}$ and this confirms the statement given in the main text.

B.2 Matrix inversion for PLRBM model with partial correlations

We start with the model (1, 6) with $\tilde{j}_{mn}=|m-n|^{-\tilde{a}}$, $j^{0}_{mn}=|m-n|^{-a_{0}}$ and compute matrix elements $\tilde{j}^{0}_{mn}+r_{mn}$ of the effective Hamiltonian (22)

[TABLE]

The first terms in the brackets of r.h.s. $\left(1+\hat{j}^{0}/E_{0}\right)^{-1}\left(E-\hat{\varepsilon}+E_{0}\right)$ scales as $|n|^{-(2-a)}$ at $a<1$ and correspond to $(1-\alpha)\hat{\tilde{j}}^{0}$ with a certain constant $\alpha$ of order of one. The calculation of it is given in Nosov et al. (2019). Further we consider the rest part

[TABLE]

In order to simplify the calculations we provide the upper bound of this term by replacing the oscillating amplitudes $h_{mn}$ of $\tilde{j}_{mn}=h_{mn}/|m-n|^{\tilde{a}}$ by their maximal absolute value $h_{mn}=1$. Within this approximation the result can be easily derived in the momentum space, since both matrices are diagonal in that basis. For purposes of clarity we reproduce main steps of similar calculations presented in Nosov et al. (2019) specifically for the case of our interest.

We start by writing down the Fourier-transformed hopping amplitude $j^{0}$ (in case of $\tilde{j}$ one needs to replace $a_{0}$ by $\tilde{a}$) in different asymptotic regimes

[TABLE]

[TABLE]

where the corresponding constants are given for $a_{0}>0$ by

[TABLE]

for $a_{0}\neq 2m+1,m\in\mathbb{N}$ and

[TABLE]

Here $\zeta_{a_{0}}$ is the Riemann zeta function.

The next step is to estimate the long-range asymptotic behavior of the effective Hamiltonian $(1-\alpha)\hat{\tilde{j}}^{0}+\hat{\tilde{j}}+\hat{r}$ given by

[TABLE]

with the zero-momentum contribution

[TABLE]

The last term in (74) can be split into three parts corresponding to different regimes of $\tilde{j}_{p}$ and $j^{0}_{p}$. However, not all of the resulting terms are equally important in the limit $1\ll n\ll N$. One can easily show that summation only over sufficiently small momenta $|p|<\alpha N$ (where $0<\alpha<1/2$) contributes to the long-range behavior of matrix elements. The rest of the summation results in the gives small contribution to the residue term $r_{n}$ and is unimportant. Thus, we focus only on the following sum

[TABLE]

In the case of our main interest ($a_{0}<1$) an additional momentum scale $p_{c}=N[E_{0}/2A_{a_{0}}+\zeta_{a_{0}}/A_{a_{0}}]^{-1/(1-a_{0})}\leq N$ emerges, and for $p<p_{c}$ the denominator is represented by the power-law contribution. Contrary, for $p>p_{c}$ a constant term is dominant and we get

[TABLE]

which for $a_{0}<1$ leads to the effective parameter governing the long-range tails of typical wave functions:

[TABLE]

in full agreement with the result (17b) mentioned in the main text. The last term in (77) contribute to the residual term as it is small compared to the second term $\sim|n|^{-\tilde{a}}$ in the considered interval $a_{0}<1$.

Appendix C Self-consistent method

In this Appendix we use the formulation of the eigenproblem in terms of exact Eqs. (26 - 31) in order to derive expressions for scaling (35) of the broadening factor $\Gamma_{n}$ in the average wave-function intensity (28).

In the first part we restrict our consideration to the first and second moments of the parameters $a_{n}$, $P_{kn}$, calculate sums $S_{1}$ (34) and $S_{2}$ (see below) in two limiting cases of (i) the completely rigid spectrum $\varepsilon_{n}-\varepsilon_{m}=(n-m)\delta$, and (ii) the Poisson statistics of $\varepsilon_{n}$, derive the mean expressions (32) for the RP model with partial correlations and the effective expression (35) for the broadening parameter $\Gamma_{n}$ depending solely on $S_{1}$. Next, we explicitly calculate $N$-scaling (35) of $\Gamma_{n}$.

C.1 Decoupling of correlations

In order to solve Eqs. (31) we assume that $P_{kn}$ and $a_{n}$ are uncorrelated from each other and from $h_{lm}$ and calculate the mean and the variance for each of them averaging over off-diagonal matrix elements $h_{kn}=h_{nk}$ assumed to be real and taking into account $\left\langle h_{lm}\right\rangle=0$, $\left\langle h_{lm}^{2}\right\rangle=1$. As a result for $k\neq n$,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Here and further we use the notation $\omega_{nl}=\varepsilon_{n}-\varepsilon_{l}$ for brevity.

This gives the following self-consistency equations

[TABLE]

with the sums

[TABLE]

While averaging the denominators $\omega_{nl}+\Gamma_{n}$ we estimate only $N$-scaling of the broadening parameter as follows

[TABLE]

and we consider the typical energy position $\varepsilon_{n}$ to lie a bit asymmetrically in the middle of the spectrum, thus the summation would be in the limits $W_{1}<\omega_{nl}<W_{2}$, with $|W_{1}-W_{2}|\sim O(1)$ and $W_{1}+W_{2}=W$.

To estimate a typical value of the sums (33) we consider two limiting cases.

(i) In the limit of the completely rigid spectrum $\varepsilon_{n}=n\delta$, with the mean bare level spacing $\delta\sim 1/[\rho(0)N]\sim W/N$ and the density of states at the Fermi level $\rho(0)$, the sums can be taken explicitly

[TABLE]

[TABLE]

and provide the following asymptotics

[TABLE]

(ii) In the opposite limit of uncorrelated eigenstates one can calculate (33) as follows

[TABLE]

[TABLE]

Unlike YS model (53) here the broadening parameter can be both smaller and larger than bare mean level spacing $\delta$, thus, during the calculation we just take into account the fact that $E_{n}-\varepsilon_{k}=\omega+\Gamma_{n}$ is off-resonant. Here the density of states is $\rho(\varepsilon)=\sum_{k}\left\langle\delta(\varepsilon-\varepsilon_{k})\right\rangle/N\ .$

In this case asymptotics read as follows

[TABLE]

Note that the expressions for sum $S_{1}$ are the same in both cases, while $S_{2}$ are different only in the non-ergodic extended phase. Let’s show that this difference do not affect the result (35) for the broadening parameter $\Gamma_{n}$. In order to prove it we consider the expression (80) in more details substituting the expressions (79) one by one.

First, let’s substitute $\sigma_{P}$

[TABLE]

As $\min\left(1,\tilde{j}\sqrt{S_{2}}\right)\leq 1$ one can neglect the whole last summand corresponding to $\sqrt{\left\langle P_{kn}^{2}\right\rangle}$ which not larger than $\left\langle a_{n}\right\rangle+\sigma_{a}+\tilde{j}$. Next, we substitute $\sigma_{a}$

[TABLE]

As $S_{1}^{2}\geq S_{2}$ in all phases of both limiting cases, one can neglect $\tilde{j}^{2}\sqrt{S_{2}}$ comparing to $\left\langle P_{nn}\right\rangle\sim\tilde{j}^{2}S_{1}$ and the whole $\sigma_{a}$ term comparing to $\left\langle a_{n}\right\rangle+\tilde{j}$ as $\left(\sqrt{S_{2}}/{S_{1}}\right)\min\left(1,j^{0}S_{1}\right)\leq 1$. As a result we come to the expression (35)

[TABLE]

depending solely on $S_{1}$. This confirms the statement given in the main text and concludes this section.

C.2 Calculation of $N$-scaling (35) of $\Gamma_{n}$

Using Eq. (83) and expressions (81, 82) for $S_{1}$, in this section we calculate the scaling (35) of $\Gamma_{n}$ in all three phases:

$\Gamma_{n}\gg W$. In this case $S_{1}=N/\Gamma_{n}$ and Eq. (83) takes the form

[TABLE]

The second term does not play any role as it is less than $\Gamma_{n}/N\ll\Gamma_{n}$, while the third term dominates over the first one and gives $\Gamma_{n}\sim N^{(1-\tilde{\gamma})/2}$ and thus $\tilde{\gamma}<1$. 2. 2.

$\delta\ll\Gamma_{n}\ll W$. In this case $S_{1}=1/\delta\sim N$ leading to

[TABLE]

As in the previous case the second term does not play any role as it is less than $\delta\ll\Gamma_{n}$, while the third term dominates over the first one and gives $\Gamma_{n}\sim N^{1-\tilde{\gamma}}$ and thus $1<\tilde{\gamma}<2$. 3. 3.

$\Gamma_{n}\ll\delta$. Here $S_{1}\sim 1/\delta$ leading to $\tilde{\gamma}>2$ and

[TABLE]

In this case the first term dominates over the third one and the concurrence of the first two terms gives the desired result $\Gamma_{n}\sim N^{-\tilde{\gamma}/2}+N^{-\gamma_{\rm YS}/2}\sim N^{-\gamma_{\rm eff}/2}$, with $\gamma_{\rm eff}$ given by (17a).