Integrable and chaotic dynamics of spins coupled to an optical cavity

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' on 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.


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 1 . In a related development, Kitaev constructed a class of quantum many-body models whose dynamics saturates this bound on chaos 2,3 and can be related to black holes through the AdS/CFT correspondence [4][5][6] . 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 7 was originally proposed as a quantum spin model with random all-to-all couplings: where S i are 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. 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. However, there is a crucial difference, already pointed out in Ref. 8, 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 [8][9][10][11] , separable all-to-all couplings are realized in numerous existing experiments [12][13][14][15] 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 spinbilinear conserved quantities, and of a novel form of integrability-labeled Integrable * in Fig. 1-with quasibilinear integrals of motion.
The class of models we consider in this paper is described by the following quantum spin Hamiltonian: where S i are SU (2)  (bottom) Schematic of the atomic sub-ensembles (red) trapped inside a single-mode optical cavity (blue). A drive field (green) at a detuning δ from the cavity resonance generates effective spin-spin interactions between the atoms (bottom right). The tunable angle θ between the spin quantization axis (along the applied magnetic field B) and the cavity's longitudinal axis leads to an anisotropy ∆ = 2 cot 2 θ. By changing the local atomic density in a region of constant coupling to the cavity mode, the effective spin size S can also be varied, allowing for the systematic exploration of the full phase diagram.
ensembles located at sites i = 1, . . . , N . The sitedependent 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 Ω i of an inhomogeneous drive field. 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 A). 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 ∆, 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 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 flex-ibility 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 α 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 ∆, 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 ∆ = 0 and ∆ = 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 ∆ = 0, 1. Specifically, we present an exact diagonalization study of the spin-1/2 model, showing the robustness of the integrable structure on the entire line 0 < ∆ < 1, with quasibilinear integrals of motion. In Section VII we simulate the classical model (S → ∞) 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 ∆ = 1, for all spin sizes S, on which the Hamiltonian from Eq. 1.2 is equivalent to a rational Gaudin model 16 . This model is quantum integrable in the mathematical sense of possessing an underlying quantum group structure 17 . 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 18,19 . 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 ∆ = 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 ∆ = 1, this integrability holds for any value of S, including the classical limit S → ∞. The integrability of the model at ∆ = 0 is connected to the existence of a non-standard class of Gaudin models [20][21][22][23][24] .
The most surprising part of the phase diagram lies in the anisotropy range 0 < ∆ < 1. In this regime, we identify a novel integrable structure that is markedly different from the type of integrability found at the two end-points, ∆ = 0 and ∆ = 1.
First, unlike the latter, integrability for 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 λ L in the classical limit (S → ∞). Nevertheless, in this latter limit, we also find that there exist modes that relax only on time scales much larger than λ −1 L . We conjecture that this is a consequence of the "quasi-integrable" nature 25 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 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 multispin terms on top of the dominant two-spin contribution is reminiscent of the quasi-local integrals of motion that characterize Many-Body Localized phases [26][27][28][29][30] .

III. PROPOSED EXPERIMENTAL SCHEME
As advertised, the full phase diagram of Fig. 1a can be accessed in experiments with atomic ensembles in singlemode 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 where each weight ξ i is set by the amplitudes of the cavity mode and drive field at the position of the i th atom. Experiments to date have realized either Ising interactions 12 H ∝ F 2 z or spin-exchange interactions 13,14 H ∝ F + F − . We now show how to extend the approach of Ref. 14 to realize generic XXZ models of the form where the anisotropy ∆ is tuned by the angle of a magnetic field. The particular experimental setup 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 B = Bẑ, which defines the quantization axis for the spins, is oriented at an angle θ to the longitudinal axisĉ 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 Ω i of the control field and vacuum Rabi frequency 2g i of the cavity, where i denotes the value for the i 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 I i , representing the local polarization and intensity of light. This effect is described by a Hamiltonian where χ is the vector ac Stark shift of a maximally coupled atom and the component of the Stokes vector along the cavity is

