Spectral Lyapunov exponents in chaotic and localized many-body quantum systems

We consider the spectral statistics of the Floquet operator for disordered, periodically driven spin chains in their quantum chaotic and many-body localized phases (MBL). The spectral statistics are characterized by the traces of powers t of the Floquet operator, and our approach hinges on the fact that, for integer t in systems with local interactions, these traces can be re-expressed in terms of products of dual transfer matrices, each representing a spatial slice of the system. We focus on properties of the dual transfer matrix products as represented by a spectrum of Lyapunov exponents, which we call spectral Lyapunov exponents. In particular, we examine the features of this spectrum that distinguish chaotic and MBL phases. The transfer matrices can be block-diagonalized using time-translation symmetry, and so the spectral Lyapunov exponents are classified according to a momentum in the time direction. For large t we argue that the leading Lyapunov exponents in each momentum sector tend to zero in the chaotic phase, while they remain finite in the MBL phase. These conclusions are based on results from three complementary types of calculation. We find exact results for the chaotic phase by considering a Floquet random quantum circuit with onsite Hilbert space dimension q in the large-q limit. In the MBL phase, we show that the spectral Lyapunov exponents remain finite by systematically analyzing models of non-interacting systems, weakly coupled systems, and local integrals of motion. Numerically, we compute the Lyapunov exponents for a Floquet random quantum circuit and for the kicked Ising model in the two phases. As an additional result, we calculate exactly the higher point spectral form factors (hpSFF) in the large-q limit, and show that the generalized Thouless time scales logarithmically in system size for all hpSFF in the large-q chaotic phase.

