Integrable and chaotic dynamics of spins coupled to an optical cavity

TL;DR

This paper explores the transition between integrable and chaotic dynamics in a class of all-to-all coupled spin models realized with atoms in an optical cavity, revealing new integrable structures and their disappearance in the classical limit.

Contribution

It introduces a novel integrable structure with higher-body conserved charges in spin models coupled to optical cavities, and analyzes the transition to chaos in the large-spin limit.

Findings

Identifies parameter regimes for integrable and chaotic dynamics.

Discovers a new integrable structure with higher-body conserved charges.

Shows the disappearance of integrability in the large-spin classical limit.

Abstract

We show that a class of random all-to-all spin models, realizable in systems of atoms coupled to an optical cavity, gives rise to a rich dynamical phase diagram due to the pairwise separable nature of the couplings. By controlling the experimental parameters, one can tune between integrable and chaotic dynamics on the one hand, and between classical and quantum regimes on the other hand. For two special values of a spin-anisotropy parameter, the model exhibits rational-Gaudin type integrability and it is characterized by an extensive set of spin-bilinear integrals of motion, independent of the spin size. More generically, we find a novel integrable structure with conserved charges that are not purely bilinear. Instead, they develop `dressing tails' of higher-body terms, reminiscent of the dressed local integrals of motion found in Many-Body Localized phases. Surprisingly, this new type of integrable dynamics found in finite-size spin-1/2 systems disappears in the large-$S$ limit, giving way to classical chaos. We identify parameter regimes for characterizing these different dynamical behaviors in realistic experiments, in light of the limitations set by cavity dissipation.

Figures

Equations

$$H=\frac{1}{\sqrt{NM}}\sum _{i,j=1}^N{U}_{ij}{\mathbf{S}}_i\cdot {\mathbf{S}}_j,$$

$$H=\frac{1}{S\sqrt{N}}\sum _{i,j=1}^N{J}_i{J}_j\left(S_i^x{S}_j^x+S_i^y{S}_j^y+\Delta S_i^z{S}_j^z\right),$$

$$\mathcal{F}\equiv \sum _i{\xi }_i{\mathbf{S}}_i,$$

$$H\propto {\mathcal{F}}_x{\mathcal{F}}_x+{\mathcal{F}}_y{\mathcal{F}}_y+\Delta {\mathcal{F}}_z{\mathcal{F}}_z,$$

$$H_I=2\chi \sum _i\left({\mathbf{I}}_i\cdot \hat{c}\right)\left({\mathbf{S}}_i\cdot \hat{c}\right),$$

$$A_{\pm ,i}=\left(\frac{{\Omega }_i{e}^{-i\delta t}}{2\sqrt{2}}+g_i{a}_{\pm }\right)/g$$

$$H_I\approx \frac{i}{2}\chi \left({\xi }_i^{\ast }v{e}^{i\delta t}-{\xi }_i{v}^{†}{e}^{-i\delta t}\right)\left({\mathbf{S}}_i\cdot \hat{c}\right),$$

$${\xi }_i=\frac{{\Omega }_i{g}_i^{\ast }}{g^2}.$$

$$H=\frac{{\chi }^2}{4\delta }\left[{\cos }^2\theta {\mathcal{F}}^z{\mathcal{F}}^z+\frac{1}{2}{\sin }^2\theta \left({\mathcal{F}}^x{\mathcal{F}}^x+{\mathcal{F}}^y{\mathcal{F}}^y\right)\right],$$

$$Q^{(n)}=\sum _{i=1}^N{q}_i^{(n)},n=2,\dots ,N,$$

$$[Q^{(m)},Q^{(n)}]=0,[Q^{(m)},H]=0.$$

$$G^{(i)}=\sum _{j=1}^N\sum _{\alpha =1}^3{w}_{ij}^{\alpha }{S}_i^{\alpha }{S}_j^{\alpha },i=1,2,\dots ,N,$$

$$[G^{(i)},G^{(j)}]=0.$$

$$H=\sum _{i=1}^N{a}_i^{(0)}{G}^{(i)},Q^{(n)}=\sum _{i=1}^N{a}_i^{(n)}{G}^{(i)},n=2,\dots ,N,$$

$$w_{ij}^{\alpha }{w}_{jk}^{\gamma }+w_{ji}^{\beta }{w}_{ik}^{\gamma }-w_{ik}^{\alpha }{w}_{jk}^{\beta }=0,$$

$$G^{(i)}(\vec{J})=\sum _{j\ne i}\frac{J_i{J}_j}{J_i-J_j}{\mathbf{S}}_i\cdot {\mathbf{S}}_j.$$

$$H=\sum _{i,j=1}^N{J}_i{J}_j{\mathbf{S}}_i\cdot {\mathbf{S}}_j=\sum _{i=1}^{N}2J_i{G}^{(i)}(\vec{J})+J_i^2{\mathbf{S}}_i\cdot {\mathbf{S}}_i.$$

$$(w_{ij}^{\beta }{w}_{ji}^{\gamma }-w_{ji}^{\beta }{w}_{ij}^{\gamma })+2w_{ji}^{\alpha }(w_{ii}^{\gamma }-w_{ii}^{\beta })=0,$$

$$(w_{ij}^{\alpha }{w}_{ji}^{\gamma }-w_{ji}^{\alpha }{w}_{ij}^{\gamma })+2w_{ji}^{\beta }(w_{ii}^{\gamma }-w_{ii}^{\alpha })=0.$$

$${\tilde{G}}^{(i)}(\vec{J})=\sum _{j\ne i}\frac{J_i{J}_j}{J_i^2-J_j^2}$$

$$+$$

$$H=\sum _{i,j=1}^N{J}_i{J}_j(S_i^x{S}_j^x+S_i^y{S}_j^y)=\sum _{i=1}^{N}2J_i^2{\tilde{G}}^{(i)}(\vec{J}).$$

$${\hat{O}}_a\equiv \frac{3}{S(S+1)}{\hat{S}}_i^{\alpha }{\hat{S}}_j^{\alpha },$$

$$\frac{1}{\mathcal{D}}\mathrm{Tr}[{\hat{O}}_a^{†}{\hat{O}}_b]={\delta }_{ab},$$

$$⟨\hat{I}⟩=⟨\hat{I}(t)⟩\equiv \sum _a{u}_a⟨{\hat{O}}_a(t)⟩=\sum _a{u}_a\overline{⟨{\hat{O}}_a⟩}.$$

$$M_{a,t}\equiv \sqrt{\frac{1}{T}}\left(⟨{\hat{O}}_a(t)⟩-\overline{⟨{\hat{O}}_a⟩}\right).$$

$$L_{a,b}$$

$$=\sum _{l=1}^{3N(N-1)/2}{\sigma }_l^2{u}_{a,l}{u}_{b,l}.$$

$${\hat{Q}}_l\equiv \sum _{a=1}^{3N(N-1)/2}{u}_{a,l}{\hat{O}}_a,l=1,\dots ,3N(N-1)/2,$$

$$\frac{1}{\mathcal{D}}\mathrm{Tr}[{\hat{Q}}_l^{†}{\hat{Q}}_k]={\delta }_{lk}.$$

$$\overline{{⟨Q_l(t)⟩}^2}-{\overline{⟨Q_l(t)⟩}}^2={\sigma }_l^2,$$

$$G_l(t)=\frac{1}{\mathcal{D}}\mathrm{Tr}\left[{\hat{Q}}_l(t){\hat{Q}}_l(0)\right],$$

Full text

I Introduction

The dynamics leading to the eventual thermalization of closed quantum systems has become a topic of intense interest over the past few years. Significant progress has been made in describing the scrambling of information through quantum chaos, which allows effectively irreversible dynamics to emerge from unitary quantum time evolution. Notably, Maldacena et al., inspired by the chaotic properties of black holes, established that quantum mechanics places an upper bound on the Lyapunov exponent that characterizes the growth of chaos Maldacena et al. (2016). In a related development, Kitaev constructed a class of quantum many-body models whose dynamics saturates this bound on chaos Kitaev (2015); Maldacena and Stanford (2016) and can be related to black holes through the AdS/CFT correspondence Sachdev (2015); Jensen (2016); Stanford and Witten (2017). The fact that these models admit controlled solutions, despite being chaotic, has further conferred on them a paradigmatic status within the field of quantum dynamics.

Finding accessible systems which realize such models is therefore highly desirable, but also, a priori, very challenging: a common feature shared by all of these maximally chaotic, holographic models is that they lack spatial locality, since they couple together an extensive number of degrees of freedom. For instance, the Sachdev-Ye (SY) model Sachdev and Ye (1993) was originally proposed as a quantum spin model with random all-to-all couplings:

[TABLE]

where $\mathbf{S}_{i}$ are $\mathrm{SU}(M)$ spin operators. A fermionic variant, the Sachdev-Ye-Kitaev (SYK) model, was subsequently introduced by Kitaev.

While infinite-range spin interactions do not occur in magnetic materials, they can be realized rather naturally in cold atomic ensembles coupled to an optical cavity mode Black et al. (2003); Majer et al. (2007); Leroux et al. (2010); van Loo et al. (2013); Barontini et al. (2015); Hosten et al. (2016); Kollár et al. (2017); Léonard et al. (2017); Norcia et al. (2018); Kroeze et al. (2018); Landini et al. (2018); Kohler et al. (2018); Guo et al. (2018); Davis et al. (2019); Braverman et al. (2019). In this setup, the delocalized cavity mode mediates infinite-range interactions between the internal states of atoms through the local coupling at each site, regardless of the distance between atoms André et al. (2002); Sørensen and Mølmer (2002); Schleier-Smith et al. (2010); Gopalakrishnan et al. (2011); Hung et al. (2016); Masson et al. (2017); Mivehvar et al. (2017); Colella et al. (2018); Mivehvar et al. (2019). However, there is a crucial difference, already pointed out in Ref. Marino and Rey, 2018, between the interactions in the SY model and the ones mediated by the cavity. The second-order process that couples the atomic degrees of freedom via the cavity mode gives a separable (rank-$1$) matrix $U_{ij}=J_{i}J_{j}$, rather than the full-rank matrix assumed in the different variants of the SY model. Although non-separable interactions are, in principle, accessible in multi-mode cavities Gopalakrishnan et al. (2011); Strack and Sachdev (2011); Swingle et al. (2016); Kollár et al. (2017); Léonard et al. (2017); Marino and Rey (2018); Guo et al. (2018), separable all-to-all couplings are realized in numerous existing experiments Leroux et al. (2010); Hosten et al. (2016); Kroeze et al. (2018); Landini et al. (2018); Norcia et al. (2018); Davis et al. (2019); Braverman et al. (2019); Kohler et al. (2018) and arise generically for interactions mediated by a single bosonic mode.

Moreover, this ostensible limitation of the cavity-QED scheme turns out to be a boon: the separability of the interaction is responsible for an even richer dynamical phase diagram (see Fig. 1), which includes regions of chaos, Gaudin-type integrability characterized by spin-bilinear conserved quantities, and of a novel form of integrability—labeled Integrable∗ in Fig. 1—with quasi-bilinear integrals of motion.

The class of models we consider in this paper is described by the following quantum spin Hamiltonian:

[TABLE]

where $\mathbf{S}_{i}$ are $\mathrm{SU}(2)$ spin-$S$ operators encoded in the magnetic sub-levels of individual atoms or atomic sub-ensembles located at sites $i=1,\dots,N$. The site-dependent coefficients $J_{i}$ are determined by the local coupling of the atoms at site $i$ to the spatially-varying cavity mode, or by the local Rabi frequency $\Omega_{i}$ of an inhomogeneous drive field. The non-uniformity of the couplings $J_{i}$ is a crucial element of the models under consideration. For perfectly uniform couplings ($J_{i}=\overline{J}$), the model is integrable and exactly solvable in terms of the macroscopic spin $\boldsymbol{\mathcal{F}}=\sum_{i}J_{i}\mathbf{S}_{i}$. The $S^{x}S^{x}$ and $S^{y}S^{y}$ terms in (1.2) describe spin-exchange interactions between pairs of atoms, mediated by virtual cavity photons, while the $S^{z}S^{z}$ terms describe state-dependent ac Stark shifts. The normalization of $H$, which is not important for the dynamical properties, ensures that the high-temperature specific heat and free energy have a proper thermodynamic limit (see Supplementary Material Section S2.1).

The dynamical phases generated by this non-local spin model, shown in Fig. 1, are accessible via two experimentally tunable parameters. The spin-anisotropy parameter $\Delta$, controlling the relative strength of the spin-exchange and $S^{z}S^{z}$ interactions, can be tuned by changing the angle of an applied magnetic field $\mathbf{B}$ (see Fig. 1). In addition, it is possible to control the strength of quantum effects by changing the spin size $S$ on each site. While the choice of internal atomic states provides some flexibility in varying $S$, a larger range of spin sizes can be achieved by varying the number of atoms trapped at each site and letting $S^{\alpha}_{i}$ represent the collective spin of the sub-ensemble at site $i$. This enables the tuning of quantum effects from semi-classical dynamics at large $S$ all the way down to a spin-$1/2$ system that is dominated by quantum fluctuations. In combination with the possibility of varying the anisotropy $\Delta$, this tunability allows for a thorough exploration of the dynamical phase diagram.

The paper and the presentation of the various regimes shown in Fig. 1 are organized as follows. We provide a brief overview of these dynamical phases in Section II and we emphasize the novel features, which constitute our main results. In Section III we describe in detail the proposed experimental scheme to realize and control the couplings of the Hamiltonian (1.2). We also describe ways of inducing perturbations that go beyond separable interactions. In Section IV we begin the derivation of the main results. We analytically construct the integrals of motion that demonstrate the integrability of the dynamics at the special points $\Delta=0$ and $\Delta=1$. In Section V we present a computational method for finding integrals of motion using numerical or experimental data. In Section VI we deploy this technique and provide numerical evidence for the existence of a novel quantum integrable regime away from the special points $\Delta=0,1$. Specifically, we present an exact diagonalization study of the spin-$1/2$ model, showing that the integrable structure persists for anisotropy values $\Delta\neq 0,1$ away from the integrable points, with quasi-bilinear integrals of motion. In Section VII we simulate the classical model ($S\rightarrow\infty$) and show that it becomes chaotic, albeit in the presence of slowly decaying modes, away from the special points. In Section VIII we discuss experimental limitations and assess the extent to which the various features of the model are accessible in the presence of dissipation. Finally, in Section IX we comment on the implications of these results before concluding.

II Overview of the phase diagram

The best understood part of the dynamical phase diagram in Fig. 1 is the line at $\Delta=1$, for all spin sizes $S$, on which the Hamiltonian from Eq. 1.2 is equivalent to a rational Gaudin model Gaudin, M. (1976). This model is quantum integrable in the mathematical sense of possessing an underlying quantum group structure Sklyanin (1989). In the context of Gaudin-type models, quantum integrability is characterized by the existence of an extensive family of commuting, bilinear conserved quantities and there exist analytical expressions for each one. Even though there is no notion of spatial locality, the conserved quantities are “$2$-local” in the complexity theory sense Kitaev et al. (2002); Aharonov and Naveh (2002). By interchanging commutators with Poisson brackets, it follows that the integrable structure persists in the classical limit.

We find that the model is integrable at $\Delta=0$ as well. We obtain analytical expressions for an extensive family of conserved quantities that are also bilinear in spin. As in the case of $\Delta=1$, this integrability holds for any value of $S$, including the classical limit $S\to\infty$. The integrability of the model at $\Delta=0$ is connected to the existence of a non-standard class of Gaudin models Balantekin et al. (2005); Skrypnyk (2005); Schmidt (2008); Skrypnyk (2009a, b, c); Lukyanenko et al. (2016); Claeys et al. (2019).

However, the most surprising part of the phase diagram occurs away from these integrable points, i.e. in the regions $\{\Delta<0\}$, $\{0<\Delta<1\}$, and $\{\Delta>1\}$. There, we find a novel integrable structure that is markedly different from the type of integrability found at the two integrable points, $\Delta=0$ and $\Delta=1$. First, unlike the latter, integrability for $\Delta\neq 0,\,1$ appears to depend crucially on the spin size $S$. We show strong evidence that the model is integrable for a spin-$1/2$ system, while it is chaotic with a finite Lyapunov exponent $\lambda_{L}$ in the classical limit ($S\to\infty$). Nevertheless, in this latter limit, we also find that there exist modes that relax only on time scales much larger than $\lambda_{L}^{-1}$. We conjecture that this is a consequence of the “quasi-integrable” nature Goldfriend and Kurchan (2019) of the classical model. The putative transition from quantum integrable to (semiclassical) chaotic dynamics, schematically shown in Fig. 1, can be probed experimentally.

Second, the integrals of motion (IOM) of the $S=1/2$ model at $\Delta\neq 0,\,1$ are not bilinear (or 2-local), but may instead be termed quasi-bilinear. We present compelling numerical evidence that each IOM has appreciable support in the space of bilinear spin operators that does not depend on the system size $N$. The fact that the integrals of motion persist while developing tails of multi-spin terms on top of the dominant two-spin contribution is reminiscent of the quasi-local integrals of motion that characterize Many-Body Localized phases Vosk and Altman (2013); Serbyn et al. (2013); Huse et al. (2014); Chandran et al. (2015); Ros et al. (2015).

III Proposed Experimental scheme

As advertised, the full phase diagram of Fig. 1a can be accessed in experiments with atomic ensembles in single-mode optical cavities. In such experiments, each spin is encoded in internal states of an individual atom. The cavity generically couples to a weighted collective spin

[TABLE]

where each weight $\xi_{i}$ is set by the amplitudes of the cavity mode and drive field at the position of the $i^{\mathrm{th}}$ atom. Experiments to date have realized either Ising interactions Leroux et al. (2010); Hosten et al. (2016); Braverman et al. (2019); Landini et al. (2018); Kroeze et al. (2018) $H\propto\mathcal{F}_{z}^{2}$ or spin-exchange interactions Norcia et al. (2018); Davis et al. (2019) $H\propto\mathcal{F}_{+}\mathcal{F}_{-}$, in the latter case directly imaging the spatial dependence of the weights $\xi_{i}$ and the resulting spin dynamics Davis et al. (2019). We now show how to extend the approach of Ref. Davis et al., 2019 to realize generic XXZ models of the form

[TABLE]

where the anisotropy $\Delta$ is tuned by the angle of a magnetic field. An alternative approach to engineering Heisenberg models has been proposed in Ref. Mivehvar et al., 2019.

The experimental setup proposed here is shown in Fig. 1b. We consider spins encoded in Zeeman states of atoms whose positions in the cavity are fixed by a deep optical lattice. A magnetic field $\mathbf{B}=B\hat{z}$, which defines the quantization axis for the spins, is oriented at an angle $\theta$ to the longitudinal axis $\hat{c}$ of the optical cavity. Driving the atoms with a control field, incident either through the cavity or from the side, allows pairs of atoms to interact by scattering photons via the cavity. The interaction strengths are governed by the spatially dependent Rabi frequency $\Omega_{i}$ of the control field and vacuum Rabi frequency $2g_{i}$ of the cavity, where $i$ denotes the value for the $i^{\mathrm{th}}$ atom.

For large detuning between the atomic and cavity resonances, the atom-cavity interaction takes the form of a Faraday effect in which each atom couples to the Stokes vector $\mathbf{I}_{i}$, representing the local polarization and intensity of light. This effect is described by a Hamiltonian

[TABLE]

where $\chi$ is the vector ac Stark shift of a maximally coupled atom and the component of the Stokes vector along the cavity is $\mathbf{I}_{i}\cdot\hat{c}\ =(A^{\dagger}_{+,i}A_{+,i}-A^{\dagger}_{-,i}A_{-,i})/2$. The field operators

[TABLE]

include the quantum field $a_{\pm}$ of the cavity for $\sigma_{\pm}$-polarized modes, weighted by the local amplitude $g_{i}$ of the cavity mode, and displaced by a classical drive field with local Rabi frequency $\Omega_{i}$. The normalization is set by the vacuum Rabi frequency $2g$ of a maximally-coupled atom. We assume that the drive field has horizontal polarization $\hat{x}=\hat{z}\times\hat{c}$ and is detuned by $\delta$ from the cavity resonance.

In the limit where the drive field is weak and far detuned, we can obtain an effective Hamiltonian for the spin dynamics by adiabatically eliminating the photon modes. To this end, we first expand $H_{I}$ to lowest order in the operators $a_{\pm}$ to obtain

[TABLE]

where $v=\left(a_{+}-a_{-}\right)/\sqrt{2}$ represents the vertically polarized cavity mode, and we have introduced the weights

[TABLE]

These weights determine the collective spin $\mathcal{F}$ defined in Eq. 3.3, which couples to the cavity mode. Then, for $\langle v^{\dagger}v\rangle\ll 1$ and for large detuning $\delta\gg\kappa,\omega_{Z}$ compared to the cavity linewidth $\kappa$ and Zeeman splitting $\omega_{Z}$, we find that the effective spin Hamiltonian is

[TABLE]

as detailed in the Supplementary Material Section S1.1. We see that Eq. 3.9 matches the Hamiltonian (1.2) with couplings $J_{i}=\chi\xi_{i}S^{1/2}N^{1/4}\sin\theta/2\sqrt{2\delta}$ and anisotropy $\Delta=2\cot^{2}\theta$. Note that arbitrary control over the set of weights $\xi_{i}$ can be obtained by designing the spatial dependence of the control field.

In addition to the coherent dynamics generated by $H$ from Eq. 3.9, the cavity-mediated interactions are subject to dissipation due to photon loss from the cavity mirrors and atomic free-space scattering. Formally, these processes can be described by a family of Lindblad operators acting within a quantum master equation (see the Supplementary Material Section S1.2). The key parameter governing the interaction-to-decay ratio is the single-atom cooperativity $\eta=4g^{2}/\kappa\Gamma$, where $\Gamma$ is the atomic excited-state linewidth. Moreover, we find that the interaction-to-decay ratio is collectively enhanced, scaling as $S\sqrt{N\eta}$ for a system of $N$ sub-ensembles consisting of $S$ atoms each.

After we discuss the various properties and measurable signatures of chaotic and integrable dynamics in (1.2), we shall return to quantifying the effects of dissipation in Section VIII. In particular, we will estimate the atom number and cooperativity $\eta$ requisite for observing these signatures in the experimental setup.

IV Integrability at $\Delta=0$ and $\Delta=1$

In this section, we demonstrate the quantum integrability of the Hamiltonian (1.2) along the two lines at $\Delta=0$ and $\Delta=1$ in the dynamical phase diagram (Fig. 1). To place our discussion in context for the non-specialist reader, we begin by recalling some key features of integrable many-body systems. Broadly speaking 111For a comprehensive discussion of the subtleties involved in achieving a generally valid definition of quantum integrability, see Ref. Caux and Mossel, 2011., such systems are characterized by an extensive number of local conservation laws that give rise to exotic transport and thermalization properties. Important examples of quantum integrable systems include the Lieb-Liniger Bose gas and the spin-$1/2$ Heisenberg chain.

To illustrate the main ideas, consider a one-dimensional, local, quantum Hamiltonian $H=\sum_{n=1}^{N}h_{n}$, on $N$ lattice sites. For this type of model, integrability means the existence of $N-1$ independent, local charges,

[TABLE]

that commute with each other and with the Hamiltonian, namely

[TABLE]

The existence of extensively many local conservation laws can be regarded as a strong constraint on the dynamics of such systems, and leads to unusual physical effects such as non-dissipative heat transportZotos et al. (1997) and equilibration to non-thermal steady-statesRigol et al. (2008); Barthel and Schollwöck (2008).

In contrast with more standard integrable systems, the Gaudin-type models that arise in the present work are somewhat unusual, since they exhibit non-local couplings and are therefore essentially zero-dimensional. To construct these models, one starts from a set of $N$ operators,

[TABLE]

that are linear combinations of spin bilinears, with real coefficients $w^{\alpha}_{ij}\in\mathbb{R}$, and satisfy the defining commutation relations:

[TABLE]

The physical Hamiltonian and the independent conserved charges are then given by linear combinations of the $G^{(i)}$, of the form

[TABLE]

where the coefficients $a^{(n)}_{i}\in\mathbb{R}$ are elements of a non-singular $N$-by-$N$ matrix. Note that by the commutation relations (4.13), the Hamiltonian $H$ and its associated charges $Q^{(n)}$ automatically satisfy the commutation relations (4.11) required for integrability. Although these operators are not local, they are sums of spin bilinears and can therefore be regarded as “2-local” in the complexity theory sense.

We now show that the Hamiltonian Eq.  (1.2) defines a Gaudin-type integrable model for $\Delta=0$ and $\Delta=1$ and all values of spin $S$. Specifically, we will demonstrate that along these lines in the dynamical phase diagram Fig. 1, there exist $N-1$ independent, conserved and mutually commuting spin bilinears. The Hamiltonian at $\Delta=1$ is related to the rational Gaudin model Gaudin, M. (1976), which is well-known to be quantum integrable in the mathematically rigorous sense of possessing an underlying quantum group structure Sklyanin (1989). Meanwhile, the Hamiltonian at $\Delta=0$ lies in a less well-known class of “non-skew” Gaudin models, which arise from Gaudin’s equations upon relaxing the constraint of antisymmetry under interchange of site indices Balantekin et al. (2005); Skrypnyk (2005); Schmidt (2008); Skrypnyk (2009a, b, c); Lukyanenko et al. (2016); Claeys et al. (2019).

It will be helpful to review the problem first studied by Gaudin Gaudin, M. (1976): under what circumstances do a set of spin bilinears, as in Eq. (4.12), define a mutually commuting set, with $[G^{(i)},G^{(j)}]=0$? If the couplings $w_{ij}^{\alpha}\in\mathbb{R}$ are taken to be antisymmetric under interchange of indices, with $w_{ij}^{\alpha}+w_{ji}^{\alpha}=0$, then the $G^{(i)}$ mutually commute if and only if the Gaudin equations

[TABLE]

hold for all pairwise distinct $\{i,j,k\}$ and $\{\alpha,\beta,\gamma\}$. The isotropic solution $w^{\alpha}_{ij}=J_{i}J_{j}/(J_{i}-J_{j})$ defines the rational Gaudin Hamiltonians

[TABLE]

The all-to-all spin model from Eq. 1.2 at $\Delta=1$ is simply a linear combination of rational Gaudin Hamiltonians and Casimirs, to wit

[TABLE]

By rotational symmetry, $H$ conserves the total spin $\mathbf{S}_{\mathrm{tot}}=\sum_{i}\mathbf{S}_{i}$, and the linear span of the $G^{(i)}(\vec{J}\;)$ includes the squared spin $\mathbf{S}_{\mathrm{tot}}\cdot\mathbf{S}_{\mathrm{tot}}=\sum_{i,j}\mathbf{S}_{i}\cdot\mathbf{S}_{j}$. The mathematical structure of traditional Gaudin models has been studied in depth Sklyanin (1989); Feigin et al. (1994).

Let us now consider relaxing the constraint of antisymmetric couplings. Then Gaudin’s equations (4.15) must be augmented by two equations constraining “on-site” couplings, which read

[TABLE]

The model from Eq. 1.2 at $\Delta=0$ arises from an “non-skew XXZ” solution $w^{1}_{ij}=w^{2}_{ij}=J_{i}J_{j}/(J_{i}^{2}-J_{j}^{2}),\,w^{3}_{ij}=J_{j}^{2}/(J_{i}^{2}-J_{j}^{2})$ to the usual Gaudin equation (4.15), augmented by onsite terms $w^{1}_{ii}=w^{2}_{ii}=1/2,\,w^{3}_{ii}=0$, which solve Eq. 4.18. The corresponding Gaudin Hamiltonians read

[TABLE]

By the Gaudin equations Eq. 4.15 and Eq. 4.18, these mutually commute and the Hamiltonian (1.2) at $\Delta=0$ can be expressed as

[TABLE]

At spin-$1/2$, this coincides with the Hamiltonian obtained in Ref. Lukyanenko et al., 2016 or the “Wishart-SYK” model Iyoda et al. (2018), and can consequently be derived as a special case of the integrable spin-$1/2$ Hamiltonians considered in the recent work Ref. Claeys et al., 2019. The integrability of (4.20) for arbitrary spin $S$ was first discussed in Refs. Skrypnyk, 2009a, b, c (see also the references therein). We conclude that there is an integrable line in the phase diagram of the model (1.2) at $\Delta=0$. By rotational symmetry about the $z$-axis, this Hamiltonian conserves $S^{z}_{\mathrm{tot}}=\sum_{i}S^{z}_{i}$ and $(S^{z}_{\mathrm{tot}})^{2}$ lies in the linear span of the $\widetilde{G}^{(i)}(\vec{J}\;)$. Finally, we note that upon replacing commutators with Poisson brackets in the derivation of the Gaudin equations, the integrable structure identified for $\Delta=0$ and $\Delta=1$ remains unaltered in the classical limit ($S\rightarrow\infty$) of the Hamiltonian.

V Extracting integrals of motion from numerical or experimental data

Having characterized the integrable structure for $\Delta=0$ and $\Delta=1$, it is natural to ask whether the integrability of (1.2) extends to other, more generic, values of the anisotropy: can we find similar extensive sets of commuting bilinear conserved charges for $\Delta\neq 0,\,1$? To tackle this question in the absence of analytical tools, such as those used in the previous section, we develop a numerical method that enables the systematic search for bilinear (2-local) integrals of motion (IOM). We emphasize that this novel technique can be applied to either numerical or experimental data.

Let us first define a set of 2-local operators $\{\hat{O}_{a}\}$:

[TABLE]

where $i>j$ and the index $a$ is a shorthand notation for $(i,j,\alpha)$. We note that this family of $3N(N-1)/2$ operators defines an orthonormal set with respect to the infinite-temperature inner product:

[TABLE]

where $\mathcal{D}\equiv\mathrm{Tr}[\mathbf{1}]=(2S+1)^{N}$ is the dimension of the Hilbert space.

Now suppose that we can measure, experimentally or numerically, the time evolution of the expectation value $\langle\hat{O}_{a}(t)\rangle\equiv\left\langle\Phi\right|\hat{O}_{a}(t)\left|\Phi\right\rangle$, where $\left|\Phi\right\rangle$ is a random initial state (i.e. far from any energy eigenstate). A bilinear integral of motion $\hat{I}$ is a special linear combination of the $\hat{O}_{a}$ that remains constant in time, to wit

[TABLE]

Here and below, the overline denotes a time average, such as $\overline{\langle\hat{O}_{a}\rangle}\equiv\int_{0}^{T}\frac{dt}{T}\langle\hat{O}_{a}(t)\rangle$ over a time interval $[0,T]$. It is useful to recast the above equation in terms of the following time series matrix:

[TABLE]

Note that $M_{a,t}$ is a rectangular matrix with $3N(N-1)/2$ rows and a continuum of columns indexed by $t\in[0,T]$, where $TJ^{2}\gg 1$. In practice, the time axis is discretized such that the number of columns in $M$ is much larger than the number of rows. We immediately see that, by Eq. 5.23, a 2-local IOM corresponds to a left zero mode of $M$, i.e. $\sum_{a}u_{a}M_{a,t}=0$ for any $t$.

Thus, to find bilinear IOMs, we want to search for zero modes of $M$. More generally, we can consider the singular value decomposition (SVD) of $M$, or equivalently, the spectrum of the real Hermitian matrix

[TABLE]

In the second line, $\sigma_{l}\geq 0$ are the corresponding singular values of $M$ and $\sigma_{l}^{2}$ are the eigenvalues of $L$; $\vec{u}_{l}$ are the left singular vectors of $M$ and eigenvectors of $L$. Equivalently, $(u_{a,l})_{a,l=1}^{3N(N-1)/2}$ is a real orthogonal matrix, defining a family of operators

[TABLE]

which are also orthonormal:

[TABLE]

As mentioned above, $Q_{l}$ is an integral of motion if and only if $\sigma_{l}=0$. Furthermore, for small $\sigma_{l}>0$, we consider $Q_{l}$ to be approximately conserved and call it a “slow mode.” The rationale for this terminology comes from the identity

[TABLE]

which means that the singular value $\sigma_{l}$ is the standard deviation of the expectation value of $Q_{l}$ over the time interval $[0,T]$. A small $\sigma_{l}$ entails that $\langle Q_{l}(t)\rangle$ exhibits small fluctuations around its time-average value.

To summarize, we propose the following procedure: compute the time series matrix $M$, perform an SVD decomposition on $M$, analyze its singular values, and identify the possible IOMs and slow modes. In the next two sections, we use this method to characterize the behavior of the model along the $S=1/2$ and $S\rightarrow\infty$ lines in the phase diagram of Fig. 1, for anisotropies $\Delta\neq 0,\,1$. In Section VI, we numerically simulate the time evolution for the quantum spin-1/2 model and we further characterize the resulting slow modes by measuring their temporal auto-correlation functions. In Section VII, we simulate the dynamics of the model (1.2) describing classical spin degrees of freedom and, upon slightly modifying the above method, we extract the behavior of the auto-correlation functions directly from the singular values.

VI Integrability∗ for $S=\frac{1}{2}$

VI.1 Identifying integrals of motion

We now focus on the spin-$1/2$ system with up to $N=14$ sites and implement the technique proposed above. We initialize the system in a random product state 222We initialize each spin $i$ in a randomly chosen (with equal probability $p=1/6$) eigenstate of either $\hat{S}_{i}^{x}$, $\hat{S}_{i}^{y}$ or $\hat{S}_{i}^{z}$. $\left|\Phi\right\rangle$ and numerically compute the time evolution of the wavefunction with the Hamiltonian (1.2) via exact diagonalization. The random fields $J_{i}$ are sampled from the normal distribution $\mathcal{N}(0,J^{2})$ and we set $J^{2}=1$. We have checked that we obtain similar results for other distributions with zero mean and unit variance. We then record the expectation values of all the operators $\hat{O}_{a}$ defined in Eq. 5.21 and construct the time series matrix $M_{a,n}$ (defined in Eq. 5.24) at each discrete time $t_{n}=n\delta t$ with $\delta t=1\;J^{-2}$, integer $n$, and up to a maximal time $T=10^{3}J^{-2}$.

Fig. 2 presents results for the singular values of $M$ obtained for two values of $\Delta$ in a fixed disorder realization. As expected, at $\Delta=1$ we find $N+1$ vanishing singular values, in agreement with the analysis of Section IV. All other singular values lie above a gap of about $0.01$, indicating that there are no other 2-local integrals of motion beyond those identified in Section IV.

The results at $\Delta=0.9$, slightly away from the integrable point, are markedly different. We find only two exactly vanishing singular values corresponding to the space spanned by the two obvious integrals of motion, $H$ and $\left(S^{z}_{\mathrm{tot}}\right)^{2}$. This behavior persists on the entire open segment $\Delta\in(0,1)$, showing unambiguously that there are no other purely bilinear integrals of motion in this range. Nonetheless, we see that the remaining set of $N-2$ nontrivial IOMs at $\Delta=1$ are transformed, upon moving to the point $\Delta=0.9$, into left singular vectors with non-zero yet small singular values. It stands to reason that these small singular values correspond to operators that exhibit a slow decay because the system is close to the $\Delta=1$ integrable point. We now test this hypothesis by directly examining the decay of these putative “slow modes.”

VI.2 Characterizing the slow operators

We have seen that the nontrivial IOMs at the points $\Delta=0$ and $\Delta=1$ transform into a set of $N-2$ “slow operators,” indicated by small singular values, away from those two points. Let us examine the dynamics of these presumed slow modes. Their decay can be studied by numerically computing the auto-correlation functions

[TABLE]

where the normalization $\mathcal{D}=\left(2S+1\right)^{N}$ ensures that $G_{l}(0)=1$. For conserved modes, we expect the auto-correlation function to remain fixed at $G_{l}(t)=1$ for all time. For generic non-conserved operators, we expect $G_{l}(t)$ to decay to values near zero as these modes thermalize.

An example of the results for a system with $N=13$ sites and $\Delta=0.75$ is shown in Fig. 3. We see that the correlation functions related to the two zero singular values, $G_{1}(t)$ and $G_{2}(t)$, are perfectly non-decaying, as they must be. Also as expected, the correlation functions $G_{l}(t)$ associated with the small non-vanishing singular values ($3\leq l\leq N$) show a slow initial decay. However, the surprise is that, at very long times, these correlation functions saturate to a non-vanishing and rather appreciable value $g_{l}$. Fig. 4 shows that this phenomenon persists when varying $\Delta$ on the segment $[-0.5,1.5]$. Moreover, we find no significant size dependence of the saturation value $g_{l}$, as shown in Fig. 5a. We have also checked that the large plateau values are not due to the overlap between the slow modes $\hat{Q}_{l}$ with higher powers of the known conservation laws $\hat{H}$ and $\hat{S}^{z}_{\text{tot}}$, such as $\hat{H}^{2},\hat{H}^{3},\dots$, nor with projectors to energy eigenstates (for details, see the Supplementary Material S2.3). In contradistinction, the operators corresponding to higher singular values ($l\gg N$) decay to a vanishing, or very small, saturation value (see the Supplementary Material S2.4).

Altogether, in addition to the obvious bilinear IOMs, $H$ and $(S^{z}_{\mathrm{tot}})^{2}$, we find $N-2$ operators whose correlation functions saturate to an appreciable non-vanishing value. This result suggests that the model remains integrable even away from the Gaudin-like points $\Delta=0$ and $\Delta=1$: the bilinear integrals of motion are transformed into quasi-bilinear ones, which retain appreciable support in the space of 2-local operators. Based on the results shown in Fig. 4, we argue that this holds everywhere away from the integrable points, namely in the regions $\{\Delta<0\}$, $\{0<\Delta<1\}$, and $\{\Delta>1\}$. In general, we can write the new integrals of motion as bilinear operators dressed by a sum over higher, $2n$-local terms:

[TABLE]

where $Z_{l}$ is the weight of the integral of motion $I_{l}$ on 2-local operators. The saturation value of the auto-correlation function of $\hat{Q}_{l}$ that we plot in Fig. 4 is, essentially, $g_{l}\sim|Z_{l}|^{2}$. It would be interesting to further characterize how the coefficients $K_{i_{1}\dots i_{2n}}^{\alpha_{1}\dots\alpha_{2n}}$, which encode the overlap of the IOMs with the different $2n$-body spin operators, decay with increasing $n$. We leave this for future work.

The structure of the integrals of motion (6.31) is, in some ways, reminiscent of the Local Integrals of Motion (LIOM) in the Many-Body Localized (MBL) state Vosk and Altman (2013); Serbyn et al. (2013); Nandkishore and Huse (2015). The latter is characterized by quasi-local integrals of motion $\tau_{i}^{z}$ that are adiabatically connected to the microscopic degrees of freedom $\sigma^{z}_{i}$. As in our case, the LIOMs are dressed versions of the microscopic bits with weight on higher $n$-body operators decaying exponentially with $n$. There are, however, crucial differences from MBL. The integrals of motion in our case are not local, but rather extensive sums of bi-local operators. Hence, the IOMs of the all-to-all spin model do not facilitate a direct-product partition of the Hilbert space into single qubit spaces. Additionally, the integrability we observe does not depend on strong disorder—in fact, we found that its signatures are more pronounced as the couplings becomes more uniform, namely as $\mathrm{std}(J_{i})\lesssim\overline{J_{i}}$.

Lastly, we also find signatures of integrability in the spectrum of $H$: the level statistics are almost perfectly Poissonian at $\Delta=0,1$ and close to Poisson (although not exactly) at intermediate $\Delta$ (see Supplementary Material section S2.5). Nonetheless, for $0<\Delta<1$ we find many level crossings and the violation of the Wigner-von Neumann non-crossing rule represents further evidence of integrability despite the fact that there seems to be some degree of correlation between the energy levels Schliemann (2010); Owusu et al. (2008); Scaramazza et al. (2016).

VI.3 Perturbing away from integrability∗

After establishing the existence of a novel integrable structure for the spin-1/2 model, characterized by quasi-2-local IOMs, it is natural to investigate its robustness to perturbations away from the class of models (1.2) with separable disorder. This question is relevant from a theoretical point of view, but also from a practical, experimental perspective.

A natural perturbation to test in this context is one that adds a non-separable, SY-like, contribution to the interaction. Specifically, we add the term

[TABLE]

where the elements $V_{ij}$ are also sampled from a normal distribution $\mathcal{N}(0,1)$.

We explicitly check that at $\epsilon>0$ and $\Delta=1$ for $H+H_{\epsilon}^{(1)}$ there are only 4 zero singular values corresponding to exactly conserved and linearly-independent 2-local quantities: the Hamiltonian, $S_{\mathrm{tot}}^{2}$, $\left(S_{\mathrm{tot}}^{x}\right)^{2}$, and $\left(S_{\mathrm{tot}}^{y}\right)^{2}$. At intermediate $0<\Delta<1$, there are only two vanishing singular values corresponding to $H+H_{\epsilon}^{(1)}$ and $\left(S_{\mathrm{tot}}^{z}\right)^{2}$. Second, we verify that the lowest bilinear modes that are not exactly conserved (i.e. either the $l=3$ one at $0<\Delta<1$ or the $l=5$ one at $\Delta=1$) decay to smaller plateau values which decrease as we increase the system size $N$, as shown in Fig. 5a. This suggests that a perturbation $H_{\epsilon}^{(1)}$, even at $\epsilon\ll 1$, can spoil the integrability for a large system $N\gg 1$.

Another type of perturbation that arises naturally in the experimental set-up, due to the driving field, is represented by random stray magnetic fields along the $z$-axis:

[TABLE]

where the fields $h_{i}$ are also sampled from $\mathcal{N}(0,1)$. Note that $H+H^{(2)}_{\epsilon}$ has a single zero singular value corresponding to $\left(S_{\mathrm{tot}}^{z}\right)^{2}$ for all $\Delta$; this is due to the fact that the full Hamiltonian is no longer purely bilinear and that $H^{(2)}_{\epsilon}$ breaks the $\mathrm{SU(2)}$ symmetry at $\Delta=1$. Aside from this effect, the behavior upon perturbing with $H^{(2)}_{\epsilon}$ is similar to that obtained by perturbing with $H^{(1)}_{\epsilon}$, as shown in Fig. 5b.

Last, we consider the effect of adding the perturbation

[TABLE]

This term appears in the model

[TABLE]

which is similar to Eq. 1.2, but differs from it by the term $H^{(3)}$, arising due to the commutator $\left[S_{i}^{+},S_{i}^{-}\right]$. As noted in Ref. Marino and Rey, 2018, the model Eq. (6.35) is also experimentally accessible in a system of cold atoms interacting with cavity photons. It is clear that the perturbation $H^{(3)}$, having a $1/\sqrt{N}$ normalization, is sub-extensive and will not matter in the thermodynamic limit. Moreover, we find that it does not qualitatively affect the integrability of our quantum model even for the small systems considered in ED (see the Supplementary Material S2.6 for the numerical results).

In sum, our numerical analysis of the response to perturbations indicates that the novel integrability of the spin-1/2 model (1.2) is not particularly robust to non-separable interactions or stray magnetic fields. Nevertheless, in a finite-size system and at finite times (see Section VIII for more details), there are signatures of proximate integrability, as shown by the finite saturation values in Fig. 5.

To recapitulate our study of the dynamical phase diagram Fig. 1 thus far, we have found that the system is integrable along three lines: at $\Delta=0,1$ for any value of the spin size $S$ (characterized by bilinear IOMs), and at $S=1/2$ for any $\Delta\neq 0,1$ (characterized by quasi-bilinear IOMs). The remaining line in the phase boundary of Fig. 1 corresponds to the classical, $S\rightarrow\infty$, limit of the model (1.2), which we now discuss.

VII Chaos for $S\to\infty$

Since Gaudin-type integrability at $\Delta=0,\,1$ persists for all values of the spin size $S$, it is natural to ask whether the integrability* structure at $S=\frac{1}{2}$, presented in the previous section, also survives for larger values of $S$. Although it is numerically challenging to extend the exact diagonalization study of the previous section to intermediate $S$, the limit $S\to\infty$ leads to classical equations of motion that are amenable to numerical simulation.

These simulations allow us to analyze another boundary in the phase diagram, namely the $S\to\infty$ line, where we find chaotic dynamics with a finite Lyapunov exponent, as explained in Section  VII.1. The presence of chaos in the infinite-$S$ limit clearly implies that the $S=\frac{1}{2}$ integrability∗ does not extend to all $S$, unlike the Gaudin-type integrability at $\Delta=0$ and $1$. Remnants of an integrability∗ structure can nevertheless be revealed by applying the SVD analysis of Section V to the classical dynamics, which we do in Section VII.2. This technique reveals the presence of a large number of slow modes, which are known to occur classically in “quasi-integrable” systems, i.e. chaotic systems in the vicinity of integrable points. We characterize these slow modes in Section VII.3.

VII.1 Classical chaos

In the infinite-$S$ limit, the model (1.2) behaves as a classical system of coupled spin degrees of freedom $S_{i}^{\alpha}$ on the unit sphere, whose Hamiltonian dynamics can be written in terms of Poisson brackets:

[TABLE]

where

[TABLE]

For our numerical investigation, we sample the random fields $J_{i}$ from the uniform distribution $\left[-J,J\right]$ and set $J=1$ (we choose a bounded distribution to avoid large $J_{i}$’s that could cause numerical instabilities). The classical spin variables $S_{i}^{\alpha}$ obey

[TABLE]

We shall probe the infinite-temperature dynamics of this classical system by direct numerical simulation.

In order to study chaos, we use the standard tangent space method Benettin et al. (1980) to study the divergence of classical trajectories and measure the leading Lyapunov exponent. Let $\mathbf{S}(t)=\left(\mathbf{S}_{1}(t),\dots,\mathbf{S}_{N}(t)\right)$ denote the $3N$-dimensional vector describing the directions of all the spins at time $t$. We initialize the system in a random infinite-temperature state $\mathbf{S}(0)$, within which each spin points in a random direction, uniformly distributed on the unit sphere $S_{2}$. We also keep track of the trajectory of the deviation vector $\delta\mathbf{S}(t)$, which lives in the tangent space of $S_{2}\times\dots\times S_{2}$ at the point $\mathbf{S}(t)$; we further set $\delta\mathbf{S}(0)$ such that $\delta\mathbf{S}_{i}(0)\perp\mathbf{S}_{i}(0)$ for all spins and $\left\|\delta\mathbf{S}(0)\right\|^{2}=\sum_{i,\alpha}\left(\delta S_{i}^{\alpha}\right)^{2}=1$.

If we define the local effective field $\mathbf{F}_{i}=\left(F_{i}^{x},F_{i}^{y},F_{i}^{z}\right)=\frac{1}{S}\left(-\frac{\partial H}{\partial S_{i}^{x}},-\frac{\partial H}{\partial S_{i}^{y}},-\frac{\partial H}{\partial S_{i}^{z}}\right)$, we see that the Hamilton equations of motion (7.36) can be written as

[TABLE]

For our model (7.37) we have $F_{i}^{x,y}=\frac{J_{i}}{\sqrt{N}}\sum_{j}J_{j}S_{j}^{x,y}$ and $F_{i}^{z}=\Delta\frac{J_{i}}{\sqrt{N}}\sum_{j}J_{j}S_{j}^{z}$.

We immediately see that the variational equations of motion for the deviation vector $\delta\mathbf{S}$ can be written as

[TABLE]

where $\delta F_{i}^{x,y}=\frac{J_{i}}{\sqrt{N}}\sum_{j}J_{j}\delta S_{j}^{x,y}$ and $\delta F_{i}^{z}=\Delta\frac{J_{i}}{\sqrt{N}}\sum_{j}J_{j}\delta S_{j}^{z}$.

We numerically integrate the coupled differential equations (7.39) and (7.40) to find the trajectory $(\mathbf{S}(t),\delta\mathbf{S}(t))$ in the tangent bundle up until a time $T=500J^{-2}$ in increments $\delta t=1\;J^{-2}$. We then compute the sensitivity, defined as $d(t)=\left\|\delta\mathbf{S}(t)\right\|^{2}$, or in full,

[TABLE]

Note that $d(0)=1$, since we have normalized the initial deviation vector. For an integrable system we expect $d(t)$ to exhibit a power-law dependence on time; the flow on invariant tori specified by the $N$ conservation laws is linear in time and, since we have defined the sensitivity as $\left\|\delta\mathbf{S}(t)\right\|^{2}$, we expect $d(t)\sim t^{2}$. In a chaotic system $d(t)$ should increase exponentially with $t$. In Fig. 6a, we average over disorder realizations $\{J_{i}\}$ and initial states $\{\mathbf{S}(0)\}$ to find $\exp\left(\left<\log d(t)\right>_{J,\mathbf{S}}\right)$. We find that the classical system exhibits chaotic dynamics and an exponential divergence of trajectories in the regions $\{\Delta<0\}$, $\{0<\Delta<1\}$, and $\{\Delta>1\}$. We also find integrable dynamics and a power law divergence of trajectories at the special points $\Delta=0,1$.

Moreover, using the multiplicative ergodic theorem, we can define the maximal Lyapunov exponent Benettin et al. (1980) as

[TABLE]

Using the normalization $\left\|\delta\mathbf{S}(0)\right\|=1$ and our definition of the sensitivity from Eq. 7.41, we see that

[TABLE]

In practice, we compute the Lyapunov exponent by fitting a line $\lambda t+b$ through the late time behavior of $\log d(t)$, as discussed in Ref. Goldhirsch et al., 1987. In Fig. 6b we plot the Lyapunov exponent $\left<\lambda\right>_{J,\mathbf{S}}$, averaged over disorder realizations $\{J_{i}\}$ and initial states $\{\mathbf{S}(0)\}$, as a function of the anisotropy $\Delta$ and find that the system exhibits the most chaotic behavior (largest Lyapunov exponent) at $\Delta=1.5$. Second, we find that $\left<\lambda\right>_{J,\mathbf{S}}$ tends to a finite value for large system sizes $N$, as shown in the inset of Fig. 6b.

VII.2 SVD analysis

Although the presence of chaos in the classical dynamics excludes proper integrability in the infinite-$S$ limit, it does not rule out the possibility of “quasi-integrability,” whereby some operators have very slow decay. We investigate this possibility by applying the SVD analysis of Section V to the classical dynamics. This allows us to determine the number of exactly conserved quantities, corresponding to zero singular values, but also to look for slow modes, corresponding to small but finite singular values.

As expected, we find an extensive number of conserved quantities at $\Delta=0,1$ and only 2 exactly conserved quantities, corresponding to the Hamiltonian $H$ and $\left(S^{z}_{\mathrm{tot}}\right)^{2}$, for all other values of the anisotropy $\Delta$. This intermediate regime, however, exhibits a large number of slow modes, which will be discussed in the next section.

Since we are now working with a classical system, a few important distinctions ought to be made from our earlier, quantum analysis. First, we consider a slightly enlarged collection of bilinear operators:

[TABLE]

where $c_{1}=\sqrt{15}/2$ and $c_{2}=\sqrt{5}/2$. As before, $a=(i,j,\alpha)$ is a composite index. In the classical case, we also include bilinears with $i=j$ (which would be trivial in the spin-$1/2$ case). Note that there are only two independent such bilinears for each $i$, and the spherical harmonics (with spin $1$) provide an orthonormal basis. Indeed, it can be checked that the bilinears $O_{a}$ defined in Eq. 7.44 satisfy the orthonormality relation

[TABLE]

where $\left<[\dots]\right>_{\mathbf{S}}$ denotes an average over the infinite temperature ensemble, while the integral $\int D\mathbf{S}_{i}$ is over the unit sphere.

Second, while a single initial state is sufficient in the quantum SVD analysis, we have to consider an ensemble of initial states in the classical setting. This is because a single classical trajectory cannot visit the whole phase space due to energy conservation (a linear superposition of configurations does not exist classically). Here, we take the infinite-temperature ensemble, namely we sample $S_{0}=\{\mathbf{S}_{1}(0),\dots,\mathbf{S}_{N}(0)\}$ as independent random points on the unit sphere. We then time evolve with (7.36) for a total time $T$, and measure the expectation value of the bilinears $O_{a}(t_{n},S_{0})$ at discrete intervals $t_{n}=n\delta t\in[0,T]$. Repeating this for a large number $\mathcal{N}$ of initial conditions $\{S_{0}\}$ in the infinite-temperature ensemble, we construct the following matrix, analogous to the one in Eq. 5.24:

[TABLE]

where $\overline{O_{a}}(S_{0})=\int_{0}^{T}\frac{dt}{T}O_{a}(t,S_{0})$ represents the time average over one trajectory. The number of rows indexed by $a$ is, according to Eq. 7.44, $3N(N-1)/2+2N$. The columns are indexed by time $t\in[0,T]$ and initial condition $S_{0}$—in practice, we discretize the time axis (with the time-step $\delta t=1J^{-2}$) and draw a large number ($10^{3}$) of samples for $S_{0}$.

Then the singular value decomposition of $M$ is equivalent to diagonalizing the real Hermitian matrix $L\equiv MM^{\dagger}$, which can be obtained by averaging over the initial conditions $S_{0}$:

[TABLE]

Note that the average $\left<\dots\right>_{\mathbf{S}}$ is with respect to $S_{0}$ in the infinite-temperature thermal ensemble and should not be confused with the quantum expectation values $\left<\dots\right>$ (i.e. without a subscript) used in Sections V and VI. Diagonalizing $L$ allows us to obtain the slow mode operators $Q_{l}$, together with their corresponding eigenvalues $\sigma_{l}^{2}$. Similarly to (5.29), we have

[TABLE]

In other words, $\sigma_{l}^{2}$ is equal to the variance, averaged over initial conditions $S_{0}$, of the fluctuations of $Q_{l}$ along a given trajectory.

The behavior of the singular values $\sigma_{l}$ (shown in Fig. 7) is similar, in several ways, to that obtained in Section VI.1 for the quantum spin-$1/2$ model 333The main difference between the classical and quantum SVD results is that the classical singular values are considerably larger than the quantum ones (compare Fig. 2 and Fig. 7). This is because the expectation value of $O_{a}$ is evaluated on a single classical configuration in the former case, while the quantum expectation value is a result of a coherent average.. At the first integrable point $\Delta=0$, we obtain $N$ zero singular values corresponding to the family of $N$ spin-bilinear conserved quantities $\widetilde{G}^{(i)}$ from Eq. 4.19. At the second integrable point $\Delta=1$, we find $N+1$ zero singular values corresponding to the conserved quantities lying in the linear span of the $G^{(i)}$s from Eq. 4.17. Lastly, as shown in Fig. 7, away from these integrable points, i.e. for $\{\Delta<0\},\{0<\Delta<1\}$, and $\{\Delta>1\}$, we find two precisely zero singular values, corresponding to the two exactly conserved spin-bilinear quantities, $H$ and $\left(S_{\mathrm{tot}}^{z}\right)^{2}$. The small magnitude of the following singular values, for $l=3,4,\dots$, signals the presence of slow modes, which will be studied in the next section.

VII.3 Decay of slow operators

The SVD analysis of the previous section revealed a large number of operators with small singular values. In principle, we could characterize the thermalization (or lack thereof) of these operators $Q_{l}$ using, in analogy to the quantum case, a two-point correlation function

[TABLE]

where $\left<\dots\right>_{\mathbf{S}}$, as before, designates an average over the $\mathcal{N}$ initial conditions $S_{0}$. As in the quantum case, the operators $Q_{l}$ are orthonormal such that $G_{l}(0)=1$ (as $\mathcal{N}\to\infty$). Yet, the accurate computation of $G_{l}(t)$ at long times is typically very demanding because it requires averaging an increasingly complex function in phase space.

Fortunately, in classical systems, the singular value $\sigma_{l}$ already informs us about the long time plateau value of $G_{l}(t)$. This can be seen from Eq. 7.48, which implies that

[TABLE]

In the second line, we used the normalization $\left<Q_{l}(t)^{2}\right>_{\mathbf{S}}=1$; in the third line, we performed a change of variables $u=|t-s|$ (recall that $G_{l}(u)=\left<Q_{l}(t)Q_{l}(t\pm u)\right>_{\mathbf{S}}$ by the invariance of the infinite-temperature ensemble under time evolution). Now, it is not hard to show that $g_{l,T}$ and $G_{l}$ have the same infinite-time limit if that exists for $G_{l}$:

[TABLE]

Thus, $g_{l,T}$ is a finite-time proxy for $g_{l}$. In the infinite-time limit, the relation (7.50) becomes

[TABLE]

Using Eq. 7.50 or 7.51, this allows us to infer the plateau values of slow modes from the data of Fig. 7. Unsurprisingly, the exactly conserved quantities have $g_{l}=1$. Away from the integrable points at $\Delta=0,1$, we find that the slowest non-conserved modes, corresponding to $l=3,4,\dots,$ (the distinction between the slow modes and the rest is less sharp here than in the quantum case, and it is suggested by the rounded cusp around $l=N$ in Fig. 7 ), have a remarkably slow decay: the plateau values $g_{l,T}$ at a finite but large time $T=10^{4}J^{-2}$ are close unity 444We have checked that this persists up to $T=10^{6}J^{-2}$. Nonetheless, we find that these plateaus eventually decay for a finite-size system, albeit after a very long time. We leave the quantitative analysis for future work., comparable to their spin-$1/2$ counterparts.

In a future companion paper, we will demonstrate that any conserved operator in the $S=1/2$ quantum case is approximately conserved in the $S\rightarrow\infty$ classical model as well, up to $1/N$ corrections in large systems. Therefore, the classical model is expected to display some signatures of integrability. In this section, we saw that such signatures cannot be found from the Lyapunov exponent, but only from the relaxation of slow modes. This is intriguing, but a similar phenomenon has previously been observed. Ref. Goldfriend and Kurchan, 2019 showed that for certain systems near integrability (called “quasi-integrable” by the authors), the relaxation time of certain operators can be significantly longer than the finite Lyapunov time $1/\lambda_{L}$. Given that our classical system is surrounded by integrable lines $\Delta=0,1$ and (arguably) $S=1/2$ in the $(\Delta,S)$ parameter plane (see Fig. 1), we conjecture that it is also quasi-integrable. From this perspective, the existence of slow modes is compatible with the finite classical Lyapunov exponent found in Section VII.1.

VIII Experimental Realities

We have provided analytical and numerical evidence for the rich dynamical phase diagram depicted in Fig. 1, including clear signatures of chaotic dynamics at large $S$ and $\Delta\neq 0,1$, along with signatures of integrability at the special points $\Delta=0,1$ for any $S$. Further, we demonstrated signatures of a novel integrability* phase at $S=1/2$ for $\Delta\neq 0,1$. We now discuss prospects for observing these signatures in the laboratory. First, what should one measure to identify the chaotic and integrable regimes of the phase diagram? Second, given the inevitable presence of dissipation in realistic experiments, what are the requirements on cavity cooperativity to access the relevant time-scales experimentally?

To identify integrals of motion, the SVD method of Section V can equivalently be implemented with experimental data. Using state-sensitive imaging of the atomic ensemble Davis et al. (2019), one may immediately extract the bilinear spin correlation functions $\langle O_{a}(t)\rangle\propto\langle S_{i}^{\alpha}(t)S_{j}^{\alpha}(t)\rangle$ defined in Eq. 5.21. As each image is obtained from a destructive measurement, one must repeat the experiment many times to obtain statistics of the spin bilinears at a fixed time $t$, and then repeat this procedure for many time-points $t$ to obtain the full matrix $M_{a,t}$. With this matrix in hand, one can then directly apply the singular-value decomposition performed above in Section V.

A caveat is that measurements of the spin bilinears can be affected by dissipation due to photon loss and atomic free-space scattering. Photon loss from the cavity mode causes a random walk in the orientation of the weighted collective spin $\boldsymbol{\mathcal{F}}$ defined in Eq. 3.3. This effect is described by Lindblad operators

[TABLE]

where the decay rate (derived in Sec. S1) is given by

[TABLE]

The collective dissipation can be suppressed by increasing the detuning $\delta$ and compensating with increased drive strength, until limited by free-space scattering.

The effect of free-space scattering is to project or flip individual spins, as described by a set of Lindblad operators

[TABLE]

where $m$ or $m^{\prime}$ indicates the spin state of an individual atom indexed by $n$, and $C_{m,m^{\prime}}$ is an order-unity branching ratio. At large detuning, the scattering rate scales as

[TABLE]

Comparing Eqs. 8.53 and 8.55 shows that the cooperativity $\eta$ will dictate an optimal detuning for minimizing the net effect of the two forms of dissipation, with higher cooperativity enabling increasingly coherent dynamics.

To determine the cooperativity required to observe the signatures of integrability, we first write down explicit equations of motion for the spin bilinears $\langle S_{i}^{\alpha}(t)S_{j}^{\alpha}(t)\rangle$ evolving under the influence of pure collective dissipation or pure single-atom decay, respectively (see Supplementary Material Section S1). We find that the spin bilinears decay exponentially at a rate $\Gamma_{\mathrm{sc}}$ due to free-space scattering and at a rate $\gamma$ due to photon loss from the cavity. Notably, the rate of spin relaxation due to photon loss is not superradiantly enhanced, thanks to the counterbalanced effects of the $L_{\pm}$ Lindblad operators. Thus, at weak to moderate cooperativity $\eta\lesssim 1$ and large detuning $\delta>\kappa$, free-space scattering dominates and the bilinears decay on a time-scale $\tau J^{2}\sim\sqrt{N}\eta S\kappa/\delta$. For strong coupling $\eta\gg 1$, where free-space scattering is suppressed relative to cavity decay, the total dissipation can be minimized at a detuning $\delta\sim\sqrt{\eta}\kappa$, leading the spin bilinears to decay on a time-scale $\tau J^{2}\sim\sqrt{N\eta}S$.

To compare the decay time $\tau$ with the characteristic time-scales for observing the signatures of integrability, we refer to the time dependence of the autocorrelation functions $G_{l}(t)$ shown in Fig. 3. To observe the slow modes, a minimum requirement is to evolve the system for a time $t\gtrsim t^{*}\approx 10\ J^{-2}$, which governs the rapid decay of all non-integrable autocorrelation functions. This time can be reached even at $S=1/2$ in a strong-coupling cavity $\eta\sim 10$ with a system of $N=10^{3}$ sites, or with weaker single-atom cooperativity at larger $S$. To observe the plateaus themselves, we must evolve the system for a significantly longer time, at least $t\approx 10^{3}\ J^{-2}$ according to Fig. 3, which places a more stringent requirement $\sqrt{N}\eta S\gtrsim 10^{3}\delta/\kappa$. This regime is challenging to access for $S=1/2$ but readily accessible with large-$S$ subensembles, e.g., at $\eta\gtrsim 1$ with $N=10^{2}$ sites each consisting of $S=10^{3}$ spin-1 atoms.

Thus, current experiments are well positioned to explore the regime of mesoscopic spin $S$, in between the quantum ($S=1/2$) and classical ($S\rightarrow\infty$) limits. This will allow for testing the prediction that the plateaus in $G_{l}(t)$ calculated for spin $S=1/2$, indicating integrability across the full range $\{\Delta<0\}$,$\{0<\Delta<1\}$, and $\{\Delta>1\}$, persist for larger spin $S$ up to $1/N$ corrections (see Sec. VII.3). Experiments with scalable spin size $S$ may furthermore shed light on the transition from quantum integrability to chaos in the classical limit, as signified by the positive Lyapunov exponent in Fig. 6.

The chaotic dynamics observed in the classical limit $S\rightarrow\infty$ can be studied experimentally via the hallmark of sensitivity to perturbations. Recent theoretical and experimental work has shown that such sensitivity is accessible in quantum systems by measuring out-of-time-order correlators (OTOCs) Larkin and Ovchinnikov (1969); Shenker and Stanford (2014); Maldacena et al. (2016); Hosur et al. (2016); Swingle et al. (2016); Gärttner et al. (2017); Vermersch et al. (2018); Lewis-Swan et al. (2019), which quantify the spread of operators in time via the commutator $C(t)=\left\langle[V(t),W(0)]^{2}\right\rangle$. The connection to classical chaos is made clear in the semi-classical limit: for operators $V=S_{i}^{z}$, $W=S_{j}^{z}$, one can show that, to lowest order in a $1/S$ expansion, $C(t)\propto(\partial S_{i}^{z}(t)/\partial\phi_{j})^{2}$ for a small rotation $\phi_{j}$ at site $j$ about the $z$-axis Cotler et al. (2018). Thus, semi-classically the OTOC $C(t)$ measures the sensitivity of the coordinate $S_{i}^{z}(t)$ to changes in initial conditions $S_{j}^{z}(0)$, and may therefore be regarded as a quantum generalization of the classical sensitivity $d(t)$ defined in Section VII.1.

One way to access out-of-time-order correlators experimentally is to “reverse the flow of time” by dynamically changing the sign of the Hamiltonian Swingle et al. (2016); Gärttner et al. (2017). In the cavity-QED system considered here Davis et al. (2019), this sign reversal is achieved by switching the sign of the laser detuning $\delta$ in Eq. 3.9. The resilience of such time-reversal protocols to experimental imperfections, including dissipation, has been analyzed theoretically in Ref. Swingle and Yunger Halpern, 2018.

To allow for probing chaos in the cavity-QED system proposed here, the rates of collective dephasing and of decoherence via single-atom decay must be small compared with the Lyapunov exponent. We thus require $\lambda\tau\gg 1$, where $\tau$ is the characteristic decay time defined above. More specifically, given the Lyapunov exponents $\lambda\leq 0.08J^{2}$ shown in Fig. 6, and the requirement of observing the system for several Lyapunov “decades” to clearly identify exponential growth (Fig. 6a), we would like to evolve the system for times $t\gtrsim 100\ J^{-2}$, which are readily accessible in the large-$S$ regime that is of interest for approaching the classical limit.

Even in this regime, the light leaking from the cavity produces a continuous weak measurement of the collective spin $\boldsymbol{\mathcal{F}}$ whose quantum back-action may have consequences for the dynamics. The interplay of measurement back-action with chaos in open quantum systems, while beyond the scope of the present work, is a subject of active inquiry Eastman et al. (2017); Xu et al. (2019) and of fundamental importance for elucidating the quantum-to-classical transition Habib et al. (1998). The proposed experimental scheme, including the possibility of tuning the strength and form of coupling to the environment, opens new prospects for exploring this interplay.

IX Discussion

We have studied a class of spin models with separable, all-to-all, random interactions and found a complex dynamical phase structure that depends on the spin size $S$ and the anisotropy $\Delta$ along the $z$-axis. We showed that our model at $\Delta=1$ is equivalent to the well-studied rational Gaudin model, and exhibits special integrable dynamics for all values of $S$. We also proved and confirmed numerically that there exists another special point at $\Delta=0$ where the model is also integrable (in the same sense), regardless of the spin size. Surprisingly, we found compelling numerical evidence that the system at $S=1/2$ is integrable for any anisotropy $\Delta\notin\{0,1\}$. In contrast to the special points $\Delta=0,1$, the integrals of motion at other values of $\Delta$ are not purely spin bilinears and develop tails on $2n$-body terms. We leave the detailed characterization of these dressing tails to future work. Lastly, we found that integrability away from $\Delta=0,1$ is a purely quantum phenomenon: by numerically solving the Hamilton equations of motion for the classical model ($S\rightarrow\infty$), we showed that its dynamics is chaotic with a non-zero Lyapunov exponent and that there exist only two exactly conserved quantities, as opposed to the extensive family of conservation laws characterizing a classically integrable system. However, even in the classical regime we find an extensive number of quasi-conserved charges, whose decay time appears to diverge in the large-$N$ limit. A more thorough study of this regime will be given in future work.

Our analysis opens up several further lines of inquiry. First, since the Hamiltonian (1.2) at the special point $\Delta=1$ (and, presumably, at $\Delta=0$ as well Balantekin et al. (2005); Skrypnyk (2005, 2009a, 2009b, 2009c)) possesses a quantum group structure, does the integrable∗ phase exhibit any algebraic structure? Is it possible to construct explicitly the dressed conserved quantities in terms of the model couplings?

Second, we have seen that even though the level statistics of the spin-1/2 system deviates from Wigner-Dyson statistics, exhibiting many level crossings (this holds also for spin-1, as shown in the Supplementary Material), its classical counterpart is chaotic with a finite Lyapunov exponent. We note that this does not contradict the Berry-Tabor conjecture Berry and Tabor (1977); Bohigas et al. (1984), which applies to the semiclassical, large-$S$, regime. In fact, the same phenomenon is known to occur in integrable quantum spin chains, such as the anisotropic Heisenberg model (or XXZ chain): its Hamiltonian $\sum_{j}\left[S^{x}_{j}S^{x}_{j+1}+S^{y}_{j}S^{y}_{j+1}+\Delta S^{z}_{j}S^{z}_{j+1}\right]$ is quantum integrable only for spin-1/2, and it is classically chaotic; its integrable higher-spin extensions have different Hamiltonians and are nontrivial to obtain Bytsko (2003). We wonder whether our integrable∗ phase admits any such extensions, which might shed light on the quasi-integrability of our classical model.

Third, we have only characterized the boundaries of the phase diagram in Fig. 1. A straightforward and interesting next step would be to study the quantum-to-classical crossover by better understanding how classical chaos (and perhaps quasi-integrability) at $S\rightarrow\infty$ emerges from the integrable∗ regime at $S=1/2$.

In fact, this putative transition between (quantum) integrability and (semiclassical) chaos may also be probed experimentally. The model (1.2) can be implemented in a near-term experiment using atomic ensembles confined in a single-mode optical cavity. This would allow for a systematic exploration of the rich physics contained in the dynamical phase diagram (Fig. 1). By changing the local atom density to increase the number of atoms in a given region of constant coupling to the cavity mode, the spin size $S$ can be varied from $S=1/2$ all the way to a semiclassical regime $S\gg 1$: this would enable the experiment to tune between quantum and classical dynamics. Meanwhile, changing the angle between the magnetic field defining the spins’ $z$-axis and the axis of the optical cavity allows for tuning of the anisotropy $\Delta$, so that both the special points $\Delta=0,1$ and the regions $\{\Delta<0\}$, $\{0<\Delta<1\}$, and $\{\Delta>1\}$ can be investigated.

Last, we emphasize that the SVD technique described in Section V can be applied directly to the experimental data, revealing the conserved quantities and slow modes. More broadly, we envision using this approach in studying a wider class of physical systems wherein the integrals of motion or their number are not a priori known.

In summary, the model (1.2) and its associated experimental setup represent a novel paradigmatic platform for studying integrability, chaos, and thermalization under closed many-body quantum dynamics.

Acknowledgements

We would like to thank Romain Vasseur, Fabian Essler, Thomas Klein Kvorning, Daniel E. Parker, Emily Davis, and Avikar Periwal for fruitful conversations. This work was supported by the DOE Office of Science, Office of High Energy Physics (GB), the grant DE-SC0019380 (XC, XLQ, MSS, and EA), the ERC synergy grant UQUAM (IDP, XC, and EA), and the Emergent Phenomena in Quantum Systems initiative of the Gordon and Betty Moore Foundation (TS). The numerical computations were carried out on the Lawrencium cluster resource provided by the IT Division at the Lawrence Berkeley National Laboratory under the DOE contract DE-AC02-05CH11231, on the Sherlock computing cluster provided by Stanford University and the Stanford Research Computing Center, and on the cluster of the Laboratoire de Physique Théorique et Modèles Statistiques (CNRS, Université Paris-Sud).

S1 Experimental Realization

S1.1 Derivation of the Effective Hamiltonian

Here we elaborate on the derivation of the effective spin Hamiltonian (3.8) from the atom-light interaction Hamiltonian (3.7), which we now repeat for completeness:

[TABLE]

To simplify the following derivation, we will assume that the weights $\xi_{i}$ are real numbers (although it is interesting to speculate whether one can access an even richer set of dynamics if the weights are allowed to have both non-uniform phases and amplitudes). In this case, we may write the full Hamiltonian as:

[TABLE]

where we have passed into a rotating frame with respect to the atomic Zeeman splitting $\omega_{Z}$.

Provided the occupation of the $v$ mode remains small, we can adiabatically eliminate it from the dynamics following the approach of Reiter and Sørenson Reiter and Sørensen (2012). This procedure is essentially a perturbation theory calculation that considers 2-photon scattering processes in which a virtual photon is scattered into the $v$ mode and reabsorbed by the atomic ensemble. Inspecting the Hamiltonian (S2), we find three distinct processes that add one photon to the $v$ mode:

[TABLE]

In addition, the Hermitian conjugates of these terms remove a photon from the mode $v$. We must consider all pairs of processes that add a photon to the mode and subsequently reabsorb it. However, only the two-photon processes that are resonant will dominate the slow, effective, ground-state dynamics. For instance, the two-photon process proportional to $vv^{\dagger}\mathcal{F}_{-}\mathcal{F}_{-}$ is an off-resonant process and is, therefore, accompanied by a rapidly rotating phase factor $e^{-2i\omega_{Z}t}$. As a result, this term quickly averages to zero on timescales $t\gg 1/\omega_{Z}$, and we are justified in ignoring it. In fact, the only resonant 2-photon processes that survive are the terms proportional to $\mathcal{F}_{z}\mathcal{F}_{z},\mathcal{F}_{-}\mathcal{F}_{+},$ and $\mathcal{F}_{+}\mathcal{F}_{-}$, since all other terms have rapidly oscillating phase factors. The result of this elimination scheme is an effective Hamiltonian for the spins:

[TABLE]

where $\delta_{\pm}=\delta\mp\omega_{Z}$ are the detunings from the two-photon resonance, $\kappa$ is the cavity linewidth, and $\zeta(\delta)=\delta/[\delta^{2}+(\kappa/2)^{2}]$. At large detuning $\delta\gg\kappa,\omega_{Z}$, Eq. S3 simplifies to Eq. 3.9, and we obtain the desired model with an anisotropy parameter $\Delta$ controlled by the angle $\theta$ of the magnetic field.

The effective Hamiltonian (3.9), however, is obtained only if the $\mathcal{F}_{+}\mathcal{F}_{-}$ and $\mathcal{F}_{-}\mathcal{F}_{+}$ terms in Eq. S3 are balanced, which occurs perfectly only in the limit $\delta\rightarrow\infty$. More generally, at finite $\delta$, the cross-terms $\mathcal{F}_{x}\mathcal{F}_{y}$ do not cancel and we obtain an effective Hamiltonian

[TABLE]

which is of the same form as Eq. 3.9, but with additional non-uniform longitudinal magnetic field terms. Such on-site terms, however, are sub-extensive relative to the spin-spin interaction terms. Moreover, we show numerically in Section S2.6 that these additional terms do not affect the integrability of the model at $S=1/2$ and finite $N$.

S1.2 Effects of Dissipation

In addition to the coherent dynamics, the driven cavity system suffers from dissipation due to photon loss from the $v$ mode and atomic free-space scattering from the excited states. These processes can be described formally in a quantum master equation by the relaxation operators

[TABLE]

and

[TABLE]

respectively. Here, $(m,m^{\prime})$ label internal states of the individual atoms indexed by $n=1,\ldots,2SN$. $C_{m,m^{\prime}}$ are branching ratios that sum to unity. $\Gamma_{\mathrm{sc},n}=(\Omega_{n}/\mathcal{D})^{2}\Gamma$ is the free-space scattering rate for the $n$th atom, where $\mathcal{D}$ is the detuning from atomic resonance and $\Gamma$ is the natural excited-state linewidth. In comparison to the scheme presented in Section III, note that the Rabi frequency $\Omega_{n}$ and detuning $\mathcal{D}$ enter likewise into the interaction strength in the Hamiltonian (3.9), where $(\chi\xi_{n})^{2}\sim(\Omega_{n}/\mathcal{D})^{2}g_{n}^{2}$. Thus, the scattering rate scales as $\Gamma_{\mathrm{sc},n}\sim(\chi\xi_{n})^{2}/(\eta\kappa)$.

The cavity dissipation described by $L_{v}$ generates collective dephasing of the atomic ensemble. Upon adiabatically eliminating the cavity mode using the Reiter-Sørenson procedure, we find that cavity dissipation leads to an effective relaxation operator

[TABLE]

Due to the rapidly oscillating phase factors $e^{\pm i\omega_{Z}t}$, this relaxation operator may be split into three effective relaxation operators

[TABLE]

where we have introduced the collective decay rate

[TABLE]

By contrast, the free-space scattering operators $L_{n,(m,m^{\prime})}$ act directly on the spins so there is no need to apply the adiabatic elimination procedure.

To what extent do these dissipative processes spoil the signatures of the integrability studied in Section VI? In the strongly quantum regime at small $S$, for instance, one might be interested in being able to evolve the system long enough to see the plateaus in the autocorrelation function $G_{l}(t)$, as observed in Fig. 5. To do so, we must ensure that the dissipation timescale is short compared to the timescale $\tau_{l}$ required for the plateaus to appear. While the full dynamics including both coherent and dissipative processes is difficult to access numerically, we can estimate the dissipative timescales by solving for the dynamics in the presence of dissipation alone. In the absence of coherent dynamics (i.e. $H=0$), one can compute the equations of motion for the expectation values of operators using the Lindblad equation:

[TABLE]

where $r$ runs over the set of relaxation operators. Plugging in the free-space scattering operators $L_{r}=L_{n,(m,m^{\prime})}$, and assuming an approximately uniform scattering rate $\Gamma_{\mathrm{sc},n}=\Gamma_{\mathrm{sc}}$, one can show that bilinears of spin-$S$ operators $S_{j}^{\alpha}$ obey the following equations of motion for $j\neq k$:

[TABLE]

where we have introduced the notation

[TABLE]

Thus, the bilinear operators decay on a timescale $\tau\sim 1/\Gamma_{\mathrm{sc}}$ due to free-space scattering.

We can perform a similar calculation for the collective relaxation operators (S8) generated by photon loss from the cavity mode. For large detuning $\delta_{\pm}\approx\delta\gg\kappa,\omega_{z}$, one can show that the spin bilinears obey the following equations of motion for $j\neq k$:

[TABLE]

Assuming that the weights $\xi_{i}$ are numbers of order one, the bilinear operators therefore decay on a timescale $\tau\sim 1/\gamma$. A key result of this analysis is that the collective dissipation rate is, somewhat surprisingly, independent of the atom number $\sim NS$, exhibiting no superradiant enhancement. The reason for this is that the average effects of the $L_{\pm}$ Lindblad operators cancel in the large-detuning limit, where $\delta_{+}\approx\delta_{-}$. In practice, at finite detuning, any imbalance in these terms can be corrected by adding a weak auxiliary drive field, with no other significant effect on the dynamics.

The relative strengths of decay via the cavity ($\gamma$) and via spontaneous emission ($\Gamma_{\mathrm{sc}}$) are controlled by the detuning $\delta$. The total decay rate $\gamma+\Gamma_{\mathrm{sc}}$, at fixed interaction strength $J^{2}$, is minimized by choosing a detuning $\delta=\sqrt{1+\eta}\kappa$. For weak to moderate cooperativity $\eta\lesssim 1$, we may instead choose a larger detuning satisfying the conditions $\delta\gg\kappa,\omega_{Z}$ assumed above.

S2 Quantum model and exact diagonalization results

S2.1 Many-body bandwidth

We now present analytical and numerical evidence that the choice of normalization for the Hamiltonian in Eq. 1.2 leads to a sensible thermodynamic limit. On the one hand, we show that the many-body bandwidth $W$ scales superlinearly, to wit: $W\sim\mathcal{O}(N^{3/2})$. An important caveat is that $W$ per se does not mean much unless it is weighed by the many-body density of states $\rho(E)$. On the other hand, the energy fluctuations $\Delta E$ (at least at high temperatures) scale as $\Delta E\sim\mathcal{O}(\sqrt{N})$ which leads to an extensive specific heat, as one would normally expect.

Let us start by rewriting the Hamiltonian as

[TABLE]

where $\boldsymbol{\mathcal{F}}\equiv\sum_{i}J_{i}\mathbf{S}_{i}$. Each operator $\mathcal{F}_{\alpha}^{2}$ has non-negative eigenvalues and we denote the largest one by $f_{0}^{2}$. Note that it is the same for all $\alpha\in\{x,y,z\}$ by isotropy. Thus, the largest eigenenergy $E_{M}$ of $H$ is bounded by

[TABLE]

We can also bound $E_{M}$ below by using any properly normalized variational state $\left|\psi\right\rangle$ since $E_{\psi}=\left\langle\psi\right|H\left|\psi\right\rangle=\sum_{n}E_{n}|\langle\psi|n\rangle|^{2}\leq E_{M}\sum_{n}|\langle\psi|n\rangle|^{2}=E_{M}$. Let us choose

[TABLE]

as a variational state. Then $\left\langle\psi\right|H\left|\psi\right\rangle=\frac{\Delta}{S\sqrt{N}}\left\langle\psi\right|\mathcal{F}_{z}^{2}\left|\psi\right\rangle$ and

[TABLE]

Moreover, this variational state is an eigenstate of $\mathcal{F}_{z}$ and has the largest eigenvalue $f_{0}=S\sum_{i}|J_{i}|$. Thus, we find that

[TABLE]

Since the fields $J_{i}$ are sampled from $\mathcal{N}(0,1)$, the random variable $|J_{i}|$ has a half-normal distribution with mean $\mu=\sqrt{2/\pi}$ and variance $\sigma^{2}=\left(1-\frac{2}{\pi}\right)$. Then $Z\equiv\sum_{i}|J_{i}|$ has mean $N\mu$ and variance $N\sigma^{2}$ so $\left<Z^{2}\right>_{J}=N^{2}\mu^{2}+N\sigma^{2}$. Thus, at large $N$, $\left<\left(\sum_{i}|J_{i}|\right)^{2}\right>_{J}\approx 2N^{2}/\pi$. After disorder averaging Eq. S17 we get

[TABLE]

which shows that the disorder averaged $\left<E_{M}\right>_{J}\sim\mathcal{O}(N^{3/2})$.

We now put bounds on the ground state energy $E_{m}$. Clearly, $E_{m}\geq 0$. More importantly, we note that

[TABLE]

is an exact eigenstate of $H$ with energy $E_{\phi}=\frac{\Delta S}{\sqrt{N}}\left(\sum_{i}J_{i}\right)^{2}$. Since $\left<\left(\sum_{i}J_{i}\right)^{2}\right>_{J}=N$, we find that

[TABLE]

which shows that $\left<{E_{m}}\right>_{J}\sim\mathcal{O}(N^{1/2})$. Since $\left<E_{M}\right>_{J}\sim\mathcal{O}(N^{3/2})$, this concludes our proof that the disorder-averaged bandwidth scales as $\left<{W}\right>_{J}\sim N^{3/2}$. We have also computed $\left<{W}\right>_{J}$ numerically and we find an excellent agreement with the above analysis, as shown in Fig. S1.

Finally, we note that although $W$ scales superlinearly, it is quite possible that the bulk of the many-body states are located around a typical energy $E_{\mathrm{typ}}\sim\mathcal{O}(N)$ and that the density of states has a long and thin tail extending all the way to $\mathcal{O}(N^{3/2})$.

S2.2 Energy fluctuations

We now show that the disorder-averaged thermal energy fluctuations ${(\Delta E)^{2}}$ at high temperatures scale as $\left<(\Delta E)^{2}\right>_{J}=cN$ for some constant $c$. Note that, by the fluctuation-dissipation theorem, this means that the specific heat per particle behaves as $\frac{C_{v}}{N}\approx\frac{c}{T^{2}}$ at high temperatures $T$.

The energy fluctuations are defined as

[TABLE]

where $\langle\dots\rangle_{T}\equiv\left(\operatorname{Tr}e^{-\beta H}\right)^{-1}\operatorname{Tr}\left[e^{-\beta H}\dots\right]$ denotes thermal averaging. For simplicity, we will focus on the $S=1/2$ and $\Delta=1$ case. At lowest order in $\beta$ we can approximate

[TABLE]

where $Z\sim\chi^{2}(N)$. Similarly, the first term becomes

[TABLE]

Putting the two results together we find that $\left(\Delta E\right)^{2}\approx\frac{3}{2N}Z^{2}$. Since $Z$ is chi-squared distributed $\chi^{2}(N)$, we know that $\mathbb{E}\left[Z^{2}\right]=N^{2}+2N$. Thus, after averaging over disorder, we get

[TABLE]

This ensures that the specific heat $C_{v}=\frac{3}{2T^{2}}N$ is extensive. We can also compute this quantity numerically using the exact spectrum of $H$ and we find an excellent agreement, as shown in Fig. S2.

S2.3 Overlaps with higher powers of the conserved quantities $\hat{H}$ and $\hat{S}_{\mathrm{tot}}^{z}$

We now address the possibility that the large plateau values (Fig. 4) observed in the quantum model at intermediate $\Delta\in(0,1)$ are caused by a significant overlap between the slow operators $\hat{Q}_{l}$ and higher powers of the conserved quantities, $\hat{H}$ and $\hat{S}_{\mathrm{tot}}^{z}$. We use a twofold approach, as detailed below.

Let us first summarize why this is a legitimate concern. Recall that at $\Delta\in(0,1)$, the SVD analysis finds two exactly conserved spin bilinears, $\hat{H}$ and $(\hat{S}_{\mathrm{tot}}^{z})^{2}$, along with another $N-2$ slow modes $Q_{l}$, where $l=3,\dots,N$. The auto-correlation function of the latter operators exhibits a large long-time plateau, which we interpreted as a consequence of the existence of further conservation laws (6.31) and of quantum integrability. However, it is possible the the $\hat{Q}_{l}$s overlap non-trivially with exactly conserved operators such as $\hat{H}^{p}(\hat{S}_{\mathrm{tot}}^{z})^{p^{\prime}}$: these higher powers of the exactly conserved quantities consist of $k$-body terms, including $k=2$, so they have non-zero support in the space of spin bilinears and, thus, can conceivably overlap with $\hat{Q}_{l}$.

Furthermore, the support in the space of $2$-local operators can be non-zero even in the thermodynamic limit $N\rightarrow\infty$ due to the all-to-all nature of the interaction. For instance, consider $\hat{\mathcal{O}}_{1}=\frac{1}{\sqrt{N}}\hat{H}$. By the previous section, we know that $\frac{1}{\operatorname{Tr}\left[\mathbf{1}\right]}\operatorname{Tr}\left[\hat{\mathcal{O}}_{1}^{2}\right]$ is of order one even in the thermodynamic limit. Similarly, let $\hat{\mathcal{O}}_{2}=\frac{1}{N}\left(\hat{S}_{\mathrm{tot}}^{z}\right)^{2}$ such that $\frac{1}{\operatorname{Tr}\left[\mathbf{1}\right]}\operatorname{Tr}\left[\hat{\mathcal{O}}_{2}^{2}\right]$ is also of order one as $N\rightarrow\infty$. It is not hard to check that the only $2$-body terms emerging from higher powers of the conserved quantities, namely operators of the form $\hat{\mathcal{O}}_{1}^{p}\hat{\mathcal{O}}_{2}^{p^{\prime}}$, that do not vanish in the thermodynamic limit are:

[TABLE]

Above, $\hat{\mathcal{O}}_{3}$ can be generated by $\hat{\mathcal{O}}_{1}\hat{\mathcal{O}}_{1}$ as long as $\Delta\in(0,1)$, and $\hat{\mathcal{O}}_{4}$ by $\hat{\mathcal{O}}_{1}\hat{\mathcal{O}}_{2}$ (in fact, one gets a prefactor $\frac{1}{N}\sum_{k}J_{k}$; it scales as $1/\sqrt{N}$ typically if $J_{k}$ has zero mean and as $O(1)$ if $\overline{J_{k}}\neq 0$; by precaution, we still include it in the following analysis). We perform a Gram-Schmidt orthonormalization of these four operators to obtain a basis set $\left\{\widetilde{\mathcal{O}}_{1},\widetilde{\mathcal{O}}_{2},\widetilde{\mathcal{O}}_{3},\widetilde{\mathcal{O}}_{4}\right\}$ wherein $\frac{1}{\operatorname{Tr}\left[\mathbf{1}\right]}\operatorname{Tr}\left[\widetilde{\mathcal{O}}_{i}\widetilde{\mathcal{O}}_{j}\right]=\delta_{ij}$. We then define a matrix $R$ whose columns vectors are the $\widetilde{\mathcal{O}}_{i}$s. This allows us to perform the SVD analysis in the operator space orthogonal to the one spanned by the $\widetilde{\mathcal{O}}_{i}$s: using the time series matrix $M$ from Eq. 5.24, we perform the SVD decomposition of $\left(\mathbf{1}-RR^{\dagger}\right)M$ instead. Trivially, we find that the first four singular values are exactly zero so we focus on the modes corresponding to $l>4$. Following the prescription from Section VI.2, we define the slow operators $\widetilde{Q}_{l}$ for $l\geq 5$ and compute their auto-correlation function $\widetilde{G}_{l}(t)=\frac{1}{\operatorname{Tr}\left[\mathbf{1}\right]}\operatorname{Tr}\left[\widetilde{Q}_{l}(t)\widetilde{Q}_{l}(0)\right]$, as shown in Fig. S3a. Thus, we see that the overlap with higher powers of $\hat{H}$ and $\hat{S}_{\mathrm{tot}}^{z}$ has a negligible effect on the plateau values of the slow modes $\hat{Q}_{l}$.

In the second approach, we numerically study the effect of $\hat{H}^{p}$ on the plateau values by analyzing the matrix elements $Q_{l,nn}\equiv\left\langle n\right|\hat{Q}_{l}\left|n\right\rangle$, where $\left|n\right\rangle$ is an energy eigenstate, $H\left|n\right\rangle=E_{n}\left|n\right\rangle$. Let us define the vector $\vec{q}_{l}$ whose elements are $q_{l}^{(n)}=\frac{1}{\sqrt{\operatorname{Tr}\left[\mathbf{1}\right]}}\left|Q_{l,nn}\right|$. It is known that the norm of this vector is just the plateau value of the $Q_{l}$ mode, whenever the latter is well-defined: $g_{l}=\left\|\vec{q}_{l}\right\|^{2}$. Another way of studying the effect of $\hat{H}^{p}$ on the plateau value is to analyze the overlap between $\vec{q}_{l}$ and the set of vectors $\left\{\vec{E}^{(p)}:E_{n}^{(p)}=\left(E_{n}\right)^{p},0\leq p\leq P\right\}$, i.e. the powers of the energy eigenvalues.

Once again, we perform a Gram-Schmidt orthonormalization of the above vectors to obtain a set $\left\{\vec{\eta}_{i}:0\leq i\leq P\right\}$ wherein $\langle\vec{\eta}_{i}|\vec{\eta}_{j}\rangle\equiv\vec{\eta}_{i}^{T}\cdot\vec{\eta}_{j}=\delta_{ij}$. We then define

[TABLE]

which is a measure of the residual plateau value for the mode $\hat{Q}_{l}$ after removing the contributions of $\hat{H}^{0},\hat{H}^{1},\dots,\hat{H}^{p}$. Naturally, for any fixed mode index $l$, $\left\|Q_{l,p}\right\|^{2}$ decreases from the plateau value $\left\|Q_{l,0}\right\|^{2}=\left\|\vec{q}_{l}\right\|^{2}=g_{l}$ as $p$ increases. Numerically, in Fig. S3b (main plot), we observe that for the slow modes, $2<l<N$, the decrease is small compared to $g_{l}$ for powers up to $p\sim N/2$, at which point $\hat{H}^{p}$ becomes an $N$-body operator. In fact, $\left\|Q_{l,p}\right\|^{2}$ depends roughly linearly on $p$ and decays approximately to zero only when $p\sim 2^{N}\gg N$. We have also checked that the result is not altered upon including low powers of $\hat{S}^{z}_{\mathrm{tot}}$ (i.e. removing contributions from $\hat{H}^{p}(\hat{S}^{z}_{\mathrm{tot}})^{r}$, for $p,r\leq N$). Therefore, the plateaus of the slow modes are not caused by simple functions of the exactly conserved quantities, $\hat{H}$ and $\hat{S}^{z}_{\mathrm{tot}}$, at $\Delta\in(0,1)$. Finally, the inset of Fig. S3b shows that slow modes $\hat{Q}_{l}$ are not a linear combination of a few projectors onto energy eigenstates.

To summarize, these analyses render highly unlikely the possibility that the slow mode plateaus at intermediate $\Delta$ are solely due to the overlap with the known conserved quantities. This reinforces the thesis of the existence of further nontrivial conserved quantities (6.31), as put forth in the main text.

S2.4 The middle of the SVD spectrum

We now check that, at $\Delta=1.0$, we do not find more than $N+1$ bilinear conserved quantities, as suggested by the analytic results of Section IV. More precisely, we numerically compute the plateau values corresponding to the $l=N+2$, $l=3N$, and $l=N^{2}$ modes obtained from the quantum SVD analysis. As shown in Fig. S4, we see that the latter two plateau values decrease with the system size $N$. However, this is not entirely clear for the $N+2$ mode as its flow with system size is inconclusive—we conjecture that this is a finite size effect rather than a sign of an additional conserved quantity.

S2.5 Level statistics

We also analyze the level statistics of the quantum Hamiltonian (1.2) by studying the ratio of adjacent energy levels $E_{n}$ defined by $r_{n}=\min\{\Delta E_{n},\Delta E_{n+1}\}/\max\{\Delta E_{n},\Delta E_{n+1}\}$, where $\Delta E_{n}=E_{n+1}-E_{n}$. We now provide details vis-à-vis the various symmetries and subsequent degeneracies we have to take into account.

At $\Delta=1$ the Hamiltonian has the full $\mathrm{SU}(2)$ symmetry and the eigenstates of $H$ are grouped in blocks corresponding to irreducible representations (irreps). Each state is simultaneously an eigenstate of $\hat{S}_{\mathrm{tot}}^{2}$ with an eigenvalue $S_{\mathrm{tot}}(S_{\mathrm{tot}}+1)$ and of $\hat{S}_{\mathrm{tot}}^{z}$ with an eingenvalue $S_{\mathrm{tot}}^{z}\in\{-S_{\mathrm{tot}},\dots,S_{\mathrm{tot}}\}$. Thus, each block contains $2S_{\mathrm{tot}}+1$ eigenstates of $H$. We keep the states in a fixed sector of both $\hat{S}_{\mathrm{tot}}^{2}$ and $\hat{S}_{\mathrm{tot}}^{z}$ by selecting the block $S_{\mathrm{tot}}$ that has the largest multiplicity in the irreps. decomposition (i.e. the most frequent block to obtain the largest number of viable states).

At $\Delta\neq 1$ the symmetry is $\mathrm{U}(1)\rtimes\mathcal{T}$: $\mathrm{U}(1)$ is generated by $U_{\phi}=\exp\left(-i\phi S_{\mathrm{tot}}^{z}\right)$ and it corresponds to the conservation of $S_{\mathrm{tot}}^{z}$, whereas $\mathcal{T}$ represents the anti-unitary time reversal symmetry. Moreover, there is a mirror symmetry corresponding to $\hat{S}_{z}\rightarrow-\hat{S}_{z}$. For the cases wherein $\mathcal{T}^{2}=-1$, we have Kramers doublets and we keep the states solely in the $S_{\mathrm{tot}}^{z}=1/2$ sector because it is the most populous. For the cases where $\mathcal{T}^{2}=+1$, we keep states solely in the $S_{\mathrm{tot}}^{z}=1$ sector in order to account for the mirror symmetry as well (the states in the $S_{\mathrm{tot}}^{z}=0$ sector are degenerate under the mirror symmetry).

Having selected the eigenstates of $H$ to account for the degeneracies due to symmetries, we compute the ratio $r_{n}$ for the adjacent levels. We then average over all $r_{n}$ and over disorder realizations to obtain $\left<r\right>_{J}$. We know that a chaotic system exhibits level repulsion and has GOE level statistics characterized by a value $\langle r\rangle_{\mathrm{GOE}}\approx 0.53$. An integrable system does not have level repulsion (its energy levels are uncorrelated) and generically has Poissonian level statistics characterized by a value $\langle r\rangle_{\mathrm{Poisson}}\approx 0.387$.

In Fig. S5a we plot $\left<r\right>_{J}(N)$ for different values of the anisotropy $\Delta$. We find that both $\Delta=0$ and $\Delta=1$ exhibit almost perfectly Poissonian level statistics. However, for intermediate $0<\Delta<1$, the system of $S=1/2$ spins (or $S=1$ for the inset) is not as close to Poisson, although it is even farther away from GOE. Moreover, we cannot extract a meaningful flow with the system size $N$ towards either GOE or Poisson. We note, in passing, that there are exceptions from the equivalence between integrability and Poissonian level statistics: as detailed in Ref. Caux and Mossel, 2011 (see also the references therein), there are well-known counterexamples Barba et al. (2008); Berry and Tabor (1977); Schliemann (2010); Scaramazza et al. (2016) of integrable systems whose level statistics are neither Poisson nor GOE.

Nonetheless, if we add the perturbation $H_{\epsilon}^{(1)}$ from Eq. 6.32 to $H$ we immediately see that the system obeys almost perfectly GOE statistics even for rather small system sizes (see the curves with square markers in Fig. S5a).

In Fig. S5b we work in a fixed disorder realization for the fields $\{J_{i}\}$ and we adiabatically change the anisotropy $\Delta$. Then we plot the exact energy levels within a fixed symmetry sector, as detailed above, and explicitly check whether there exist level crossings or if they are avoided (level repulsion). We find that there are numerous real crossings at intermediate $\Delta$, which provides further qualitative evidence that $H$ is integrable for $S=1/2$.

S2.6 Stability to the commutator perturbation

We conclude by studying the effect of the perturbation (6.34) that arises upon recasting the Hamiltonian into the “flip-flop” form of Eq. 6.35. The term $J_{i}J_{j}S_{i}^{+}S_{j}^{-}$ from (6.35) in the main text can be written as

[TABLE]

which leads an overall term $\frac{1}{S\sqrt{N}}\sum_{i}J_{i}^{2}S_{i}^{z}$ in addition to the Hamiltonian $H$ defined in Eq. 1.2. We have argued that this is an irrelevant perturbation in a thermodynamic system with $N\gg 1$ degrees of freedom.

We now explicitly study its impact on the plateau values of the auto-correlation function $G_{l}(t)$. As shown in Fig. S6, we find that the effect of this perturbation is only quantitative: the modes reach plateau values that are lower than those occuring in the unperturbed model, but they do not decay as the system size $N$ increases. We note that the quantitative shift is also due to the fact that the slow modes we compute are bilinear operators even though the perturbation $\frac{1}{S\sqrt{N}}\sum_{i}J_{i}^{2}S_{i}^{z}$ is a linear combination of $1$-body terms.