Survival probability in Generalized Rosenzweig-Porter random matrix ensemble

We study analytically and numerically the dynamics of the generalized Rosenzweig-Porter model, which is known to possess three distinct phases: ergodic, multifractal and localized phases. Our focus is on the survival probability $R(t)$, the probability of finding the initial state after time $t$. In particular, if the system is initially prepared in a highly-excited non-stationary state (wave packet) confined in space and containing a fixed fraction of all eigenstates, we show that $R(t)$ can be used as a dynamical indicator to distinguish these three phases. Three main aspects are identified in different phases. The ergodic phase is characterized by the standard power-law decay of $R(t)$ with periodic oscillations in time, surviving in the thermodynamic limit, with frequency equals to the energy bandwidth of the wave packet. In multifractal extended phase the survival probability shows an exponential decay but the decay rate vanishes in the thermodynamic limit in a non-trivial manner determined by the fractal dimension of wave functions. Localized phase is characterized by the saturation value of $R(t\to\infty)=k$, finite in the thermodynamic limit $N\rightarrow\infty$, which approaches $k=R(t\to 0)$ in this limit.


Abstract.
A recent generalization of the Rosenzweig-Porter random matrix model offers a simple but remarkable example of a system possessing a whole phase of multifractal extended states. This phase is separated from the localized and ergodic phases by the Anderson localization and the so-called ergodic transition, respectively. Here, we study the dynamics of this system initially prepared in a highly-excited non-stationary state (wave packet) confined in space and containing a fixed fraction of all eigenstates. Our focus is on the probability for finding the initial state later in time, the so-called survival probability R(t). Three main aspects are identified in different phases. The ergodic phase is characterized by the presence of oscillations in R(t) with the frequency equal to the band width of the wave packet. These oscillations surviving in the thermodynamic limit is a signature of the ergodic phase. In multifractal extended phase the survival probability shows an exponential decay which rate is determined by the fractal dimension of wave functions. Localized phase is characterized by the saturation value R(t → ∞) which is non-zero in the thermodynamic limit N → ∞. In the case when the initial state is localized on a single site R(t → ∞) → 1 in the limit N → ∞.

Introduction
A rapidly expanding field of Many-Body Localization (MBL) pioneered in Refs. [1,2,3] includes, among other directions, a search for non-ergodic extended (NEE) phase nicknamed as "bad metal" [4,5]. In such a phase many-body wave functions are neither localized nor occupying all the available Hilbert space either, similar to multifractal one-particle wave functions at the transition point of an ordinary Anderson localization problem [6]. One of the consequences of the lack of ergodicity in such a phase is the behavior of the survival probability R(t) at large times t → ∞ in a system of finite dimension N of the Hilbert space. In an interacting disordered system an initially prepared product state Ψ(t = 0) is evolving with time such that the quantity | Ψ(t)|Ψ(0) | 2 averaged over disorder realizations, known as "survival probability" R(t), tends to a time-independent limit: Ergodicity implies an asymptotic equipartition over all available many-body states, and thus D = 1. Many-body localization means that only finite number of states is involved, and therefore D = 0. The NEE phase is characterized by 0 < D < 1, so that the volume in the Hilbert space occupied at t → ∞ is infinite N D → ∞ in the thermodynamic limit, yet it constitutes zero fraction N D−1 → 0 of all states. A very simple system where all three phases: ergodic, localized and non-ergodic extended phases exist in the corresponding range of parameters, is the Generalized Rosenzweig-Porter (GRP) random matrix ensemble [7,8]. The existence of the nonergodic extended state and a transition to the ergodic state in this model has been recently suggested in Ref. [8], and further confirmed in Ref. [9,10] with a rigorous prove in Ref. [11]. The role of initial product state Ψ 0 = Ψ(t = 0) in this model is played by the basis vector in which only one component is 1 and all the others are equal to 0. The final state Ψ(t → ∞) is a random N -vector in which N D elements are of the order of ∼ N −D and all the others are much smaller.
In this paper we consider the time dependence of the survival probability: where the N × N matrix Hamiltonian is drawn from the GRP ensemble [8] characterized by independent identically distributed fluctuating entries H nm , such that: Here A denotes ensemble average of A, N is the matrix size, and γ, λ are real parameters. As was shown in Ref. [8], it is the parameter γ that determines to which of the three above phases the system belongs. The "fractal dimension" D in the NEE phase has been shown [8] to be equal to: The other parameter, λ, is largely redundant and is important only at the point where the transitions from the NEE to the localized (D = 0) and ergodic (D = 1) phases occurs. These transitions will be referred to as the localization and the ergodic transitions, respectively. The main results of this paper can be formulated as follows, see • In the ergodic phase 0 < γ < 1 the survival probability exhibits oscillations which decay algebraically with time and persist even if only a small fraction f of all states is involved in dynamics; the frequency of oscillations is equal to the spectral width of the initially prepared wave packet and thus is inversely proportional to f . The presence of such oscillations in the thermodynamic limit N → ∞ is the signature of the ergodic phase. It reflects the rigidity of system of energy levels.
• In the NEE phase the survival probability decays exponentially R(t) = exp[−Γ(N ) t] with time as long as R(t) N −D . The characteristic energy Γ(N ) ∝ N D−1 . Some oscillations in R(t) may exist, but only as a finite-size effect.
• In the localized phase R(t) tends to a constant R(∞) as t → ∞, where R(∞) remains finite and tends to its initial value R(0) as N → ∞.