One of the fundamental goals of quantum statistical mechanics is to understand the basic hallmarks of chaotic dynamics. From a practical perspective, the presence of chaos is associated with memoryless evolution, so that the thermodynamic description is well-justified on the basis of the ergodic hypothesis. However, the classical notion of chaos does not extend directly to the quantum world, as Schrodinger evolution is linear and unitary and cannot admit diverging trajectories in Hilbert space 1 . Nevertheless, a large class of interacting many-body systems are believed to show quantum chaotic behaviour, as embodied in the eigenstate thermalization hypothesis (ETH) [2][3][4] . By contrast, many-body localization provides a generic mechanism which prevents the onset of chaos in arXiv:2012.05295v1 [cond-mat.stat-mech] 9 Dec 2020 FIG. 1. Summary of results: Behaviour is indicated using red lines for the quantum chaotic phase and blue lines for the MBL phase. (a): Schematic dependence of the SFF on t for many-body chaotic and MBL Floquet systems with system size L. The dashed red line is the RMT circular unitary ensemble (CUE) behavior. The generic behavior for many-body chaotic systems is characterized by two time scales: (i) the Thouless time t Th , which marks the onset of RMT behaviour in the SFF; (ii) the Heisenberg time tHei, which is of the order of the inverse of mean level spacing and so scales exponentially with L. For MBL systems with localization length ξ, the SFF grows quickly with t, reaching a plateau at a time O(q ξ ) that is independent of L. (b): Schematic dependence of the leading Lyapunov exponent λ> [see Eqns. (5) and (10)] on t in the two phases. At large t, λ> tends to zero in the chaotic phase and to a finite value in the MBL phase. (c): Schematic dependence at late times of the leading Lyapunov exponent λ0(k) in each momentum sector on the momentum eigenvalue k [see Eq. (12)]. In the chaotic phase λ0(k) converges to zero for all k. In the MBL phase λ0(k) converges to a smooth function of k for k = 0, with a distinct, larger value at k = 0. (d): Schematic dependence at late times of the scaled cumulant generating function Ft(α) on α [see in Eq. (3)]. The gradient of Ft(α) near α = 0 captures the L-dependence of fluctuations in the SFF. We find that Ft(α) has zero gradient in the chaotic phase and finite gradient in the MBL phase, reflecting fluctuations of the SFF that grow rapidly with L in the second case.
Random matrices 8 have long played a key role in providing minimal prototypes for properties of quantum chaotic systems. One important outcome is that spectral correlations have been identified as an indicator of chaotic behavior: as originally conjectured by Bohigas, Giannoni and Schmidt 9 , chaotic quantum systems exhibit the same spectral correlations as those of random matrices in the appropriate symmetry class. In particular, a distinctive fingerprint of quantum chaos is the presence of level repulsion between energy eigenvalues 10,11 .
Spectral fluctuations can be conveniently characterized via the Fourier transform of the two-point correlator of eigenvalues, known as the spectral form factor (SFF): Here W is the generator of the time evolution and W (t) denotes its t-power, while {θ m } are the spectral levels of the system under consideration (energies for systems with a time-independent Hamiltonian, or eigenphases of the Floquet operator for a periodically driven system).
The analysis of the SFF in many-body systems has recently been spurred on by the development of two novel approaches to Floquet models, where W generates the time evolution for a single period. Both in a long-range version of the kicked Ising model, 12,13 and in Floquet random circuits [14][15][16][17] in the limit of large local Hilbert space dimension, the average SFF K(t) ≡ K(t) (where . . . denotes the average over an ensemble of statistically similar systems) was shown to reproduce the RMT result for times t larger than a scale t Th known as the Thouless time (see Fig. 1a). It has been reported on the basis of analytical and numerical calculations that t Th diverges with the system size L in generic quantum systems [14][15][16][17][18] , with the exception of specific fine-tuned models in the absence of conservation laws 19 . For this reason, it is important to understand which features control the behavior of K(t) for intermediate times 1 t t Th , which can nevertheless be arbitrarily large in the thermodynamic limit. A simple argument suggests that, in this time regime, K(t) is typically exponentially large in L: because of locality of interactions, different portions of the system for t t Th have not had time to generate correlations of their eigenphases; as a consequence, the trace in (1) can be factorized into contributions from the Hilbert space of each decoupled region 15 . Moreover, in this regime, since the SFF is not self-averaging 20 and for many-body systems has exponentially large fluctuations in the system size, its average may not be sufficient to characterize its behavior.
In this paper, we study signatures of spectral statistics of quantum many-body systems with local interactions by using the fact that Tr[W (t)] can be expressed as a product of dual transfer matrices, each associated with a spatial slice of the system. The dual transfer matrix product is characterized by a set of Lyapunov exponents, which we dub the spectral Lyapunov exponents, and by an associated cumulant generating function. The dual transfer matrix product grows exponentially with system size: average growth rates are given by Lyapunov exponents and sample-to-sample fluctuations in growth rate are described by the cumulant generating function. There are several motivations for this approach. Knowledge of the Lyapunov exponents allows one to investigate both the spectral statistics of quantum many-body systems in the thermodynamic limit and spectral statistics at times earlier than t Th . In addition, knowledge of the cumulant generating function allows one to study fluctuations of the SFF. Finally, the study spectral Lyapunov exponents provide a different way of characterising localized systems already in the thermodynamic limit.
A summary of our results is as follows. At fixed time t, the spectrum of Lyapunov exponents can be organized into t momentum sectors, associated with the invariance under discrete time translations of the evolution operator. We characterize the behavior of the leading Lyapunov exponent in each sector, showing that there is a clear distinction at large time ruled by the ergodicity properties of the dynamics (see Fig. 1b and c): For chaotic systems, the largest Lyapunov exponent at each momentum sector converges to zero at large time, signaling the absence of exponential growth of K(t) with system size and the emergence of random matrix behavior in the spectral correlations. For many-body localized systems, the Lyapunov exponents remain non-zero at large time, with a limiting but non-universal form of their spectrum. We also discuss the fluctuations of the leading Lyapunov in the zero-momentum sector (see Fig. 1d) by introducing a (scaled) cumulant generating function. We argue that higher cumulants are not important except in some non-generic settings. These results are justified by considering two different models and a combination of analytical and numerical analyses.
As a side result, in the chaotic phase, we compute exactly the higher point spectral form factors (hpSFF) in the limit of large local Hilbert space dimension q → ∞ and thermodynamic limits. The hpSFF is closely related to other diagnostics of chaos. As an example, the out-oftime-order correlator 21-26 is known to be related to the hpSFF for local operators at late times 27 , and for global operators 28,29 . We also define the generalized Thouless times as the time after which hpSFF behaviour of a quantum many-body system reduces to the RMT result. We show that the generalized Thouless times derived from all hpSFF scale logarithmically in system size in the large q limit.
Our calculations complement earlier work that has been concerned with the averaged SFF and its relation to chaos and localization 15,30,31 . In Ref. [30], the growth rate of the ensemble-averaged SFF at fixed time was studied specifically for the kicked Ising model (see Sec. III B) across the many-body localization transition. The authors introduce an appropriate ensemble-averaged transfer matrix, study its symmetries, and discuss the role of the time-momentum operator. More recently a general picture was presented in Ref. [18] for the long-time behaviour of the ensemble-averaged transfer matrix, together with numerical results for a random quantum circuit in the ergodic phase. Behaviour of the ensembleaveraged transfer matrix across the MBL transition is discussed in Ref. [32]. In contrast to this previous work, our focus here is on the average of the log of SFF rather than of the SFF itself, and on the notion of spectral Lyapunov exponents and fluctuations in the growth rate of the dual transfer matrix product.
The remainder of this paper is organized as follows. In Sec. II we introduce the spectral Lyapunov exponents and cumulant generating function. In Sec. III we define two quantum circuit models which each display both a quantum chaotic phase and an MBL phase as a coupling parameter is varied. In Sec. IV, we compute exactly the Lyapunov exponents and the generating function for a random circuit model in the large-q limit. The results demonstrate that the leading Lyapunov exponent in the chaotic phase tends toward zero at large times. In Sec. V, we discuss the Lyapunov exponents for models of non-interacting systems, systems with small coupling and systems with local integrals of motion. In this way we argue that the leading Lyapunov exponent remains finite at large t in the MBL phase. In Sec. VI, we present numerical results for the two quantum circuit models. Within the limitations imposed by the maximum computationally accessible values of t, results are consistent with distinct types of behaviour in each phase as described above. Finally, we conclude and discuss the outlook in Sec. VII.

