Thermalization Universality-Class Transition Induced by Anderson Localization

We study the disorder-induced crossover between the two recently discovered thermalization slowing-down universality classes -- characterized by long- and short-range coupling -- in classical unitary circuits maps close to integrability. We compute Lyapunov spectra, which display qualitatively distinct features depending on whether the proximity to the integrable limit is short or long ranged. For sufficiently small nonlinearity, translationally invariant systems fall into the long-range class. Adding disorder to such a system triggers a transition to the short-range class -- implying a breaking of this invariance -- and in the very limit of vanishing non-linearity Anderson localization emerges. The crossover from long- to short-range class is attained by tuning the localization length, $\xi$, from $\xi \approx N$ to $\xi \ll N$, where $N$ is the system size. As a consequence, the Lyapunov spectrum becomes exponentially suppressed, depending on how strongly its translational invariance is destroyed. We expect that this disorder-induced crossover will lead to prethermalized phases and, following quantization, to many-body localization.

Thermalization is a universal property of nonintegrable many-body systems, with characteristic timescales that diverge upon approaching integrability [1][2][3][4][5][6][7].In the nearintegrable regime, one can interpret the system as a perturbation of an integrable one, with the perturbation's overall effect being to couple the action-angle variables of the unperturbed system in a nonlinear manner [4][5][6][7].Thermalization slowing-down was shown to strongly depend on the perturbation's coupling range [5,7,8], which can be classified according to two universality classes: short-range network (SRN) and long-range network (LRN) [5].[9][10][11][12].To quantify these slowing-down processes, one can study finite time averages of observables.These, however need to be selected with care [5,11,13] to ensure that one obtains proper ergodization timescales [9][10][11][12].The ambiguity in the choice of suitable observables led to the use of a novel method based on the Lyapunov spectrum (LS) to distinguish the two universality classes [11].In this approach, which can also be employed to diagnose phase transitions [14][15][16][17], significantly more information than the typical Lyapunov time (given by the inverse of the largest Lyapunov exponent) is available, since each Lyapunov exponent (LE) in the spectrum carries its own characteristic timescale.
The numerical challenges of dealing with weakly nonintegrable many-body Hamiltonian systems led to the study of one-dimensional unitary circuits maps, for which time-evolution is exact.This removes time discretization errors and allows for substantially larger evolution times, and thus for higher resolution in the LS [11,13].The resulting universal scaling properties of the LS rendered possible an unambiguous identification of the different SRN and LRN universality classes.These were also observed in a recent study of multidimensional Josephsonjunction networks across all possible lattice dimensions [18].The predictive power of unitary circuits maps was thereby confirmed, and above all the LS has established itself as an invaluable tool in the study of the thermalization of many-body systems [13].
Our objective is to gain insight into the interplay of disorder with nonintegrability, which might reveal a connection with the celebrated phenomenon of many-body localization [19].For this, we employ unitary circuits maps to investigate how the disorder impacts the two thermalization slowing-down universality classes.We use tailored disorder which leads to Anderson-localized states [20] in the integrable limit of linear maps.The tunable localization length, ξ is universal for all eigenmodes and is determined solely by the hopping-like parameter associated with the unitary circuits map [21].We then demonstrate that the system's thermalization universality class changes from the LRN to SRN as the localization length is tuned from ξ ≈ N to ξ ≪ N , where N is the system size.Our findings intertwine the fields of many-body localization and thermalization of weakly non-integrable systems and open a new venue for connecting the slowing down of classical many-body dynamics in the presence of disorder and localization physics in quantum many-body systems like the ones proposed in Refs.[22][23][24][25][26][27].
We employ a modification of the classical unitary circuits maps introduced in [11,13].These consist of a one-dimensional periodic chain of N complex numbers, ψ n , with n = 1, 2, . . ., N denoting the sites (we assume that N is even).An initial vector ⃗ Ψ with these complex scalar components ψ n is evolved by applying iteratively the unitary map Here, we interpret the number of iterations as a discretetime, i.e., ⃗ Ψ(t + 1) = U ⃗ Ψ(t).A pictorial description of U is provided in Fig. 1.The operators C (even/odd) are linear transformations.In matrix representation they are block diagonal with each block consisting of a 2×2 unitary matrix C n,n+1 coupling the components ψ n (t) and ψ n+1 (t), parameterized by a hopping-like angle parameter θ, All remaining matrix elements are zero and, due to periodicity, ψ N +1 = ψ 1 .Successive applications of C (odd)  and C (even) intertwine neighbors to the left and to the right of site n (green boxes in Fig. 1).In matrix representation the operator G corresponds to a diagonal matrix whose elements G n,n are nonlinear on-site potentials acting as Similar to an anharmonic oscillator whose frequency depends on its level of excitation, the imposed phase shift (during one given time step) depends on the amplitude |ψ n |.This term is responsible for the chaoticity of the dynamics (yellow boxes).Finally, D acts on a given site, that is, in matrix representation again is a diagonal matrix with elements where the ϵ n are site-dependent disorder potentials that are uniformly distributed in [−π, π] (blue boxes).The unitarity of U implies that the total squared-norm | ⃗ Ψ(t)| 2 is a conserved quantity.Accordingly, in order to allow all possible typical scenarios for the temporal behavior with equal probability, we chose the initial values of the squared moduli of the rescaled components η n = N |ψ n | 2 uniformly spread over the N sphere with the jointprobability distribution P ({η n }) ∝ δ N − N n=1 η n .This yields for the probability distribution of the η n P (η) For the computation of the LS, Λ, comprising the LEs, Λ i , with we follow the calculation procedure of Ref. [13], outlined in more detail in the appnedix.Note that we are considering the part of the LS which is composed of non-negative LEs in the spectrum.We can do so because, similarly to time-independent Hamiltonian systems, the N positive LEs come in pairs with negative ones of the same modulus, Λ i = −Λ 2N −i+1 [13].The final simulation time is denoted by T s .
In the absence of hopping, achieved by setting θ = 0, the sites decouple and the squared-norms |ψ n | 2 are conserved.This results in N functionally independent conserved quantities.Since the system has N degrees of freedom implying a 2N dimensional phase space, this limit is integrable.The squared-norms then provide a discrete analog to the actions in continuous Hamiltonian systems when θ = 0 [11,13].A second integrable limit is realized by setting g = 0.In this case, the sites remain coupled but there is no source of nonlinearity, and the squared-norms of the corresponding normal modes, i.e. of the eigenvectors of Û (g = 0), are conserved [11,13].
Let us consider first the disorder-free case.When approaching the decoupled integrable limit, i.e., when θ ≪ 1, each action is only weakly coupled to its nearest neighbors, resulting in a SRN [11,13].In the regime g ≪ 1, on the other hand, all normal mode actions are weakly coupled to each other, such that all-to-all linear interactions persist and implicate a LRN [11,13].Furthermore, it was demonstrated that, in the absence of disorder, the Lyapunov spectra for LRN and SRN systems behave significantly differently as the integrable limit is approached: While most of the LEs in the LRN class are of the same order of magnitude at any distance from the integrable limit, the LEs in the SRN class are damped and end up spanning several orders of magnitude, with the span increasing upon approaching the integrable limit (see Fig. 2).We derived an ansatz for the rescaled Lyapunov spectrum Λ of the LRN and SRN.This rescaled spectrum is obtained by dividing the original LEs by the maximal LE (mLE), Λi = Λ i /Λ 1 .Starting point of the derivation is the analytical result for the number of Lyapunov exponents N (Λ) below Λ = Λ i , given by N (Λ i ) = N + 1 − i, i = 1, 2, . . ., N .Away from the limiting values, 0 < Λi < 1, it is well described by the integrated Wigner semicircle law [29][30][31] yielding an inverse semicircle law for Λ versus ρ i = i/N .The damping rate of the spectrum in the SRN regime, on the other hand, is exponential.An approximation of the resulting expres- sion, which for ρ > 0 reflects these features, is provided by the ansatz The parameters α, β, γ are determined from the fit of this ansatz to the rescaled LS.Yet, the qualitative behavior of Λ(ρ) is already well captured with γ = 1 [18].In addition, we compute the rescaled Kolmogorov-Sinai (KS) entropy κ = 1 Λ dρ, which is a very useful quantifier of the different regimes.In the SRN regime, κ tends to 0 when approaching the integrable limit, while in the LRN regime, it saturates at some value, as shown in Fig. 3.The insets of Fig. 3 show that the fitting parameter α does not vary significantly for all considered cases.In contrast, the exponent β grows with a power law in the SRN upon lowering θ, while it saturates similarly to α in the LRN with decreasing g.As a consequence, the rescaled LS Λ(ρ) saturates on some analytic invariant curve in the LRN for g → 0, and knowing one LE (e.g. the mLE) results in knowing all others as well, whereas the SRN is characterized by an exponential damping of LEs compared to the mLE.5) evaluated with the coefficients from the fits.
Our next objective is to get insight into the effect of disorder on the LRN with the aim to realize a SRN for weak, but nonzero nonlinearity g ≪ 1.For this we now move to the disordered case to connect to Anderson localization [20,32,33].For g = 0, all normal mode eigenstates of Û are exponentially localized [21].The localization length, ξ depends only on the hopping-like parameter θ, Let us expand an arbitrary ⃗ Ψ(t) in the eigenmode basis at g = 0, Û (g = 0) ⃗ Here, the eigenmodes ⃗ Ψ k are Anderson-localized with components |ψ n k | ∼ e −|n|/ξ .The quasienergies ω k are real, and the expansion coefficients c k (t) are complex.The actions {|c k | 2 } are the constants of motion.Away from the integrable limit, i.e. for 0 < g ≪ 1, the expansion coefficients are coupled as where For small g and | π 2 −θ| the length ξ becomes larger than the system size, implying that the normal modes extend over the entire system, and Eq. ( 7) turns into a LRN with essentially all-to-all action couplings.For small g and θ ≪ 1 the length ξ tends to zero, the eigenmodes are strongly localized, and Eq. ( 7) turns into a SRN with essentially nearest neighbor action couplings.
We, therefore, predict that the universality class of a finite disordered system will depend on the ratio ξ/N .In Fig. 4 we show the LS dependence on g for (a) ξ/N = 0.2 and (b) ξ/N = 0.01.We clearly observe that case (a) To gain further insight, we extrapolate the rescaled LS down to g = 0 as a function of θ and fit the ansatz Eq. ( 5) to the outcome; see the appendix.The results are summarized in Fig. 5 and for α in Fig. A5 of appendix.Panel (a) shows the dependence ξ(θ) for convenience.Figure 5(b) exhibits the dependence of the exponent β on θ (filled squares), i.e. on the localization length ξ with a growth over several decades as expected when crossing over from a LRN to a SRN class.In Fig. 5(c) the dependence of the asymptotic rescaled KS entropy κ a versus θ is plotted (filled circles) showcasing the vanishing of κ a for small θ.For reference, we plot in panels (b) and (c) the corresponding results for the ordered case using open symbols.In this case no sizable change of β(θ) and κ a (θ) is observed.In addition to inducing the transition from a LRN to a SRN by varying the localization length ξ, this can be achieved by increasing the system size N .To demonstrate this, we have calculated the quantity κ a for different system sizes, all with the same localization length ξ = 2, exhibited in Fig. A6 of the appendix, yielding that κ a indeed decreases with increasing system size N .The Lyapunov spectrum (or better the inverse of its Lyapunov exponents) captures the timescale of thermalization.These times diverge in a qualitatively different way upon approaching integrable limits for short and long range network classes.Previous studies used the limit of weak lattice coupling to realize a SRN, while weak nonlinearities resulted in a LRN.Here, we show that the addition of disorder allows us to realize a SRN in the limit of weak nonlinearities, which in many situations correspond to weak two-body interactions.To achieve our goals we compute the rescaled Lyapunov spectrum and fit the ansatz Eq. ( 5) to it to extract an exponent whose divergence signals the presence of a SRN, and measure the rescaled KS entropy of the rescaled LS as another highly useful quantifier to tell SRN and LRN regimes apart.We also extrapolate the rescaled LS down to vanishing nonlinearity strength.These procedures allow us to unambiguously identify the crossover from a long range network to a short range network in a disordered system upon reducing the localization length.
Let us discuss and speculate about some consequences of the observed crossover.Approaching the integrable limit in the regime of a short range network implies that at any time the dynamics of the system is mostly regular, with rare local spots of nonlinear resonances leading to weak chaos (which is probably associated to the largest Lyapunov exponent).The density of these rare spots diminishes the closer the system is tuned to the integrable limit.Large but finite systems are therefore expected to display prethermalization features, where certain parts of the system show thermal properties which will fluctuate from part to part.Additional quantization of the considered systems might lead to a suppression of the chaotic resonances and the prethermalization, and ultimately to many-body localization related features.It remains to be studied whether the well defined short range network regimes of the studied classical systems will indeed result in many-body localization for the corresponding quantum systems, at variance to the long range network classes.
Our studies were confined to finite system sizes.We expect that the thermodynamic limit N → ∞ will result in a short range network class for all cases of finite localization length; see Sect.VII of the appendix.It remains to be studied whether large but finite values of the localization length will or will not result in qualitative changes of the Lyapunov spectrum scaling in the limit of weak nonlinearities and infinite system size.state ⃗ Ψ(t) = ⃗ x(t) + ⃗ w(t) follows the equations of motion Eq. (A1).
Because of the linearity property we can expand the norm in the exponential term in Eq. (A1).We only keep the first order of ⃗ w, Expanding the nonlinear part in a Taylor series and retaining only the first order of ⃗ w(t), we obtain Subtracting the contribution from the linear trajectories yields the growth of the perturbation  As we approach an integrable limit, the divergent timescales become inversely proportional to the LS.However, due to computational limitations, we are required to terminate our evolution at a specific number of time steps, typically ranging from 10 8 to 10 9 in this study.Despite this limitation, even at these time steps, certain parts of the LS may not have reached saturation.To quantify the saturation level, we calculate the standard deviation σ t of the data between the final time, 10 ns , and its previous "generation", 10 ns−1 .This interval typically covers more than 90% of the total evolution time.
To control the accuracy of our simulations, we utilize the parameter L = σt Λmax , where Λ max is the maximum Lyapunov exponent.We found out that the requirement L ≤ L max = 0.11 is sufficient to obtain reliable results for the Lyapunov exponents Λ i .

