Spectral statistics and many-body quantum chaos with conserved charge

TL;DR

This paper analytically and numerically studies spectral statistics in chaotic many-body quantum systems with a conserved charge, revealing diffusion-driven scaling and identifying the Thouless time using Bethe Ansatz techniques.

Contribution

It introduces an analytical approach to spectral form factors in systems with conserved charge using a minimal Floquet model and Bethe Ansatz, connecting spectral statistics to spin chain physics.

Findings

Spectral form factor exhibits diffusion-driven scaling for times less than the Thouless time.

Bethe Ansatz allows extraction of the Thouless time in the model.

Numerical results support the analytical predictions.

Abstract

We investigate spectral statistics in spatially extended, chaotic many-body quantum systems with a conserved charge. We compute the spectral form factor $K(t)$ analytically for a minimal Floquet circuit model that has a $U(1)$ symmetry encoded via auxiliary spin-$1/2$ degrees of freedom. Averaging over an ensemble of realizations, we relate $K(t)$ to a partition function for the spins, given by a Trotterization of the spin-$1/2$ Heisenberg ferromagnet. Using Bethe Ansatz techniques, we extract the 'Thouless time' $t^{\vphantom{*}}_{\rm Th}$ demarcating the extent of random matrix behavior, and find scaling behavior governed by diffusion for $K(t)$ at $t\lesssim t^{\vphantom{*}}_{\rm Th}$. We also report numerical results for $K(t)$ in a generic Floquet spin model, which are consistent with these analytic predictions.

Equations

$$K(t)\equiv \sum _{m,n=1}^{\mathcal{D}}⟨e^{ı\left({\theta }_m-{\theta }_n\right)t}⟩=⟨|\mathrm{Tr}[\hat{W}(t)]{|}^2⟩.$$

$$K(t,s)\equiv ⟨\underset{s}{\mathrm{Tr}}[\hat{W}(t)]\underset{s}{\mathrm{Tr}}[{\hat{W}}^{†}(t)]⟩$$

$$\underset{q\to \infty }{\lim }K(t,s)=\left|t\right|\underset{s}{\mathrm{Tr}}[{\hat{M}}^t],$$

$${\hat{\mathbb{T}}}_{j,j^{\prime }}^{}=\frac{1}{2}({\hat{1}}_{j,j^{\prime }}^{}+{\hat{\mathbb{P}}}_{j,j^{\prime }}^{}),$$

$$|S,k⟩=\frac{1}{\sqrt{L}}\sum _{j=1}^L{e}^{ıjk}{\hat{\sigma }}_j^{-}|S+1⟩,k=\frac{2\pi p}{L},$$

$$\underset{q\to \infty }{\lim }K(t,s)\underset{t\gg L^2}{=}\left|t\right|\left(1+e^{-\frac{4{\pi }^{2}t}{L^2}}+\dots \right).$$

$$\varphi (t,s)=-\underset{L\to \infty }{\lim }{L}^{-1}\ln \left[K\left(t,s\right)/|t|\right],$$

$$\varphi (t,s)=-\frac{C}{\sqrt{t}}+\frac{1}{2(2s+1)t}+\dots ,$$

$$\underset{t,L\to \infty }{\lim }\ln \left[K(t,s)/t\right]=\kappa (x,s).$$

$$\kappa (x,s)\underset{x\ll 1}{\sim }C{x}^{-1/2}\mathrm{and}\kappa (x,s)\underset{x\gg 1}{\sim }{e}^{-4{\pi }^{2}x}.$$

$$\kappa \left(x,s\right)=-\sum _{n\ne 0}\ln \left[1-e^{-xD{\left(2\pi n\right)}^2}\right],$$

$$\hat{W}=e^{-i{t}_4{\hat{H}}_4}{e}^{-i{t}_3{\hat{H}}_3}{e}^{-i{t}_2{\hat{H}}_2}{e}^{-i{t}_1{\hat{H}}_1},$$

$${\hat{H}}_1$$

$${\hat{H}}_3$$