II. SPECTRAL LYAPUNOV EXPONENTS AND GENERATING FUNCTION
Consider an ensemble of disordered systems, each associated with a Floquet operator W which we assume spatially inhomogeneous due to the presence of local disorder. Let {θ m } be the eigenphases of W . We introduce the higher point spectral form factor (hpSFF) 28,33 as where · is the ensemble average, and the subscript denotes the system size L with periodic boundary conditions. For α = 1, we have the standard SFF, K L (t) = K L (t) . To study fluctuations of the hpSFF (which are exponentially large in L) in the thermodynamic limit, we introduce the scaled cumulant generating function As we will see below, the function F t (α) captures the large-L scaling of all cumulants of the SFF. Knowledge of it gives access to the large-deviation distribution of K(t). By definition, F t (α) is a convex function. The behavior of F t (α) can be analysed by considering a dual picture 34,35 , using a 90-degree rotation which exchanges space and time. To be more concrete without losing generality, we can represent W as a matrix-product operator, where the vertical bonds have the physical dimension q and the auxiliary horizontal ones have dimensionq. Then we can rewrite its trace in the dual picture as where the operators V i areq t ×q t matrices defined implicitly by the diagram in Fig. 2. To avoid confusion, we have written explicitly the Hilbert space where the trace is taken as a subscript and we set V (L) ≡ V L . . . V 1 . Since the operator W is inhomogeneous in space, the matrices V i are different one from the other and randomly distributed due to the presence of local disorder.
At this stage, one can proceed in two ways. One possibility is to perform the average over the disorder by considering 2α layers By using Eq. (4), this amounts in practice to computing the disorder average of 2α replicas of the single-slice transfer matrix V i for integer α. The resulting transfer matrix leads directly to K(t) α . Additionally, after averaging, the resulting transfer matrix is invariant under spatial translations and so it is sufficient to study a single slice, and its leading eigenvalues and associated eigenvectors. This approach was employed recently in several studies [13, 15, 18, 19, 30, 32, 36, and 37]. We will use this method to compute analytically F t (α) in the limit of large local Hilbert space dimension.
Another possibility is to consider the transfer matrix for a single layer W (t). This has the advantage for numerical calculations that its size (q t ×q t ) is smaller and independent of α. However, since there is no sense in averaging W (t), we have to study this transfer matrix for individual samples. That means at large L, V (L) is the product of many random matrices. The natural quantities that characterise this product are the Lyapunov exponents. More precisely, in order to define them, we note that the trace in Eq. (4) enforces periodic boundary conditions and homogeneity in time ensures that the matrices V i are invariant under translations in the time directions. There is therefore a momentum quantum number associated with the time direction. The spectral decomposition of V (L) can thus be organised into the different momentum sectors k = 2πj/t, with j = 0, . . . t − 1, in the form where λ . . are growth rates which have sample-to-sample fluctuations for finite L but converge with probability one to the spectral Lyapunov exponents with momentum k. We refer below to these growth rates as finite-size spectral Lyapunov exponents. The φ (L) a (k)'s are the corresponding phases, while | a and |r a are respectively the left and right eigenvectors, which are biorthogonal and normalized such that r a | b = δ ab . We find that the largest Lyapunov exponent always lies in the zero-momentum sector, so for convenience we denote Furthermore, for any finite t, there is always a gap ∆λ (L) between λ (L) > and the other Lyapunovs, so that at large L and that derivatives of F t (α) at α = 0 provide cumulants of the largest finite-size Lyapunov exponent: In particular, when L → ∞, we extract the average and variance Therefore, provided F t (α = 0) is not divergent, in the limit L → ∞, the distribution of λ > is concentrated on its mean λ > almost surely and the function F t (α) encodes its large deviations. In contrast with K L (t), the Lyapunov exponents are thus self-averaging. In general we will denote which defines the spectral Lyapunov exponents. In the following, we will study the k-dependence of the leading Lyapunov exponents λ 0 (k) and the fluctuations of λ (L) > as encoded by the generating function F (α) for chaotic and MBL systems.

III. MODELS
For our analytical and numerical analysis we will consider two main models: the random phase model (RPM) 15 and the kicked Ising model (KIM) 19,30 . Below, we summarise their definitions and main features.

A. Random Phase Model (RPM)
The RPM consists of q-state 'spins' arranged with nearest-neighbour coupling on a one-dimensional lattice. We use site labels n = 1 . . . L and orbital labels a n = 1 . . . q on the n-th site. The q L × q L Floquet operator W = W 2 · W 1 is a product of two factors.
generates rotations at each site n, with q × q unitary matrices U n chosen randomly and independently from the circular unitary ensemble (CUE). W 2 couples neighbouring sites and is diagonal in the basis of site orbitals. The phase of the diagonal elements is a sum of terms depending on the quantum states of adjacent sites, so that We take each ϕ an,an+1 to be an independent Gaussian random variable with mean zero and standard deviation , which effectively controls the coupling between neighbouring spins. For fixed q, the model exhibits a many-body localization transition as a function of 15,38 , with a critical value c separating an MBL ( < c ) from a chaotic phase ( > c ). We will employ this model for exact analytic calculations within the chaotic phase, in the limit q → ∞. Note that accessing the MBL phase in this limit is problematic as c → 0 when q → ∞. We will therefore complement the analysis with numerical studies at q = 3, for which the model has c ≈ 0.25 15 .

B. Kicked Ising Model (KIM)
The kicked Ising Model (KIM) is a Floquet Ising spin-1/2 chain defined by the time evolution operator W = W 2 · W 1 with with h j , J and b real parameters. Similarly to the RPM, this model has a many-body localization transition at a critical coupling strength J c = 0.23 [30], so that it exhibits a MBL phase for J < J c and a chaotic phase for J c < J ≤ π/4. This model has recently received a lot of attention, because of the existence of a "self-dual point" in the parameter space: |J| = |b| = π/4 and arbitrary local longitudinal fields h j . For these special values of the parameters, not only the evolution operator W , but also its duals V j (see Eq. (4) and Appendix B for the detailed definition) acting in the space direction, can be chosen to be unitary and with the same form of Eq. (15). In this case, at all times t, not only the average SFF 19 , but also its higher moments are in perfect agreement with the prediction of an appropriate random matrix ensemble which takes care of all the symmetries 13 . Indeed, unlike the RPM, this model is time-reversal invariant, and consequently, the behaviour of SFF for t Hei > t t Th is expected to follow the circular orthogonal ensemble (COE), which is given in the limit of large random matrices for t t Hei by At the self-dual point, additional discrete symmetries have been identified for the dynamics induced by Eq. (15), but they become irrelevant at large t 13,19,30 . Although solvable, the behavior at the self-dual point is not generic 39 , as it implies for instance that t Th does not diverge with the system size but remains O(1). Here, we will mainly use this model for numerical analysis without restricting to the self-dual point, taking advantage of its particularly small finite-time corrections near the selfdual point.