Asymptotic curve of rescaled Lyapunov spectrum
Fig. A2 depicts the extrapolation process used to obtain the results plotted in Fig. 5 of the main text.The values of the smallest and second smallest perturbation strength g were determined based on the required accuracy L ≤ L max .To obtain information on the Λ i when g approaches zero, we perform linear extrapolation of these two points.The red stars in Fig. A2 exemplary indicate the values of the Lyapunov exponents (LEs) that are used in the plots of asymptotic curves resulting from such an extrapolation procedure.These LEs capture the behavior of the system as g approaches zero and provide valuable insights into the system properties in the nearintegrable regime.
Using this extrapolation procedure, we have obtained the asymptotic curves for the rescaled Lyapunov spectra of both ordered and disordered cases, considering different values of θ.These curves are depicted in Fig. A3.The rescaled Lyapunov spectra provide valuable insights into the thermalization universality-class transition induced by disorder and the impact of localization on the system's behavior in the near-integrable regime.

Appendix E: Fitting of Lyapunov Spectrum
To further illustrate the distinct characteristics of the short-and long-range networks, we employ Eq. ( 5) from the main text to fit the rescaled LS.
Fig. A4 presents one example for the fit of the proposed ansatz Eq. ( 5) of the main text to the rescaled LS for the ordered case and one for the disordered case.This fitting procedure allows us to better understand the behavior of the system and the underlying thermalization universality-class transition induced by Anderson localization.The fitted curves provide valuable insights into the relationship between the localization length ξ and the thermalization behavior in the near-integrable regime.In addition to inducing the transition from a longrange network to a short-range network by varying the localization length ξ, we can also achieve this transition by increasing the system size N .We have calculated the quantity κ a for different system sizes, all with the same localization length ξ = 2, as presented in Fig. A6.
We observe that as the system size N increases, κ a decreases, providing further evidence for the long-to shortrange network transition.Notably, the error bars deduced from ensemble averages are of comparable magni- tude to the average value.This observation clearly indicates non-ergodicity and suggests that the maximum achievable time T s = 10 8 may not be sufficient for ergodization in the system.Note also that the ergodization time diverges much more rapidly in the short-range network regime [5,6], which will be explored in forthcoming works.