The field operators
include the quantum field a ± of the cavity for σ ±polarized modes, weighted by the local amplitude g i of the cavity mode, and displaced by a classical drive field with local Rabi frequency Ω i . The normalization is set by the vacuum Rabi frequency 2g of a maximally-coupled atom. We assume that the drive field has horizontal polarizationx =ẑ ×ĉ and is detuned by δ 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 ± to obtain where v = (a + − a − ) / √ 2 represents the vertically polarized cavity mode, and we have introduced the weights These weights determine the collective spin F defined in Eq. 3.3, which couples to the cavity mode. Then, for v † v 1 and for large detuning δ κ, ω Z compared to the cavity linewidth κ and Zeeman splitting ω Z , we find that the effective spin Hamiltonian is 9) as detailed in the Supplementary Material Section S1 A. We see that Eq. 3.9 matches the Hamiltonian (1.2) with couplings J i = χξ i S 1/2 N 1/4 sin θ/2 √ 2δ and anisotropy ∆ = 2 cot 2 θ. Note that arbitrary control over the set of weights ξ 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 B). The key parameter governing the interaction-to-decay ratio is the single-atom cooperativity η = 4g 2 /κΓ, where Γ is the atomic excited-state linewidth. Moreover, we find that the interaction-to-decay ratio is collectively enhanced, scaling as S √ N η 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 η requisite for observing these signatures in the experimental setup.
We now discuss the two lines at ∆ = 0 and ∆ = 1 in the dynamical phase diagram, Fig. 1. In particular, we demonstrate that for these values of ∆ and for all values of the spin S, the Hamiltonian (1.2) is quantum integrable, in the sense of possessing N independent, commuting bilinears. The Hamiltonian at ∆ = 1 is related to the rational Gaudin model 16 , which is well-known to be quantum integrable in the mathematically rigorous sense of possessing an underlying quantum group structure 17 . Meanwhile, the Hamiltonian at ∆ = 0 lies in a less wellknown class of "non-skew" Gaudin models, which arise from Gaudin's equations upon relaxing the constraint of antisymmetry under interchange of site indices [20][21][22][23][24] .
We begin by reviewing the problem first studied by Gaudin 16 : under what circumstances do the quantum spin Hamiltonians define an extensive set of commuting operators, [G (i) , G j) ] = 0? Here S α i denote quantum mechanical spin-S degrees of freedom, satisfying standard commutation relations [S α i , S β j ] = iδ ij αβγ S γ j . If the couplings w α ij are taken to be antisymmetric, with w α ij + w α ji = 0, then the G (i) mutually commute if and only if the Gaudin equations hold for all pairwise distinct {i, j, k} and {α, β, γ}. The isotropic solution w α ij = J i J j /(J i −J j ) defines the rational Gaudin Hamiltonians The all-to-all spin model from Eq. 1.2 at ∆ = 1 is simply a linear combination of rational Gaudin Hamiltonians and Casimirs, to wit By rotational symmetry, H conserves the total spin S tot = i S i , and the linear span of the G (i) ( J ) includes the squared spin S tot · S tot = i,j S i · S j . The mathematical structure of traditional Gaudin models has been studied in depth 17,31 .
Let us now consider relaxing the constraint of antisymmetric couplings. Then Gaudin's equations (4.11) must be augmented by two equations constraining "on-site" couplings, which read (4.14) The model from Eq. 1.2 at ∆ = 0 arises from an "non- to the usual Gaudin equation (4.11), augmented by onsite terms w 1 ii = w 2 ii = 1/2, w 3 ii = 0, which solve Eq. 4.14. The corresponding Gaudin Hamiltonians read By the Gaudin equations Eq. 4.11 and Eq. 4.14, these mutually commute and the Hamiltonian (1.2) at ∆ = 0 can be expressed as At spin-1/2, this coincides with a Hamiltonian obtained in Ref. 23, and can consequently be derived as a special case of the integrable spin-1/2 Hamiltonians considered in the recent work Ref. 24. However, the integrability of this model for arbitrary spin S does not appear to have been discussed explicitly in the literature before. We conclude that there is an integrable line in the phase diagram of the model (1.2) at ∆ = 0. By rotational symmetry about the z-axis, this Hamiltonian conserves S z tot = i S z i and (S z tot ) 2 lies in the linear span of the G (i) ( J ). Finally, we note that upon replacing commutators with Poisson brackets in the derivation of the Gaudin equations, the integrable structure identified for ∆ = 0 and ∆ = 1 remains unaltered in the classical limit (S → ∞) of the Hamiltonian.