A. General behavior
We start by focusing on systems belonging to the CUE symmetry class and on the case F t (α = 1), which is simply related to the usual average of the spectral form factor K(t) . We make use of the defining property Eq. (3) to estimate the behavior of F t (α = 1) at large t in the chaotic phase. As observed in [15][16][17][18], for systems in the CUE symmetry class, the SFF approaches the random matrix prediction The specific details controlling the behavior t Th (L) are not yet fully understood, but in different set-ups 15-17 one expects t Th ∝ L ν , with ν > 0 40 . As a consequence, as already stated, the Thouless time t Th (L) → ∞ when L → ∞. Although F t (α) is formally defined only in the limit L → ∞, we expect it to capture well the finite-L behavior of K L (t) when t ∼ t Th (L). From Eq. (3), we can write and this suggests that, in order for the exponential growth in L of K L (t) to be suppressed, we must have F t Th (α = 1) 1/L. We thus deduce the scaling This argument can be extended to other values of α and we reach the conclusion that the chaotic phases must be characterized by In the next subsection, we will quantitatively justify this statement by computing explicitly F t (α) for the RPM in the limit q → ∞.
Additionally, we see that not only the leading Lyapunov in the zero-momentum sector λ > , but t Lyapunov exponents have to vanish in the large-t limit in order to reproduce the linear growth in time in (17). The most natural assumption is that the t vanishing Lyapunov exponents correspond to the different λ 0 (k) in the t momentum sectors. Similarly, for systems belonging to the COE symmetry class, in order to fulfill Eq. (16), we expect two vanishing Lyapunov exponents for t → ∞ in each momentum sector. We support these conjectures with numerical simulations in Sec. VI B. Ft(α) for RPM at q → ∞ As a solvable model of the chaotic phase in a spatially extended many-body quantum system we consider the RPM 15 , and compute analytically F t (α) and λ > in the large-q limit. We first of all consider integer values of α = n. We map the computation of K(t) n to the partition function of a one-dimensional statistical mechanical problem with nearest-neighbour interactions. As explained above, we introduce a transfer matrix in the space direction which allows the exact computation of K(t) n in the limit of large q: the value of F t (n) corresponds to the leading eigenvalue of the transfer matrix for the statistical mechanics problem, in a way that generalises the approach described in [15]. Lastly, we analytically continue F t (α) to non-integer α and obtain λ > .
To derive the transfer matrix for K(t) n , we construct the associated Hilbert space by performing the Haaraverage over W 1 for each site independently, as illustrated in Fig. 3. This independent averaging over W 1 is legitimate because the 1-gates are drawn independently across different sites, and because W 2 consists of diagonal 2gates only. Using the procedure explained in [14] and [15], we find a total of n! t n diagrams at each site in the limit of large-q. To each diagram we associate a state in the Hilbert space, labelled by a vector v in Z n t and by p = (σ(1), σ(2) . . . σ(n)) where σ belongs to the permutation group S n of n elements. Fig. 3a is the diagrammatic representation of K(t) α . Fig. 3b is the diagrammatic representation of a given site where each Haar-random 1-gate is represented by a single dot. 41 Upon averaging, the j-th loop (out of n loops) on the left is paired with the p j -th loop on the right in Fig. 3b. Furthermore, the pairing of j-th loop will have 1 out of t possible configurations, labelled by v j . Fig. 3c and d are two examples.
The average over W 2 in the large-q limit gives the matrix elements of the transfer matrix T which is constructed by counting the unmatched configurations and pairings between configuration (p, v) and (p , v ), since each unmatched configuration gives a factor of exp(− nt). As an example, the matrix element between the states in Fig. 3c and d is exp(−6 ), since n = 3, t = 2 and none of the pairings or configurations match. In summary, we have shown that the evaluation of K(t) n can be mapped to a one-dimensional statistical mechanical model where each site has n!t n states and where the interaction is defined by Eq. (21).
In Appendix A, we compute the leading eigenvalue E t (α) of T and analytically continue the result from in-teger n to arbitrary α ≥ 0 to obtain where Γ(a, b) denotes the incomplete Gamma function and we parameterize x = e −t . As a consistency check, at = 0, E t (α) = α!t α as expected since all the entries of the transfer matrix are unity. In general, we have the relation At large times (x 1), we obtain the expansion (24) Using the replica trick, the Lyapunov exponent can be computed as These analytical solutions are plotted in Figs. 4 and 5. Fig. 4 shows that F t (α) becomes flat as t (main panel) and (inset) increase, which implies λ > tends to zero for increasing and t. This behaviour of λ > in time is shown more explicitly in Fig. 5. Note that for small t (in particular t = 1) λ > is negative. We will see that this short-time feature also appears in the finite-q numerics. Th associated with the hpSFF as the time after which hpSFF behaviour (of a quantum many-body system) coincides with the RMT result. For the CUE in the large-q limit, the hpSFF is exactly α!t α due to the same diagrammatic approach explained in Fig. 3. The transfer matrix (21) becomes the identity matrix in the limit of large-t, and its trace gives the hpSFF CUE result as expected. To compute the t

V. THE MBL PHASE
In this section, we discuss the general features of F t (α) and of the Lyapunov spectrum λ 0 (k) in the MBL phase. In order to obtain some intuition, we first treat the case of uncoupled sites by analyzing the RPM at = 0. Then we consider systems with small coupling using a perturbative analysis applicable to both the RPM and KIM. Lastly, we analyze the leading Lyapunov exponent for an effective model of MBL in terms of local integrals of motion (LIOM).