FIG. 1 .
FIG. 1. Schematic representation of the disordered unitary circuits map, where violet and brown circles indicate sites.Large green blocks represent unitary matrices Ĉ(θ), small yellow blocks indicate local nonlinear unitary operators Ĝ(g), and small blue blocks indicate local disorder unitary operators D({ϵn}).One time step contains four steps (unitary transformations).

FIG. 3 .
FIG. 3. Ordered case.Rescaled KS entropy κ in the LRN regime (red circles, top, θ = 1.13) versus g, and SRN regime (blue circles, bottom, g = 1) versus θ.Solid lines connect the data points and guide the eye.Here Ts = 10 8 , N = 200.The insets show the fit coefficients α (triangles) and β (squares) of Eq. (5) in the LRN (top left) and SRN (bottom right) regimes.The error bars represent the standard deviation deduced from an ensemble of 12 trajectories.The dashed curves in the main plot show the integral of Eq. (5) evaluated with the coefficients from the fits.

FIG. 5 .
FIG. 5. Disordered case after extrapolating down to g = 0. Ordered case results are added for comparison.(a) Localization length ξ versus θ.The black line is the analytical result of Eq. (6), and the blue circles mark the θ values used for the computations.(b) Fit coefficient β of the asymptotic rescaled LS Λ(ρ; g → 0) versus θ, filled blue squares.Open red squares are obtained from the corresponding ordered case data; see appendix.(c) Asymptotic rescaled KS entropy κa versus θ, filled blue circles.Open red circles are obtained from the corresponding ordered case data; see appendix.The error bars represent the time (σt, black) and ensemble (σg, blue and red) standard deviations; see appendix.Lines guide the eye.

