Scale-invariant critical dynamics at eigenstate transitions

The notion of scale-invariant dynamics is well established at late times in quantum chaotic systems, as illustrated by the emergence of a ramp in the spectral form factor (SFF). Building on the results of the preceding Letter [Phys. Rev. Lett. 131, 060404 (2023)], we explore features of scale-invariant dynamics of survival probability and SFF at criticality, i.e., at eigenstate transitions from quantum chaos to localization. We show that, in contrast to the quantum chaotic regime, the quantum dynamics at criticality do not only exhibit scale invariance at late times, but also at much shorter times that we refer to as mid-time dynamics. Our results apply to both quadratic and interacting models. Specifically, we study Anderson models in dimensions three to five and power-law random banded matrices for the former, and the quantum sun model and the ultrametric model for the latter, as well as the Rosenzweig-Porter model.


I. INTRODUCTION
Since the pioneering paper by M. Berry [1], the spectral form factor (SFF) is considered as one of the key diagnostics of quantum chaos.Specifically, one says that the dynamics of a quantum system is chaotic if the SFF after sufficiently long time follows the so-called ramp, i.e., a linear increase in time that is predicted by the random matrix theory (RMT) [2].The SFF is commonly normalized such that different quantum chaotic systems follow the same ramp, irrespective of the system size [3][4][5].Emergence of such a scale-invariant dynamical behavior represents a convenient tool to detect universality of chaotic dynamics, as shown by a great amount of recent studies .The concept of SFF can be extended to survival probabilities [26] of initial nonequilibrium wavefunctions, which, at least within the RMT, also exhibit scale-invariant dynamics at long times [2,[27][28][29].Of particular importance is the time of the onset of scale-invariant chaotic behavior in the SFF, denoted as the Thouless time t Th .However, t Th may be very long and it typically increases with system size.For example, in systems of linear size L that exhibit diffusion, one expects t Th ∝ L 2 [3,8,10,[30][31][32].
The main focus of this paper are the dynamical transitions from chaos to localization, for which critical behavior in finite systems is associated with an abrupt transition in eigenstate properties.The eigenstate transitions are relatively well understood in systems described by quadratic Hamiltonians, with Anderson localization being a paradigmatic example [33][34][35][36][37]. Recently, a new class of interacting models based on the avalanche theory [38] was introduced, with the potential to provide a stepping stone into a more detailed understanding of critical behavior at ergodicity breaking phase transitions [5,[38][39][40][41][42][43][44][45][46][47][48][49][50].
Motivated by the emergence of scale-invariant dynamics in quantum chaotic systems, we here ask two questions.First, does some notion of scale-invariant dynamics persist at criticality?If the answer is affirmative, in which time regimes may one then expect scale-invariant behavior?In particular, does one need to wait for very long times, e.g., to the Thouless time, to observe scale invariance, or does it emerge already at much shorter times?
The main result of this paper and the preceding Letter [49] is that for both quadratic systems exhibiting Anderson localization, and interacting systems exhibiting ergodicity breaking phase transition, the answer to the first question is indeed affirmative, and that scale invariance at criticality emerges already in mid-time dynamics, as sketched in Fig. 1.The central measures of scale invariance are the SFF k, which exhibits a broad plateau in the mid-time dynamics, and the survival probability p of site-localized states, which exhibits a power-law decay.Importantly, both quantities are measured in units of the typical Heisenberg time (proportional to the typical inverse level spacing), and are normalized such that they also exhibit the late-time scale invariance in the chaotic regime.They hence represent useful indicators of both chaotic dynamics deep in the chaotic regime as well as the critical dynamics at the boundary of quantum chaos.
We note that the main results of this paper and the preceding Letter [49] is to establish the link to the critical dynamics of interacting systems at the ergodicity breaking transition point.In the context of quadratic systems, many studies contributed to the identification of certain properties that are responsible for the emergence of critical dynamics.These studies range from the observation of scale invariance in eigenstate correlations [51] and in short-range spectral statistics [52], to studies of longrange spectral statistics, in particular their spectral rigidity [53][54][55][56].The long-range spectral correlations were further argued to be related to the eigenfunction correlations at criticality in Refs.[55][56][57][58], where the mid-time features of the SFF and survival probability were also discussed [58].The later insights were sharpened by noting that the emergence of a broad scale-invariant plateau in the SFF of the Anderson models is a strong indicator of criticality [4,59].An important connection to interacting models was recently made by observing a very similar plateau in the SFF of the interacting quantum sun (QS) 1. Concept of mid-time and late-time dynamics at eigenstate transition, for the dynamics of the SFF k and survival probability p of initially site-localized states.We measure time t in units of typical Heisenberg time t typ H , such that τ = t/t typ H , and we normalize k and p such that they approach k = p = 1 at τ > 1.The late-time dynamics is defined as the dynamics in between the scaled Thouless time τ Th = t Th /t typ H and the scaled typical Heisenberg time τ typ H = 1.At the transition, the late-time dynamics is preceded by the scale-invariant mid-time dynamics, where k exhibits a broad plateau, k ≈ kχ, and p decays with power-law dependence p ≈ a0 τ −β .While the late-time dynamics are universal in quantum chaotic sense, the properties of midtime dynamics, such as the values of Kχ and β, depend on the properties of transition point.model at the ergodicity breaking transition point [5].
Here we follow the perspective that the SFF is a survival probability of a special initial state [60], which is an equal superposition of Hamiltonian eigenstates.It is then natural to seek for scale invariance also in survival probabilities of other initial states.In quadratic systems, several studies contributed to the understanding that the survival probability of the initial site-localized states exhibits a power-law decay in the vicinity of Anderson transition [61][62][63].This observation paved way towards scale invariance in mid-time dynamics, as well as towards making the connection of the power-law exponent with the wavefunction fractal dimension [63][64][65].The first contribution of our Letter [49] was to introduce the scaled survival probability p mentioned above, which exhibits scale invariance in both mid-and latetime dynamics at criticality.The later quantity allows for establishing the connection between the power-law exponent and the wavefunction fractal dimension also in other quadratic models such as the Aubry-André model.The main contribution of [49] was then to generalize these findings to the QS model of the avalanche theory, which exhibits scale-invariant mid-time and late-time dynamics at the ergodicity breaking transition point.In this paper we first provide a comprehensive overview of the subject, and then go beyond the Letter [49] by studying a broad range of different models, ranging from the quadratic models such as the Anderson models in dimensions three to five and power-law random banded matrices, to the interacting models such as the QS model and the ultrametric model defined within the avalanche picture, as well as the Rosenzweig-Porter model that we consider as a separate class of models.Moreover, we also explore to which degree scale invariance in survival probabilities emerges for other initial states, such as plane waves, and we find that it indeed exhibits certain fingerprints of scale invariance, in particular at weak multifractality.
The paper is organized as follows.In Sec.II we define the models studied in this paper.In Sec.III we introduce the measures of quantum dynamics, i.e., the SFF and survival probabilities, and we comment on their similarities.In Sec.IV we discuss the key properties of scaleinvariant dynamics in the chaotic regime and at criticality, and we introduce the scaled survival probabilities that are the central object of investigation.Numerical results for the quadratic models, interacting models and the Rosenzweig-Porter model are presented in Secs.V, VI and VII, respectively.In particular, we establish similarity between the quadratic and interacting models, and we highlight subtle differences present in the Rosenzweig-Porter model.We summarize our findings in Sec.VIII, and we conclude in Sec.IX.

II. MODELS
In this section we introduce the models that exhibit eigenstate transitions and whose critical dynamics are studied in this paper.We split them into two main categories, quadratic models and interacting models, and finally we also consider the Rosenzweig-Porter model as a separate category of models.

A. Quadratic models
We start with the well-known Anderson model [33] on a d-dimensional hypercubic lattice of linear size L, where ĉ † j (ĉ j ) are the fermionic creation (annihilation) operators at site j, J is the hopping matrix element between nearest neighbor sites, ni = ĉ † i ĉi is the site occupation operator, and ϵ i are the on-site energies.The later are independent and identically distributed random variables, drawn from a box distribution ϵ i ∈ [−W/2, W/2].We consider periodic boundary conditions.Since the model is quadratic, the number of lattice sites coincides with the single-particle Hilbert space dimension, D = L d .
Theoretical arguments suggest the transition in this model to occur in dimensions d > 2 [34].In three di-mensions (3D), numerical studies of transport properties of single-particle eigenstates in the center of the energy band [66][67][68], based on the transfer-matrix technique, showed that the system is insulating at W > W c ≈ 16.54J [69], and at W < W c it becomes diffusive [63,70,71].The value of the transition point W c grows with the dimension d and was accurately calculated using numerical techniques [4,31,[71][72][73][74][75][76].At the transition, the model exhibits subdiffusion [63] and multifractal single-particle eigenfunctions [37,73,77].The transition point is energy dependent, i.e., at W > W c all single-particle states are localized in site-occupation basis, while at W < W c the system exhibits a mobility edge [78].
Next, we introduce the ensemble of power-law random banded (PLRB) matrices.We define the corresponding PLRB Hamiltonian as a quadratic Hamiltonian, whose matrix elements in the single-particle Hilbert space of size D = L are given by (2) i.e., the matrix elements decrease as power laws with distance |i − j| [37,55,56,79].In Eq. ( 2), µ i,j are independent and identically distributed random variables drawn from a box distribution, µ i,j ∈ [−1, 1] [80], but other distributions can also be considered [55,56,79].
The PLRB model is parameterized by two parameters a and b, and it exhibits an eigenstate transition at a = 1 for all b.The single-particle eigenstates are delocalized at a < 1 and localized in site-occupation basis at a > 1.The model was designed as a toy model to understand the features typical for the Anderson model close to and at criticality.At the critical value a = 1 the system exhibits multifractality [55,56,79,81] and spectral statistics are intermediate between the Wigner-Dyson and Poisson statistics [55,56,79].These properties are tuned by the parameter b, ranging from strong multifractality at b ≪ 1 to weak multifractality at b ≫ 1.Some of the dynamical aspects of critical ensembles were considered both analytically [82][83][84] and numerically [80].

B. Interacting models
In the domain of interacting models we consider two representatives, which are both expected to describe the avalanche mechanism of ergodicity breaking phase transitions.
The first model, dubbed the quantum sun (QS) model [5,49,50,85], shares many similarities with the initially proposed toy model of quantum avalanches [38,40].The model consists of N + L spin-1/2 degrees of freedom in a Fock space of dimension D = 2 N +L .It is divided into a quantum dot with N spins and a remaining subsystem with L spins outside the dot, described by the Hamilto- The interactions within the dot, denoted by R, are all-toall and they exclusively act on the dot subspace.They are represented by a 2 N × 2 N random matrix drawn from the Gaussian orthogonal ensemble (GOE) [86].The spins outside the dot are subject to local magnetic fields h i ∈ [0.5, 1.5] that are drawn from a box distribution.Each of the spins outside the dot is coupled to one spin in the dot, with the interaction strength α ui .For a chosen spin i outside the dot, we randomly select an in-dot spin n i .
The coupling to the first spin outside the dot (i = 0) is set to one since u 0 = 0, while at i ≥ 1, is drawn from a box distribution.
The second model is the ultrametric (UM) model.This model was initially introduced in the single-particle picture to study Anderson localization on ultrametric lattices [87][88][89][90][91][92][93][94][95].However, recent study has highlighted similarities between the UM model and the QS model when defined for interacting systems in a Fock space of dimension D [85].Following the convention in the QS model, we define a Fock space of N + L spin-1/2 particles, with D = 2 N +L .The model Hamiltonian is constructed as a sum of block-diagonal random matrices Ĥk with k = 0, 1, . . ., L. At each k, the matrix structure of Ĥk consists of 2 L−k diagonal blocks of size 2 N +k × 2 N +k .We sample each random block independently from the GOE distribution and use normalization such that see [85] for details.Here, the superscript i denotes the ith random block.We sample its matrix elements in analogy to the QS model, hence The full Hamiltonian of the UM model then reads The first term Ĥ0 of size 2 N × 2 N describes the initial quantum dot that consists of N spin-1/2 particles in the absence of coupling to the external particles.The sum in the second term mimics the exponentially decaying coupling between the dot and the k-th localized spin-1/2 particle through the exponentially decaying values of α k .Additionally, we have also included the parameter J c tuning the overall perturbation strength, which carries certain analogies with the parameter g 0 in the QS model in Eq. ( 3).The UM model exhibits a sharp transition at α c = 1/ √ 2 and the properties of the transition point can be tuned by J c .However, here, we fix N = 1, J c = 1 and vary α only.

C. Rosenzweig-Porter model
Finally we consider the Rosenzweig-Porter (RP) random matrix ensemble [96] and its generalized form [97].We introduce the corresponding RP model in a form of a quadratic model given by Eq. ( 2), for which the Hamiltonian matrix elements are given by even though later we consider the RP model independently from other quadratic and interacting models.The parameters λ and c in Eq. ( 6) are real, and µ i,j are independent and identically distributed random variables drawn from a box distribution, Previous studies of the RP model reported emergence in three different regimes [97][98][99][100][101]: at c < 1 the system is fully ergodic, at 1 < c < 2 the eigenstates are fractal and the system is in a nonergodic extended regime, and at c > 2 the eigenstates are localized in the computational basis.Numerical investigation of the short-range spectral statistics [102] indicates that there is a change of the nearest-neighbour level statistics at c = 2 from the GOE to the Poisson statistics.We refer to this transition as the eigenstate transition in the RP model.The statistics of the nearest-neighbour gaps at the transition can be tuned by the parameter λ.We note that the system at the transition, c = 2, is not multifractal but fully localized [97,102].This is an important difference when compared to the PLRB model, which, as we argue in this paper, is also manifested in the dynamical properties of the RP model [80,103].

III. MEASURES OF QUANTUM DYNAMICS
A. Spectral form factor (SFF) We study quantum dynamics of initial states evolved under the Hamiltonian Ĥ.The Hamiltonian eigenspectrum is characterized by the SFF where E ν are eigenenergies of Ĥ.The SFF is further averaged over different Hamiltonian realizations, denoted by ⟨...⟩ H , which gives rise to At long times, K (t) saturates to the corresponding longtime value which equals the time-averaged value over a sufficiently large time interval.The SFF is a Fourier transform of energy level distances and it gives access to both shortrange and long-range spectral correlations.Here we argue that it may also represent a useful tool for understanding the behavior of survival probabilities in the quench dynamics studied below.

B. Survival probability
We are interested in quantum quenches from the eigenstates {|m⟩} of the initial Hamiltonian Ĥ0 to the final Hamiltonian Ĥ with eigenstates {|ν⟩}.The survival probability of an eigenstate |m⟩ is, for a given realization of the Hamiltonian Ĥ, defined as where we set ℏ ≡ 1, D is the Hilbert-space dimension, c νm = ⟨ν|m⟩ is the overlap of |m⟩ with |ν⟩, and E ν is an eigenenergy of Ĥ.The averaged survival probability is defined as where ⟨...⟩ m denotes the average over all eigenstates |m⟩ of the initial Hamiltonian Ĥ0 , and ⟨...⟩ H denotes the average over different realizations of the final Hamiltonian Ĥ.
At long times, P (t) approaches P , which is equal to the average inverse participation ratio of eigenstates of Ĥ in the eigenbasis of Ĥ0 , P = ⟨⟨ ν |c νm | 4 ⟩ m ⟩ H .We express P as i.e., as a sum of the nonzero asymptotic value and a part that vanishes in the thermodynamic limit D → ∞ as ∝ D −γ , where γ > 0 is the fractal dimension.We distinguish between two scenarios that give rise to nonzero P ∞ > 0: the case of a mobility edge, when a fraction of eigenstates |m⟩ are localized in the eigenbasis spanned by |ν⟩ and the remaining fraction is not localized, and the case of complete localization when all eigenstates |m⟩ are localized.When using Eq. ( 12), the definition of γ as the fractal dimension is only meaningful in the former case.

C. Relationship between the SFF and survival probability
A common interpretation of the SFF ( 7) is that it is only a measure of the Hamiltonian spectrum {E ν }.An alternative perspective that we pursue here is that the SFF can be viewed as the survival probability (10) of a special initial state |T ⟩, namely the infinite-temperature pure state such that K H (t) → P H T (t), and the averaging is then only performed over different Hamiltonian realizations.This relationship is valid for any Hamiltonian Ĥ that dictates the dynamics.
The main goal of this paper is to study survival probabilities from different types of initial states, including the state ( 14) that provides the link to the SFF, with the focus on the dynamics at criticality.Specifically, the initial states are of three types.(a) The one from Eq. ( 14), for which the survival probability corresponds to the SFF.(b) Product states |m⟩ = |I⟩ in site occupation (computational) basis for quadratic models, and in Fock space for interacting models.We refer to this type of states as site-localized (Fock-space-localized) states.(c) Product states in quasimomentum occupation basis, , which correspond to conventional plane waves for quadratic models, and to plane waves in the Fock space for interacting models.For brevity we refer to them as plane waves for all models under consideration.
We note that in the special case of the GOE matrices, there is an exact asymptotic relationship between the SFF and survival probabilities from other types of initial states such as site-localized (Fock-space-localized) states and plane waves.We further elaborate on this relationship in Sec.IV B in the context of Eq. (23).
The choice of the initial states imply the value of the time averaged P from Eq. ( 12).In the case (a), it follows from Eq. ( 9) that P ∞ = 0 and γ = 1.In the case (b) one gets P ∞ = 0 in the fully delocalized regime of Ĥ, while P ∞ > 0 in the localized regime or in the regime with a mobility edge.If the initial wave function at the transition exhibits (multi)fractal properties in the eigenbasis of Ĥ, one expects γ < 1.In the case (c) one again gets P ∞ = 0 and γ ≈ 1 for all cases under consideration here.
We note, however, that the latter is not the case in the Aubry-André model [49].
We study survival probabilities in the units of scaled time τ , where t typ H is the typical Heisenberg time, δE typ is the typical level spacing, and ⟨...⟩ ν denotes the average over all pairs of nearest levels.We note that in Hamiltonian systems, the SFF is commonly studied after spectral unfolding that removes the impact of the density of states.However, if one seeks to establish the connection of the SFF with the survival probabilities, which are measured either as a function of raw time t or the scaled time τ , no unfolding is carried out.In Appendix A we compare the SFF with and without unfolding, and we show that the main features of our results are independent of this choice.

IV. SCALE-INVARIANT DYNAMICS
As discussed in introduction, a well-established universality of the SFF occurs in the quantum chaotic regime, when the SFF follows the so-called ramp, predicted by the RMT [2,[27][28][29].Specifically, the quantity that actually exhibits the scale-invariant (i.e., L-independent) behavior in the quantum chaotic regime is the normalized SFF k(τ ), such that at long times one gets k(τ ≫ 1) = 1.The normalized k(τ ) from Eq. ( 16) has been applied in recent studies of quantum chaos in interacting systems, see, e.g., Refs.[3-5, 14, 25, 31, 49].For brevity we refer to the normalized SFF as the SFF.
In what follows, our main goal is to explore whether a similar form of scale-invariant behavior can also be observed at eigenstate transitions from quantum chaos to localization.

A. From SFF to survival probability
Inspired by the scale-invariant behavior of the SFF in the quantum chaotic regime (16), and the relationship between the SFF and survival probabilities from Sec. III C, we now explore features of scale invariance in survival probabilities.In Ref. [49] we introduced the scaled survival probability p(τ ), henceforth survival probability, which should be considered as an analog to the SFF defined in Eq. ( 16).In particular, the survival probability saturates at long times to p(τ ≫ 1) = 1, and the The GOE results for p GOE (τ ; cp) from Eq. ( 20) at D = 2500.Main panel: k GOE (τ ), p loc GOE (τ ) and p pw GOE (τ ), see text for definitions.(Inset) k GOE (τ ), 3p loc GOE (τ ) − 2 and 2p pw GOE (τ ) − 1.The curves exhibit a perfect collapse, as expected from Eq. ( 23).
subtraction with P ∞ is used to make sure that in the case of mobility edges in the spectrum, P (τ ) decays to zero in the thermodynamic limit D → ∞.As discussed above, for certain initial states such as those that make the connection to the SFF, one gets P ∞ = 0 and hence no subtraction is needed in Eq. ( 16).Also, we will show in Sec.VII that at the eigenstate transition in the RP model, at which the eigenstates are fully localized and hence P = P ∞ for site-localized initial states in finite systems, again no subtraction is needed in Eq. ( 17).
Below we illustrate scale invariance of survival probabilities for the studied initial states in the case of GOE Hamiltonians, and in the remainder of the paper, see Secs.V-VII, we numerically test scale-invariant properties at criticality for various other models.

B. Scale invariance of GOE Hamiltonians
The case in which the Hamiltonian is represented by a GOE matrix serves as a benchmark for quantum chaotic dynamics as it allows for establishing analytical predictions [2].Here we provide exact expressions for survival probabilities from all the initial states studied in this paper.More generally, these results will also serve as a motivation in the subsequent sections in the search for scale-invariant dynamics at criticality.
For concreteness, we define the GOE Hamiltonian in analogy to Eq. ( 2), where the matrix elements h i,j are i.i.d.random numbers drawn from a Gaussian distribution with zero mean and variance 2 (1) for the diagonal (off-diagonal) matrix elements.Analytical predictions for the survival probability from initially site-localized states, evolving under the GOE Hamiltonian, were de-  16) (green) compared to the analytical prediction for the SFF k GOE (τ ) from Eq. (21) (dashed black) at D = 2500.We additionally perform the running average to obtain τ -averaged k(τ ) (dashed-dotted blue).(Inset) Details in the comparison of the GOE ramp from Eq. ( 22) to the numerical results that show slight differences in the slope of the ramp.
rived in [27][28][29].While Refs.[27][28][29] focused on the particular initial state, we here adapt the formula therein to describe the survival probabilities of all initial states under consideration.We first express time in units of the corresponding Heisenberg time τ = t/t typ H from Eq. ( 15), which is t typ H = 2 √ D for the GOE Hamiltonian.In this units, Eq. ( 8) from Ref. [29] can be written as where J 1 is the Bessel function of the first kind, and the two-point function b 2 (τ ) [2] is given by The first term on the right-hand side of Eq. ( 18) quickly vanishes at τ > 1 and hence the limit P GOE → P is reached at long times τ > 1.
For the GOE Hamiltonian, the time-averaged survival probabilities P (12) have P ∞ = 0 and γ = 1 for all initial states of interest.Hence, the dependence of P on P can be expressed as the dependence on the coefficient c P .For the site-localized states one expects P = 3/D [27][28][29], while for the infinite-temperature state in Eq. (14) [i.e., the SFF] one expects P = 1/D and for the plane waves one expects P = 2/D.According to Eq. ( 12), this gives: c P = 3 for the site-localized states, c P = 1 for the initial state in Eq. ( 14), and c P = 2 for the plane waves.Then, following Eq.( 17) we rescale the survival probability from Eq. ( 18) to obtain where we assumed (1 The term in Eq. ( 20) that contains the Bessel function is the only term that contains the information about the Hilbert space D, however, it quickly decays to zero since its envelope function is given by 1/(c p 8πD 2 τ 3 ).It is then clear from Eq. ( 20) that at sufficiently long times p GOE becomes scale-invariant, i.e., it only depends on τ for all system sizes.The crossover time to the scaleinvariant dynamics is nowadays known as the Thouless time, which in the case of GOE Hamiltonians scales as . Scale invariance of p GOE (τ ) in the GOE case can be understood as the extension of the scale invariance of the SFF k GOE (τ ) from Eq. ( 16).Specifically, Eq. ( 20) suggest that and hence k GOE (τ ) = p GOE (τ ; c P = 1), giving rise to the well-known expression for the GOE ramp in the SFF, Beyond the emergence of scale invariance at τ > τ T h GOE it is also interesting to consider the relationships between different survival probabilities, which are valid for arbitrary time τ .Let us define p loc GOE (τ ) = p GOE (τ ; c p = 3) for the initial site-localized states and p pw GOE (τ ) = p GOE (τ ; c p = 2) for the initial plane waves.Equation (21) then implies This relation will also be tested in the dynamics at criticality in the next sections.
The results of the above discussion are illustrated in Fig. 2 at D = 2500.In the main panel of Fig. 2 we show k GOE (τ ), p loc GOE (τ ), and p pw GOE (τ ).The SFF k GOE (τ ) exhibits a sharp decay until it hits the GOE ramp at the Thouless time τ GOE Th .In contrast, the other survival probabilities exhibit a broad minimum, which is also referred to as the correlation hole [29].From these results only, one may not be able to exclude the possibility that the Thouless time τ GOE Th depends on the specific type of initial states.However, when one plots the rescaled survival probabilities 3p loc GOE (τ ) − 2 and 2p pw GOE (τ ) − 1 from Eq. (23), see the inset of Fig. 2, one observes a perfect collapse of all survival probabilities suggesting that, as expected, τ GOE Th is independent of the choice of the initial state.
While the results in Fig. 2 are obtained from the analytical expressions, we also test them numerically and compute k(τ ), p pw (τ ) and p loc (τ ) of the GOE Hamiltonians.We first verified (not shown) that the relationship between these quantities, given by Eq. ( 23), also holds for the numerical results.Thus it is sufficient to only focus on the SFF k(τ ). Figure 3 compares the numerical and analytical results k(τ ) and k GOE (τ ), respectively.The numerical results follow the behavior of the analytical result, and we smoothen the noise by carrying out the running averages in time.The inset of Fig. 3, however, highlights that the numerical result differs from the analytical result, since the numerical result lies higher than the analytical one and it follows a slightly different slope of the ramp.We attribute this deviation to the lack of spectral unfolding and filtering in the numerical results for the SFF.In Appendix A, see Fig. 20, we compare the SFFs with and without unfolding and filtering, and we show that in the former case the numerical results fit the analytical predictions.In the reminder of the paper, when the results are compared to the GOE predictions, they are always compared to the numerical results for the GOE Hamiltonians.In Appendix B we also provide further details about the numerical averaging over Hamiltonian realizations and the extraction of the Thouless time τ Th .

C. Arguments for scale invariance at criticality
We now turn our attention to criticality and give arguments for the emergence of scale invariance in quadratic models.Perhaps one of the earliest papers addressing scaling properties at criticality is Ref. [51], where the twoparticle spectral function was investigated.The study revealed a novel scaling form of the two-particle spectral function emerging at the critical point (and coexisting with diffusive behavior), which was conjectured to be a consequence of the fractal structure of eigenstates at criticality.The ansatz of Ref. [51] was then argued to be consistent with the emergence of power-law behavior seen in the survival probabilities from initially localized states, i.e., P loc (t) ∝ t −β [62,64,65,82,83].The studies confirmed that the power-law exponent β is related to the multifractal dimension of the eigenstates at criticality [62,64,65,82,83] even though they did not investigate systematically the system size dependence of P loc (t).While the ansatz of Ref. [51] predicts scale-invariant properties, we note that it is based on properties of eigenstates a within microcanonical window around critical states.General initial states may ultimately be projected to all eigenstates and the scale invariance might be hindered in P loc (t) by the effect of mobility edges.
Our ansatz for the scaled survival probability, introduced in Ref. [49] and in Eq. ( 17) of this paper, generalizes the above considerations.For the initial sitelocalized states it predicts the scale-invariant power-law decay at criticality, where a 0 and β are fitting parameters.In particular, the ansatz from Eq. ( 17) includes two refinements: First, the survival probability is subtracted by P ∞ , which is necessary in order to observe scale invariance in models with mobility edges, such as the Anderson models, and second, time is measured in units of the typical Heisenberg time t typ H .The later is crucial in models such as the Aubry-André model [49], in which the typical and the average Heisenberg time scale differently.These refinements give rise to two important advantages of the ansatz: first, the scale invariance is not only present in the power-law regime in the mid-time dynamics, but also in the late-time dynamics (as in the SFF), and second, it allows for a generalized relationship between the fractal dimension γ and the power-law exponent β, where n determines the scaling of the typical Heisenberg time with the Hilbert-space dimension D. We note that the ansatz of scale-invariant power-law decay as in Eq. ( 24) seems to work also for driven critical systems [104].
We here add two remarks.The first is that that the scale-invariant power-law decay of p loc (τ ) in Eq. ( 24) implies scale invariance of the power-law decay of the subtracted, but unscaled survival probability P loc (t) − P ∞ in the mid-time dynamics.This statement is illustrated by rewriting Eqs. ( 17) and ( 12) as which, using the relationship from Eq. ( 25), yields The second remark is that the reverse of the first remark is not true.This can be illustrated by the power-law decay P loc (t) − P ∞ ∝ t −β , which may emerge in the localized regime in finite systems [49,62].However, the latter do not imply scale invariance of p loc (τ ) in the midtime dynamics [49].Here, the dynamics is dominated by the projection to localized eigenstates, for which the ansatz from Ref. [51] is lost.A seemingly different, but not unrelated aspect is the scale invariance of the SFF at criticality.Early study addressed the level statistics at criticality focusing on short-range statistics [52], which was shown to be neither of Wigner surmise nor of Poisson statistics.The longrange spectral correlation [53,54], measured by spectral rigidity, were further argued to be related to the eigenfunction correlations at criticality [55,57,58].Remarkably, in Ref. [58] the connection between the power-law decay of P loc (t) and the plateau in the SFF K(t) was conjectured.Then, it may not be unexpected, at least for quadratic models, that the scale-invariant power-law decay of p loc (τ ) is accompanied by the scale-invariant plateau in the SFF k(τ ).Recent studies [4,5] investigated system-size dependence of the plateau in the SFF k(τ ) and indeed observed a scale-invariant behavior.

V. RESULTS FOR QUADRATIC MODELS A. From single-particle quantum chaos to localization
We start our analysis of the dynamics with two quadratic models, the 3D Anderson model ( 1) and the PLRB model (2).We first study the dynamics across the critical point at a fixed system size L, and we focus on the unscaled SFF K(τ ) from Eq. ( 8) and the survival probability P (τ ) from Eq. (11).
Figure 4 shows K(τ ) and P (τ ) in the PLRB model at different values of the parameter a while keeping the parameter b fixed.The limit of a = 0 is, up to normalization, the limit of the GOE Hamiltonians, and the numerical simulations fit the corresponding analytical predictions rather well, see Fig. 4(a).Increasing the value of a but keeping it below the transition, see Fig. 4(b), produces the most significant effect to P loc (τ ) while K(τ ) and P pw (τ ) change only mildly.At the transition point a = 1, see Fig. 4(c), P loc (τ ) develops a power-law decay in the mid-time dynamics and K(τ ) develops a plateau in approximately the same time regime.We also note a change in P pw (τ ), which starts deviating from the analytical result at a = 1.In the localized regime, see Fig. 4(d), the deviations from the analytical results further increase.
In Fig. 5 we compare the dynamics in the PLRB model to those in the Anderson model, in which we vary the disorder W .At small W/J = 1 the long-time limit of P pw (τ ) does not approach the limit of GOE but has a higher value, see Fig. 5(a).This can be interpreted as a proximity to the translationally invariant limit where the single-particle eigenstates are plane waves.Increasing the value of W to W/J = 10 within the chaotic regime, see Fig. 5(b), the influence of the translationally invariant limit is weaker and the behavior is similar to the delocalized site of the PLRB model, compare Fig. 4(b) and Fig. 5(b).At the transition point W c , see Fig. 5(c), there is even higher similarity with the PLRB model, compare Fig. 4(c) and Fig. 5(c).P loc (τ ) develops a power-law decay in the mid-time dynamics and K(τ ) develops a plateau in approximately the same time regime.The situation is similar in the localized regime, see Fig. 5(d).
To summarize, while the small a and W limits of both models are distinct, their behavior close to criticality is almost identical.

B. Nearest level spacing statistics
We now focus on the scale-invariant properties at criticality.We first demonstrate them for short-range spectral statistics, i.e., we calculate the ratio of consecutive level spacings of the single-particle spectrum, where δE ν = E ν+1 − E ν is the nearest level spacing.We define the average nearest level spacing ratio (shortly, the average gap ratio) as r = ⟨⟨r ν ⟩ ν ⟩ H , with ⟨...⟩ ν denoting the average over all pairs of nearest level spacings and ⟨...⟩ H denoting the average over 100 Hamiltonian realizations.
While the change of the parameters a and W determines how close the systems are to the critical point, the parameters b (in case of the PLRB model) and the dimensionality d of the hypercubic lattice (in case of the Anderson models) are expected to change the spectral properties and the eigenstate properties at the critical point [55,56].While results for the average gap ratios r of the Anderson models were reported elsewhere (see, e.g., Refs.[4,76]), we here focus on the PLRB models.
The results for r as a function of a are shown in Fig. 6 for several values of b and L. For each value of b there exist a scale-invariant value r at the critical point a = 1, see the crosses in Fig. 6(b).By increasing b the scaleinvariant value of r increases from the Poisson-like value r P = 2 ln 2 − 1 ≈ 0.3863 [105] to the GOE-like value r GOE ≈ 0.5307 [106].The dependence of r on b is summarized in Fig. 19(b) (see below).The change of the spectral properties at criticality is expected to be accompanied with the change of the multifractality of states [79].In particular, the GOE-like value of r is associated with weak multifractality where γ ≈ 1, and the Poisson-like value of r is associated with strong multifractality where γ ≈ 0. We next study an alternative perspective on the degree of multifractality of the critical wavefunctions, ob-  Survival probability p loc (τ ) from Eq. ( 17) in the PLRB model at the critical point a = 1, for different system sizes L. Different panels correspond to different values of the parameter b, ranging from small b = 0.03 in (a) to large b = 10 in (i).The dotted lines are the fits ∝ τ −β from Eq. ( 24), and in each panel we provide the value of β obtained from the fit.tained from the mid-time dynamics.

C. Critical dynamics in the PLRB model
As discussed in Sec.IV C, the survival probability p loc (τ ) from Eq. ( 17) is expected to exhibit a scaleinvariant behavior at criticality, starting from the initial site-localized states.In particular, p loc (τ ) is expected to follow a power-law decay from Eq. ( 24) at midtimes, which is connected to the fractal dimension γ via Eq.(25).
While Ref. [49] demonstrated these properties for the 3D Anderson model, we complement them here by the study of the PLRB model.Figure 7 shows results for p loc (τ ) at the critical point a = 1 for different system sizes L and different values of the parameter b. Results confirm the two expectations raised above: (a) emergence of scale-invariant dynamics at both mid-times and latetimes, and (b) a power-law decay p loc (τ ) ∝ τ −β in the mid-time dynamics, described by Eq. (24).
Different panels of Fig. 7 correspond to different values of the parameter b, and in each panel we list the value of the power-law exponent β obtained from the fit.The later is connected to the fractal dimension γ via γ = nβ (25), where n ≈ 1 for the PLRB model.We observe that β, and hence γ, increases with increasing b, i.e., weaker multifractality is associated with larger b.The dependence of β on b is summarized in Fig. 19

(a) (see below).
We next study survival probabilities from other initial states, specifically, the SFF k(τ ) and the survival prob-ability from initial plane waves p pw (τ ).In Fig. 4(c) we studied the corresponding unscaled quantities K(τ ) and P pw (τ ) at criticality, which were shown to approach at late times (τ → 1) the values K → 1/D and P pw → 2/D.Their values match the predictions of the GOE Hamiltonians studied in Sec.IV B. For the quantities k(τ ) and p pw (τ ), which are of interest here, the GOE case implies the relationship k = 2p pw − 1, as established in Eq. (23).Therefore, in Figs. 8 and 9 we explore to which degree the similarity between k(τ ) and 2p pw (τ ) − 1 extends to the dynamics at criticality.
Figure 8 shows the results at weak multifractality, b > 1, and Fig. 9 shows the results at strong multifractality, b < 1.We observe emergence of a broad plateau in k(τ ), which is marked by the horizontal dotted lines in Figs.8(a) and 9(a).The plateau emerges in the midtime dynamics, in which p loc (τ ), see Fig. 7, exhibits a power-law decay.This was first observed in [49] for the 3D Anderson model and is here complemented by the analysis in the PLRB model.
We extract the two features of k(τ ) from Figs. 8(a) and 9(a), which we then further analyze in Fig. 19 of Sec.VIII: the value of k at the plateau k χ , where χ refers to spectral compressibility [55,56] responsible for the plateau, and the value of k at the Thouless time k(τ T h ).The former characterizes the mid-time dynamics, see also the discussion in Sec.VIII, while the latter characterizes the late-time dynamics.Both values are b-dependent, and they decrease with increasing b.
It is interesting to observe that 2p pw (τ ) − 1 at criticality indeed shares many similarities with k(τ ).This similarity is in particular apparent at weak multifractality, see Fig. 8, for which it persists both in the mid-time and late-time dynamics.This similarity can be interpreted as being a consequence of proximity of the weak multifractal regime to the chaotic behavior.Still, it contains several nontrivial implications, e.g., the extraction of Thouless time can be reliably carried out from 2p pw (τ ) − 1.At strong multifractality, 2p pw (τ )−1 exhibits additional features in the mid-time dynamics, as shown in the insets in Figs.8(b) and 9(b).In particular, we observe a dip below the plateau value k χ , which is also associated with the breakdown of scale-invariant behavior.

D. Critical dynamics in the Anderson models
We complement results for the PLRB model with those for the Anderson models in hypercubes with dimensions d = 3, 4, 5. Figure 10 shows the dynamics of the survival probability from the initial site-localized states p loc (τ ).While the results in d = 3 were already reported in Ref. [49], the results in d = 4 and 5, see Figs. 10(b) and 10(c), respectively, establish generality of scale-invariant mid-time and late-time dynamics at criticality.The exponent β of the power-law decay in the mid-time dynamics, p loc (τ ) ∝ τ −β , decreases with increasing the dimensionality d, which implies stronger    multifractality.The shift towards stronger multifractality in higher-dimensional Anderson models was previously studied from the perspective of level statistics [76] and optical conductivity [71].
In Fig. 11 we then study the survival probabilities from other initial states, specifically, the SFF k(τ ) and the survival probability from the initial plane waves p pw (τ ).The SFF k(τ ), see Fig. 11(a), exhibits a plateau in the midtime dynamics, which is consistent with the results for the PLRB model in Figs.8(a) and 9(a), and with the results for the 3D Anderson model in Ref. [4].We note, however, that scale invariance of the plateau in Fig. 11(a) is reasonably good but not perfect, which we attribute to the absence of spectral unfolding and filtering in the calculation of the SFF k(τ ).In Fig. 21 of Appendix A we show results for the SFF after spectral unfolding and filtering, for which scale invariance of the plateau is further improved.In Fig. 19 (see below) we summarize the behavior of β, k χ in the mid-time dynamics, and r, k(τ T h ) in the late-time dynamics, for the 3D, 4D, and 5D Anderson models.
Finally, motivated by the observation of similarities between k(τ ) and 2p loc (τ ) − 1 in the PLRB model in Figs. 8 and 9, we ask to which extent this similarity translates to the Anderson models.Results for the latter are shown in Fig. 11(b).They share certain similarities with the PLRB model at strong multifractality, see Fig. 9(b).Nevertheless, they do not show convincing evidence neither of scale invariance nor of building a plateau in the mid-time dynamics, at least for the system sizes of investigation.Still, based on the similarity with the PLRB models, for which the latter properties appear to be present at weak multifractality, see Fig. 8(b), we can not exclude that they also appear at strong multifractality in system sizes much larger than those considered here.

VI. RESULTS FOR INTERACTING MODELS A. From many-body quantum chaos to Fock space localization
We now study the dynamics of two interacting models, the QS model ( 3) and the UM model (5).The initial states are many-body states, which belong to the Fock space of dimension D = 2 L+N .Previous study has established that both models exhibit an ergodicity breaking phase transition from many-body quantum chaos at α > α c to Fock space localization at α < α c , with the critical point exhibiting multifractality [85].
We first study the dynamics across the critical point at a fixed system size L, focusing on the unscaled SFF K(τ ) from Eq. ( 8) and the survival probability P (τ ) from Eq. (11). Figure 12 compares the dynamics in both models.The results exhibit strong similarities in all phases, in particular in the mid-time and late-time dynamics.These similarities can be interpreted as the extension of the established similarity between the two models in Ref. [85] to time domain.In the ergodic phase, see Figs. 12(a) and 12(b), K(τ ) and P pw (τ ) are close to the GOE predictions, while P loc (τ ) is not.At criticality, see Figs. 12(c) and 12(d), P loc (τ ) develops a power-law decay and K(τ ) exhibits a plateau.In the localized phase, see Figs. 12(e) and 12(f), properties of the critical point gradually fade away in both models.

B. Similarity of critical dynamics in QS and UM models
We next focus our analysis on the dynamics at criticality, α = α c .As before, we study survival probabilities introduced in Eq. ( 17) from different initial states, p loc (τ ), k(τ ), and p pw (τ ).
The survival probability from the initial Fock-spacelocalized states p loc (τ ) exhibits a scale-invariant powerlaw decay in the mid-time dynamics of both models, shown in Figs.13(a) and 14(a).The exponent β of the power-law decay, however, depends on the specific properties of the models.For the model parameters set in Sec.II, the systems are at moderate to strong multifractality.
The SFF k(τ ) exhibits a scale-invariant plateau in both models, see Figs. 13(b) and 14(b).Scale invariance of the plateau is well converged in the UM model in Fig. 14(b).On the other hand, in the QS model the tendency towards scale invariance is apparent but the convergence is not yet optimal.In Ref. [5] the SFF was studied using spectral unfolding and filtering, giving rise to clearer convergence towards scale-invariant behavior.
For the survival probability from the initial plane waves p pw (τ ) we study its rescaled version 2p pw (τ ) − 1, which is identical to the SFF k(τ ) in the GOE Hamiltonians (23).This relationship also holds for both models at criticality in the late-time dynamics, see Figs. 13(c) and 14(c).As a consequence, the Thouless times, which are very similar in both models, can be extracted either from k(τ ) or from 2p pw (τ )−1.However, at mid-times the validity of the relationship between k(τ ) and 2p pw (τ ) − 1 is less obvious due to the emergence of dips in 2p pw (τ )−1 below the plateau values.

C. Similarity with quadratic models
Apart from the similarity in the mid-time and late-time dynamics in the interacting QS and UM models, there is also an apparent similarity between the interacting models considered here and the quadratic models (the PLRB and Anderson models) studied in Sec.V.The emergence of this similarity is nontrivial provided that in the latter case one considers single-particle states and singleparticle spectrum, while in the former case all quantities are of genuinely many-body origin.
The most prominent features that exhibit similarity in the mid-time dynamics are the scale-invariant powerlaw decay of p loc (τ ) ∝ τ −β and the plateau value k χ of the SFF k(τ ).We characterize them by studying the behavior of β and 1 − k χ in Fig. 19(a) (see below) as a function of the degree of multifractality (see Sec. VIII for details).The results for the interacting and quadratic models exhibit, for each of the measures β and 1 − k χ , a nearly perfect collapse to a single function.
While the above properties hint at the universality in mid-time dynamics, it is interesting note that quantitatively similar behavior can also be observed in the late-time dynamics.In particular, in Fig. 19(b) (see below) we plot 1 − k(τ T h ), i.e., the distance of the SFF k(τ ) at the onset of quantum chaos to its long-time average, and again observe a nearly perfect collapse.This suggests that fingerprints of criticality may be encoded in both the mid-time as well as the late-time dynamics.24), with β = 0 in all panels.Results are shown for different linear system sizes L = 1000, 2500, 5000, 10 000, 20 000, which perfectly overlap.els.Here we show that the RP model exhibits a slightly different behavior that, although consistent internally, features distinct mid-time critical dynamics at the eigenstate transition.

A. Transitions in the RP model
Even though we defined the RP model in Sec.II C as a quadratic model with random on-site potentials and all-to-all hoppings, the model is different from the other models considered so far since it is expected to exhibit two transitions, at c = 1 and c = 2.In the context of this paper we argue that only the second transition (at c = 2) can be considered as the eigenstate transition with scale-invariant critical dynamics.At this transition, the short-range spectral correlations change from the GOE statistics to the Poisson statistics at any λ [102].Here we are mostly concerned with the properties of long-range spectral correlations [24,97,103,107].
We first study the dynamics at and in the vicinity of the first transition at c = 1.In Fig. 15, we show an example of the unscaled survival probability P (τ ) from Eq. (11) as well as the SFF K(τ ).We vary the parameter c while keeping the parameter λ fixed.Figure 15(a) shows that the limit of c = 0 is, up to normalization, the limit of the GOE Hamiltonians, and the numerical simulations fit the corresponding analytical predictions very well.Increasing the value of c within the range c ≤ 1, see Figs. 15(b) and 15(c), gives rise to the visible differences in P loc (τ ), while K(τ ) and P pw (τ ) remain virtually identical.Above the transition point to the non-ergodic extended phase, i.e., at c > 1 shown in Fig. 4(d), the deviations of K(τ ) and P pw (τ ) from the analytical results become more apparent.
Next we study the dynamics at and in the vicinity of the second transition at c = 2.By increasing c toward the value c = 2, the results for K(τ ) and P pw (τ ) further deviate from the GOE predictions, see Figs. 16(a)-16(c).In particular, the agreement between K(τ ) and the GOE predictions is only preserved at long times close to the Heisenberg time τ H = 1.It is interesting to note that by increasing c, P pw (τ ) approaches the time dependence of K(τ ), and at c > 2, see Fig. 16(d), they indeed become very similar.This excludes the validity of the relationship k(τ ) = 2p pw (τ ) − 1 from Eq. ( 23), which is further confirmed in the next section where we study scale invariance of the dynamics at c = 2.
Another important aspect concerns the mid-time dynamics of K(τ ). Figure 16(c) shows that at c = 2, which we refer to here as the eigenstate transition, K(τ ) develops a plateau, similarly to the critical points for eigenstate transitions in other models under considerations.However, the plateau value is maximal in the sense that it equals the long-time average.As discussed in the next section, this is consistent with the absence of the decay of P loc (τ ) in the mid-time dynamics, which is indeed observed in Fig. 16(c).

B. Critical dynamics in the RP model
We now focus on the dynamics at the critical point c = 2, and we study the survival probabilities p loc (τ ) in Fig. 17, and k(τ ), p pw (τ ) in Fig. 18.The most important property of the critical point at c = 2 is that Hamiltonian eigenstates are localized [97,103], in contrast to all other critical points considered before.This has implications both for the time dependence in the mid-time dynamics, as well as on the scale invariance since it requires a refinement of the definition of the survival probability p loc (τ ).
As a consequence of localization at the critical point one needs to consider Eq. ( 12) with some care.For the initial site-localized states, we observe P loc = P loc ∞ for any finite system under consideration.Then, there is no need to apply Eq. ( 17), but a simpler rescaling, p loc (τ ) = P loc (τ ) which is analogous to the one applied to the SFF in Eq. ( 16), is sufficient to observe scale invariance.The resulting p loc (τ ) is plotted in Fig. 17 at c = 2 and several values of λ, exhibiting perfect scale invariance.Another important feature of Fig. 17 is the absence of any decay of p loc (τ ) in the mid-time dynamics for all values of λ.This can be interpreted as a power-law decay with β = 0, and hence the fractal dimension is γ = 0, which is consistent with localization.A manifestation of this property is the emergence of a plateau in k(τ ) in the mid-time dynamics, shown in Fig. 18(a), which equals k χ = 1.We note that the plateau is more clearly visible in the unfolded and filtered version of the SFF k(τ ) shown in Fig. 22 of Appendix A. A certain tendency towards scale invariance is also manifested in the dynamics of 2p pw (τ ) − 1 shown in Fig. 18(b), however, as already in the context of Fig. 16(c), the time evolution of the later quantity can not be related to k(τ ).
The main properties of the RP model are summarized in the next section, see Fig. 19, and compared to other models.For the long-time dynamics, one can show that the RP model complies with the universal features observed in other models, see the next section for details.However, the mid-time dynamics of the RP model is not consistent with the universal properties of other models.We attribute localization of the Hamiltonian eigenfunctions as the main source of deviations, since all the other models exhibit multifractality at criticality.

VIII. MID-TIME VS LATE-TIME DYNAMICS
We summarize our findings by combining the two perspectives on the properties of the critical points, one obtained from the mid-time dynamics and another from the late-time dynamics.The mid-time dynamics is characterized by the power-law exponent β and the distance of the SFF long-time average to the plateau, 1 − k χ , while the long-time dynamics is characterized by the normalized average gap ratio at criticality, and the distance of the SFF at the onset of quantum chaos (i.e., at τ = τ Th ) to its long time average, 1 − k(τ Th ).The goal is to study all these quantities as a function of a single parameter that determines the degree of multifractality of the critical wavefunctions.We choose this parameter to be b from the PLRB model, see Eq. ( 2).For the other models, we determine the corresponding value of b by requesting the gap ratio r N of that model to match the gap ratio of the PLRB model.Results are shown in Fig. 19.For the late-time dynamics, see Fig. 19(b), the curve r N versus b is smooth by construction.It is interesting to observe that also the distance of the SFF long time average to the SFF at the onset of quantum chaos, 1 − k(τ Th ), appears to approach a well-defined function of the parameter b.This observation is in particular important since we consider both quadratic and interacting models within the same framework.
For the mid-time dynamics, see Fig. 19(a), we observe that both the power-law exponent β and the distance of the SFF to the plateau 1−k χ appear to be monotonously increasing functions of b.Even though the quantities are extracted purely from the dynamics, the results for quadratic models are fully consistent with the findings of Refs.[55,56], where similar quantities (fractal dimension and spectral rigidity) were extracted from the analysis of energy eigenspectra and properties of eigenvectors.Again, our main result is that also the critical properties of interacting models may belong to the same singleparameter function.
The exception to the observed universal properties of the mid-time and late-time dynamics is the RP model.Figure 19(a) shows that both β and 1 − k χ deviate from the single-parameter function, since β = 1 − k χ = 0 in all cases under consideration.In fact, since the wavefunctions at criticality are localized, they should rather correspond to the b = 0 limit of the PLRB model.This does not match with the properties of other models for which the average gap ratio is given by the Poisson statistics if critical point is localized.The observed deviations only hint at a lack of the correspondence between the mid-time and late-time dynamics in the RP model.

IX. CONCLUSIONS
In this paper we studied quantum dynamics of survival probabilities from different initial states, including those that define the SFF, and we provided answers to the two questions posed in the Introduction.The most important statement is that there exist the notion of scale invariance in quantum dynamics at criticality.Remarkably, the scale invariance is observed in quantities that bear direct analogies with the well-established measures of quantum chaos.Another important statement is that scale invariance at criticality may be observed at times that are much shorter than those required to detect quantum chaotic dynamics.These properties open possibilities to detect fingerprints of criticality at experimentally accessible time scales (in the so-called mid-time dynam-ics), which are much shorter than the longest time scales of finite systems, characterized within the notion of the Heisenberg time.
While universality of the chaotic dynamics such as the emergence of a ramp in the SFF is expected for all quantum chaotic models, understanding of the universal features at criticality is far from understood.The survival probabilities introduced here and in the preceding Letter [49] allow for establishing similarity of quantum dynamics in different models.In particular, we established the link within two classes of models: for quadratic models, between the Anderson model and the PLRB model, for which connections have often been explored in the past [37], and for interacting models, between the QS model and the UM model, for which connections have been noted only recently [85].Most intriguingly, however, we also established similarity in quantum dynamics between the quadratic and interacting models.One of the possible outcomes of this result is that ergodicity breaking phase transitions in the interacting models under consideration may be described as Anderson transitions in the corresponding Hilbert spaces.
A common feature of all the models for which similarity in the critical dynamics has been established is the wavefunction structure, which is multifractal in the corresponding Hilbert spaces.This motivated us to quantify features of both mid-time and late-time dynamics, and express them in terms of a single parameter that quantifies the degree of multifractality.Results shown in Fig. 19 suggest that indeed a universal, single-parameter description of the main dynamical features is possible.In particular, it allows for understanding the mid-time dynamics via the late-time dynamics, and vice versa.We also showed that the RP model represents an exception from the latter property, and we attributed this by the lack of multifractal structure in the RP model at the critical point.22), while k(τ ) better agrees with the numerical results for the GOE Hamiltonians.At criticality, see Figs. 20(c) and 20(d), k(τ ) exhibits a clearer plateau in the mid-time dynamics, however, it also exhibits a dip below the plateau at short times.We also note that the plateau value of k(τ ) emerges at slightly higher value when compared to the plateau value of k(τ ).On the localized side, see Figs. 20(e) and 20(f), the results are again similar, apart from the dips at short times that emerge in k(τ ).
Figure 21 extends the analysis of k(τ ) at criticality to the 4D and 5D Anderson models.This result should be compared to Fig. 11(a), in which we show k(τ ).The tendency for the emergence of a scale-invariant plateau in the mid-time dynamics is clearer in k(τ ) in Fig. 21.However, k(τ ) also exhibits some deviations when compared to the analytical prediction for the GOE ramp at late times.
We also study the impact of unfolding and filtering on the SFF in the RP model.Figure 18(a) shows k(τ ) at criticality, which exhibits a tendency toward forming a plateau in the mid-time dynamics at k χ = 1. Figure 22 shows that the emergence of a plateau can be demonstrated more convincingly in k(τ ).On the other hand, the agreement between k(τ ) and the analytical prediction for the GOE ramp at late times τ ≈ 1 is again less accurate.
FIG. 1. Concept of mid-time and late-time dynamics at eigenstate transition, for the dynamics of the SFF k and survival probability p of initially site-localized states.We measure time t in units of typical Heisenberg time t typ H , such that τ = t/t typ H , and we normalize k and p such that they approach k = p = 1 at τ > 1.The late-time dynamics is defined as the dynamics in between the scaled Thouless time τ Th = t Th /t typ H and the scaled typical Heisenberg time τ typ

FIG. 3 .
FIG.3.Numerical results of the SFF k(τ ) from Eq. (16) (green) compared to the analytical prediction for the SFF k GOE (τ ) from Eq. (21) (dashed black) at D = 2500.We additionally perform the running average to obtain τ -averaged k(τ ) (dashed-dotted blue).(Inset) Details in the comparison of the GOE ramp from Eq. (22) to the numerical results that show slight differences in the slope of the ramp.

5 FIG. 4 .
FIG. 4. K(τ ), P loc (τ ) and P pw (τ ) in the PLRB model at b = 0.5 and D = L = 20 000.(a) a = 0, (b) a = 0.75, (c) the critical point a = 1, and (d) a = 1.5.Green lines are results before time averaging.The corresponding results for GOE Hamiltonians from Sec. IV B (see Fig. 2) are shown as dashed black lines.The dotted lines in (c) indicate the power-law decay for P loc (τ ) and the plateau in K(τ ).

20 FIG. 5 .
FIG. 5. K(τ ), P loc (τ ) and P pw (τ ) in the 3D Anderson model at D = L 3 = 32 768.(a) W/J = 1, (b) W/J = 10, (c) the critical point W/J = 16.54, and (d) W/J = 20.Green lines are results before time averaging.The corresponding results for GOE Hamiltonians from Sec. IV B (see Fig. 2) are shown as dashed black lines.The dotted lines in (c) indicate the power-law decay for P loc (τ ) to a nonzero constant and the plateau in K(τ ).

FIG. 6 .
FIG. 6.Average gap ratio r in the PLRB model as a function of a, for several values of b and L. The results in (b) are identical to those in (a), but zoomed into a narrow window around the critical point a = 1.The scale-invariant values of r at the critical point in (b) are denoted by the crosses (from bottom to top: b = 0.01, 0.05, 0.1, 0.25, 0.5, 1).
FIG. 7.Survival probability p loc (τ ) from Eq. (17) in the PLRB model at the critical point a = 1, for different system sizes L. Different panels correspond to different values of the parameter b, ranging from small b = 0.03 in (a) to large b = 10 in (i).The dotted lines are the fits ∝ τ −β from Eq. (24), and in each panel we provide the value of β obtained from the fit.

2 FIG. 8 .
FIG. 8. (a) The SFF k(τ ), and (b) the survival probability from the initial plane waves 2p pw (τ ) − 1, in the PLRB model at criticality, a = 1.Results are shown for several values of b > 1 at weak multifractality, and different system sizes L, as indicated in the legend [the same parameters as in Figs.7(f)-7(i)].The horizontal dotted lines denote the plateau values k χ , and the dashed-dotted lines denote the linear ramp of the numerically evaluated GOE Hamiltonians, cf.Fig. 3.The circles correspond to the extracted Thouless times τ Th , at which k(τ ) approach the numerically evaluated GOE values.The inset in (b) shows a detail of the additional short-time feature of 2p pw (τ ) − 1 compared to k(τ ).

3 FIG. 9 .
FIG. 9. (a) The SFF k(τ ), and (b) the survival probability from the initial plane waves 2p pw (τ ) − 1, in the PLRB model at criticality, a = 1.Results are shown for several values of b < 1 at strong multifractality, and different system sizes L, as indicated in the legend [the same parameters as in Figs.7(c)-7(e)].All other features are identical to those in Fig. 8.

FIG. 10 .
FIG. 10.Survival probability p loc (τ ) from Eq. (17) in the Anderson models at the critical points, for different linear system sizes L. Different panels correspond to different dimensionalities d of the hypercubic lattice: (a) d = 3 (the critical point W/J = 16.54),(b) d = 4 (the critical point W/J = 34.5),and (c) d = 5 (the critical point W/J = 57.5).The dotted lines are the fits ∝ τ −β from Eq. (24), and in each panel we provide the value of β obtained from the fit.

k 5 L = 4 , 5 , 6 , 7 , 8 L = 6 , 8 , 10 , 12 , 14 L = 12 4 FIG. 11 .
FIG. 11.(a) The SFF k(τ ), and (b) the survival probability from the initial plane waves 2p pw (τ )−1, in the Anderson models at the critical points, for different linear system sizes L. Results are shown for different dimensionalities d of the hypercubic lattice: d = 3 (the critical point W/J = 16.54),d = 4 (the critical point W/J = 34.5),and d = 5 (the critical point W/J = 57.5).The horizontal dotted lines denote the plateau values k χ , and the dashed-dotted lines denote the linear ramp of the numerically evaluated GOE Hamiltonians, cf.Fig. 3.The circles correspond to the extracted Thouless times τ Th , at which k(τ ) approaches the numerically evaluated GOE values.The inset in (b) shows a detail of the additional short-time feature of 2p pw (τ ) − 1 compared to k(τ ), which approaches a plateau.

11 FIG. 13 . 14 FIG. 14 .
FIG. 13.Survival probabilities in the QS model at criticality, α = αc = 0.715, for several system sizes L with the Fock-space dimension D = 2 N +L and N = 5.(a) p loc (τ ) from the initially Fock-space-localized states, (b) the SFF k(τ ), and (c) 2p pw − 1 from the initial plane waves.The dotted line in (a) denotes the power-law fit from Eq. (24) and the horizontal dotted lines in (b) and (c) denote the plateau values k χ .The thick dashed-dotted lines in (b) and (c) denote the numerical results for the SFF of the GOE Hamiltonians, cf.Fig. 3.The black circles correspond to the extracted Thouless times τ Th .

FIG. 18 .
FIG. 18.(a) The SFF k(τ ), and (b) the survival probability from the initial plane waves 2p pw (τ ) − 1, in the RP model at the critical point c = 2, for different linear system sizes L. Results are shown for different values of λ, as indicated in the legend.The horizontal dotted line in (a) denotes the plateau value k χ = 1, and the horizontal dotted lines in (b) denote the estimates of the plateau values for 2p pw (τ ) − 1.The dashed-dotted lines in both panels denote the linear ramp of the numerically evaluated GOE matrices, cf.Fig.3.The circles correspond to the extracted Thouless times τ Th , at which the quantities approach the numerically evaluated GOE value.

FIG. 19 .
FIG. 19.Characterization of scale-invariant properties at criticality in the (a) mid-time dynamics, and (b) late-time dynamics, see also Fig. 1.(a) The power-law decay exponent β (red) and the distance of the SFF long-time average to the plateau, 1 − kχ (blue).(b) The distance of the SFF long-time average to the SFF at the onset of quantum chaos, 1 − k(τ Th ) (blue), and the normalized average gap ratio rN from Eq. (30) (black).The crosses (×) are results for the PLRB model, as a function of the parameter b.For all other models, we assign the corresponding value of b by matching the values of rN in the studied model with those in the PLRB model.This yields the rN versus b curve smooth by construction.Results for the 3D, 4D and 5D Anderson models are marked with circles (•), for the QS model with squares (□), for the UM model with right triangles (▷), and for the RP model with left triangles (◁).

FIG. 21 .
FIG.21.The SFF k(τ ) after spectral unfolding and filtering, see Eq. (A1), at the critical points of the Anderson models on 3D, 4D and 5D hypercubic lattices, shown at different system sizes L. The dotted lines denote the plateau values of k(τ ).The dashed line denotes the analytical expression for the GOE ramp from Eq. (22).

Figure 20
Figure20compares the unfolded and filtered k(τ ) in the 3D Anderson model, studied in[4], with k(τ ) studied in the main text.Results are qualitatively very similar.In the chaotic regime, see Figs.20(a) and 20(b), we observe that in the late-time dynamics after the Thouless time, k(τ ) follows the analytical prediction of the GOE from Eq. (22), while k(τ ) better agrees with the numerical results for the GOE Hamiltonians.At criticality, see Figs.20(c) and 20(d), k(τ ) exhibits a clearer plateau in the mid-time dynamics, however, it also exhibits a dip below the plateau at short times.We also note that the plateau value of k(τ ) emerges at slightly higher value when compared to the plateau value of k(τ ).On the localized side, see Figs.20(e) and 20(f), the results are again similar, apart from the dips at short times that emerge in k(τ ).Figure21extends the analysis of k(τ ) at criticality to the 4D and 5D Anderson models.This result should be compared to Fig.11(a), in which we show k(τ ).The tendency for the emergence of a scale-invariant plateau in the mid-time dynamics is clearer in k(τ ) in Fig.21.However, k(τ ) also exhibits some deviations when compared to the analytical prediction for the GOE ramp at late times.We also study the impact of unfolding and filtering on

FIG. 22 .
FIG.22.The SFFs k(τ ) from Eq. (16) (solid blue line) and k(τ ) from Eq. (A1) (dashed blue line) in the RP model at c = 2, λ = 1 and L = D = 20000.The dotted lines denote the plateau values of k(τ ) and k(τ ), which are both equal to one.The dashed black line denotes the analytical expression for the GOE ramp from Eq. (22), while the dashed-dotted black line denotes the numerical result for the ramp in the GOE Hamiltonians, cf.Fig.3.