$${\hat{H}}_2$$

$${\hat{\mathcal{R}}}_{a,b}^{}(\lambda )=\frac{\lambda }{\lambda +2ı}{\hat{1}}_{ab}^{}+\frac{2ı}{\lambda +2ı}{\hat{\mathbb{P}}}_{ab}^{},$$

$$\hat{\mathcal{T}}(\lambda )\equiv {\hat{\mathcal{R}}}_{1,a}^{}(\lambda -{\xi }_1){\hat{\mathcal{R}}}_{2,a}^{}(\lambda -{\xi }_2)\dots {\hat{\mathcal{R}}}_{L,a}^{}(\lambda -{\xi }_L),$$

$$\displaystyle\hat{\mathcal{T}}(\lambda)=\left(\begin{array}[]{ll}\hat{A}(\lambda)&\hat{B}(\lambda)\\ \hat{C}(\lambda)&\hat{D}(\lambda)\end{array}\right)\;.$$

$$\hat{M}=\underset{\delta \to ı}{\lim }\hat{F}{\left(ı-\delta \right)}^{-1}\hat{F}\left(ı+\delta \right),$$

$$\hat{M}{|{\lambda }_1,\dots ,{\lambda }_N⟩}_S=e^{-{\sum }_j\epsilon ({\lambda }_j)}|{\lambda }_1,\dots ,{\lambda }_N⟩,$$

$${\left(\frac{{\lambda }_j+2ı}{{\lambda }_j-2ı}\right)}^{L/2}=\prod _{\begin{array}{c}{j}^{\prime }=1\\ j^{\prime }\ne j\end{array}}^N\left(\frac{{\lambda }_j-{\lambda }_{j^{\prime }}+2ı}{{\lambda }_j-{\lambda }_{j^{\prime }}-2ı}\right),$$

$$\epsilon (\lambda )\equiv -2\ln \cos k(\lambda ).$$

Full text

Introduction. Statistical mechanics is a fundamental tool in understanding condensed matter systems, allowing for their description in terms of a few state variables, rather than thermodynamically many degrees of freedom. For quantum or classical systems in equilibrium with their environment, thermodynamics arises naturally from the exchange of conserved quantities with a thermal reservoir. For generic isolated systems, thermalization is not guaranteed, but rather must emerge dynamically. In recent years, there has been substantial theoretical D’Alessio et al. (2016) and experimental effort Bloch et al. (2008, 2012) to understand how many-body quantum systems in isolation equilibrate under their own dynamics to reproduce the familiar results of statistical mechanics. The Eigenstate Thermalization Hypothesis (ETH) Deutsch (1991); Srednicki (1994); Rigol et al. (2008) provides a universal mechanism for establishing ergodicity of isolated quantum systems: in its simplest formulation, it is based on the convergence of the expectation values of local observables in nearby energy eigenstates when the thermodynamic limit is considered. Indeed, this assumption is enough to recover thermal behavior from the long-time dynamics of many-body quantum systems. The notion of quantum ergodicity associated with ETH is intertwined with random matrix theory (RMT): first, quantum chaotic systems are characterized by an RMT eigenvalue distribution Bohigas et al. (1984); Guhr et al. (1998); second, their eigenfunctions can be understood as random vectors Berry (1977); Deutsch (1991); Srednicki (1994); Borgonovi et al. (2016). One consequence is that quantum thermalization is always associated with level repulsion between energy eigenvalues. In practice, this spectral rigidity has often been used as an efficient means to pinpoint quantum ergodicity breaking Rigol (2009); Santos and Rigol (2010); Biroli et al. (2010); Pal and Huse (2010); Alba (2015).

The validity and possible regimes of violation of ETH have been scrutinized in different types of chaotic systems Rigol and Srednicki (2012); Luitz and Lev (2016); Dymarsky (2018); Foini and Kurchan (2019a, b); Chan et al. (2018); however, numerical tests of ETH are challenging as they require the diagonalization of Hamiltonians whose size grows exponentially with the number of microscopic degrees of freedom Brandino et al. (2011); Beugeling et al. (2015); Kim et al. (2014); Luitz and Lev (2016). An obvious limitation is that, while RMT captures several aspects of quantum chaos, replacing the microscopic time evolution by a random matrix overlooks a key facet of the former, namely locality.