Survival probability in terms of eigenfunctions and eigenvalues
The survival probability R(t) can be expressed in terms of eigenvalues ε n and eigenfunctions ψ n (r) of the state n: Indeed, the initial wave function Ψ(r, t = 0) = δ r,0 which is non-zero only on one site r = 0, can be represented as follows using the completeness of the set of eigenfunctions: The Schrödinger dynamics ψ n (r, t) = ψ n e i εn t then leads to: which immediately results in Eq. (6) for R(t) = |Ψ(0, t)| 2 . The summation in Eq. (6) goes over all states n. If, however, the summation in Eqs. (7), (6) is restricted to a finite number of states constituting a fraction f < 1 of all the states with energies ε n close to the center of the band, ε = 0, the corresponding initial density |Ψ(r, t = 0)| is spread over a number of sites close to r = 0. We will consider evolution of such a wave packet for different fractions f < 1 which is not only convenient for computational purposes but also leads to an important test for ergodicity.

Ansatz for eigenvectors of GRP model
The key point of our analytical consideration is the following ansatz (see also [10,12,13]) for the amplitude of a random wave function ψ n (r) of GRP model which maximum is on a site n: For NEE states (1 < γ < 2) [8] the characteristic energy Γ(N ) has a meaning of the width of a mini-band [5,14] in the local spectrum formed by N D hybridized quasiresonance sites with the mean local level spacing δ mb : The typical |ε n − ε r | ∼ 1 corresponds to the typical amplitude |ψ n (r) The level spacing δ mb (N ) = C mb δ(N ) ∝ N −1 inside a mini-band scales with N like the global mean level spacing δ(N ) ∝ N −1 , albeit with a coefficient C mb < 1 which depends on γ and fluctuates from mini-band to mini-band.
For ergodic states γ < 1 the global density of states (DoS) has a semi-ellipse form [15,8]: so that in the center of the band δ(N ) = (ρ(0) N ) −1 ∝ N −(1+γ)/2 [8]. The mini-bands are absent in this case, and grows with N like the global band-width f N δ. As Γ 1 δ(N ), the typical and the maximal amplitudes are of the same order |ψ n (r) For localized states, is just the mean splitting |H nm | 2 of two resonance levels [8]. In this case the peak amplitude at n = r is the maximal |ψ n (r)| 2 of a typical wave function at r = n is of the order of: while the typical scales in the same way with γ as in NEE phase [8].
Let us define the eigenvector overlap function: where . . . off denotes averaging over the off-diagonal matrix elements H nm only. Plugging Eq. (9) into Eq. (17) one obtains for delocalized phases with Γ δ: In deriving this equation we assumed that the eigenvalue fluctuations are independent from those of the off-diagonal matrix elements H nm . Then we switched from summation over r to integration over ε r and used the fact that the convolution of two Cauchy functions is also the Cauchy function with the double width. The first, most crucial assumption is based on the large number of different paths contributing to a given quantum state.
In the localized phase, where Γ δ, the main contribution to Eq. (17) is done by the two terms with r = n and r = m ‡. Taking into account that |ψ n (n)| 2 ≈ 1 one obtains:

R(t) and the eigenvalue distribution
Finally we obtain the expression for R(t) in terms of K(ε n − ε m ): where . . . ε denotes averaging over the eigenvalue distribution. In this equation we also took into account a possibility that only a fraction of all states f < 1 close to zero energy is taken into account. Eq. (20) can be rewritten (21) ‡ In the delocalized phases, Γ δ, the contribution of r = n and r = m are also of the same order as K(ω) itself.
using the DoS correlation function:

Limiting cases
Let us first assume that the energy levels are uncorrelated C(ε, ε ) =ρ(ε)ρ(ε ). Then close to ε = 0 we obtain from Eq. (21): This immediately results in an exponential behavior for both delocalized phases: with a finite-size saturation value Another limiting case (which is not literally realizing in disordered systems) is the case of absolutely rigid energy spectrum ε n − ε n+ = δ. In this case it is convenient to use Eq. (20) and do summation first over (n + m)/2 and then over = n − m. We obtain: Eq. (25) can be evaluated using the Poisson summation formula: For delocalized states, the function F (x) takes the form: where A = 2π fρ(0) 2 λ 4 N 2(1−γ) /(Γδ),ρ(0) is the mean DoS and g = 2Γ/δ. Evaluating the integral: with y = xδt/(2πk), k = k m,± = tΓ π ± 2(Γ/δ)m, and the ratio of widths of the band and the mini-band we can express the r.h.s. of Eq. (26) in terms of the cosine integral (ci(x)) and the sine integral (si(x)) functions of complex arguments c ± i. The first two lines in Eq. (28) correspond to the first term in (1 − |x|/c), and the last two lines follow from the term −|x|/c.
In the limit 2π k √ c + 1 1 this expression is very well approximated by: Note that the term ∝ sin(2πkc)/k from the first two lines of Eq. (28) cancels out by the corresponding term from the last two lines. At t 2π/δ it is sufficient to take only one term with m = 0 in the sum Eq. (26), so that |k| = tΓ/π, and Y (k, c) is proportional to R(t) − R(∞) (see Fig. 2): where A/g ∼ 1/f is an N -independent constant for extended states (0 < γ < 2) and A/g ∝ N 1−γ/2 vanishes in the thermodynamic limit for localized states (γ > 2). Thus we conclude from Eqs. (30),(31) that the r.h.s. of Eq. (25) oscillates with the period inversely proportional to the band-width: This period is N -independent for γ > 1 and it is proportional to N −(1−γ)/2 for γ < 1. According to Eqs. (30),(31) the amplitude of oscillations decays as 1/t 2 . For ergodic phase, γ < 1, in agreement with Eqs. (30, 31) the survival probability, Eq. (25), coincides with the Fourier transform of the global DoS, Eq. (22), of the considered fraction f < 1 of states. These oscillations have the same origin as those considered in [16,17] for the whole band-width with f = 1 in the ergodic phase. However, because of the sharp edge of the DoS for f < 1 they decay as t −2 , but not as t −3 as in the case of a semi-elliptic DoS [16,17]. Another advantage of considering f < 1 is that in this case the oscillations are a direct evidence of the spectrum rigidity rather than an indirect one [16,17] through the sharp edge of the semi-ellipse DoS.

Oscillations in R(t) and the spectral rigidity
The physical reason for oscillations in R(t) is the fixed number f N of states participating in R(t) and the assumed regularity of spectrum ε = δ. For random Hamiltonians the spectrum is also random but exhibits the spectral rigidity in the ergodic phase. It implies the variance var(n) of the number of levels in an energy interval containing n levels on the average, to be much smaller than for statistically independent (Poissonian) fluctuation of levels typical for localized states: var loc (n) = n.
This level rigidity ensures that the behavior of R(t) is practically indistinguishable from the one in the simple example of the regular spectrum considered above. Indeed, the fluctuation of the width of the energy strip containing levels results in smearing of the edge x = ±f N of the θ-function in Eq. (27) by w ∼ δ ln( ). This makes the amplitude of oscillations acquiring an exponentially decaying factor exp(−w t) in addition to a power-law t −2 decay due to edge singularity. However, the oscillations will not be affected by this exponential decay in the time interval ∆t t w −1 . The existence of such an interval is guaranteed in the ergodic phase, as it requires only a condition: The lack of level rigidity in the non-ergodic phase by itself does not kill oscillations, provided that Γ w ∼ δ edge ∼ N − 1 2 §. In this case the condition similar to Eq. (36) reads: which is easily fulfilled. Much stronger effect on oscillations comes from the reduced size of the jump in F (x) in Eq. (27) rather than from its smearing. This is an effect of eigenvector fluctuations and it is quantified by a large value of the parameter c in Eq. (29) when N δ ∼ 1/ρ(0) ∼ 1 and Γ ∼ N −(1−γ) 1. In Fig. 2 we plot the function Y (k, c) which determines R(t) − R(∞) via Eq. (31) in two cases: (a) small value of c typical for the ergodic phase at a small fraction f 1 of states considered and (b) for large value of c typical for the non-ergodic extended phase at large enough system size N . One can see that oscillations disappear at large c, since the main contribution to Eq. (30) in its region of validity is given by the first, exponentially decreasing but non-oscillatory term.