V. EXTRACTING INTEGRALS OF MOTION FROM NUMERICAL OR EXPERIMENTAL DATA
Having characterized the integrable structure for ∆ = 0 and ∆ = 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 ∆ ∈ (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 {Ô a }: where i > j and the index a is a shorthand notation for (i, j, α). We note that this family of 3N (N − 1)/2 operators defines an orthonormal set with respect to the infinite-temperature inner product: where D ≡ Tr[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 Ô a (t) ≡ Φ|Ô a (t) |Φ , where |Φ is a random initial state (i.e. far from any energy eigenstate). A bilinear integral of motionÎ is a special linear combination of thê O a that remains constant in time, to wit (5. 19) Here and below, the overline denotes a time average, such as Ô a ≡ T 0 dt T Ô a (t) over a time interval [0, T ]. It is useful to recast the above equation in terms of the following time series matrix: Note that M a,t is a rectangular matrix with 3N (N −1)/2 rows and a continuum of columns indexed by t ∈ [0, T ], where T J 2 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.19, a 2-local IOM corresponds to a left zero mode of M , i.e. 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 which are also orthonormal: As mentioned above, Q l is an integral of motion if and only if σ l = 0. Furthermore, for small σ 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 25) which means that the singular value σ l is the standard deviation of the expectation value of Q l over the time interval [0, T ]. A small σ l entails that Q l (t) 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 S = 1/2 and S → ∞ for ∆ ∈ (0, 1) lines in the phase diagram of Fig. 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 autocorrelation functions directly from the singular values.

A. Identifying integrals of motion
We now focus on the spin-1/2 system with up to N = 14 sites and implement the technique proposed At ∆ = 1, we see N + 1 zeros corresponding to the N + 1 conserved charges that can be written as a sum over bilinear operators; these zeros are separated from the rest of the singular values by a "spectral" gap. At ∆ = 0.9, we see two zeros corresponding to the conservation of H and (S z tot ) 2 . We also see the liftoff of N − 2 singular values corresponding to the previously conserved bilinear charges at ∆ = 1. Note that they, too, are separated from the rest by a "spectral" gap.
above. We initialize the system in a random product state 54 |Φ 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 N (0, J 2 ) and we set J 2 = 1. We then record the expectation values of all the operatorsÔ a defined in Eq. 5.17 and construct the time series matrix M a,n (defined in Eq. 5.20) at each discrete time t n = nδt with δt = 1 J −2 , integer n, and up to a maximal time T = 10 3 J −2 .
Plot of the auto-correlation function G l (t) in a given disorder realization for N = 13 spins at ∆ = 0.75. The solid curves represent the numerically computed G l (t): the black curve corresponds to either of the two exactly conserved bilinear quantities; the red, blue, green, and magenta curves correspond to the next four modes (arranged by increasing singular value); the yellow curve corresponds to a mode in the middle of the singular value "spectrum." The dashed curves represent fits of the form G l (t) = ζ l exp (−t/τ l ) + g l through the data. that these small singular values correspond to operators that exhibit a slow decay because the system is close to the ∆ = 1 integrable point. We now test this hypothesis by directly examining the decay of these putative "slow modes." B. Characterizing the slow operators for 0 < ∆ < 1 We have seen that the nontrivial IOMs at the points ∆ = 0 and ∆ = 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 where the normalization D = (2S + 1) N ensures that G l (0) = 1. For conserved modes, we expect the autocorrelation 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 ∆ = 0.75 is shown in Fig. 3. We see that The brackets . . . J denote an average over 10 3 disorder realizations for the {Ji}. Different colors correspond to different modes: black corresponds to the two lowest and exactly conserved modes; red, blue, green, and magenta correspond to the next four modes; yellow corresponds to a mode in the middle of the singular value "spectrum". We find no strong dependence on the system size N : see Fig. 5 for a plot of the plateau value g l J as a function of the system size N for the l = 3 (red) mode at ∆ = 0.5. 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 ≤ l ≤ 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 ∆ on the entire segment [0, 1]. 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 modesQ l with higher powers of the known conservation lawsĤ andŜ z tot , such asĤ 2 ,Ĥ 3 , . . . , nor with projectors to energy eigenstates (for details, see the Supplementary Material S2 C). In contradistinction, the operators corresponding to higher singular values (l N ) decay to a vanishing, or very small, saturation value (see the Supplementary Material S2 D).
Altogether, in addition to the obvious bilinear IOMs, H and (S z 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 ∆ = 0 and ∆ = 1: the bilinear integrals of motion are transformed into quasi-bilinear ones, which retain appreciable support in the space of 2-local operators. In general, we can write the new integrals of motion as bilinear operators dressed by a sum over higher, 2n-local terms: i2n , (6.27) where Z l is the weight of the integral of motion I l on 2-local operators. The saturation value of the autocorrelation function ofQ l that we plot in Fig. 4 is, essentially, g l ∼ |Z l | 2 . It would be interesting to further characterize how the coefficients K α1...α2n i1...i2n , 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.27) is, in some ways, reminiscent of the Local Integrals of Motion (LIOM) in the Many-Body Localized (MBL) state 26,27,32 . The latter is characterized by quasi-local integrals of motion τ z i that are adiabatically connected to the microscopic degrees of freedom σ 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 directproduct 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 std(J i ) J i .
Lastly, we also find signatures of integrability in the spectrum of H: the level statistics are almost perfectly Poissonian at ∆ = 0, 1 and close to Poisson (although not exactly) at intermediate ∆ (see Supplementary Material section S2 E). Nonetheless, for 0 < ∆ < 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 33,34 .
C. 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  (2) (right panel) from Eq. 6.28 and Eq. 6.29, respectively: the round markers correspond to the unperturbed Hamiltonian H (1.2); the triangular and square markers correspond to = 0.01 and = 0.1, respectively. The error bars related to disorder averaging . . . J are included, but they are smaller than the size of the markers. We see that the plateau value for the unperturbed H is independent of the system size. Converserly, upon adding even a small perturbation 1, the plateau value decreases with N , suggesting that the autocorrelation function G l (t → ∞) vanishes for a thermodynamic system (N → ∞).
interaction. Specifically, we add the term 28) where the elements V ij are also sampled from a normal distribution N (0, 1).
We explicitly check that at > 0 and ∆ = 1 for H + H (1) there are only 4 zero singular values corresponding to exactly conserved and linearly-independent 2-local quantities: the Hamiltonian, S 2 tot , (S x tot ) 2 , and vanishing singular values corresponding to H + H (1) and (S z tot ) 2 . Second, we verify that the lowest bilinear modes that are not exactly conserved (i.e. either the l = 3 one at 0 < ∆ < 1 or the l = 5 one at ∆ = 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 (1) , even at 1, can spoil the integrability for a large system N 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: where the fields h i are also sampled from N (0, 1). Note that H +H (2) has a single zero singular value corresponding to (S z tot ) 2 for all ∆; this is due to the fact that the full Hamiltonian is no longer purely bilinear and that H (2) breaks the SU(2) symmetry at ∆ = 1. Aside from this effect, the behavior upon perturbing with H (2) is similar to that obtained by perturbing with H (1) , as shown in Fig. 5b.
Last, we consider the effect of adding the perturbation This term appears in the model which is similar to Eq. 1.2, but differs from it by the term H (3) , arising due to the commutator S + i , S − i . As noted in Ref. 8, the model Eq. (6.31) 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/ √ 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 F 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) for 0 < ∆ < 1 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 ∆ = 0, 1 for any value of the spin size S (characterized by bilinear IOMs), and at S = 1/2 for any ∆ ∈ (0, 1) (characterized by quasibilinear IOMs). The remaining line in the phase boundary of Fig. 1 corresponds to the classical, S → ∞, limit of the model (1.2) for ∆ ∈ (0, 1), which we now discuss.
Since Gaudin-type integrability at ∆ = 0, 1 persists for all values of the spin size S, it is natural to ask whether the integrability* structure at 0 < ∆ < 1 and S = 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 → ∞ leads to classical equations of motion that are amenable to numerical simulation.
These simulations allow us to analyze another boundary in the phase diagram, 0 < ∆ < 1 and S → ∞, where we find chaotic dynamics with a finite Lyapunov exponent, as explained in Section VII A. The presence of chaos in the infinite-S limit clearly implies that the S = 1 2 integrability * does not extend to all S, unlike the Gaudin-type integrability at ∆ = 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 B. 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 C.