Recent efforts have endeavored to establish and improve upon minimal models of chaotic many body quantum systems starting from RMT, and enforcing locality via random local unitary gates to form ‘circuits’ Nahum et al. (2017); von Keyserlingk et al. (2018); Chan et al. (2018a, b, 2018); Rakovszky et al. (2018); Khemani et al. (2018); Sünderhauf et al. (2018); Kos et al. (2018); Gharibyan et al. (2018). Such models generally display quantum chaos, as characterized by entanglement entropy, the decay of local observables, and out-of-time-ordered correlation functions Nahum et al. (2017); von Keyserlingk et al. (2018); Chan et al. (2018a). By considering Floquet random circuits, it has been possible to derive analytically RMT spectral rigidity, in the limit of large local Hilbert space dimension Chan et al. (2018a, b) or at fine-tuned solvable points Kos et al. (2018); Bertini et al. (2018). More precisely, these works demonstrated that RMT behavior only appears for eigenvalue separations small on the scale of the inverse of the Thouless time, $t^{\vphantom{\dagger}}_{\rm Th}$, named in analogy with single-particle disordered conductors Thouless (1977); Altshuler and Shklovskii (1986). The value of $t^{\vphantom{\dagger}}_{\rm Th}$ depends on the linear system size $L$ and characterizes the time scale for the onset of quantum chaos. It remains an open question to understand which mechanisms control the scaling of $t^{\vphantom{\dagger}}_{\rm Th}$ with $L$.

In this work, we investigate the effect of a local conserved quantity $\hat{Q}$ on the behavior of $t^{\vphantom{\dagger}}_{\rm Th}$ and, more generally, on the spectral properties of the evolution operator, $\hat{W}$, of a Floquet system. In particular, we consider the two-point spectral form factor (SFF) $K(t)$, defined as

[TABLE]

Here $\{\theta_{m}\}$ are the eigenphases of $\hat{W}$, $\mathcal{D}$ is the Hilbert space dimension, $\hat{W}(t)$ indicates the $t$-th power of $\hat{W}$, and $\langle\ldots\rangle$ denotes the average over an ensemble of statistically similar systems. The SFF is the Fourier transform of the two-point correlation function of eigenphases. For uncorrelated eigenphases, $K(t)=\mathcal{D}$, while for random matrices belonging to the Circular Unitary Ensemble (CUE), $K_{\rm CUE}(t)=\left|t\right|$ until the Heisenberg time $t^{\vphantom{\dagger}}_{\rm Heis}=\mathcal{D}$, after which $K_{\rm CUE}(t)=\mathcal{D}$. The linear ramp is thus a fingerprint of level repulsion. For spatially-extended one-dimensional (1$d$) systems without a conserved density Chan et al. (2018b); Braun et al. (2019), $K(t)\simeq t^{L/\xi(t)}$ for $t\ll t^{\vphantom{\dagger}}_{\rm Th}$: the system can be seen as partitioned into $L/\xi(t)$ chaotic blocks, with a length $\xi(t)$ that grows with $t$. RMT behavior is recovered for $t\gtrsim t^{\vphantom{\dagger}}_{\rm Th}$ with $\xi(t=t^{\vphantom{\dagger}}_{\rm Th})\sim L$. In the presence of a conserved quantity with diffusive transport, it is natural to expect $t^{\vphantom{\dagger}}_{\rm Th}\sim L^{2}/D$, where $D$ is the diffusion constant. The idea that the timescale $L^{2}/D$ controls the onset of RMT spectral correlations was proposed on a heuristic basis in Gharibyan et al. (2018), with support from a variety of estimates and numerical studies. Here we establish this result in an exact treatment of a minimal model. We also set out the scaling behavior of $K(t)$ in the time interval $1\ll t\lesssim t^{\vphantom{\dagger}}_{\rm Th}$ and show that this holds in a computational study.