A. Uncoupled sites
We use the RPM at = 0 as a toy model for the MBL phase. In this case W 2 is simply the identity and the model reduces to L non-interacting spins, each independently evolving with a random CUE matrix. From Eq. (2), we obtain for all moments where the average is performed within the CUE from which U is drawn. Given the trivial dependence of (26) on the system size L, we see from (7) that except for λ 0 (k = 0), all the other Lyapunov exponents (thus including all λ 0 (k = 0)) are degenerate with the value −∞: this is a general feature of models with uncoupled sites. From (3), we obtain an expression for F t (α) in terms of average of a single CUE matrix. In particular, for q = 2, we obtain the explicit formula lim t→∞ F t (α) = ln 4 α Γ(α + 1/2) Γ(α + 1) , q = 2 .
Note that the large time limit washes away many microscopic details and this expression holds more generally for non-interacting spins 1/2 with an arbitrary distribution of random fields, thus including the KIM at J = 0, as well as disordered free fermions in one dimension, i.e. the Anderson model 42 . For q > 2, one cannot get an analytic expression; nevertheless at large q but t q, one can use that Tr[U (t)] behaves as a gaussian-distributed complex random variable with zero average and variance q, leading to Note that in realistic models, the limit of large time is reached quite quickly, whenever t is larger than the single-spin Heisenberg time, i.e. t q = O(1). By contrasting (27) and (28) with (20), we observe a first indication of the different behaviour in a non-ergodic phase: F t (α) converges to a non-zero function at large t.
In the next sections we will see that this feature also characterises the MBL phase.
While accessing numerically the whole function F t (α) can be problematic, we will show in Sec. VI that the neighbourhood of α = 0 can be studied efficiently. Indeed, with the exception of special cases (e.g. for noninteracting spin 1/2, Eq. (27) leads to F t (0) t→∞ −→ 0), the behaviour F t (0) = 0 at large times provides a sufficient indication of a non-ergodic phase.

B. Perturbative analysis at small coupling
The dual transfer matrix provides an interesting framework in which to perform a perturbative expansion at small coupling between sites. The technique can be applied to both the KIM and the RPM, but we focus on the first. In Appendix B and C, we show that the transfer matrix V i corresponding to the two-layer structure introduced in Eqs. (15) can be written as and f (a) = arctanh(e −2ıa ). Note the resemblance with Eqs. (15) whose unitary form is recovered at the selfdual unitary point |b| = |J| = π/4 [19]. Here, we focus on b = π/4 and small J. The operator V 1 is easily diagonalised, and at small J the leading eigenstate is |0 = |+ . . . + , with σ x |± = ± |± . Every spin flip σ z j |0 is suppressed by a power of J. At the leading order in J, we thus truncate the Hilbert space of the trace in (4) to states only involving up to one spin flip σ z µ |0 , where we use Greek letters µ = 0, . . . , t − 1 to parameterise the position in time in the dual Hilbert space. Additionally, we employ the translational invariance in the time direction so that, within this truncation, we have a single magnon in each momentum sector k We can thus obtain an expression for λ 0 (k) for every k = 0, which takes the form (see Appendix C for the full derivation) where P (h) is the probability distribution of the random fields h i . A comparison between Eq. (31) and numerically exact results is shown in Fig. 7. Note that at first order in J, the time variable does not appear explicitly in Eq. (31). We can thus take the t → ∞ limit, where k becomes a continuous variable k ∈ [−π, π]. We leave for further investigation the study of the convergence of higher order corrections, but quite interestingly Eq. (31) provides an explicit result in the limits of both large times and large system sizes. The case k = 0 needs a different treatment because even at the leading order, O(J), the zero-momentum sector is two-dimensional, containing |0 and the zeromomentum magnon |k = 0 in Eq. (30). This fact is at the origin of the discontinuity observed in the spectrum at k = 0 (see Fig. 7). The resulting Lyapunov exponents λ 0 (0) and λ 1 (0) cannot be written analytically but can easily be computed numerically (see Appendix C).

C. Local Integrals of Motion
To describe the general behavior of the (fully) MBL phase, we consider an effective model based on the hypothesis that the MBL phase is characterised by an extensive number of LIOM with exponentially decaying interactions 43,44 , FIG. 6. The leading (solid lines) and sub-leading (dashed lines) Lyapunov exponents vs time for different ratios between the variances J 2 2,1 and J 2 1 . For any finite ratio J 2 2,1 /J 2 1 > 0, the leading exponent is finite at sufficiently large t. Note that the sub-leading Lyapunov at J 2 2,1 /J 2 1 = 0 converges to a large negative number and is not shown in the plot.
length, and || · || the operator norm. The τ z i provide an extensive set of integrals of motion that do not relax. Relaxation for real spins σ z i operators is thus induced by the accumulating random phases between different components of the system. This dephasing dynamics in MBL is the origin of logarithmic growth of entanglement 45,46 and power-law relaxation of local observables 47,48 .
For simplicity, we focus on the two-body model where J Furthermore, we will consider the simplest non-trivial LIOM in the main text where J 2,1 = 0 and J 2,r = 0 for all r > 1.
To analyze the behaviour of Lyapunov exponents in LIOM, we construct a 2 × 2 transfer matrix for all time t, . We numerically compute the two Lyapunov exponents using the method of QR decomposition described in Sec. VI. We see that for any finite ratio J 2 2,1 /J 2 1 > 0, the leading exponent λ > ≡ λ 0 and sub-leading exponent λ 1 converge to positive and negative finite values respectively, as shown in Fig. 6.
Moreover, we analyze K(t) α for the LIOM (32) with two-body terms for integer α in Appendix D. We map K(t) α to the partition function of stacked spin chains with two-body interactions, which can be written in terms of another transfer matrix, whose size increases as α increases and as we include longer range two-body terms in (33). We numerically diagonalize the transfer matrix and show that the F t (α) for integer α are qualitatively consistent with the Lyapunov exponents calculation above, and with the form of F t (α) computed for the RPM and KIM in MBL regime, as discussed below in Sec. VI.
We have used the LIOM picture to show that the leading Lyapunov exponent converges to a finite value as a function of time. We expect the existence of a positive finite λ > to persist for general LIOM with exponentially decaying support (see examples in Appendix D) and higher-body interaction terms. As one includes interaction terms of larger supports in the analysis, the size of the transfer matrix (34) and, consequently, the number of Lyapunov exponents increases. However, intriguingly, there is not a notion of time-momentum sectors for the Hamiltonian in Eq. (32) once expressed in the LIOM basis. This seems to indicate the possibility of a further structure for the LIOM effective Hamiltonian which would retain the notion of a time-momentum quantum number. We will leave the analysis of Lyapunov exponents for Hamiltonian systems for future studies.