A. Classical chaos
In the infinite-S limit, the model (1.2) behaves as a classical system of coupled spin degrees of freedom S α i on the unit sphere, whose Hamiltonian dynamics can be written in terms of Poisson brackets: where For our numerical investigation, we sample the random fields J i from the uniform distribution [−J, J] 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 obey 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 35 to study the divergence of classical trajectories and measure the leading Lyapunov exponent. Let S(t) = (S 1 (t), . . . , S N (t)) denote the 3Ndimensional vector describing the directions of all the spins at time t. We initialize the system in a random infinite-temperature state 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 δS(t), which lives in the tangent space of S 2 × · · · × S 2 at the point S(t); we further set δS(0) such that δS i (0) ⊥ S i (0) for all spins and δS(0) If we define the local effective field , we see that the Hamilton equations of motion (7.32) can be written as We immediately see that the variational equations of motion for the deviation vector δS can be written as where δF x,y We numerically integrate the coupled differential equations (7.35) and (7.36) to find the trajectory (S(t), δS(t)) in the tangent bundle up until a time T = 500J −2 in increments δt = 1 J −2 . We then compute the sensitivity, defined as d(t) = δS(t) 2 , or in full, on invariant tori specified by the N conservation laws is linear in time and, since we have defined the sensitivity as δS(t) 2 , we expect d(t) ∼ 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 {S(0)} to find exp log d(t) J,S . We find that the classical system exhibits chaotic dynamics and an exponential divergence of trajectories at intermediate 0 < ∆ < 1, and integrable dynamics at the special points ∆ = 0, 1.
Moreover, using the multiplicative ergodic theorem, we can define the maximal Lyapunov exponent 35 as Using the normalization δS(0) = 1 and our definition of the sensitivity from Eq. 7.37, we see that In practice, we compute the Lyapunov exponent by fitting a line λt+b through the late time behavior of log d(t), as discussed in Ref. 36. In Fig. 6b we plot the Lyapunov exponent λ J,S , averaged over disorder realizations {J i } and initial states {S(0)}, as a function of the anisotropy ∆ and find that the system exhibits the most chaotic behavior (largest Lyapunov exponent) at ∆ = 0.5. Second, we find that λ J,S tends to a finite value for large system sizes N , as shown in the inset of Fig. 6b.