To probe $K(t)$ we build on a Floquet circuit model introduced in Chan et al. (2018a), consisting of a chain with $q$-state ‘spins’ at each site. The model has a time-evolution operator $\hat{W}$ constructed from unitary gates that act on neighboring pairs of sites. Gates are randomly selected in space but repeated periodically in time. Using a diagrammatic method to average over the individual matrices Brouwer and Beenakker (1996), to leading order at large $q$, one finds $K(t)=K_{\rm CUE}(t)$ for any $t\neq 0$, so that $t^{\vphantom{\dagger}}_{\rm Th}\to 0$ as $q\to\infty$. In the following, we formulate and characterize an extension of this model that hosts a $U(1)$ symmetry corresponding to a local, conserved operator $\hat{Q}$ that commutes with $\hat{W}$. In this way, the limit $q\to\infty$ has a twofold convenience: first, it allows controlled diagrammatic calculations; second, it washes out any effect on $K(t)$ not due to $\hat{Q}$.

Circuit model. The minimal model is a Floquet random unitary circuit (FRUC) defined on a chain of $L$ sites with local Hilbert space $\mathcal{H}_{\rm loc}\equiv\mathbb{C}^{q}\otimes\mathbb{C}^{2}$ – the tensor product of a $q$-dimensional color and a spin-$1/2$. The former facilitates Haar averaging Chan et al. (2018a, b); Khemani et al. (2018), and we encode a $U(1)$ symmetry in the latter, following Khemani et al. (2018), corresponding (in standard notation) to conservation of $\hat{Q}=\hat{S}^{z}=\frac{1}{2}\sum_{j=1}^{L}\hat{\sigma}^{z}_{j}$.

The single-period – or Floquet – evolution operator $\hat{W}\equiv\hat{W}^{\vphantom{*}}_{2}\cdot\hat{W}^{\vphantom{*}}_{1}$ is a depth-two circuit comprised of local two-site gates: assuming even $L$, the two layers correspond respectively to odd and even bonds, with $\hat{W}^{\vphantom{*}}_{1}=\hat{U}^{\vphantom{*}}_{1,2}\otimes\hat{U}^{\vphantom{*}}_{3,4}\otimes\dots$ and $\hat{W}^{\vphantom{*}}_{2}=\hat{U}^{\vphantom{*}}_{2,3}\otimes\hat{U}^{\vphantom{*}}_{4,5}\otimes\dots\hat{U}^{\vphantom{*}}_{L,1}$. We require that each $\hat{U}^{\vphantom{*}}_{j,j+1}$ preserves the local magnetization $S^{z}_{j,j+1}=\frac{1}{2}\left(\hat{\sigma}^{z}_{j}+\hat{\sigma}^{z}_{j+1}\right)$ Rakovszky et al. (2018); Khemani et al. (2018). Thus $\hat{U}^{\vphantom{*}}_{j,j+1}$ is a $4q^{2}\times 4q^{2}$ block diagonal matrix, acting as a $q^{2}\times q^{2}$ matrix in each of the $\uparrow\uparrow$ and $\downarrow\downarrow$ subspaces, and as a $2q^{2}\times 2q^{2}$ matrix in the $\uparrow\downarrow\,,\,\downarrow\uparrow$ subspace, with all three blocks independently drawn Haar random unitaries.

To characterize spectral correlations in this quantum circuit, we compute the SFF (1). Since $[\hat{W},\hat{S}^{z}]=0$, $\hat{W}$ is block-diagonal, and levels from different $\hat{S}^{z}$ sectors do not repel. Thus, we define

[TABLE]

where ‘$s$’ indicates restriction to the subspace $\hat{S}^{z}=S=L\,s$ 111For finite $L$, only a discrete set of values of $s$ are allowed. and $\left\langle\cdots\right\rangle$ denotes Haar averaging.

