{'fraction_non_alphanumeric': 0.05067888472928424, 'fraction_numerical': 0.02039447952170755, 'mean_word_length': 4.303111111111111, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 16, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Universal dynamics of density correlations at the transition to many-body localized state

Within one-dimensional disordered models of interacting fermions we perform a numerical study of several dynamical density correlations, which can serve as hallmarks of the transition to the many-body localized state. Results confirm that density-wave correlations exhibit quite abrupt change with the increasing disorder, with nonvanishing long-time value characteristic for nonergodic phase. In addition, our results reveal in a wide time-window a logarithmic variation of correlations in time, which we can bring in connection with the anomalous behavior of the dynamical conductivity near the transition. Our result support the view that transition to many-body localization can be characterized by universal dynamical exponents.


I. INTRODUCTION
The idea of many-body localization (MBL) emerged from the well understood Anderson localization of non-interacting (NI) fermions, [1][2][3][4] by taking into account the many-body interaction. 5,6 The basic claim that localization can, at large disorder, persist in the whole spectrum and consequently at all temperatures, has been by now supported by numerous studies on one-dimensional (1D) disordered models. There are essential features of the MBL phase 7-17 which have been confirmed numerically, predominantly for interacting spinless fermions: the change of the level statistics from a Wigner-Dyson to a Poisson-like in the nonergodic phase, 18 the vanishing of d.c. transport at any T 19-25 even beyond the regime of the linear response 26 , the logarithmic growth of the entanglement entropy in the MBL phase, [27][28][29][30][31][32][33] the nonergodic behavior of correlation functions related also with the existence of local conserved quantities. 30,31,[34][35][36][37][38][39][40][41][42][43][44][45][46][47][48][49] Experimental support for the MBL comes so far from studies of cold-atom systems. 37,[50][51][52] Absence of d.c. transport under constant force 50 and the nonergodic evolution of initially quenched state 37,51 have been used as the experimental criteria. In particular, in the latter studies the long-time remainder of the charge imbalance has been used as a practical hallmark of the MBL which becomes finite within the nonergodic phase. Such a quantity is a natural counterpart of the d.c. transport quantities, as e.g., the d.c. conductivity σ 0 being nonzero in the "normal" ergodic phase (including also the possibility of subdiffusive transport). [19][20][21][22][23] However, the detailed theoretical and numerical analysis of such indicators of the MBL is still missing. 31,37 On the other hand, several numerical calculations of dynamical conductivity 21,24,25 as well as more general studies of dynamics using the renormalizationgroup approach 22,31,53,54 indicate that the transition is primarily characterized by dynamical critical exponent, e.g. of the dynamical conductivity σ(ω).
In this paper, we study dynamical density correlations within the prototype 1D model of the MBL system. We concentrate on the aspect how the Anderson localization, established for NI fermions in 1D for any disorder strength W , is destroyed by a modest repulsive interaction V . In particular, we study the time-dependent density-wave (DW) correlation functions C(t), closely related to the charge imbalance. 31,37,51 We show that C(t = ∞) reveal quite sharp transition at large disorder W ∼ W c , hence they can serve as hallmarks of the MBL phase. In addition, our results show so far novel and universal logarithmic time-dependence of C(t) in a very wide time-range being particularly extended near the MBL transition. Using the memory-function analysis we show that this anomalous density dynamics is closely related to the scaling of the dynamical conductivity, confirming that such universal dynamical scaling appears to be more fundamental hallmark of the MBL than the stationary value C(t = ∞) itself. Finally, we find evidence that such anomalous dynamics is characteristic also for other disordered systems with no apparent connection with NI Anderson localization.

II. MODEL AND NUMERICAL METHODS.
We consider the prototype model for the MBL, i.e. the 1D model of interacting spinless fermions with random local potentials, where n i = n i − 1/2 = c † i c i − 1/2. We take quenched disorder h i with a uniform distribution −W < h i < W . The model (1) is in 1D equivalent to the anisotropic Heisenberg model with random fields. Taking in the following t 0 = 1 as the the energy unit, we mostly consider cases of modest interaction V /t 0 = 1 being closer to the NI case V = 0. We also fix the density of fermions to half-filling, i.e. having N = L/2 fermions on a system with L sites and with periodic boundary conditions. Since the MBL at larger W occurs at any temperature, we study the limiting case T → ∞ being the optimal one for numerical calculations. However, in the Appendix A we show also results for finite temperatures.
The dynamical quantities are obtained mainly from exact diagonalization (ED) of the Hamiltonian, (1), where one can reach L = 16. Further on we also perform calculations of the DW correlations using the microcanonical Lanczos method (MCLM), 55,56 very suitable for the study of dynamical quantities at T ≫ 1 and allowing for substantially larger system, e.g., L = 24. In the latter method, the main restriction is finite frequency resolution of dynamical spectra, with typically δω ∼ 10 −3 (for ∼ 10 4 Lanczos steps) restricting the reachable times t < 1000.

III. DENSITY-WAVE CORRELATIONS.
As a very practical tool to investigate nonergodicty in disordered system, we employ the staggered DW operator O s , and its normalized autocorrelation function C s (t) can be directly related to the time-dependent imbalance measured in the cold-atom systems. 37,51 Since C s (t → ∞) > 0 marks a nonergodic behavior, the stiffness C s (∞) can be used as an indicator for the MBL transition. Here, ... represents average over ∼ 10 2 ÷ 10 3 configurations of h i . The decay of C s (t) is at short t ∼ 1 masked by oscillations with frequency ω = 2 which emerge from H 0 only. 26,31,37 These oscillations are clearly visible in Fig. 1. As discussed later on, they can also be observed from the Fourier transform C s (ω) = 1/2π ∞ −∞ dt exp(iωt) C s (t). It is convenient to consider also a modified DW operator as well as the corresponding correlation function C l (t). Index l enumerates sorted energies of the single-particle Hamiltonian, H − H V = l ǫ l ϕ † l ϕ l . Consequently, labelling the position of the localized states by l is to some extent arbitrary. The correlation function C s (t) describes the decay of the initial staggered density wave with the wave-vector q = π. Contrary to this, C l (t) describes how/whether the system retains the information about a random density distribution. Here, we search for common properties of both correlation functions which should be generic for most of the spatial particles distributions. Note also that the decay of C l (t) is solely due to many-body interactions, while C l (t) = 1 for V = 0.
In order to further demonstrate universality of the longtime behavior we study also a system with a homogeneous single-particle Hamiltonian (h i = 0) but with a disordered many-body interaction 57 where V i ≥ 0 (positive to avoid localization due to bound states) are uniformly distributed variables, 0 ≤ V i ≤ 2W . We study the autocorrelation function C V (t) for the related DW operator The DW decay function C s (t) is calculated numerically using both methods described above. C l (t) and C V (t) are obtained from ED. In the Appendix B we show for moderate time-window that C l (t) can also be obtained in a reduced basis for larger systems than accessible for ED calculations. C s (t) and C l (t) are obtained for the Hamiltonian (1) with V = 1 , while C V (t) for interaction (6) with V ′ = 1. Whenever our discussion applies to all correlation functions, we omit subscripts s, l, V denoting C = C (s,l,V ) and O = O (s,l,V ) . To recognize the essential features corresponding to t ≫ 1 (or ω ≪ 1) it is convenient to present the integrated spectra where H|α = E α |α . When carrying out the ED calculations we first obtain I(τ ) from Eq. (9) and then C(ω) by differentiating Eq. (8). Neglecting accidental degeneracies, C(t → ∞) can formally be obtained from Eqs. (8) and (9) by restricting summation to diagonal terms α = α ′ . In Figure 1 we compare the real-time correlation function C s (t) with the integrated spectrum I s (τ ). Note that I(τ ) represent to some extent the real-t evolution, since I(t) = C(t) in both limits t → 0 and t → ∞. Since I(τ ) doesn't show the transient oscillations, it is more convenient to use the latter quantity for the studies of the long-time dynamics which turns out to be very slow. Is(τ )   (1). It is interesting that roughly the same W c is obtained for a system with the disordered interaction (6). For W > W c , the extrapolation L → ∞ leads to finite stiffnesses I * , which can be used as convenient indicator for the MBL phase. 37 (c) Most remarkable, all I(τ ) reveal for W ≥ 3 a very wide time-window beyond τ > 1 with a slow, logarithmic-like decay, e.g. I(τ ) ∼ a − b log(τ ). As shown in Figs. 2c and 2d, this behavior is particularly clear at the transition W ≃ W c . In the latter case, deviations from the logarithmic time-dependence diminish when the system size increases, hence these deviations seem to represent the finitesize effects. Such a decay extends typically to τ * ∼ 1000 for largest systems L = 16 available for ED, before saturating at C(∞). The decay continues apparently to the largest τ * ∼ 1000 available by MCLM (L = 24). The comparison of results for different L confirms that τ * can extend at least for one decade when increasing the system from L = 12 to L = 16.
While the observations (a) and (b) have been at least partly reported before, the universality of the slow (logarithmic) variation 58,59 appears to be a novel one. It is evident that I(τ ) ∝ log(τ ) has to emerge from an anomalous ωdependence. Taking the ED results at finite L as the input, it requires Here, A = C(t → ∞) is the stiffness, while C reg (ω) represents the regular part. We note that at finite L the latter is meaningful only for ω > δ 0 where δ 0 is typical level spacing (with δ 0 ∼ 10 −2 , 10 −3 for L = 12, 16, respectively). We display in Fig. 4 the ED results for C reg (ω) on a log scale, which clearly reveal that ζ ≤ 1 in a wide range of W in the vicinity of W ∼ W c . On the other hand, the saturation with ∆ ∼ 10 −2 > δ 0 is among results well resolved only in the case of weak disorder, i.e. W = 2. The question remains whether ∆ is finite also, e.g. for W=3. Still, our results in Fig. 4c strongly suggest that ∆ vanishes within the MBL phase, i.e. for W > W c , while more detailed behavior of ∆ in the regime W W c remains an open problem. A strict logarithmic dependence of C(t) implies that ζ → 1, whereas Fig. 4 shows that ζ is close to but smaller than 1 leading to time dependence as t −ε , with ε = 1 − ζ > 0. We cannot judge whether this tiny deviations are real or show up as numerical artifacts. However, since t −ε = 1 − ε log(t) + O[ε 2 log(t) 2 ] such deviations may become relevant first for long times t ∼ exp(1/ε).

IV. THE RELATION TO DYNAMICAL CONDUCTIVITY
In order to relate such a behavior to other dynamical observables, it is convenient to analyze the full complex response function so that C ′′ (ω) = C(ω). Since C(ω) is analytical function of ω (for ℑω > 0), it can be represented in terms of the complex In particular, for the DW correlations C s (ω), M (ω) is related to (an effective) dynamical conductivity at the same q = π, i.e., M ′′ (ω) ∼ sin 2 (q/2) σ ′ (q, ω). Note, however, that σ(q, ω) is the current response function only in the limit q → 0. [60][61][62] In general, the representation of C(ω), Eq. (13), in terms of M (ω) has a clear advantage that instead of diverging C(ω → 0), we are dealing with a regular Γ(ω) = M ′′ (ω) > 0, representing the DW relaxation-rate function.
Results obtained for all correlation functions are compared in Fig. 5.
For W > W c the nonergodic contribution A > 0 is unavoidable. The latter leads to Γ(ω) ∼ ω α with α = 2 − ζ. It is evident that the marginal case corresponds to ζ = 1 and α = 1, which is another criterion for the MBL transition. Our results (not shown) indicate that for W ≫ W c the critical exponent is increasing, α > 1, nevertheless we find in a broad range of W nearly constant α ∼ 1.

V. CONCLUSIONS
We presented numerical results for several dynamical density correlations. We have shown that the correlation functions display universal long-time behavior. It holds true for the prototype MBL Hamiltonian as well as for a system with homogenous single-particle Hamiltonian but with disordered many-body interaction. 57 All quantities reveal a nonergodic behavior, well visible in the stiffnesses, C(t → ∞) > 0, which remain finite even after the extrapolation of finite-size results to L → ∞. In this sense, the extrapolated values C * can be used as indicators of the nonergodic MBL phase, in direct correspondence to the imbalance stiffness in the coldatom experiments. 31,37 Still, the main message of our study is that all the correlation functions C(t) exhibit near the MBL transition anomalously slow relaxation towards the presumable t → ∞ limit. Such a logarithmic time-dependence is visible over several decades in the window 1 < t < t * where t * > 1000 seems to be limited at W ≥ W c only by finite size restrictions of our numerical methods. Although we are dealing in our study with DW at wave-vectors q ≫ 0, our analysis reveals a close similarity of the relaxation functions Γ(ω) to the behavior of the optical conductivity σ(ω) which is the q → 0 property. The observed low-ω behavior Γ ∼ Γ 0 + gω α indicates that the MBL transition is best characterized by the critical exponent α = 1, consistent with several other numerical and renormalization-group analysis of dynamical quantities. 21,22,24,25,54 The extremely slow relaxation together with finite time-windows, both in experiments as well in numerical studies, suggest that identifying MBL from the stiffness (i.e. from the saturated correlation functions) might be very challenging. Hence a more proper definition of the MBL should be related just to the critical dynamics, since this dynamics can be established within much shorter timewindows.
There are nevertheless clear open questions, e.g., whether Γ 0 > 0 for arbitrary W < W c (as well a d.c. σ 0 > 0) 20,24,25 or there exists intermediate regime with subdifussive behavior 21,22,63 consistent with Γ 0 = 0. Beyond the fundamental importance, even more relevant is the question, whether experimental results on the imbalance relaxation could also confirm logarithmic character of relaxation towards the limiting nonergodic stiffnesses. 37  All the correlation functions discussed in the main text have beed obtained at infinite temperature. While our numerical methods are inappropriate for the studies of the lowtemperature regime, these methods provide reliable results also for large but finite temperatures. Here, we demonstrate that our main conclusion concerning the dynamics of correlation functions remain valid also in the latter regime, i.e., for the state ρ = exp(−βH)/Z, Z = Tr[exp(−βH)] with the inverse temperature β ≪ 1. Deep in the MBL regime one may formally assume the system to be in a thermal state. However such assumption may be unphysical simply because systems do not thermalize in the latter regime. In order to avoid this possible inconsistency we restrict our studies to the case W = 4, i.e., we stay close to the MBL transition.
At finite temperature, the real-time correlation function becomes while the integrated spectra can be calculated from In Fig. 6 we show integrated correlation function, I s (τ ), obtained from the exact diagonalization for L = 16. The increase of β from 0 to 0.5 only weakly affects I s (τ ), as shown in figure 6a. As expected, the stifness C(∞) = I(τ → ∞) slightly increases with β. Unfortunately, the limited accuracy of the finite-size scaling of the disorder-averaged data does not allow us to judge whether the critical disorder for the MBL transition depends on β. However, results in Fig. 6 clearly show that our main claim concerning the extremely slow quasi-logarithmic decay of correlation functions remains valid also for large but finite temperatures. then written in terms of localized states as Considering the many-particle states within such localized basis  (1). The goal is to use RBA with N r ≪ 2 L basis states and evaluate within it the dynamical quantity via a direct time evolution. In this way one may bypass and even monitor the question of MB resonances, the well known problem within the theory of localization 1,5,54,64,65 . We typically take N r ≃ 4 · 10 4 . In Fig. 7 we compare results for the integrated correlation functions obtained from various methods, see equations (8) and (9) in the main text. For τ ≤ 10 2 results obtained from exact diagonalization (L ≤ 16) nicely overlap with the data from RBA for L=20. While, the ED results saturate for larger τ , the logarithmic decay continues apparently even further for L = 20 up to the largest τ ∼ 2000 available for the RBA. The comparison of results for different L confirms that the range of logarithmic decay can extend at least for one decade when increasing the system from L = 10 to L = 20. It confirms also our main result that a more proper/practical definition of the MBL should be related to the critical dynamics rather than to stiffness. The latter quantity becomes available first after correlation functions saturate, while former one can be measured already during the relaxation.