B. 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 ∆ = 0, 1 and only 2 exactly conserved quantities, corresponding to the Hamiltonian H and (S z tot ) 2 , in the intermediate regime 0 < ∆ < 1. 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 where c 1 = √ 15/2 and c 2 = √ 5/2. As before, a = (i, j, α) 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.40 satisfy the orthonormality relation where [. . . ] S denotes an average over the infinite temperature ensemble, while the integral DS 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 = {S 1 (0), . . . , S N (0)} as independent random points on the unit sphere. We then time evolve with (7.32) 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δt ∈ [0, T ]. Repeating this for a large number N of initial conditions {S 0 } in the infinite-temperature ensemble, we construct the following matrix, analogous to the one in Eq. 5.20: Note that the average . . . S is with respect to S 0 in the infinite-temperature thermal ensemble and should not be confused with the quantum expectation values . . . (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 σ 2 l . Similarly to (5.25), we have In other words, σ 2 l 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 σ l (as shown in Fig. 7a) is similar, in several ways, to that obtained in Section VI A for the quantum spin-1/2 model 55 . At the first integrable point ∆ = 0, we obtain N zero singular values corresponding to the family of N spin-bilinear conserved quantities G (i) from Eq. 4.15. At the second integrable point ∆ = 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.13. Lastly, as shown in Fig. 7a, for intermediate 0 < ∆ < 1 we find two precisely zero singular values, corresponding to the two exactly conserved spin-bilinear quantities, H and (S z tot ) 2 . The small magnitude of the following singular values, for l = 3, 4, . . . , signals the presence of slow modes, which will be studied in the next section.

C. 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 where . . . S , as before, designates an average over the N initial conditions S 0 . As in the quantum case, the operators Q l are orthonormal such that G l (0) = 1 (as N → ∞). 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 σ l already informs us about the long time plateau value of G l (t). This can be seen from Eq. 7.44, which implies that In the second line, we used the normalization Q l (t) 2 S = 1; in the third line, we performed a change of variables u = |t − s| (recall that G l (u) = Q l (t)Q l (t ± u) 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 : Thus, g l,T is a finite-time proxy for g l . In the infinitetime limit, the relation (7.46) becomes Using Eq. 7.46 or 7.47, 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.
At intermediate 0 < ∆ < 1, we find that the slowest non-conserved modes, corresponding to l = 3, 4, . . . , (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 56 , 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 → ∞ 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. 25 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/λ L . Given that our classical system is surrounded by integrable lines ∆ = 0, 1 and (arguably) S = 1/2 in the (∆, S) parameter plane (see Fig. 1), we conjecture that it is also quasiintegrable. From this perspective, the existence of slow modes is compatible with the finite classical Lyapunov exponent found in Section VII A.

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 intermediate 0 < ∆ < 1, along with signatures of integrability at the special points ∆ = 0, 1 for any S and integrability at 0 < ∆ < 1 for S = 1/2. 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 14 , one may immediately extract the bilinear spin correlation functions O a (t) ∝ S α i (t)S α j (t) defined in Eq. 5.17. 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 timepoints 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 F defined in Eq. 3.3. This effect is described by Lindblad operators where the decay rate (derived in Sec. S1) is given by The collective dissipation can be suppressed by increasing the detuning δ 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 L n,(m,m ) = C m,m Γ sc |m m | n , (8.50) where m or m indicates the spin state of an individual atom indexed by n, and C m,m is an order-unity branching ratio. At large detuning, the scattering rate scales as Comparing Eqs. 8.49 and 8.51 shows that the cooperativity η 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 S α i (t)S α j (t) 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 Γ sc due to free-space scattering and at a rate γ 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 ± Lindblad operators. Thus, at weak to moderate cooperativity η 1 and large detuning δ > κ, free-space scattering dominates and the bilinears decay on a time-scale τ J 2 ∼ √ N ηSκ/δ. For strong coupling η 1, where free-space scattering is suppressed relative to cavity decay, the total dissipation can be minimized at a detuning δ ∼ √ ηκ, leading the spin bilinears to decay on a time-scale τ J 2 ∼ √ N ηS.
To compare the decay time τ 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 t * ≈ 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 η ∼ 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 ≈ 10 3 J −2 according to Fig. 3, which places a more stringent requirement √ N ηS 10 3 δ/κ. This regime is challenging to access for S = 1/2 but readily accessible with large-S subensembles, e.g., at η 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 → ∞) 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 0 < ∆ < 1, persist for larger spin S up to 1/N corrections (see Sec. VII C). 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 for 0 < ∆ < 1 in Fig. 6.
The chaotic dynamics observed in the classical limit S → ∞ 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 outof-time-order correlators (OTOCs) 1,11,[37][38][39][40][41][42] , which quantify the spread of operators in time via the commutator C(t) = [V (t), W (0)] 2 . The connection to classical chaos is made clear in the semi-classical limit: for operators V = S z i , W = S z j , one can show that, to lowest order in a 1/S expansion, C(t) ∝ (∂S z i (t)/∂φ j ) 2 for a small rotation φ j at site j about the z-axis 43 . Thus, semi-classically the OTOC C(t) measures the sensitivity of the coordinate S z i (t) to changes in initial conditions S z j (0), and may therefore be regarded as a quantum generalization of the classical sensitivity d(t) defined in Section VII A.
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 11,40 . In the cavity-QED system considered here 14 , this sign reversal is achieved by switching the sign of the laser detuning δ in Eq. 3.9. The resilience of such time-reversal protocols to experimental imperfections, including dissipation, has been analyzed theoretically in Ref. 44.
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 λτ 1, where τ is the characteristic decay time defined above. More specifically, given the Lyapunov exponents λ ≤ 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 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 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 45,46 and of fundamental importance for elucidating the quantum-to-classical transition 47 . 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 ∆ along the z-axis. We showed that our model at ∆ = 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 ∆ = 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 along the entire line 0 < ∆ < 1. In contrast to the special points ∆ = 0, 1, the integrals of motion at intermediate ∆ 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 at 0 < ∆ < 1 is a purely quantum phenomenon: by numerically solving the Hamilton equations of motion for the classical model (S → ∞), 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 ∆ = 1 (and, presumably, at ∆ = 0 as well 20,21 ) 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 48,49 , 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 j S x j S x j+1 + S y j S y j+1 + ∆S z j S z j+1 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 50 . 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-toclassical crossover by better understanding how classical chaos (and perhaps quasi-integrability) at S → ∞ 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 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 ∆, so that both the special points ∆ = 0, 1 and the intermediate range 0 < ∆ < 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. probability p = 1/6) eigenstate of eitherŜ x i ,Ŝ y i orŜ z i . 55 The 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 Oa is evaluated on a single classical configuration in the former case, while the quantum expectation value is a result of a coherent average. 56 We 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.

A. 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: To simplify the following derivation, we will assume that the weights ξ 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: where we have passed into a rotating frame with respect to the atomic Zeeman splitting ω 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 51 . 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: 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 † F − F − is an off-resonant process and is, therefore, accompanied by a rapidly rotating phase factor e −2iω Z t . As a result, this term quickly averages to zero on timescales t 1/ω Z , and we are justified in ignoring it. In fact, the only resonant 2-photon processes that survive are the terms proportional to F z F z , F − F + , and F + F − , since all other terms have rapidly oscillating phase factors. The result of this elimination scheme is an effective Hamiltonian for the spins: where δ ± = δ ∓ ω Z are the detunings from the two-photon resonance, κ is the cavity linewidth, and ζ(δ) = δ/[δ 2 + (κ/2) 2 ]. At large detuning δ κ, ω Z , Eq. S3 simplifies to Eq. 3.9, and we obtain the desired model with an anisotropy parameter ∆ controlled by the angle θ of the magnetic field.
The effective Hamiltonian (3.9), however, is obtained only if the F + F − and F − F + terms in Eq. S3 are balanced, which occurs perfectly only in the limit δ → ∞. More generally, at finite δ, the cross-terms F x F y do not cancel and we obtain an effective Hamiltonian which is of the same form as Eq. 3.9, but with additional non-uniform longitudinal magnetic field terms. Such onsite terms, however, are sub-extensive relative to the spin-spin interaction terms. Moreover, we show numerically in Section S2 F that these additional terms do not affect the integrability of the model at S = 1/2 and finite N .

B. 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 respectively. Here, (m, m ) label internal states of the individual atoms indexed by n = 1, . . . , 2SN . C m,m are branching ratios that sum to unity. Γ sc,n = (Ω n /D) 2 Γ is the free-space scattering rate for the nth atom, where D is the detuning from atomic resonance and Γ is the natural excited-state linewidth. In comparison to the scheme presented in Section III, note that the Rabi frequency Ω n and detuning D enter likewise into the interaction strength in the Hamiltonian (3.9), where (χξ n ) 2 ∼ (Ω n /D) 2 g 2 n . Thus, the scattering rate scales as Γ sc,n ∼ (χξ n ) 2 /(ηκ). 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 Due to the rapidly oscillating phase factors e ±iω Z t , this relaxation operator may be split into three effective relaxation operators where we have introduced the collective decay rate By contrast, the free-space scattering operators L n,(m,m ) 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 τ 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: where r runs over the set of relaxation operators. Plugging in the free-space scattering operators L r = L n,(m,m ) , and assuming an approximately uniform scattering rate Γ sc,n = Γ sc , one can show that bilinears of spin-S operators S α j obey the following equations of motion for j = k: where we have introduced the notation Thus, the bilinear operators decay on a timescale τ ∼ 1/Γ 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 δ ± ≈ δ κ, ω z , one can show that the spin bilinears obey the following equations of motion for j = k: Assuming that the weights ξ i are numbers of order one, the bilinear operators therefore decay on a timescale τ ∼ 1/γ.
A key result of this analysis is that the collective dissipation rate is, somewhat surprisingly, independent of the atom number ∼ N S, exhibiting no superradiant enhancement. The reason for this is that the average effects of the L ± Lindblad operators cancel in the large-detuning limit, where δ + ≈ δ − . 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 (γ) and via spontaneous emission (Γ sc ) are controlled by the detuning δ. The total decay rate γ +Γ sc , at fixed interaction strength J 2 , is minimized by choosing a detuning δ = √ 1 + ηκ. For weak to moderate cooperativity η 1, we may instead choose a larger detuning satisfying the conditions δ κ, ω Z assumed above.

A. 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 ∼ 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 ρ(E). On the other hand, the energy fluctuations ∆E (at least at high temperatures) scale as ∆E ∼ O( √ N ) which leads to an extensive specific heat, as one would normally expect. Let us start by rewriting the Hamiltonian as where F ≡ i J i S i . Each operator F 2 α has non-negative eigenvalues and we denote the largest one by f 2 0 . Note that it is the same for all α ∈ {x, y, z} by isotropy. Thus, the largest eigenenergy E M of H is bounded by We can also bound E M below by using any properly normalized variational state |ψ since as a variational state. Then ψ| H |ψ = ∆ S √ N ψ| F 2 z |ψ and (S16) that (S17) Since the fields J i are sampled from N (0, 1), the random variable |J i | has a half-normal distribution with mean µ = 2/π and variance σ 2 = 1 − 2 π . Then Z ≡ i |J i | has mean N µ and variance N σ 2 so Z 2 J = N 2 µ 2 + N σ 2 . Thus, at large N , ( i |J i |) 2 J ≈ 2N 2 /π. After disorder averaging Eq. S17 we get which shows that the disorder averaged E M J ∼ O(N 3/2 ). We now put bounds on the ground state energy E m . Clearly, E m ≥ 0. More importantly, we note that which shows that E m J ∼ O(N 1/2 ). Since E M J ∼ O(N 3/2 ), this concludes our proof that the disorder-averaged bandwidth scales as W J ∼ N 3/2 . We have also computed W 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 typ ∼ O(N ) and that the density of states has a long and thin tail extending all the way to O(N 3/2 ).

B. Energy fluctuations
We now show that the disorder-averaged thermal energy fluctuations (∆E) 2 at high temperatures scale as (∆E) 2 J = cN for some constant c. Note that, by the fluctuation-dissipation theorem, this means that the specific heat per particle behaves as Cv N ≈ c T 2 at high temperatures T . The energy fluctuations are defined as where . . . T ≡ Tr e −βH −1 Tr e −βH . . . denotes thermal averaging. For simplicity, we will focus on the S = 1/2 and ∆ = 1 case. At lowest order in β we can approximate   These results correspond to a fixed disorder realization for the {Ji} in a system of N = 9 spins (S = 1/2) with ∆ = 0.9. a) The correlation function G l (t) for the slow modes Q l with l > 4 in a space orthogonal to the one spanned by the operators defined in Eq. S25. We use a color gradient between l = 5 (blue curve) and l = 14 (red curve). b) The plateau values after we substract the contributions fromĤ p as a function of the power p (see the definition in Eq. S26). We use a color gradient for the solid curves in which l = 3 corresponds to the blue curve and l = 9 corresponds to the red curve; the red, blue, green dashed curves correspond to l = 10, 11, 12, respectively. (inset) The diagonal matrix elements n|Q3 |n 2 for the mode l = 3 as a function of the energy eigenvalue En. This example provides further evidence that the slow modesQ l are far away from being projectors on many-body energy eigenstates.
Let us first summarize why this is a legitimate concern. Recall that at ∆ ∈ (0, 1), the SVD analysis finds two exactly conserved spin bilinears,Ĥ and (Ŝ z tot ) 2 , along with another N − 2 slow modes Q l , where l = 3, . . . , 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.27) and of quantum integrability. However, it is possible the theQ l s overlap non-trivially with exactly conserved operators such asĤ p (Ŝ z tot ) p : 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 withQ l .
Furthermore, the support in the space of 2-local operators can be non-zero even in the thermodynamic limit N → ∞ due to the all-to-all nature of the interaction. For instance, considerÔ 1 = 1 √ NĤ . By the previous section, we know that 1 Tr [1] Tr Ô 2 1 is of order one even in the thermodynamic limit. Similarly, letÔ 2 = 1 N Ŝ z tot 2 such that 1 Tr [1] Tr Ô 2 2 is also of order one as N → ∞. 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Ô p 1Ô p 2 , that do not vanish in the thermodynamic limit are:Ô Above,Ô 3 can be generated byÔ 1Ô1 as long as ∆ ∈ (0, 1), andÔ 4 byÔ 1Ô2 (in fact, one gets a prefactor 1 N k J k ; it scales as 1/ √ N typically if J k has zero mean and as O(1) if J k = 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 Tr [1] Tr O i O j = δ ij . We then define a matrix R whose columns vectors are the O i s. This allows us to perform the SVD analysis in the operator space orthogonal to the one spanned by the O i s: using the time series matrix M from Eq. 5.20, we perform the SVD decomposition of 1 − RR † 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 B, we define the slow operators Q l for l ≥ 5 and compute their auto-correlation function G l (t) = 1 Tr [1] Tr Q l (t) Q l (0) , as shown in Fig. S3a. Thus, we see that the overlap with higher powers ofĤ andŜ z tot has a negligible effect on the plateau values of the slow modesQ l .
In the second approach, we numerically study the effect ofĤ p on the plateau values by analyzing the matrix elements Q l,nn ≡ n|Q l |n , where |n is an energy eigenstate, H |n = E n |n . Let us define the vector q l whose elements are q |Q l,nn |. 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 = q l 2 . Another way of studying the effect ofĤ p on the plateau value is to analyze the overlap between q l and the set of vectors E (p) : E (p) n = (E n ) p , 0 ≤ p ≤ P , i.e. the powers of the energy eigenvalues.
Once again, we perform a Gram-Schmidt orthonormalization of the above vectors to obtain a set { η i : which is a measure of the residual plateau value for the modeQ l after removing the contributions ofĤ 0 ,Ĥ 1 , . . . ,Ĥ p . Naturally, for any fixed mode index l, Q l,p 2 decreases from the plateau value Q l,0 2 = q l 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 ∼ N/2, at which pointĤ p becomes an N -body operator. In fact, Q l,p 2 depends roughly linearly on p and decays approximately to zero only when p ∼ 2 N N . We have also checked that the result is not altered upon including low powers ofŜ z tot (i.e. removing contributions fromĤ p (Ŝ z tot ) r , for p, r ≤ N ). Therefore, the plateaus of the slow modes are not caused by simple functions of the exactly conserved quantities,Ĥ andŜ z tot , at ∆ ∈ (0, 1). Finally, the inset of Fig. S3b shows that slow modesQ 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 ∆ are solely due to the overlap with the known conserved quantities. This reinforces the thesis of the existence of further nontrivial conserved quantities (6.27), as put forth in the main text.

D. The middle of the SVD spectrum
We now check that, at ∆ = 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.

E. 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{∆E n , ∆E n+1 }/ max{∆E n , ∆E n+1 }, where ∆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 ∆ = 1 the Hamiltonian has the full SU(2) symmetry and the eigenstates of H are grouped in blocks corresponding to irreducible representations (irreps). Each state is simultaneously an eigenstate ofŜ 2 tot with an eigenvalue S tot (S tot + = 1.0 : Semilog plot of of the disorder-averaged plateau values g l J = G l (t → ∞) J as a function of N for three modes in the middle of the SVD spectrum at ∆ = 1.0: the N + 2 mode which is expected to be the first decaying mode, the 3N mode, and the 3N 2 mode. We see that the latter two decay as N increases, but this is not completely clear for the N + 2 mode.
1) and ofŜ z tot with an eingenvalue S z tot ∈ {−S tot , . . . , S tot }. Thus, each block contains 2S tot + 1 eigenstates of H. We keep the states in a fixed sector of bothŜ 2 tot andŜ z tot by selecting the block S 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 ∆ = 1 the symmetry is U(1) T : U(1) is generated by U φ = exp (−iφS z tot ) and it corresponds to the conservation of S z tot , whereas T represents the anti-unitary time reversal symmetry. Moreover, there is a mirror symmetry corresponding toŜ z → −Ŝ z . For the cases wherein T 2 = −1, we have Kramers doublets and we keep the states solely in the S z tot = 1/2 sector because it is the most populous. For the cases where T 2 = +1, we keep states solely in the S z tot = 1 sector in order to account for the mirror symmetry as well (the states in the S z tot = 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 r J . We know that a chaotic system exhibits level repulsion and has GOE level statistics characterized by a value r GOE ≈ 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 r Poisson ≈ 0.387.
In Fig. S5a we plot r J (N ) for different values of the anisotropy ∆. We find that both ∆ = 0 and ∆ = 1 exhibit almost perfectly Poissonian level statistics. However, for intermediate 0 < ∆ < 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. 52 (see also the references therein), there are well-known counterexamples 34,48,53 of integrable systems whose level statistics are neither Poisson nor GOE.
Nonetheless, if we add the perturbation H (1) from Eq. 6.28 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 ∆. 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 ∆, which provides further qualitative evidence that H is integrable for S = 1/2.
=0.1 from Eq. 6.28 renders the level statistics of H +H (1) close to GOE (the square markers), signaling level repulsion and a chaotic behavior. (inset) We obtain qualitatively similar results for a Hamiltonian (1.2) consisting of S = 1 degrees of freedom. b) Energy levels of H within a fixed S z tot = 1/2 sector as a function of ∆ in a single disorder realization {Ji} for a system of N = 11 spins (S = 1/2). We see many "real" (i.e. not avoided) level crossings, as shown in the inset which is a zoom on two such events. This violation of the Wigner-von Neumann non-crossing rule provides further evidence of quantum integrability at intermediate 0 < ∆ < 1.

F. Stability to the commutator perturbation
We conclude by studying the effect of the perturbation (6.30) that arises upon recasting the Hamiltonian into the "flip-flop" form of Eq. 6.31. The term J i J j S + i S − j from (6.31) in the main text can be written as which leads an overall term 1 S √ N i J 2 i S z i 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 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 1 S √ N i J 2 i S z i is a linear combination of 1-body terms. The overlapping circular markers correspond to the lowest four modes of the unperturbed model (1.2), which has N + 1 bilinear conserved quantities. The square markers correspond to adding a perturbation 1 S √ N i J 2 i S z i and the different colors correspond to the lowest four modes (that are not exactly conserved) in increasing order of their singular values: red, blue, green, and magenta. While the perturbation decreases the plateau value, we note that there is no flow with system size: the perturbation is sub-extensive (suppressed by 1/ √ N ).