Effective spin-$1/2$ model.— The ensemble averaging in (2) maps $K(t,s)$ to the partition function of a Trotterized Heisenberg ferromagnet. Evaluating this average amounts to generating all diagrams Chan et al. (2018a) by pairing unitaries $\hat{U}^{\vphantom{*}}_{j,j+1}$ with their complex conjugates $\hat{U}^{\dagger}_{j,j+1}$ at each bond 222See supplementary material at [url].. As $q\to\infty$, the leading contributions come from $t$ diagrams, each of which has an identical ‘cyclical’ pairing at all sites Chan et al. (2018a). These diagrams can be expressed algebraically as (see Note (2) for details)

[TABLE]

where the factor of $\left|t\right|$ comes from there being $t$ such leading diagrams. The trace over the effective spin-$1/2$ evolution operator, $\hat{M}$, accounts for the sum over the color and spin degrees of freedom in a given leading diagram. Like $\hat{W}$, $\hat{M}\equiv{\hat{M}}^{\vphantom{*}}_{2}\cdot{\hat{M}}^{\vphantom{*}}_{1}$ consists of two layers: ${\hat{M}}^{\vphantom{*}}_{1}={\hat{\mathbb{T}}}^{\vphantom{*}}_{1,2}\otimes{\hat{\mathbb{T}}}^{\vphantom{*}}_{3,4}\otimes\dots$ and ${\hat{M}}^{\vphantom{*}}_{2}={\hat{\mathbb{T}}}^{\vphantom{*}}_{2,3}\otimes{\hat{\mathbb{T}}}^{\vphantom{*}}_{4,5}\otimes\dots$. $\hat{M}$ is hermitian, owing to contraction of a unitary and its conjugate, and is invariant under a shift by two sites due to ensemble averaging. The matrix $\hat{\mathbb{T}}^{\vphantom{*}}_{j,j^{\prime}}$ acts only on sites $j$, $j^{\prime}$ as

[TABLE]

where $\hat{\mathbb{P}}^{\vphantom{*}}_{j,j^{\prime}}=\frac{1}{2}(\hat{\mathbb{1}}^{\vphantom{*}}_{j,j^{\prime}}+\vec{\sigma}^{\vphantom{*}}_{j}\cdot\vec{\sigma}^{\vphantom{*}}_{j^{\prime}})$ is the ‘swap operator’. We note that $\hat{M}$ describes a discrete-time symmetric simple exclusion process (SSEP) for a classical lattice gas Schütz (2001). Although our original FRUC featured a $U(1)$ symmetry, after Haar averaging and taking $q\to\infty$, $K(t,s)$ exhibits an enlarged $SU(2)$ invariance in the remaining spin-$1/2$ variables; we believe this is specific to the large-$q$ limit. Additionally, as we clarify below, $\hat{M}$ belongs to a family of commuting transfer matrices, unveiling an emergent integrability, and the possibility of computing $K(t,s)$ exactly Vanicat et al. (2018).

This model leads to a Thouless time which scales diffusively. To see this, note $\hat{\mathbb{T}}_{j,j+1}\equiv\hat{\mathbb{1}}_{j,j+1}-\hat{\mathbb{H}}_{j,j+1}$, where $\hat{\mathbb{H}}_{j,j+1}=-\frac{1}{4}(\vec{\sigma}_{j}\cdot\vec{\sigma}_{j+1}-\hat{\mathbb{1}}_{j,j+1})$ describes the spin-$1/2$ Heisenberg ferromagnet. Thus, we can interpret ${\rm Tr}[\hat{M}^{t}]$ in (3) as a Trotterization of the partition function at inverse temperature $\beta=t$, i.e. ${\rm Tr}_{s}\,[\hat{M}^{t}]\simeq{\rm Tr}_{s}\,[e^{-tH_{\rm XXX}}]$, with $H_{\rm XXX}=\sum_{j}\hat{\mathbb{H}}_{j,j+1}$. Hence, the behavior of $K(t,s)$ at late times reflects the low-temperature properties of the Heisenberg ferromagnet, $H_{\rm XXX}$, which has $(L+1)$-fold degenerate ground states with vanishing energy. Each $S^{z}$ sector has a unique ground state, ${|S\rangle}\equiv(\hat{S}^{-})^{N_{\downarrow}}{|\uparrow\ldots\uparrow\rangle}$ with $\hat{S}^{\pm}\equiv\sum_{j}\hat{\sigma}^{\pm}_{j}$ and $N_{\downarrow}=L/2-S=L\left(1/2-s\right)$. Low-lying excitations above each ${|S\rangle}$ are magnons, i.e. plane-wave superpositions of spin flips