VI. NUMERICS
The advantage of the dual formulation is that the Lyapunov exponents can be computed efficiently via an iterative procedure at arbitrarily large L. Indeed, by using the QR decomposition, we can write where Q i is an orthogonal matrix and R i is an upper triangular matrix. An estimate of the a-th Lyapunov exponent in the momentum sector k is then where we define for convenience By iteratively acting with the matrices V i and projecting onto the momentum sector k, we can generate a large number L of µ a,i . In this way we can obtain the behavior of F t (α) for α in the neighbourhood of 0. However, in order to access larger values of α 1, it is necessary to access values of λ where is chosen such that the spatial correlation between µ 0 and µ is sufficiently small. Our data suggest that, for both the RPM and KIM simulations, it is sufficient to have = 10, which we will take hereafter. We then define an effective cumulant generating function that approximates Eq. (3) as where · denotes the average over all realizations of consecutive µ's in (38). We can perform this numerical procedure exactly and the main limitation is represented by the exponential growth in the size of the matrices V i with t. Alternatively, one can adopt some approximate scheme based on matrix-product states (MPS) and the density-matrix renormalization group (DMRG) algorithm. However, we will see below that this is effective only deep in the MBL phase.
Using these methods, we compute the leading Lyapunov spectrum λ 0 (k), focusing in particular on two main representative cases λ > ≡ λ 0 (k = 0) and λ 0 (k = π) as functions of time t. We also extract the cumulant generating function F t, (α) in the chaotic and MBL phases. At late time in the MBL phase, we expect λ 0 (k) to have a non-uniform shape as a function of k with a positive finite λ > in the k = 0 momentum sector. In the chaotic phase, we expect the leading Lyapunov λ 0 (k = 0) to approach zero at late time, and we further conjecture that the largest Lyapunov exponents in the other momentum sector approach zero as well, so that λ 0 (k) is flat in the chaotic phase. Finally, we expect F t, (α) to have a finite positive gradient in the MBL phase, and to have zero gradient in the chaotic phase.
We summarize the result of numerics as follows: For the KIM, the data are in agreement with the theoretical expectations above. Note that, exactly at the self-dual point of the KIM, V i is unitary. Consequently, the SFF does not grow exponentially in space, and λ > is identically zero at the self-dual point. For this reason, even away from the self-dual point, the finite-time corrections are small. For RPM with on-site dimension q = 3, the data are compatible with the theoretical expectations, but agreement is not conclusive due to the limited times that are accessible within our numerics.
In Fig. 7 and 8, we show the largest Lyapunov exponents λ 0 (k) in each momentum sector k for the KIM and RPM respectively. For the KIM in the chaotic phase, λ 0 (k) is very small for all k. On the other hand, in the MBL phase, λ 0 (k = 0) is positive (except for very small J, see below), and λ 0 (k = 0) is negative. In Fig. 7 we include data for J as small as 0.001 and show that, for k = 0, it agrees well with the result from perturbation theory given in Eq. (31) (full equation in (C12)), and that for k = 0 it agrees with the result from degenerate perturbation theory evaluated numerically. Note that in Fig. 7 we observe a peculiarity in λ 0 (k = 0) for J = 0.001, 0.01, where the λ 0 (k = 0) have small negative values. We find that the window of J where λ 0 (k = 0) < 0 gets smaller as t gets larger and we expect this to be only a finite-time effect. For RPM, the data shown in Fig. 8 are limited by finite-t effects, but they are compatible with and seem to tend towards the expected behaviours. Next, in order to characterize the t-dependence of the spectral Lyapunov exponents, we focus on two distinctive cases: k = 0, π. In Fig. 9 and 10, we show λ > ≡ λ 0 (k = 0) and λ 0 (k = π) respectively as a function of t for the KIM. Consistently with our picture, in the chaotic phase both λ 0 (k = 0) and λ 0 (k = π) are small at large t. In the MBL phase, λ 0 (k = 0) converges towards a positive value while λ 0 (k = π) tends towards a finite negative value as t increases. Note that there are decay-ing oscillations in time with a periodicity of 4 for small J which are still visible at the accessible time with exact matrix multiplication (t ∼ 20). In order to access larger values of t, we employ a variation of the DMRG algorithm: after the application of each transfer matrix, we re-project the dual Hilbert space onto a matrix product state at fixed bond dimension χ. With this method, we can access much larger times (t ∼ 40) and confirm that the oscillations are suppressed in t, as shown in Fig. 11. However, the accessible values of ξ are limited by the necessity of using periodic boundary conditions in the time direction, and the non-unitarity of the dual transfer matrix. In the chaotic phase, the DMRG algorithm applied in the dual picture cannot be exploited for large t since the Lyapunov exponents obtained in this way do not converge for accessible values of χ.
In Fig. 12 and 13, we show λ 0 (k) against t for the RPM for k = 0 and π respectively. In the MBL phase, λ 0 (k) behaves as expected for both momentum sectors. However, the behaviour of λ 0 (k) in the chaotic phase is affected by the finite time effects. While λ 0 (k = π) for the chaotic phase tends towards zero and is small relative to the corresponding Lyapunov exponents in the MBL phase, λ 0 (k = 0) remains finite for the accessible values of t.
In Fig. 14 and 15, we show the cumulant generating function (39) computed for the KIM and RPM respectively. Recall that the first and second cumulants of λ > are the first and second derivatives of the cumulant generating function at α = 0. For the KIM, F t, (α) shows obviously distinctive behaviours in the chaotic and MBL phases. In particular, F t, (α) has zero derivative in the former phase, which is consistent with the expectation that λ > = 0, discussed in earlier sections. However, again, for RPM, F t, (α) does not show such a clear difference in behaviour between the two phases for the accessible t (Fig. 15).
Finally, we recall the different symmetry classes of the KIM and RPM, namely COE and CUE respectively. The former symmetry class has K COE ≈ 2t. Therefore, for the KIM, it is natural to expect in the chaotic phase at large times that there are 2t (not just t) zero Lyapunov exponents contributing to K(t) ∼ k,a e λa(k)L , two from each of the t momentum sector. In order to check this, we compute the gaps ∆λ a (k) ≡ λ a (k) − λ a+1 (k) in Appendix. E, and verify that ∆λ 0 (k) is indeed small at large t in the chaotic phase for the KIM. In RPM, the corresponding computation shows that the gap ∆λ a (k) is much larger.