Fig. A1 exemplifies
Fig. A1 exemplifies the evolution of the positive transient LS for one specific case.
Appendix D: Saturation criterion and Asymptotic curve of rescaled Lyapunov spectrum 1. Saturation criterion of Lyapunov spectrum FIG. A2.The colored lines show Λi versus g.The black dashed lines show the extrapolation process.The red stars show the points that are used in the plots of the asymptotic curves.
FIG. A3.Asymptotic curves for rescaled LS versus rescaled index ρ for the (a) ordered cases and (b) disordered cases in log scale.The angle θ varies from 1.43 (blue) to 0.38 (red).The error bars represent the time (σt black) and ensemble (σg colored) standard deviations.For all cases, the system size is N = 200.Unphysical error bars that would generate negative exponents are removed in panel (b).
Fig.A5displays the fit coefficients α obtained from Eq. (5) of the main text as a function of θ for the asymptotic curves shown in Fig.A3.In the plot, the red triangles connected by red dashed lines represent the fit coefficients for ordered systems, while the blue triangles connected by blue straight lines represent the fit coefficients for disordered systems.
Appendix G: Variation of κa with Different System Sizes FIG. A4.Asymptotic curves (red) for the rescaled LS versus rescaled index ρ for the (a) ordered and (b) disordered cases in log scale.The angle θ is equal to 0.38.The error bars represent the ensemble (σg) standard deviations.The blue dashed curves are the fittings.