[TABLE]

characterized by a quadratic dispersion relation $\varepsilon(k)\propto k^{2}$ at small $k$. Expanding in $t\gg L^{2}$, one expects only the lowest energy magnon contributes and

[TABLE]

This suggests diffusive scaling of the Thouless time, $t^{\vphantom{\dagger}}_{\rm Th}\propto L^{2}$; a similar correspondence with $H_{\rm XXX}$ was established in Gharibyan et al. (2018) for random unitary circuits with a conserved density, lending support for the generality of this result. However, to investigate the regime $1\ll t\ll L^{2}$, we must consider states with extensive numbers of magnons, and many-body effects.

Scaling form. We define the function

[TABLE]

which can be computed exactly for any integer $t$, either by solving an infinite set of coupled integral equations—i.e. the Thermodynamic Bethe Ansatz (TBA) Takahashi (2005); Yang and Yang (1969)—or, perhaps more efficiently, via the ‘quantum transfer matrix method’ Suzuki (2003); Klümper (2004), which requires the solution of an algebraic equation in $2\left|t\right|$ variables Note (2). While the latter is better suited to calculating $K(t)$ at a particular time $t$, the former affords analytic insight into behavior at large times. Since the limit $L\to\infty$ implies $t^{\vphantom{\dagger}}_{\rm Heis}\to\infty$, we expand (7) about large $t$ using TBA,

[TABLE]

where the constant $C=\zeta(3/2)/\sqrt{4\pi}$ [$\zeta(z)$ is the Riemann Zeta function] and $q\to\infty$ is taken implicitly. Ignoring the Trotterized structure of $\hat{M}$ and taking $\hat{M}\sim e^{-H_{\rm XXX}}$ relates (8) to the low-temperature expansion of the specific heat close to the ferromagnetic ground state of $H_{\rm XXX}$ Schlottmann (1985); Takahashi (1986).

The form of (8) implies diffusive scaling even for $t\ll DL^{2}$. From the behavior $\phi(t,s)\sim(Dt)^{-1/2}$ it is apparent that the value of $D$ is independent of $s$, a consequence of the emergent $SU(2)$ symmetry at $q\to\infty$. However, the scaling limit relevant for $K(t)$ in the regime $1\ll t\lesssim t^{\vphantom{\dagger}}_{\rm Th}$ is distinct from that recovered from TBA (8): the former requires $t,L\to\infty$ with $x\equiv t/L^{2}$ fixed, while the latter requires the thermodynamic limit $L\to\infty$ at fixed $t$. Nevertheless, these results suggest a scaling form

[TABLE]

Despite the inherent integrability, exact calculation of $\kappa(x,s)$ is a challenging task. Nevertheless, its asymptotic behavior can be read off from (6) and (8): for early times ($x\ll 1$) one inserts (7) into (8); for late times ($x\gg 1$) one expands the log of the right side of (6). Thus

[TABLE]

By treating the magnons as non-interacting bosons, using ${{\rm Tr}}[{\hat{M}}^{\,t}\,]\approx{\rm Tr}[e^{-tH}]$, we recover Note (2)

[TABLE]

which precisely agrees with (10) if one uses the diffusion constant $D=1$ associated to the true dispersion (22) at small $k$. Although these predictions are obtained for $q\to\infty$, we expect their qualitative features to be valid for generic chaotic many-body systems with conserved charges. Numerical simulation. We now turn to numerical simulation to test (9) in chaotic quantum systems at finite $q$. At $q=1$, the FRUC considered above exhibits a numerically small diffusion constant that makes it difficult to avoid finite-size effects at the accessible values of $L$. Instead, we use a model adapted from Rakovszky et al. (2018), defined by