VII. CONCLUDING REMARKS
We have proposed a new set of physical quantities, the spectral Lyapunov exponents, which allow us to explore the fluctuations and the generic behaviour of the SFF in the thermodynamic limit. We have shown that the spectral Lyapunov exponents have distinct long-time behaviours in the chaotic and MBL phases: For chaotic systems, the largest Lyapunov exponent in each momentum sector k converges to zero at large time, implying the absence of exponential growth of K(t) with system size  Our results for behaviour of the spectral Lyapunov exponents in each phase are complementary to and consistent with recent studies based on a transfer matrix that generates the average SFF 18,32 .
There are many interesting directions to pursue in the future. First, it would be exciting to look at the behavior of spectral Lyapunov spectrum when the MBL-ETH transition is approached and where universality is expected and could manifest itself both in the fluctuations F t (α) and the spectrum λ 0 (k). Second, it remains to understand how the existence of conserved quantities affects the behavior of the spectral Lyapunov exponents. One possible extension would be the inclusion of a U (1) charge conservation 16 . More generally, one could look at the behavior of Hamiltonian systems for which the energy provides a natural conserved quantity. In such cases, the time variable in the dual picture is continuous and the time momentum operator becomes a local conserved quantity in contrast to the Floquet case. This should be at the origin of the different scaling expected for the Thouless time in these systems.

VIII. ACKNOWLEDGEMENT
AC is supported by fellowships from the Croucher foundation and the PCTS at Princeton University. JTC is supported in part by EPSRC Grants EP/N01930X/1 and EP/S020527/1.
In this Appendix, we compute K(t) α for the RPM in the limit of large q and large L by obtaining the leading eigenvalue of the transfer matrix (21). Furthermore, we analytically continue the results to compute F t (α) and λ > (t) in the same limits.
To obtain the leading eigenvector of the transfer matrix (21) with integer α = n, note that all of its matrix elements are non-negative. So there is a unique largest real eigenvalue and a corresponding eigenvector with nonnegative components due to the Perron-Frobenius theorem. Furthermore, due to the symmetry of the diagrams, the eigenvector must be invariant under permutation, and hence we find (1, . . . , 1) T as the leading eigenvector.
To find the leading eigenvalue E 1 , we sum over any given row of T , and obtain where P (n, d) is the number of elements in S n with distance d from any given reference permutation 49 , say the identity p = (1, 2, . . . , n); C(t, , n, d) = v (1, . . . , n), (1, . . . , 1)| T |p, v is the sum of t n matrix elements at fixed p. From the leading eigenvalue in (A1a), we can then recover F t (n) = log E 1 . Above we derived an expression for K(t) α at integer α = n. Now we re-express Eq. (A1) in a different form where the dependence on α can be easily analytically continued to real values. First of all, we can rewrite the sum in (A1b) as where e is the Neper number. Plugging (A2) in (A1b), we can exchange the order of sums in (A1a) and perform the sum over d. After some manipulations, the final result takes the compact form valid for arbitrary α ≥ 0 and lim q→∞ F t (α) = log E 1 .
Leaving the large-q limit implicit, we can now compute the leading Lyapunov exponent The derivative of the incomplete Gamma function can be evaluated as ∂ ∂α Γ(α + 1, y) α→0 = e −y log(y) + Γ(0, y) .