Numerical test for ergodicity
In Fig. 3 we plot R(t) obtained by numerical diagonalization of GRP matrix Hamiltonians of the size N up to N = 2 14 . The results, which should be compared with Fig. 2, are quantitatively described by the model of the rigid system of energy levels discussed in the previous sections (33).
The oscillations in R(t) can be used as the numerical test for ergodicity of wave functions not only in GRP but also in other models [16,17,18].

Extracting fractal dimension D from R(t) in the NEE phase
As it follows from Eq. (29) the parameter c ∝ N 1−D in Eq. (30) which determines R(t) via Eq. (31), is large in the non-ergodic extended phase characterizing by D < 1. dimension D through measuring the decrement 2Γ of the exponential decay. In Fig. 4 (which should be compared with Fig. 2(b)) we present the plots of R(t) obtained by exact diagonalization of matrix Hamiltonian from GRP ensemble and show an example of extraction of D(γ) at γ = 1.5.

Residual oscillations in R(t) in the NEE phase
One should also take care of the residual oscillations in R(t) with the same universal period (32) which are present at a finite N in the NEE phase close to ergodic transition at γ = 1 when c ∝ λ −2 N γ−1 is not large enough. They can be suppressed either by increasing the system size or by decreasing a value of the constant λ in Eq. (3), cf. Figs. 4 and 6.
In Fig. 6 we demonstrate how with increasing the system size in the NEE phase, γ > 1, the oscillations die out giving way to the exponential decay. We also demonstrate that the exponential decay has a different scaling with N than the oscillating part, the period of oscillations being N -independent.

R(t) in the localized phase
In Fig. 1 we give an overview of the evolution of the form of R(t) with increasing γ above the ergodic transition. The localization transition at γ = 2 is marked by the scale-independent R(t) that has a finite limit R(∞) ∝ I 2 at N → ∞. In the localized phase γ > 2 the limiting value R(∞) remains finite.
For finite N the asymptotic value R(∞) depends on the constant λ in Eq. (3) ,  Fig. 7. At a fixed λ and increasing N the value R(∞) stays constant in the critical point γ = 2 of the localization transition but it moves towards its initial value R(0) = 1 in the insulating phase (see Fig. 8 and Fig. 9). This is related with the special property [8] of the localized phase in GRP ensemble that there is a gap between the peak value (14) of |ψ n (n)| 2 ≈ 1 and the typical maximal value of |ψ n (r = n)| 2 ∼ N −(γ−2) 1. This implies that the Lyapunov exponent (which does not coincide with the inverse localization length in GRP model) is divergent in the thermodynamic limit, as for the Bethe lattice with infinite connectivity [24].

R(t) in the Yuzbashyan-Shastry models
Surprisingly, the GRP model is very different from the model (referred to as Yuzbasgyan-Shastry (YS) model) with deterministic off-diagonal elements H nm = λ N −γ/2 [21,19,20,22], which naively could be thought of as the closest relative to GRP model. However, the YS model belongs to a different universality class of models [23] where long-range (and even distance-independent like in YS model) hopping does not lead to delocalization. The YS model is the simplest particular case of exactly solvable random matrix models in which H n =m = g n g * m [21,22]. In such models R(∞) is finite in N → ∞ limit no matter what the nature of g n : g n can be deterministic and n-independent (with a possible N -dependence g n ∼ N −γ/4 as in YS model), or, alternatively, a random gaussian variable with zero mean. The majority of eigenstates (except for measure zero number of states at the edge of spectrum) are like in the RP model for γ ≥ 2: they are either critical (as for γ = 2 in GRP) or localized (as for γ > 2 in GRP) but never extended [13]. This is the reason why R(t) in such models always has a finite limit as t → ∞ and N → ∞.

Conclusions
In summary, the detailed analysis of the survival probability R(t) in the generalized Rosenzweig-Porter ensemble [8] demonstrates three distinctly different types of behavior in different phases. The ergodic phase shows robust oscillations in R(t) with a wide region of polynomial decay ∼ t −2 of their amplitude with time and a universal frequency equal to the band width. These oscillations have a universal nature related with spectral rigidity and they can be used to identify the ergodic extended phase in different models. An exponential decay of the survival probability in multifractal extended phase is nearly free of the finite size effects and provides an accurate way to extract the multifractal dimension of eigenstates from the decay rate. Localized phase is characterized by the universal size-independent saturation value R(t → ∞) = R(∞) coinciding with the initial value R(t → 0) = R(0) in the thermodynamic limit N → ∞. For all γ > 2 in the insulating phase R(∞)/R(0) increases with increasing the system size and should approach 1 in the thermodynamic limit. At the localization transition γ = 2 the ratio R(∞)/R(0) is saturated for N > 10000 at a λ-dependent value smaller than 1.