[TABLE]

where

[TABLE]

with periodic boundary conditions. We take $J^{1}_{z}=\left(\sqrt{3}+5\right)/6$, $J^{2}_{z}=\sqrt{5}/2$, and $J^{\vphantom{*}}_{xy}=\left(2\sqrt{3}+3\right)/7$, with $h^{1,2}_{n}$ drawn independently from the uniform distribution $\left[-1.0,1.0\right]$ for ensemble averaging. We choose $t_{1}=0.4$, $t_{2}=0.1$, $t_{3}=0.3$, and $t_{4}=0.2$ to avoid time-reversal symmetry around any instant in the period (we check that nearest-neighbor level statistics are CUE). In contrast to a recent study Šuntajs et al. (2019), we did not investigate $K(t)$ at strong disorder, as our concern is with the behavior of ergodic systems.

We restrict to half-filling, and measure $K(t)$ for sizes $L=12,14,16$, each averaged over $\gtrsim 10^{4}$ disorder configurations Note (2). Fig. 1 (lower panel) shows the general behavior of $K(t)$ in a system with a conserved charge, with values much larger than RMT during an initial interval, followed by a linear ramp regime in which $K(t)=\left|t\right|$, and finally a plateau for $t\geq t^{\vphantom{\dagger}}_{\rm Heis}$. Fig. 1 (upper panel) shows scaling collapse of $\ln\left[K\left(t\right)/t\right]$ versus $x=t/L^{2}$ for different $L$, following (10), with a comparison to (11) with $D$ a fitting parameter, here taken to be $0.05$.

Bethe Ansatz solution. We conclude by sketching the analysis of (3) using Bethe Ansatz. The computation of ${\rm Tr}_{s}\,[{\hat{M}}^{t}]$ would be simplified by knowledge of the full eigenspectrum of $\hat{M}$. Following the coordinate Bethe Ansatz for $H_{\rm XXX}$ Gaudin (2014); Franchini , one seeks multi-magnon eigenfunctions, i.e. plane-waves, along with a scattering matrix describing the exchange of the excitations’ momenta. However, this approach suffers from technical complications due to the circuit construction of $\hat{M}$. A more direct approach is instead based on the equivalent algebraic Bethe Ansatz formulation Korepin et al. (1997). Thus, we introduce the $R$-matrix

[TABLE]

which acts on spins-$1/2$ labelled $a$ and $b$. We next introduce the transfer matrix, which acts jointly on the $L$ physical spins and an auxiliary spin, $a$:

[TABLE]

where the rapidity $\lambda$ and inhomogeneities $\xi$’s are arbitrary complex numbers. Note that the subscripts in (15) label Hilbert spaces, as is customary in the Bethe-Ansatz literature (see also Fig. 4 in Note (2)). Schematically, $\lambda$ parameterizes the quasi-momentum $k(\lambda)=\operatorname{arccot}(\lambda/2)\in(-\pi/2,\pi/2)$ carried by the auxiliary particle while traversing the chain. Alternatively, in the auxiliary space, $\hat{\mathcal{T}}(\lambda)$ can be written as a $2\times 2$ matrix of operators acting on the physical spins,

[TABLE]

This construction is useful because the $R$-matrix in (14) satisfies the Yang-Baxter relation, implying a set of algebraic relations between the coefficients in (18), computed at the same inhomogeneities Faddeev (1996); in particular, setting $\hat{F}\left(\lambda\right)\equiv\hat{A}\left(\lambda\right)+\hat{D}\left(\lambda\right)$, one has $[\hat{F}\left(\lambda\right),\hat{F}\left(\lambda^{\prime}\right)]=0\,,~{}\forall~{}\lambda,\lambda^{\prime}$. The presence of a one-parameter family of commuting quantities establishes integrability for any choice of $\left\{\xi_{j}\right\}$, but only particular choices give rise to interesting local models. For instance, the isotropic Heisenberg spin chain is recovered for the homogeneous case, $\xi_{j}=\imath/2$. However, the brick-wall geometry relevant to $\hat{M}$ is realized via $\xi_{j}=\imath\,(1+\left(-1\right)^{j})$, from which it follows Note (2) that