(A6)
After some straightforward manipulations, we arrive Eq. (25), reproduced below, Appendix B: Explicit form of the dual circuit Here, we derive an explicit form for the dual transfer matrix for the two models introduced in Sec. III. Both models are composed of a layer W 1 of single-site unitaries and a layer W 2 of 2-site unitaries diagonal in the computational basis. We will therefore treat them both at once. To be more specific, we use the notation introduced in Sec. III A for the RPM in Eqs. (13,14), i.e.
For the RPM, the unitary matrices U j are drawn from the CUE and the phases φ an,an+1 are Gaussian variables with zero average and standard deviation . With the same notation, the KIM can be recovered setting q = 2, with U j = e ıhj σ z j e ıbσ x j and ϕ (n) an,an+1 = Je ıπ(an+an+1) (a n = 1, 2). In order to deduce the form of the transfer matrix in the space direction we write explicitly the trace in (4). We introduce a compact notation for the indices a = (a 1 , . . . , a L ) and we have We now introduce a dual Hilbert spaceH = ⊗ t µ=1 C q with dimensionÑ = q t and the computational basis b = {b 1 , . . . , b t } with each b µ = 1, . . . , q. Then, defining the j-dependent dual layers and V j = V 2,j V 1,j , we have that Note that in the dual formulation the 1-body unitary matrices in W 1 are converted into 2-body diagonal matrices in V 2 , while the 2-body phases in W 2 are converted into the 1-body V 1 .

Appendix C: Weakly coupled spins
In this Appendix we provide the details of the calculation of the Lyapunov spectrum in the limit where different sites are weakly coupled. This corresponds to J → 0/ → 0 respectively for the KIM/RPM. For the sake of clarity, we will focus on the KIM, although the discussion can be easily adapted to the RPM. From Eqs. (B3), we have where in the last equalities we used the matrix identity holding for any operator O 2 = 1 and f (a) = arctanh(e −2ıa ). Setting σ x |± = ± |± , we define At small J, the largest eigenvalue is associated with the vacuum ferromagnetic state |0 and spin flips are suppressed with powers of tan(J). At the leading order in J, we can restrict our Hilbert space to a single spin flip (M = 1 in (C4)). In order to compute the trace in Eq. (B4) in this limit, we need the matrix elements of V (j) 2 between pairs of single spin-flip states. They can be written explicitly by going back to the original time direction as where the trace is performed in the Hilbert space of a single spin. Additionally we can make use of the translational invariance in the time direction to decompose the trace in (B4) in momentum sectors. We thus define a spin wave with momentum k as The trace in the single spin flip of momentum k = 0 can then be written as L ] = (2 cos J) tL (ı tan J) L k|V 1 |k . . . k|V (2) L |k (C7) We deduce Setting we can rewrite which can be easily diagonalized and we arrive at the final expression At large t, we can make the replacement inside the integral ln | sin(tθ h )| 2 → −2 log 2 and for b = π/4, we get the final expression ∼ 2t ln | cos J| + 2 ln | tan J| + dhP (h) ln cos(h) 2 (2 − cos(h) 2 )(cos k + sin(h) 2 ) 2 .
(C12) with up to nextto-nearest-neighbour 2-body interactions, which is mapped to spin chains with additional next-to-nearest-neighbour (green) 2-body interactions. The associated Hilbert space has dimension 16. (c) Representation of K(t) 2 with nearest-neighbour 2-body interaction, which is mapped to four spin chains with 1-body (red) and nearest-neighbour terms (blue). Again, the Hilbert space dimension is 16.
For the zero momentum sector, instead two states can contribute to the trace, i.e. the vacuum |0 and the zeromomentum magnon |k = 0 . The trace in this sector can then be rewritten as In this Appendix we analyze K(t) α for the LIOM model (32) with 2-body nearest-neighbour terms for integer α. We map K(t) α to the partition function of stacked spin chains with 2-body interactions, which can be written in terms of a transfer matrix 50 . We numerically diagonalize the transfer matrix constructed from the LIOM and show that the results are qualitatively compatible with the numerical results from the RP and KIM model in MBL regime.
It is instructive to construct the transfer matrix for K(t) for (32) with nearest-neighbour 2-body terms, and then generalize the procedure for general 2-body terms and hpSFF. Before averaging, the argument of the (1st point) SFF is exp ıtJ (1) k (m k − n k ) + ıtJ (2) k,k+1 (m k m k+1 − n k n k+1 ) , With periodic boundary condition, K(t) | J2=0 = (2 + 2e −2J 2 1 t 2 ) L → 2 L at large t as expected. The evaluation of K(t) can be generalized to LIOM (32) with general (not just nearest-neighbour) 2-body terms. We take the variance of 2-body coupling between spins separated by r sites to be (J (2) i,i+r ) 2 = J 2 2 e −2(r−1)/ξ ≡ J 2 2,r . Using the same approach, the ensemble average becomes (D9) This is the partition function of a stack of two spin chains with 2-body interactions up to a distance of r max . Consequently, the Hilbert space associated with the transfer matrix is a tensor product of r max copies of on-site Hilbert spaces, and contains degrees of freedom labelled by (m k , n k , . . . , m k+rmax−1 , , n k+rmax−1 ), where n k , m k · · · = ±1. The cases of r max = 1 and 2 are illustrated in Fig. 16 a and b. The resulting transfer matrix has 4 rmax eigenvalues: a genuine MBL phase has an infinite number of non-trivial Lyapunov exponents which are recovered in the limit r max → ∞.
We can further generalize this approach to the evaluation of K n (t) with integer exponent n and with only 2-body nearest-neighbour terms. In this case we have K n (t) = m (1) ,n (1) ,m (2) ,n (2) ...
which is the partition function of 2n copies of spin chains with 2-body nearest-neighbour interaction, as illustrated in Fig. 16 c, so that the transfer matrix Hilbert space size grows as 4 n . We numerically diagonalize the transfer matrix, and plot the value of F t (α) in Fig. 17 for integer α up to α = 5. Although this approach does not allow analytical continuation of K n (t) , we see that the form of F t (n) is compatible with the expectation that λ > = F t (α = 0) is finite, as discussed in Sec. V.