Sub-diffusive Thouless time scaling in the Anderson model on random regular graphs

The scaling of the Thouless time with system size is of fundamental importance to characterize dynamical properties in quantum systems. In this work, we study the scaling of the Thouless time in the Anderson model on random regular graphs with on-site disorder. We determine the Thouless time from two main quantities: the spectral form factor and the power spectrum. Both quantities probe the long-range spectral correlations in the system and allow us to determine the Thouless time as the time scale after which the system is well described by random matrix theory. We find that the scaling of the Thouless time is consistent with the existence of a sub-diffusive regime anticipating the localized phase. Furthermore, to reduce finite-size effects, we break energy conservation by introducing a Floquet version of the model and show that it hosts a similar sub-diffusive regime.

Although MBL has been extensively studied, many of its aspects remain puzzling and they are still under intense debate. For instance, the mechanism and nature of the transition is still unclear (see also the discussion L = 6 FIG. 1. Sketch of a random regular graph (RRG) with diameter L = 6 and N = 2 L = 64 nodes and branching number K = 2, i.e. each node has K + 1 = 3 neighbors. The particle can hop along the edges (black lines) of the graph with hopping amplitude t = 1 and is subject to on-site disorder µi (color intensity) on the nodes of the graph (red dots of different intensity).
in [45]). This is for instance exemplified by early scaling attempts which yield critical exponents in conflict with general constraints, i.e., so-called Harris bounds [46][47][48][49][50][51][52][53][54]. The Fock-space structure and the multifractal nature of the eigenstates close to the MBL transition within the ergodic side have been at the center of a recent spur of research activity [24][25][26][27][28]30]. In the ergodic phase, numerical simulation has found sub-diffusive transport at intermediate time scales [32,33,[55][56][57][58][59][60][61][62][63][64], contrary to the expected diffusive behaviors of a metal. The sub-diffusion has been argued to stem from the existence of rare Griffith regions which suppress transport. However, the existence of these regions and of the corresponding subdiffusive behavior as a phase, separated from the diffusive arXiv:2201.04673v2 [cond-mat.dis-nn] 14 Sep 2023 one, is still under intense investigation and it is not clear whether a phase transition between diffusion and subdiffusion exists. Furthermore, several recent works have been questioning the existence of a genuine MBL transition [65][66][67]. In these works, using exact diagonalization (ED), a systematic flow towards the ergodic regime has been observed with increasing system sizes, which would imply that MBL is not a stable phase of matter. Subsequent works have pointed out however that this result might be due to the limitation in system sizes reached using ED [68][69][70], or by a large finite-size crossover regime before the MBL phase [45,71,72]. For instance, also systems having a firmly established metal-insulator transition, such as the Anderson model on random regular graphs (RRGs) (locally tree-like graphs without boundaries, see Fig. 1) [68,70], present similar finite-size corrections compared to the MBL problem in one dimension.
Settling this controversy is hindered by the exponential scaling of the Hilbert space with the size of the system and other approaches to make progress are therefore needed. Given these difficulties, a promising route is to focus on more tractable models, that reproduce the main features of MBL systems. Following the original idea of mapping a disordered quantum dot to a localization problem in the Fock space [73], the problem of Anderson localization in hierarchical structures [74], such as RRGs, has been suggested to be a useful proxy to describe MBL systems [75,76]. In this spirit, the sites of the RRG are interpreted as Fock space basis states of the non-interacting model, the on-site energies as renormalized random fields/potentials, and the hopping between sites as an interaction connecting non-interacting states. The reduction of the (one dimensional) MBL problem to the Anderson model on the RRG is a simplification on several accounts: (i) the structure of the Fock space is a hypercube with extensive connectivity of nodes, while in RRG models the connectivity is fixed and intensive; (ii) the disorder in real space translates to correlated disorder in Fock space, while in the RRG model the on-site energies are uncorrelated; (iii) the hypercube has many short-range loops unlike RRG, but for both structures the typical loop size is a finite fraction of the graph diameter.
Furthermore, the Anderson model on the RRG is not just a proxy for the interacting problem but exhibits a rich and interesting phenomenology on its own. Recently, it has been argued that the Anderson model on the RRG might host a new intermediate phase in between the ergodic and the localized one. This phase, dubbed nonergodic extended (NEE) phase [75][76][77][78][79], might be composed by critical/multifractal states, i. e., states having strong space fluctuations. Thus, unlike the Anderson model on the hypercubic lattice Z d with d > 2, where multifractal states appear only at the localization transition point, in the RRG an entire phase composed by critical states might exist [99]. However, several works have pointed out that this intermediate phase in RRG might be merely a finite-size effect and therefore ergodicity would be restored in thermodynamic limit [84,92,[94][95][96]98]. Another intriguing suggestion is the possible existence of a sub-diffusive phase. By inspecting the spread of a particle initially localized in one of the sites of the RRG, Refs. 88-91 have provided evidence that the transport at finite time scales could be sub-diffusive. This sub-diffusive propagation has been found for a range of disorder strengths within the extended phase and longtime scales. This propagation should be compared with the behavior of the Anderson model on the hypercubic Z d lattice, where sub-diffusion happens only at the critical point.
Our work aims to shed light on the existence of subdiffusion in the Anderson model on the RRG solely from the spectral perspective. Our main probes are the spectral form factor and the power spectrum. Both, the spectral form factor [65,68,70] and the power spectrum [100][101][102][103][104] are measures of long-range correlations between eigenvalues of the Hamiltonian. Importantly, these two measures are efficient probes for the dynamics in the system. The Thouless time, t Th , defined as a time scale beyond which the system dynamics is universal and described by the random matrix theory of Gaussian ensembles [105] (see Fig. 2) can be extracted from the spectral form factor and the power spectrum. The Thouless time can be interpreted as the time that a particle needs to diffuse throughout the system. As a result, depending on the scaling of t Th ∼ L 1/β with system size L (diameter of the graph in the case of RRG), one can define different kinds of propagation based on the scaling exponent β, ranging from super-diffusion, β > 1, and diffusion, β = 1 to sub-diffusion, 0 < β < 1, and localization, β → 0. In particular, the Thouless time should be compared to the Heisenberg time, t H ∼ 2 L , which is the largest meaningful time scale and is defined as the inverse level spacing of the energy spectrum. The system is localized if t Th /t H goes to a finite constant in thermodynamic limit. By analyzing the scaling of t Th with system size L, we confirm the existence of a sub-diffusive regime and discuss its finite-size effects. To better pin down the sub-diffusive behavior by reducing finite-size effects, we allow energy fluctuations by introducing a Floquet version of the Anderson model on the RRG. In agreement with the static model, we find also for the Floquet model sub-diffusive transport for a wide range of parameters anticipating the Anderson transition.
The rest of the work is organized as follows. In Sec. II we define the static and Floquet models and discuss their phase diagrams. In Sec. III, we introduce the two dynamical probes, the spectral form factor, and the power spectrum. The main results of our work are presented in Sec. IV and Sec. V contains concluding remarks.

II. MODELS
We consider a single particle hopping on a RRG with N sites, subject to on-site disorder {µ i }. The Hamiltonian is given by where {µ i } are independent random variables with box distribution [−W/2, W/2] and the sum in the second term runs over the edges {i, j} of the graph, i.e. the set of connections between sites in a given realization of the RRG (cf. Fig. 1 for an example, where edges are indicated by lines). Each site of an RRG has a fixed number of neighbors K + 1, where K is the branching number, while the precise configuration of edges in the graph is subject to random sampling. In this work, we focus on the smallest nontrivial branching number K = 2 different from a onedimensional problem, such that each site has 3 neighbors. The total number of nodes in the graph is N = 2 L but we refer to it via the diameter of the graph L ∼ ln N in order to make the analogy with the Hilbert space of 1/2-spin chains of length L. This model has an Anderson localization transition at W AT ≈ 18.1 [76,79,81,96], which is blurred by a finite size crossover regime and appears at lower disorder in finite graphs [45,84,92,[94][95][96]98]. For instance the average level spacing for L = 17 has crossings close to W ≈ 16 (see Fig. 3).
The Hamiltonian in Eq. (1) hosts a single-particle mobility edge, which separates extended from localized states as a function of energy and its density of states has a non-trivial shape peaked close to the edges of the spectrum (see Appendix A and cf. [106]). In order to overcome finite-size effects from the localized energy bands and from the energy dependence of the density of states, we introduce a Floquet version of the Anderson model on RRG where H 1 and H 2 are hopping Hamiltonians as in Eq. (1) on the same RRG but with different diagonal disorder {µ i }. 2T is the driving period and we set T = 2 for the rest of the work. The results are not qualitatively affected by the concrete value of T , see Appendix C for details. Both models, static and periodic, have transitions from Wigner-Dyson to Poisson level statistics as can be seen in the average gap ratio r in Fig. 3. The Anderson localization transition in the static case, Eq. (1), has been extensively studied [74-79, 81, 83-86, 91-98] and we expect the transition in the driven model to have similar features. The advantage of introducing a Floquet drive is the flat density of states on the unit circle, which allows us to perform an unbiased study of long-range spectral correlations.
In the following we refer the eigenvalues {E n } of the static model, Eq. (1), as "spectrum" and {θ n = arg(ω n )} of the driven model, Eq. (2), as the eigenphases or quasienergies, where {ω n = e iθn } are the eigenvalues of U , which lie on the complex unit circle. We expect that the Floquet model (2) yields cleaner results for correlations in the spectrum at large phase differences between the eigenvalues, which is confirmed by our numerical results.

III. METHODS
In this section, we introduce the quantities used to study spectral correlations and detail the computation of the Thouless time.

A. Level spacing
We start with a well known probe to detect short range spectral correlations, which are captured by the statistics of spacings between two adjacent eigenvalues The spacing ratio, defined as has been found to be a useful resource to separate a delocalized phase from a localized one [9,107]. In an ergodic phase, the spacings between eigenvalues are distributed like the ones of a Gaussian random matrix. In contrast, in a localized phase energy levels tend to cross each other as a function of a parameter with little interaction between eigenvalues and the spacings are thus Poisson distributed. In particular, for the Gaussian orthogonal ensemble (GOE) the average spacing ratio is given by the Wigner-Dyson value ⟨r⟩ GOE ≈ 0.53590, while for Poisson spectra it is ⟨r⟩ P oisson = 2 ln 2 − 1 ≈ 0.38629 [9,107] where ⟨...⟩ indicates averaging over disorder realization of graph and diagonal disorder as well as eigenstates. The averaged energy spacing is the smallest energy scale in the problem and it defines the Heisenberg time, Consequently, the Heisenberg time t H is the largest meaningful time scale in a finite system, on which discrete energy levels can be resolved. It is proportional to the dimension of the Hilbert space t H ∝ 2 L .

B. Spectral unfolding
The r-statistics defined in Eq. (4) provides only limited insight into the dynamics of the system because it only involves the computation of small energy scales (asymptotically large times). As a result, the r-statistics can separate only a delocalized phase from a localized one, without providing information about the transport properties.
To inspect finite time scales relevant for the transport, we have to consider the long-range spectral correlations. However, before defining the two main measures to analyze energy correlation, we first need to introduce the notion of spectrum unfolding. To remove the influence of the non-uniform density of states in the spectra of the static Hamiltonian, Eq. (1), we perform an unfolding of the spectrum, which maps the eigenvalues {E n } of the Hamiltonian with non-uniform density of states to the "unfolded" eigenvalues {ε n }, which have a homogeneous density of states.
Instead of working directly with the density of states, we use the cumulative distribution function (CDF), which is defined by the fraction of eigenvalues smaller than E: where for a given energy E, #(E n < E) ∈ N is the number of eigenvalues smaller than E. For each disorder realization, this function is the empirical CDF, a stepwise function with steps of size 1/N at the positions of the eigenvalues {E n }.
In order to obtain the average CDF over n real disorder realizations, we combine the spectra {E Each eigenvalue E (i) n of a realization (i) can then be mapped to its unfolded and normalized counterpart where the density of states of the unfolded ε (i) n is constant. The procedure is illustrated in Fig. 8 in Appendix A. For measuring the Thouless time in units of the mean level spacing, it is more convenient to work with an equidistant spectrum with unit mean level spacing, thus we work with the unfolded spectrum: Note that in the case of the Floquet model defined in Eq.
(2), spectral unfolding is not needed and the simple rescaling ε n /2π yields a uniform spectrum with unit level spacing similar to ε in Eq. (8). Henceforth the realization index (i) will be dropped from the unfolded spectrum {ε (i) n } and we use only {ε n } to denote the n-th level of the spectrum, making the disorder index implicit.

C. Spectral form factor
The first probe that we introduce is the spectral form factor, defined as [105] K(τ ) = 1 where {ε n } is the Hamiltonian unfolded spectrum or the Floquet quasienergy spectrum. The function is a filter for mitigating spectral-edge effects in the static unfolded spectra.ε is the center of the spectrum for a given disorder realization, σ is the width of the spectrum and η is a parameter which controls the relative width of the filter. To avoid finite size-effects from the localized energy-bands we focus on only 60% of the states centered at the middle of the spectrumε. In addition, we apply the Gaussian filter mentioned above and set η = 0.3 through the entire work, except where specified otherwise. The definition of spectral form factor given in Eq. (9) guarantees K(τ ) = 1 at asymptotic large times τ ≈ τ H = 1.
The spectral form factor for GOE spectra is described by a linear growth, K GOE (τ ) ≈ 2τ − τ ln(1 + 2τ ) [108] for τ < 1, while for a Poisson spectrum, K Poisson (τ ) ≡ 1. The Thouless time τ Th identifies the time scale after which the dynamics is described by random matrix theory. In particular, for the times τ < τ Th the quantum dynamics is governed by the locality of the system, quantum correlations spread dynamically before reaching the boundary of the system. Thus, in the way how τ Th scales with L, it is possible to probe different classes of the propagation of correlations, e. g., diffusion, sub-diffusion or localized states. At later times, τ > τ Th , quantum correlations have been scrambled to a large extent across all length scales, and thus at such late times the system becomes indistinguishable from a random matrix model, leading to the GOE behavior of the spectral form factor. Hence, the Thouless time τ Th corresponds to the time beyond which K(τ ) ramps linearly and can be approximated with the behavior of K GOE (τ ), i.e. within a certain threshold (K(τ > τ Th ) ≈ K GOE (τ )). Figure 2 shows the typical behavior of K(τ ) for a fixed system size and several disorder strengths. As one can observe in Fig. 2, it is possible to identify a time scale τ Th , from which K(τ ) follows the GOE curve (dashed line). In Appendix B the procedure for extracting the Thouless time is explained in detail.
To analyze the scaling of the Thouless time with system size, we rescale it with the actual Heisenberg time in Eq. (5): As we will discuss later in Sec. III E, for a diffusive system in a tree-like structure we expect t Th ∝ L like in the classical diffusion problems on Bethe lattice [109] or RRG [110], while in the localized phase, the spectrum is Poissonian and thus K Poisson (τ ) ≡ 1 and t Th = t H ∝ 2 L .

D. Power Spectrum
A complementary measure to detect long-range level correlations is the spectral power spectrum. It can be defined through the spectral statistic {δ n } given by with n ∈ {1, ..., N } and the Hilbert space dimension N = 2 L . For both models, static and Floquet, ⟨ε n ⟩ = n. Consequently, {δ n } measures the distance of the n-th level from a rigid equidistant spectrum. The sequence {δ n } can be interpreted as a discrete "time" series with zero mean and the index n as "timepoint". A useful way to analyze this time series {δ n } is to consider its Fourier transform, the power spectrum: As for the spectral form factor, we have introduced the filter function (cf. Eq. (10)) to reduce finite-size energy edge effects for the Hamiltonian model. N = n g 2 n is the normalization constant of the Fourier transform in the filtered spectrum. Like for the spectral form factor,ε and σ 2 are the mean and variance of the unfolded spectrum for each disorder realization. The strength of the filtering is set to η = 0.3. For reducing edges effects even further, the edges of the spectrum are cut off and only 60% of the spectrum centered at the middle of the spectrumε. 1 ≤ k ≤ N y is integer, with N y = N/2 being the largest possible meaningful frequency in the system, called the Nyquist frequency.
For Poisson spectra, when there are no correlations at any range in the spectrum this time series is similar to a sample of a Brownian motion with displacement δ n . In the limit k/N ≪ 1 and N ≫ 1, the asymptotic form of the power spectrum is given by P Poisson (k) = 1/f 2 , where f is the Fourier frequency f = 2πk/N [111]. In the GOE case, the asymptotic form is P GOE (k) = N 2πk = 1/f . This 1/f noise has been argued to be an unique characterization for quantum chaotic systems [111][112][113].
The variable k in the power spectrum does not have units of inverse time, it is rather a dimensionless "energy momentum" alluding to the fact that it comes from the argument of a Fourier transform of an energy coordinate. In the same spirit of the spectral form factor, the dimensionless Thouless energy momentum k Th is interpreted as the smallest momentum for which P (k) ∼ P GOE (k) ∝ 1/k withk = k/N (see Fig. 2). l = 1/k can be interpreted as an average energy distance in the spectrum, henceforth at distances l > l Th , with l Th = 1/k Th , the levels are uncorrelated whilst at l < l Th they are correlated and well described by random-matrix theory. Figure 2 shows the power spectrum as a function of k for L = 17 and several disorder strengths. The dashed line in Fig. 2 is the GOE behavior, P (k) ∝ 1/f , from where one can read the Thouless momentum k Th . In analogy with the spectral form factor, the Thouless time is given in units of the Heisenberg time: where we use the same expression for t H as in Eq. (11). In Appendix B the procedure for computing k Th is explained in detail.

E. (Sub-)diffusion and Thouless time
Having defined the spectral form factor and the power spectrum, our main probes we will use to detect longrange correlations in the energy spectrum, we now discuss the relations between (sub-) diffusion and the scaling of the Thouless time with L.
Let us start by pointing out the differences between the propagation in hypercubic, Z d -like, lattices and hierarchical structures, e. g., RRG or the Bethe lattice. Diffusive dynamics is usually defined through the spread of an initially localized wave packet. In the Z d lattice, diffusion is quantified by the Gaussian profile of a propagating wave-packet. As a result, in the Z d lattice the survival probability decays as Π(t) ∼ t −d/2 and the mean squared displacement ∆X 2 (t) ∼ t 2β with β = 1/2. Sub-diffusive propagation is consequently given by β < 1/2. In particular, the related Thouless time scales with the linear system size as t Th ∼ L 1/β , which is the time to cross the system. On the other hand, in diffusive processes on tree-like structures the survival probability decays exponentially fast ∼ e −Ω(K)t , with a certain decay rate Ω(K), depending on the branching number K, and the mean square displacement grows as ∼ t 2β with β = 1 [109,110]. As a result, linear growth of the Thouless time with the diameter L (t T h ∼ L) implies diffusion and sub-diffusion takes place if ∆X 2 (t) ∼ t 2β (t T h ∼ L 1/β ) with β < 1.
In a localized phase, where the degrees of freedom are frozen at long times, the Thouless time is comparable with the Heisenberg time. Indeed, the ratio t Th /t H can be used as a probe to detect an delocalization-localization transition. In the extended phase t Th /t H → 0, while in a localized one t Th /t H ∼ O(1) in thermodynamic limit.

IV. RESULTS
In the following, we present results obtained by solving the models in Eqs. (1) and (2) using exact diagonalization techniques and computing the spectral form factor and the power spectrum to extract the Thouless time as a function of system size. Unlike the static Hamiltonian case, the critical disorder, W AT , for the Anderson transition in the Floquet model, Eq. (2), is unknown. Thus, before starting our investigation on the long-range correlations, we start with the short-range ones. This preamble will allow us to obtain a rough estimate of the critical value for the Floquet model and, as well as to determine the range of disorder strengths in which finite-size corrections are negligible. Figure 3 shows the r-statistics as a function of disorder strength W and several system sizes L for both models. The critical strength for the static Hamiltonian model in Eq. (1) has been estimated analytically at W AT ≈ 18.1 [76,79,81,96]. As expected at weak disorder the r-statistics is GOE and at strong disorder the spectrum is Poisson. For the available system sizes, the r-statistics in Eq. (4) shows crossings close to W ≈ 16 < W AT with an apparent drift towards larger W for growing graph diameters L, as one can observe in Fig. 3. Therefore, at W ≳ 16, when our system sizes are smaller than the correlation volume, the physics will be largely dominated by finite-size effects, and the system looks already localized. Since we are interested in the transport properties of the system, we will hence mainly focus on disorder strengths W < 16. For our Floquet model of Eq. (2) there is no exact estimate of the critical disorder. The finite-size gap ratio crossover between the ergodic value and the Poisson one happens around W/T ≈ 18. This value, W/T ≈ 18, provides us with a lower bound of the exact critical value and we consider disorders W/T > 18 to be already within the localized phase. A more precise estimate of the critical disorder of the Floquet model is not needed here and beyond the scope of this work. Since we consider transport properties, we focus on the delocalized part of the phase diagram at a safe distance from the critical regime. We therefore restrict the discussion to disorder strengths W < 16 for the Hamiltonian model and W/T < 18 for the Floquet model.
Having ascertained a reliable delocalized regime for our simulations, we first look at the raw data for the spectral form factor K(τ ) and the power spectrum P (k) for the Hamiltonian and Floquet model. The full data, for several system sizes and relevant disorder strengths, can be found in Fig. 4 and 5, respectively. As one can observe in Fig. 4, at weak disorder W = 4 and W/T = 6.0, deep in the ergodic phase, the GOE "ramp" at late time τ , K(τ ) ≈ 2τ , gets longer with the increasing system size. When approaching the Anderson transition, the GOE ramp gets shorter until it is barely visible at disorder strengths W = 14.0 and W/T = 16.0. These disorder amplitudes are not yet within the localized phase accord-ing to the mean level spacing in Fig. 3, although in a finite system they show very slow dynamics signaled by very long Thouless times [90]. In Fig. 5 we show the power spectrum for the same disorder values. Similarly to K(τ ), at low disorder W = 4.0 and W/T = 4.0 the power spectrum falls on top of the GOE prediction (P (k) ∼ 1/k), which is shown by the dashed line. Upon increasing disorder the power spectrum meets the GOE curve at larger k. At W = 12.0 and W/T = 14.0 the GOE scaling region gets shorter but still increases with system size, denoting a flow towards delocalization. On the other hand, at W = 14.0 and W/T = 16.0 all accessible system sizes seem to follow the Poisson behavior (P Poisson (k) ∼ 1/k 2 ), although at W = 14.0 the power spectrum has a small visible offset from the Poisson value. This offset is also seen in the Floquet model at W/T = 10.0 for small system sizes, hence that might be a sign of flow towards delocalization at larger system sizes, compatible with the previously estimated location of the Anderson transition in the static model at W AT = 18. At very strong disorder, W/T = 22, where the Floquet model is likely in the localized phase (cf. Fig. 3), we find excellent agreement of the spectral form factor with Poisson statistics.
Based on this data, we now turn to the main aim of this work, the scaling of the Thouless time with system size t Th ∼ L 1/β(W ) . As discussed in Sec. III E, the exponent β is connected to the transport properties. The Thouless time can be understood as the time for a particle to propagate throughout a graph of diameter L. In particular, β = 1 corresponds to diffusion and 0 < β < 1 to sub-diffusion. In Fig. 6 we show the Thouless time, t Th , extracted from K(τ ) and P (k) as a function of L on a doubly logarithmic scale for several disorder strengths W . For comparison in Fig. 6 the Heisenberg time (dashed line) and the diffusive scaling of the Thouless time with β = 1 (solid grey line) are also reported. The Heisenberg time signals the localization behavior (t Th ∼ t H ). At weak disorder, the scaling of the Thouless time is indeed compatible with a power law L 1/β , while at strong disorder we observe significant curvature, stemming from the exponential scaling of the Heisenberg time and t Th ∼ t H in the localized phase.
From the apparent power law scaling of the Thouless time t Th ∼ L 1/β , we extract the dynamical exponent β, which is reported in Fig. 7. We find that even at weak disorder, in both models β < 1, corresponding to subdiffusion propagation [90]. Indeed, for disorders W = 6.0 and W/T = 6.0 the dynamical exponent β ≈ 0.5 and, then decreases upon increasing disorder as expected. At stronger disorder, W = 15.0 and W/T = 15.0, β stops decreasing and becomes difficult to determine due to the crossover to the exponential scaling, t T h ∼ t H ∼ 2 L , characteristic of a localized phase. Remarkably, in the Floquet case, both the spectral form factor and the power spectrum yield very similar estimates of the Thouless time for any disorder value and system size, confirming the expectation that the Floquet model exhibits less finite-size effects compared with the Hamiltonian case. In the static case, only 60 % of the centered at the middle are taken. The "energy momentum" k is normalized by the effective size of the spectrumÑ = 0.62 L , thus, the arguments of PGOE(k) and PP oisson(k) are also rescaled byÑ . In Floquet case, k is rescaled by the matrix size N = 2 L such that the Nyquist frequency becomes 0.5.  On the other hand, we observe a discrepancy at weak disorder W ≲ 8 between the extracted Thouless time from the spectral form factor and the power spectrum in the Hamiltonian model in Fig. 7(a)). The result from the spectral form factor seems to have a much smaller derivative with respect to W between W = 4.0 and W = 6.0 compared to the dynamical exponent extracted from the power spectrum in the same range. One possible cause of this discrepancy could be due to the spectral cutoff used for mitigating unwanted edge effects (see Appendix B for more details). Discarding 20% of the states at both edges implies a reduction of the energy bandwidth. This effective bandwidth ∆E sets the smallest time scale we have access to (∼ 1/∆E). In the static Hamiltonian model the bandwidth is proportional to W and independent of L, therefore at small W s and Ls the Thouless time might be out of reach because it is smaller than the shortest accessible time scale. When the Thouless time is getting close to the inverse bandwidth, both spectral measures will be strongly affected by edge effects even in the presence of the filter Eq. (10). This effect can be observed in Fig. 6 at L = 8 where t T h is almost equal for W = 4.0 and W = 6.0. As a result, the spectral form factor and the power spectrum might behave differently due to edge effects yielding different scaling when the matrix size is not big enough for bringing the Thouless time within the bulk of the spectrum. The latter strongly limits the study of the Thouless time scaling at low disorder in the current system sizes. This effect is absent at the intermediate disorder because the Thouless time is larger than the inverse energy bandwidth.

V. CONCLUSIONS
In this work, we have studied the long-range spectral correlations in the Anderson model on the random regular graph in terms of two probes: the spectral form factor and the power spectrum. We provide numerical evidence that the sub-diffusive phase in the Anderson model on the random regular graph, claimed in [88][89][90][91], can be probed by the global spectral statistics which does not include any information about the eigenstate structure. To reduce finite-size effects due to mobility edges and energy dependence of the density of states, we remove the constraint of energy conservation by introducing a Floquet version of the Anderson model on the RRG. The Floquet model has the advantage of having a structureless density of states and the absence of mobility edges.
In this setting, we extract the Thouless energy from both the spectral form factor and the power spectrum which agree with each other in the reliable range of Hamiltonian parameters and show good algebraic scaling with the diameter of RRG. The above scaling of the extracted Thouless energy with the graph size allows us to unambiguously observe the sub-diffusive character of this dependence in the entire delocalized phase and find the exponent β < 1 of the sub-diffusion. This finding is in a good agreement with a recent random-matrix consideration [94] of the RRG where the dynamical characteristics of the model, like the return probability, has been claimed to show the sub-diffusion up to zero disorder W = 0. As the Anderson model on the RRG provides a proxy for the dynamics of more realistic manybody disordered systems, our consideration provides the ground to the tentative sub-diffusive phase in the finite size many-body localized regime [45] close to the true MBL transition and to its observation with the global spectral probes like the spectral form factor [65,68,70].

VI. ACKNOWLEDGMENTS
states of both folded and unfolded spectra are shown. By design, the density of states of the unfolded spectrum is constant. The folded spectra are rescaled to the interval [0,1] for perspective purposes. In this section we explain the procedure of the extracting of the Thouless time from the equidistant spectra with unit mean level spacing ε ∈ [0,N ]. In the Floquet case (when spectral edge effects are not present), the spectral form factor and the power spectrum are computed without filtering and then compared to the GOE functions K GOE (τ ) = 2τ [108] and P GOE (k) = N/2π 2 k [111] mentioned in sections III C and III D, respectively. Henceforth we refer only to the spectral form factor, however the procedure is the same for both quantities with the caveat that k/N is used instead of k. At early τ we have K > K GOE , consequently we identify the Thouless time as the smallest τ for which K = K GOE , we instead use log(K/K GOE ) as a measure of the distance between the two functions. In practice, it is useful to set a threshold like log(K(τ Th )/K GOE ) = 0.4 and take the time that satisfies this relation as τ Th . This is illustrated in Fig. 9 and 10. It can be seen that a finite but not too small threshold captures the system size scaling of the curves which encodes the Thouless time system-size scaling as well. We have verified that our results are not affected by the choice of the threshold beyond the noise fluctuations, see Figs. 9 and 10.
In the Hamiltonian case, states at the edges of the spectrum are usually more localized and so have strong unwanted influence on the spectral measures. In order to remove such spurious effects, the spectrum is first cut off from the edges. Throughout this work 20% of the states at each edge have been discarded leaving 60% of the total spectrum centered at the middle of it. At the end, one works with a spectrum of reduced dimensioñ N = ⌊0.6N ⌋. This cut is made after the unfolding, meaning that the unfolding is carried out on the whole spec- trum before cutting off the edges. Edge effects are further mitigated by multiplying a Gaussian function centered at the middle of the spectrum and reduced variance compare to the effective spectral width. In other words one does exp(−i2πετ ) → g(ε) exp(−i2πετ ) with g(ε) given in Eq. (10). The same applies to the spectral statistic δ n = ε n − n where its Fourier transform carries a weight g(ε n ) given by Eq. (14). The parameter η controls the variance of the Gaussian filter, in Fig. 9 the function log(K/K GOE ) is shown with η = 0.3 and η = 0.4. The small difference in filter width does not have major effects. We have set η = 0.3 throughout the main text. Next, the spectral form factor, the power spectrum, and Thouless time are computed as explained above.   Fig. 11. The same fitting procedure explained in caption of Fig. 7 of the main text was used.
Appendix C: Sub-diffusion at different driving period In order to check the sensitivity of sub-diffusion in the Floquet-driven RRG to changes in the half driving period T we have computed the spectral form factor for T = 1.5, 1.8, 2.2, 2.5. Recall that in the main text T = 2.0 is fixed. First of all, the gap ratio shows the crossing of the curves for different system sizes at roughly the same disorder value for all T (not shown), therefore we still see a localization transition within the same range of disorder. In Fig. 11 the Thouless time extracted from the spectral form factor is plotted. We can see that at low disorder the scaling the Thouless time is polynomial in L regardless of the half period T . Indeed Fig. 11 is quite similar to panels (c) and (d) of Fig. 6.
In order to confirm a sub-diffusive character of the transport for different T , we estimate the dynamical exponent fitting t T h ∼ L 1/β for each curve in Fig. 11. The resulting dynamical exponent β is shown in Fig. 12 as a function of disorder strength W . We can see that around W ≳ 25 the dynamical exponents saturate pointing out the localization threshold for all periods. At weak disorder strengths β < 1 independently of T . The latter suggests that the sub-diffusion regime does not arise due to fine-tuned choice of driving period in the model of Eq. 2, but it is quite robust and generic.