[TABLE]

where the limit is needed to account for the non-invertibility of $\hat{\mathbb{T}}^{\vphantom{*}}_{j,j^{\prime}}$ in (4). The common eigenstates of the conserved quantities $\hat{F}(\lambda)$ (and thereby $\hat{M}$) can be obtained via algebraic properties, from which one can interpret $\hat{B}(\lambda)$ as an effective magnon creation operator on the vacuum ${|S\rangle}$, which decreases $\hat{S}^{z}$ by one. Thus $\hat{M}$ has eigenstates

[TABLE]

where ${|\lambda_{1},\ldots,\lambda_{N}\rangle}_{S}=\hat{B}(\lambda_{1})\ldots\hat{B}(\lambda_{N}){|S\rangle}$. The integer $N$ encodes the magnetization eigenvalue via $\hat{S}^{z}{|\lambda_{1},\ldots,\lambda_{N}\rangle}_{S}=(S-N){|\lambda_{1},\ldots,\lambda_{N}\rangle}$. Due to interactions, the parameters $\lambda_{1},\ldots,\lambda_{N}$ are not free, but satisfy

[TABLE]

the solution of which provides the full spectrum of $\hat{M}$.

The dispersion relation in (20) is given by

[TABLE]

At small $k$, a quadratic dispersion relation is recovered, as is expected since the discrepancies between $\hat{M}^{t}$ and $e^{-tH_{\rm XXX}}$ become irrelevant at long wavelengths. Although (20,21) simplify substantially the evaluation of (3), the main bottleneck remains the exponential growth in $L$ of the Hilbert space dimension. However, in the thermodynamic limit $L,S\to\infty$ with $s=S/L$ fixed, the solutions of (21) acquire a simple structure Note (2), and the function $\phi(t,s)$ in (7) can be computed analogously to thermodynamic quantities in integrable spin chains.

Discussion. We have presented analytical and numerical evidence showing the significance of the Thouless time $t^{\vphantom{\dagger}}_{\rm Th}=L^{2}/D$ for spectral correlations in systems with a conserved charge. These results are consistent with a scaling form $K(t)\sim\left|t\right|\exp\left[\kappa(t/L^{2})\right]$. We believe this form to be generic in describing the onset of chaos in quantum systems with a conserved quantity. We also provide a form (11) of the scaling function $\kappa$ (9) by neglecting the interactions between magnons in our FRUC, which is in good agreement with our numerical simulations. The question of the exactness and universality of (11) is an interesting topic for future study. Another interesting perspective for future work is the study of systems supporting non-Abelian symmetries, where the interplay of different conserved quantities can give rise to anomalous transport Das et al. (2018); Gopalakrishnan and Vasseur (2019); De Nardis et al. (2019). Finally, while this manuscript was in preparation, related numerical results were presented in Ref. 49 for the scaling of $t^{\vphantom{\dagger}}_{\rm Th}$ with $L^{-2}$ in a Hamiltonian (rather than Floquet) model.

Acknowledgments. We thank A. Nahum, T. Rakovszky, and B. Ware for useful discussions. This research was supported in part by the National Science Foundation via Grants DGE-1321846 (Graduate Research Fellowship Program) and DMR-1455366 (AJF), by EPSRC Grant No. EP/N01930X/1 (JTC) and the European Unions Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie Grant Agreement No. 794750 (A.D.L.). We also acknowledge the hospitality of the Kavli Institute for Theoretical Physics at the University of California, Santa Barbara (supported by NSF Grant PHY-1748958), and of the University of Massachusetts Amherst (AJF), where parts of this work were completed.