Entanglement and replica symmetry breaking in a driven-dissipative quantum spin glass

We describe simulations of the quantum dynamics of a confocal cavity QED system that realizes an intrinsically driven-dissipative spin glass. A close connection between open quantum dynamics and replica symmetry breaking is established, in which individual quantum trajectories are the replicas. We observe that entanglement plays an important role in the emergence of replica symmetry breaking in a fully connected, frustrated spin network of up to fifteen spin-1/2 particles. Quantum trajectories of entangled spins reach steady-state spin configurations of lower energy than that of semiclassical trajectories. Cavity emission allows monitoring of the continuous stochastic evolution of spin configurations, while backaction from this projects entangled states into states of broken Ising and replica symmetry. The emergence of spin glass order manifests itself through the simultaneous absence of magnetization and the presence of nontrivial spin overlap density distributions among replicas. Moreover, these overlaps reveal incipient ultrametric order, in line with the Parisi RSB solution ansatz for the Sherrington-Kirkpatrick model. A nonthermal Parisi order parameter distribution, however, highlights the driven-dissipative nature of this quantum optical spin glass. This practicable system could serve as a testbed for exploring how quantum effects enrich the physics of spin glasses.


I. INTRODUCTION
In a spin glass, quenched disorder and the resulting frustration of spin-spin interactions generate a rugged free energy landscape with many minima.This means that in some cases, below a critical temperature, the single paramagnetic thermodynamic state fractures into a multitude of distinct possible thermodynamic states [1].The number of such states is exponential in the system size.A consequence of this is that exact copiesreplicas-of such a system may cool into distinct thermodynamic states.This is replica symmetry breaking (RSB), which Parisi invoked [2,3] to solve the Sherrington-Kirkpatrick (SK) model [4].The SK model describes a network of spins with all-to-all couplings with random signs.The Parisi solution showed how RSB arises by studying the distribution of spin overlaps between different replicas, as captured by the Parisi order parameter and the ultrametric, clustered tree-like structure of the distances between replicas.Since these features depend on details of the different thermodynamic states, they cannot be identified purely by looking at averaged properties [5].
Replica symmetry breaking is one example of the idea of ergodicity breaking.In an ergodic system, the dynamics of the system explores all allowed states, such that time averages are equivalent to configuration-space averages; this equivalence between time-and configurationaverages fails in states with RSB.This has important consequences when considering the relation between in-dividual quantum trajectories of the system and the trajectory-averaged density matrix.Theoretically, studying individual trajectories corresponds to stochastic unraveling of the density matrix equation of motion [6][7][8].Physically, unraveling the dynamics into trajectories corresponds to treating the system-environment interaction as a generalized measurement of the system by the environment.Each measurement projects the system into a specific state, conditional on the measurement outcome.The sequence of measurement outcomes (and associated states) is called a quantum trajectory.
We show that quantum trajectories can act as replicas to directly probe RSB.To see this link, we first discuss a simpler case, that of symmetry breaking in a standard second-order phase transition.When a perfectly isolated quantum system undergoes spontaneous symmetry breaking, it enters a macroscopic superposition, or 'cat state,' of the symmetry-broken states.This cat state is extremely fragile: Any interaction between the system and the environment allows the environment to learn which state the system chose.Backaction from measuring the environment stochastically collapses the system into one of the symmetry-broken states.Thus, each run of the experiment yields a symmetry-broken state, although the ensemble of states, averaged over experimental runs, remains symmetric.
The above simple picture also extends to the case where, rather than a small number of symmetry-broken states, one has many complex ergodicity-breaking patterns, as in a spin glass.In the case of a cavity QED system, each thermodynamic state emits a characteristic pattern of photons into the environment.On each run of the experiment, the backaction from observing this field collapses the system into a distinct thermodynamic state.This corresponds to the notion of a 'weak' symmetry that is broken in individual experimental runs but not in the ensemble-averaged density matrix [9].Thus, because there is a one-to-one correspondence between thermodynamic states and emission patterns, the overlap distribution is accessible through the correlations between the photon measurement records on distinct runs of the experiment.
We investigate the emergence of RSB in the open quantum system dynamics of confocal cavity QED, an experimentally practicable setting [10][11][12][13].In this system, atoms represent individual spins, while the cavity provides an all-to-all but sign-changing random interaction, dependent on the position of the atoms.This position dependence means it is possible to achieve random but repeatable interactions by controlling the placement of the atoms.By monitoring the spatiotemporal correlations of the light leaking out of the cavity, one can reconstruct the dynamics along individual trajectories.Because monitoring provides access to these correlations, the cavity QED setting gives us a powerful way to study RSB.In the RSB phase, the dynamics along each trajectory reaches a specific nonergodic state, so the spin configuration (and hence that of the emitted light) is stable over time on that trajectory.We quantify the distribution of overlaps among the patterns from different trajectories and the resulting Parisi order parameter.Together these show the distinctive features predicted by the Parisi ansatz for the SK spin glass.
We consider a realization where spins correspond to single atoms, giving spin-1/2 (and thus quantum) degrees of freedom, allowing entanglement to play a role in the emergent spin organization.We show that the quantum dynamics are distinct from the semiclassical limit, in which a semiclassical energy barrier severely inhibits passage to a low-energy manifold of states.By transitioning via entangled states, the trajectory avoids semiclassical energy barriers that would otherwise bar access to the low-energy spin manifold where RSB occurs.
Previously, we considered this same setup, but in the limit far above threshold where semiclassical approaches are valid, and considered the process of memory recall in such a device [14].Here, we address a very different question, focusing on the (necessarily) quantum dynamics near threshold, and the resulting distribution of lowenergy states found.A key point of this paper is that the final states at the end of the pumping sequence are classical, yet the ability to recover them relies on quantum dynamics.
The theoretical possibility of a spin glass phase in a multimode cavity QED system was suggested in Refs.[15][16][17].We note that the driven-dissipative nature of the confocal cavity QED quantum spin glass discussed differs from previous theoretical investigations of trans-FIG.1.(a) Sketch of the confocal cavity QED system.Transverse pump lasers (red) illuminate a network of atomic spin ensembles (orange) and scatter light into the cavity.The atomic spin ensembles at each node create either a collective spin or an effective spin-1/2 via Rydberg-blockade.Either way, the atomic spin ensembles are held in place by optical dipole 'tweezer' lasers (not shown).The confocal cavity field is composed of a local field (blue) at each spin ensemble and a nonlocal field (green) that mediates interactions between spin ensembles.The spin states are read out by imaging the cavity emission on a camera via spatial heterodyne detection [11].(b) Example simulated detection traces of integrated cavity emission for five spins; for clarity, the y-axis is normalized to the maximum signal magnitude after 4 ms.Each spin organizes into one of two orientations above the semiclassical superradiant threshold, indicated by a dashed line.(c) Plot of the ramp schedules versus time.'Normal' (N) and 'superradiant' (SR) regimes are to either side of the semiclassical threshold (dashed vertical line).(d) Pumping scheme for a 87 Rb atom.Balanced Raman transitions realize a pseudospin-1/2 degree of freedom.See text for details.
verse SK models in the closed-system context [18][19][20][21][22]. Experimental observations of RSB in physical settings have been reported in the spectra of semiclassical systems such as random lasers [23,24] and nonlinear wave propagation [25].Recent experimental results indicate that RSB in a confocal cavity QED system has been realized using a network with many XY spins per node [13].Recent theoretical work has noted that there can be phase transitions in the entanglement and correlations along individual quantum trajectories [26,27], even when such transitions are absent in the trajectory-averaged density matrix [28].
The paper is organized as follows.The next section describes the physical system we aim to simulate, its model Hamiltonian, and the Lindbladian dynamics to be unraveled by the trajectory simulations.Section III presents the quantum trajectory simulation method we employ, followed by results from individual trajectories in Sec.IV.Also discussed in Sec.IV are the evolution of entanglement entropy per spin and the difference in the lowest reachable configuration energy between entangled and semiclassical trajectories.The connection between quantum trajectories and replicas is presented in Sec.V. Evidence for RSB using the spin overlap order parameter is shown in Sec.VI.Section VII discusses the nonequilibrium nature of the system's overlap distribution via comparisons to equilibrium distributions.Section VIII provides evidence for spin glass, ferromagnetic, and paramagnetic phases.The emergence of an ultrametric structure between replicas is presented in Sec.IX.A summary and discussion of broader implications are in Sec.X.Eight appendices provide information on: A, the form of the effective spin connectivity matrices in confocal cavity QED; B, the derivation of the semiclassical critical coupling strength; C, the derivation of the atomonly theory; D, the stochastic unraveling of the master equation; E, the semiclassical limit of the spin dynamics; F, the Parisi distribution in terms of quantum trajectories; G, error bar estimation via bootstrap analysis; and H, the full set of overlap distributions used in forming the Parisi overlap order parameter.

II. THE CONFOCAL CAVITY QED SPIN SYSTEM
We now provide a description of a practicable system upon which we base this numerical study.All parameters have been experimentally realized or are plausible using existing technology [29].Photon-mediated spin interactions in the confocal cavity QED system were previously discussed in [14] and experimentally explored in [11,13,30].The system is depicted in Fig. 1.Two transversely orientated lasers of pump strength Ω ± scatter photons off a network of N spin ensembles, each with M spins per ensemble, where M can vary between one (realizing the spin-1/2 quantum limit) and 10 5 (describing current experiments [10,11,13,29,30]); see Sec.X and Ref. [31] for description of a practicable spin-1/2 scheme.These experiments employ 87 Rb in a 1-cm-long confocal cavity.To realize the network, each spin may be trapped at a position r i in the midplane of the cavity using an array of optical tweezers [32] (or optical dipole traps of larger waist [11,13]).
Similar to recent experimental work in which spins couple to a transversely pumped cavity [33][34][35], the (pseudo)spin states considered here correspond to 87 Rb hyperfine states The two pumps scatter light into the cavity via Raman transitions [36].We will consider an atomic detuning of ∆ ± ≈ −2π×100 GHz from the 5 2 P 3/2 atomic excited state and a controllable two-photon de-tuning ω z ≈ 2π×10 kHz.The maximum single-atom, single-mode light-matter coupling strength can reach a magnitude of g 0 = 2π × 1.5 MHz; see Sec.X for more details.We consider a confocal cavity of even symmetry under reflection in the cavity center axis [37].This restricts the possible cavity modes to the set of Hermite-Gauss TEM lm modes with indices l + m = 0 mod 2 [30].An even confocal cavity retains modes of sine and cosine longitudinal character, and trapping the atoms in one of these two longitudinal quadratures with optical tweezers further restricts the set of participating modes to l + m = 0 mod 4.This results in the effective Ising coupling we consider here; see Ref. [14] and App.C. We denote the remaining mode functions by Ξ µ (r), indexed by µ for brevity.A total of N m modes participate in the near-degenerate family of confocal-cavity modes to which the atoms couple [38].The modes are detuned from the mean pump frequency by ∆ µ = ω p −ω µ ≈ −2π×80 MHz.
The Hamiltonian in the rotating frame of the pump, after adiabatic elimination of the atomic excited state, is a multimode Hepp-Lieb-Dicke model [39][40][41]: The cavity modes are described by the bosonic operators a µ , while the spin ensembles are described by the collective spin operators S x/y/z i to facilitate generalization to the M > 1 case.In the spin-1/2 limit, S x/y/z i = σ x/y/z i /2.The effective coupling strength g = √ 3g 0 Ω ± /12∆ ± is the same for each pump laser, which can be achieved by controlling their pump intensities [14,33,34].Dissipation of the field of each cavity mode is incorporated using the Lindblad master equation where D[a] = 2aρa † − {a † a, ρ}.We consider a uniform cavity loss rate κ = 2π×260 kHz, similar to recent experiments [29].
The cavity-mediated interaction J ij between spin ensembles i and j may be derived using a polaron transformation of the Hamiltonian [14,30].In the ideal confocal limit, it takes the form where ∆ C is the detuning from the fundamental mode.
The first two terms are local and mirror-image interactions.The local interactions are shown in blue in Fig. 1(a); for clarity, the mirror images are not shown.They arise from the constructive interference of cavity modes at the positions of the spins and their image across the cavity axis (due to the modes' even parity).The finite spatial extent of the spin ensemble and cavity imperfections regularize the delta functions to form short, but finite-range interactions with tunable length scale > ∼ 2 µm in realistic cavities [10,29].The interaction range is much smaller than the waist of the fundamental mode w 0 = 35 µm.We provide in App.A formulas for the J matrix that incorporates finite N m and size effects.
The third, nonlocal term generates all-to-all, signchanging interactions [10,11].The nonlocal field is depicted in green in Fig. 1(a).By choosing positions r i either close to or far from the cavity center, the nonlocal interaction can yield J matrices that interpolate between ferromagnetic (all J ij > 0) and spin glass regimes [14].The glass regime results from Js with randomly signed off-diagonal elements.That is, these have approximately independent and identically distributed off-diagonal elements, roughly equal parts positive and negative.When atoms are distributed over a sufficiently large area, the confocal cavity-induced J connectivity matrices exhibit eigenvalue statistics that approximately those of a SK spin glass, the Gaussian orthogonal ensemble [14].Glassiness might also be achievable in a confocal cavity without position disorder [42].
The transversely pumped system realizes a nonequilibrium Hepp-Lieb-Dicke phase transition [41].At t = 0, the system is in the 'normal' state with all cavity modes in the vacuum state and the spins pointed along S z i in |↓⟩.As discussed below, we consider a protocol where the coupling strength is ramped up as a function of time.At threshold, the system transitions into to a 'superradiant' phase characterized by macroscopically populated cavity modes and spins that spontaneously break the global Z 2 symmetry to align along ±S x i .A sharp increase in cavity emission heralds the superradiant transition.Concomitantly, the spins order and the phase of the light emitted from each spin ensemble is locked to the spin orientation.Thus, the spins can be imaged in real time by spatially measuring the emitted phase of the cavity light.Holographic (spatial heterodyne) imaging has already been demonstrated [11].
In App.B, we provide a derivation of a general expression for the critical coupling strength g c at which the superradiant threshold is reached in the semiclassical limit.This is given by Note that it depends on the J matrix through its largest eigenvalue λ max : The form of the J matrix determines both the threshold and the character of the ordered phase-e.g., ferromagnetic versus spin glass [14].
(4) The ramp is the lowest-order polynomial that smoothly interpolates between two points with vanishing first derivative; the results that follow are insensitive to the precise functional form.The final pump intensity is 5× the critical value in the semiclassical limit.Given the experimentally relevant parameters employed, the timescales chosen are sufficient to allow the spins to reach an organized steady state before the onset of spontaneous emission, which occurs on approximately the 10-ms timescale.In addition, we choose to simultaneously ramp down the transverse field as ω z (t) = 2π×10[1−f (t)] kHz.This could be accomplished by changing the two-photon detuning of the pumps [33].Ramping ω z to zero turns off spin flips between different S x states because the Hamiltonian becomes diagonal in the S x basis.The ramps for g(t) 2 and ω z (t) are plotted in Fig. 1(c).

III. QUANTUM TRAJECTORY SIMULATIONS
Exact numerical simulations of (open) quantum manybody dynamics can be computationally expensive, especially in the confocal system due to the large number of modes in play.To explore the spin dynamics throughout the superradiant transition, we simplify the full dynamics to an atom-only Lindblad master equation whose derivation is a multimode generalization of the method of Jäger et al. [43]; see App.C for derivation.The atomonly Hamiltonian has the form where J ij is the same matrix as in Eq. (2).The other coefficients are where we restrict the treatment to the case of a completely degenerate cavity with uniform detuning ∆ µ = ∆ C for all modes.This is not an unreasonable approximation in the far-detuned regime |∆ C | ≫ κ, ω z [29].
In this limit, The Hamiltonian thus resembles a transverse-field Ising model with an additional term S x i S y j that is sufficiently small to play little role in the present simulations.The full atom-only master equation has the Lindblad form ρ  given by Here, v k i is the i'th element of the k'th eigenvector of the J matrix; all eigenvalues λ k ≥ 0. Each of the N collapse operators represents an orthogonal superposition of spin operators.Appendix C presents its derivation.
As noted above, the experimental protocol we consider involves ramping ω z to zero at late times.In this limit, α − goes to zero, so the final Hamiltonian has the simple Ising form H ∝ − ij J ij S x i S x j .Likewise, the collapse operators contain only S x i operators.As such, any S x i eigenstate may become a steady state above threshold (though some are energetically preferred; see Sec.IV).Quantum trajectory simulations of the atom-only master equation provide a continuous record of the state of each spin.These trajectories arise from simulating a sequence of balanced homodyne measurements of the field emitted from each spin ensemble; see App.D for details.This mimics experimentally practicable heterodyne measurements: Emitted cavity light is interfered on a camera with local oscillator (LO) light derived from the pumps to provide a phase reference [33].This procedure enables holography of the spin states; see Fig. 1(a) for illustration.Figure 1(b) shows what such data would look like as the homodyne signal for each spin is integrated in time.Each signal is dominated by noise below the superradiance threshold.Above threshold, the homodyne signals become phase-locked to the spins and undergo a bifurcation.The sign of each homodyne signal thus serves as a measurement of the corresponding superradiant spin state.
A single quantum trajectory evolves under the non-Hermitian Hamiltonian and is interrupted by quantum jumps with displaced collapse operators (C k ± iβ)/ √ 2; see App.D for details.These operators represent the two quadratures of the balanced homodyne scheme.The real number β is proportional to √ κ multiplied by the coherent state amplitude of the LO and the spatial overlap of the LO and emitted cavity light.Each collapse operator can induce quantum jumps.To simulate the detection of cavity emission, these are stochastically generated to independently occur at rates ||(C k ± iβ) |ψ⟩ || 2 /2.The trajectories of each spin are derived from these simulated detections, as shown in App.D. While the quantum state diffusion method is simpler to define, the quantum jump method with high LO strength (large β) provides similar results with greater numerical stability at late times.While β influences the timescale over which the global Z 2 symmetry is broken, we find it does not affect the ensemble of steady-state spin configurations that are found.

IV. QUANTUM SPIN DYNAMICS, ENTANGLEMENT, AND ENERGY BARRIERS
We now explore the dynamics of the spin trajectories and elucidate the role of quantum entanglement therein.Figure 2(a,b) plots two independent quantum trajectories for a network of N = 15, M = 1 (i.e., spin-1/2) particles that share the same J matrix.A glassy J matrix is selected by assigning the spins to random positions in the cavity midplane according to a 2D Gaussian distribution with a 2w 0 -wide standard deviation [14].
We observe that a spin-aligned σ x i steady-state configuration emerges within a few milliseconds of crossing the semiclassical transition threshold.Beforehand, a combination of unitary quantum dynamics and stochastic projections from the continuous measurement drive the spins away from their initial ⟨σ z i ⟩ = −S configuration.Measurement acts to break the spin's Z 2 symmetry along σ x i .The rate at which this happens is proportional to κg 2 /∆ 2 C and is time-dependent through g.The timescale is approximately 5 ms at threshold and decreases to approximately 200 µs at the end of the ramp schedule.However, the organization timescale also depends on the structure of eigenvectors v k of the J matrix.
We also observe that collective spin-flip events between different low-energy states occur beyond this timescale.A diverse range of collective spin behavior occurs.For example, Fig. 2(a) exhibits a group of three spins approaching a ⟨σ x i ⟩ = −1 steady state before collectively flipping toward ⟨σ x i ⟩ = 1 at around 750 µs.By contrast, Fig. 2(b) shows another behavior in which a group of four spins undergo an extended period of unbroken Z 2 symmetry before rapid organization into a steady-state configuration.
The spin trajectory in Fig. 2(a) may be visualized using the Bloch sphere representation in Fig. 2(c).As the many-body quantum state remains pure within a single trajectory, paths through the interior of the Bloch sphere indicate entanglement between spins.We see that the quantum spins first take a non-classical trajectory of unbroken global Z 2 symmetry through the interior of the Bloch sphere.After initially moving upward towards the center of the Bloch sphere, the continuous measurement breaks spin-flip symmetry.The spins then emerge from the interior of the Bloch sphere to reach a steady-state spin configuration.Figure 2(d) shows the average of the paths the spins take.
Entanglement is present during both the initial organization near threshold and during subsequent spin-flip events.We consider the entanglement entropy for each spin, given by − tr[ρ i log(ρ i )], where ρ i is the reduced density matrix for spin i.In general, the entanglement entropy can be nonzero for either entangled states or mixed states.Thus, the entropy would not be a good measure of entanglement when applied to the density matrix of the system.However, the entropy does provide a good measure of entanglement when applied at the level of individual quantum trajectories because each trajec- tory remains globally pure at all times.We choose the entropy over other measures of entanglement, such as the negativity [44], as it is more computationally tractable while faithfully capturing entanglement.
The entanglement entropy per spin is shown in Figs.2(e,f) for the trajectories in panels (a) and (b), respectively.An initial increase near threshold can result from unbroken global Z 2 symmetry and transitions to other low-energy states.At first, the superposition of low-energy states largely preserves the global Z 2 symmetry.Measurement then begins to lead to superposition projection around 0.5 ms, resulting in decreasing entanglement and the breaking of Z 2 symmetry.Subsequent spikes in the entanglement accompany the spin-flip events.Last, the entanglement decays to zero as the spins reach a steady-state configuration that corresponds to a classical state.Figure 2(g) plots the entanglement entropy for each spin averaged over 200 trajectories.The initial peak slowly decays to zero, reflecting the occurrence of later spin-flip events.
Figure 3 contrasts these quantum spin trajectories with semiclassical trajectories for the same J matrix used in Fig. 2. The spin-1/2 degrees of freedom are now replaced with semiclassical, collective spins, each comprised of M = 10 5 spin-1/2 atoms.The semiclassical equations of motion are derived in App.E. We see that, in contrast to the quantum dynamics, the semiclassical trajectories exhibit a rapid organization at the semiclassical transition threshold but are confined to the surface of the Bloch sphere, indicating the lack of entanglement.Unlike the quantum limit, large oscillations are observed around the x-axis on the Bloch sphere.
To investigate the role entanglement might play in the evolution toward low-energy, steady-state spin configurations, we plot in Fig. 4 the energy of 20 quantum and semiclassical trajectories for the same J matrix and ramp schedule considered in Fig. 2. The shaded region is inaccessible to any unentangled spin state constrained to the surface of the Bloch sphere.We identified the boundaries of this semiclassically forbidden region through Monte-Carlo sampling of semiclassical spin states (i.e., those states constrained to the surface of the Bloch sphere) followed by gradient descent to the lowest possible energy state.We find that entanglement enables the quantum spins to follow trajectories (through the Bloch sphere interior) that bypass this semiclassical energy barrier.This allows the quantum spin network to reach lower-energy steady-state configurations.Slower ramps could allow the semiclassical trajectories to follow a more adiabatic path back downward to similarly low-energy states.However, we find that such ramps must be at least an orderof-magnitude slower than those considered here.The addition of noise to the initial state can also yield a ∼25% decrease in the semiclassical steady-state energies, but this remains an order-of-magnitude higher compared to quantum trajectories.
The steady-state energy of the quantum trajectories seems to be primarily controlled by the ramp rate.Evolution through the superradiance transition has the form of a many-body Landau-Zener problem with many-body gaps controlling the adiabatic timescale of the transition.The many-body gap near the transition is on the order of ω z for a ferromagnetic J. Thus, ω −1 z sets the timescale for adiabatic evolution through the transition to either of the two ferromagnetic ground states.By contrast, spin glasses are characterized by nearly degenerate spin configurations that become exponentially numerous with N .This results in much smaller gaps near the transition and nonadiabatic evolution is more likely to occur, as we see in Fig. 4(b).The chosen ramp rate is slow enough to prevent nonadiabatic transitions to highly excited states, but not enough to prevent transitions to the nearly degenerate local minima states.Unitary evolution through the transition then produces an entangled superposition of low-energy states, as seen in Fig. 2, before projection into a single spin state occurs.The final energy of the trajectories is thus controlled by the nonadiabatic transitions experienced during the ramp as well as the measurement projection before ω z is ramped to zero.

V. THE OVERLAP ORDER PARAMETER
We now establish the link between quantum trajectories and the spin-glass order parameter.Order in glassy systems can be identified through correlations between the many symmetry-broken thermodynamic states.This is captured by the replica overlap [3], defined classically as q αβ = N i ⟨s α i ⟩⟨s β i ⟩/N where s α,β are replica spin states and brackets denote a time average.The overlap distribution is given by where n R is the number of replicas.Each P J (q) can have structure that varies depending on the disorder re- alization J, even in the thermodynamic limit; this is the lack of self-averaging inherent in spin glass [5].Disorderdependent fluctuations are averaged out in the Parisi distribution P (q) ≡ [P J (q)] J , where [•] J denotes an average over disorder realizations.We discuss the central features of P (q) in Sec.VI.The overlap distribution for glassy quantum systems is defined similarly after performing a Suzuki-Trotter mapping to an equivalent classical system [45][46][47]; see App.F for details.
To connect the overlap distribution to quantum trajectories, we first cast the overlap into a particular form that applies directly at the quantum level.In App.F we show that the overlap distribution is closely related to the operator where the σ x i are Pauli operators for each site.We refer to the above as the overlap operator given its close correspondence to the classical overlap q αβ .It acts on the doubled Hilbert space of ρ J ⊗ ρ J , where ρ J is the density matrix for a given J realization.The statistical moments of the Parisi distribution are shown to be given by q (k) = [⟨O k ⟩] J .This close relation allows for a simple expression of the overlap distribution in terms of O.To do so, we use the eigenstate representation O = q qP q , where the sum over q includes all N + 1 overlap values linearly spaced in [−1, 1].The operators P q are projections onto the space of spin states with overlap q.The overlap is then given by The connection to quantum trajectories is now established using the pure state representation ρ J = n T α=1 |ψ α J ⟩⟨ψ α J | /n T , where n T is the total number of pure states |ψ α J ⟩.Each quantum trajectory is one of these pure states.Inserting this form into Eq.( 10) yields (11) To summarize, trajectories of the same disorder realization can find different symmetry-broken states of the glassy landscape due to stochastic evolution induced by the environment.Each pair of symmetry-broken states |ψ α J ⟩ and |ψ β J ⟩ then contribute to P J (q) through their projection onto the subspace with overlap q.
The classical expression for P J (q) is recovered when the trajectories |ψ α J ⟩ are in spin eigenstates.Each state then corresponds to a classical spin vector with elements 11) then reduces to Eq. ( 8) with the replica overlap given by We refer to q αβ above as the mean-field overlap, as it corresponds to the overlap operator O in the mean-field limit-i.e., when trajectories are spin eigenstates and thus the wavefunction factorizes between sites.The fundamental difference between the classical and quantum overlap is how entanglement between spins can allow for superpositions of spin states.This allows a pairs of quantum states |ψ α J ⟩ and |ψ β J ⟩ to contribute to multiple values of the overlap distribution at once, which does not occur in the classical limit.

VI. REPLICA SYMMETRY BREAKING
We now explore the emergence of RSB as the system is pumped through the transverse Ising transition.To do so, we first analyze the correlations between independent quantum trajectories of a system with the same quenched disorder for all trajectories.That is, a frustrated spin system with the same J matrix for all trajectory simulations.Because the trajectories ultimately reach classical steady-state spin configurations, despite entangled quantum dynamics at intermediate times, we study the meanfield overlap q αβ in Eq. ( 12).This q αβ directly yields to the overlap distribution predicted by replica theory once steady-state spin configurations are reached, which occurs after ∼2 ms; it provides a mean-field estimate at earlier times.This q αβ depends on only first-order expectation values, and thus has the significant advantage of being directly observable from the trajectory measurement record.
Once a steady-state configuration is found, the overlap takes on one of N + 1 possible values ∈ [−1, 1].The overlap distribution is always symmetric about 0 due to the global Z 2 symmetry, in the absence of a longitudinal field.An ordered phase will exhibit an overlap distribution containing 'goalpost' peaks at q = ±q EA , where q EA = q αα is the Edwards-Anderson order parameter, also known as the self-overlap.(Paramagnets do not have such peaks, but ferromagnets and spin glasses do.)Peaks may also arise associated with overlaps q αβ between replicas that settle into different spin states.These additional non-vanishing peaks between the goalposts indicate RSB and arise from the smaller overlap between distinct, lowenergy states [5].
We first consider the overlap distribution for the same J matrix considered in Figs.2-4. Figure 5(a-e) shows the time evolution of q αβ as it approaches steady state, around 4 ms.To construct the overlap distribution of a fixed J matrix, we consider 200 quantum trajectories with identical initial conditions and the same ramp schedule as shown in Fig. 1(c).We then compute the overlap q αβ between every pair of the 200 replicas and bin the results as a function of time.We exclude the self-overlaps q αα because, while they have vanishing weight in the limit of an infinite number of replicas, they provide an asymmetric bias to only the positive q = q EA peak of the overlap for finite sample sizes.At t = 0, the system is in the normal (paramagnetic) state and the overlap between any two replicas is zero because ⟨σ x iα ⟩ = 0 for all spins in the initial σ z i state.A nonzero overlap emerges as the spins transition to the superradiant regime and align along ±σ x i .The final overlap distribution shows goalpost peaks at ±q EA .We find q EA ≈ 1 in the parameter regime of our simulation, but note that q EA is commonly overestimated due to finite-size effects [48].Interior peaks indicate RSB.These interior peaks arise from correlations between distinct spin configurations that are local minima of the Ising energy E = − ij J ij s i s j , where all s i = ±1.By local minimum, we refer to spin states for which flipping any single spin raises the total energy.We note that the steady-state overlap distribution is independent of the exact measurement scheme; for any LO strength β > 0 the same ensemble of spin states are found, leading to the same overlap distribution.
Figure 5(f-j) shows the full overlap matrix q αβ versus time before we bin into the above histograms.In the thermodynamic limit, the Parisi ansatz for the solution to the SK model predicts an ultrametric overlap structure emerging from RSB.Specifically, the ansatz predicts a nested block-diagonal structure where the overlap magnitudes are larger in the diagonal blocks than in the offdiagonal blocks.The diagonal blocks are then expected to further divide into smaller diagonal blocks with larger overlap and off-diagonal blocks with smaller overlap, and so on.The ansatz thus predicts a self-similar overlap matrix in the limit N → ∞, while for finite-size systems the self-similarity truncates at a finite depth.For the N = 15 case in Fig. 5(j), we find evidence for up to three levels of the RSB block structure.A primary 2×2 block structure emerges that approximately separates replicas 1-100 from 101-200.The primary block of replicas 1-100 is further subdivided into 2×2 blocks separated near replica 30, distinguishing regions of higher overlap from lower.Evidence for tertiary block structure may be found in the sub-block containing replicas 30-100; a final subdivision may be seen near replica 55.We leave to future work quantitative analyses of self-similarity, the depth of RSB, and how RSB scales with N .However, we later quantify the degree of ultrametricity in the overlap distribution in Sec.IX.
To provide further insight, we delve into the structure of the steady-state overlap distribution produced by the J matrix considered in Figs.2-5.This J matrix induces a rugged Ising energy landscape that contains six local minima not related by the global Z 2 symmetry.These were found by numerically enumerating all spin states.Of the 200 quantum trajectories in Fig. 5, 66% reached one of these six local minima in steady state.An additional 20% were within one spin flip of a local minimum, while the remaining 14% were between two-to-four spin flips away.To show the relation of each minima's oc-currence probability to its energy, we plot these together in Fig. 6(a), binning each trajectory by its nearest local minimum.The result shows a clear anticorrelation: The trajectories with lowest-energy, steady-state spin configurations are observed most frequently.Though the system is not in thermal equilibrium, this tendency to lowenergy states is reminiscent of a low-temperature system; see Sec.VII for a discussion of system temperature.
The overlap matrix between local minima is plotted in Fig. 6(b).The numbers in the matrix entries are the absolute value of the overlap values between the indicated minima.The diagonal entries correspond to the self-overlap, which is always unity.The overlap matrix allows us to pinpoint the pairs of spin configurations that create each peak in the overlap distribution of Fig. 5(j).This plot is reproduced in Fig. 6(c).Every peak in the distribution can be understood by considering the overlaps between the first 3 local minima in Fig. 6(b).Each peak in the distribution is annotated with the pair of minima x:y that produce that value of the overlap.The remaining local minima were found too infrequently to produce any distinct peaks in the overlap distribution.
Each J matrix produces a different set of local minima, and thus different overlap distributions.This is evident in an ensemble of 100 confocal J matrices produced by assigning spins to different random locations in the cavity midplane with standard deviation 2w 0 , which lies in the spin glass regime [14].The overlap distribution for each J matrix is constructed from 200 quantum trajectories as in Fig. 5. Resulting steady-state overlap distributions for three representative J's are shown in Fig. 7. Appendix H has plots of all 100 overlap distributions.
The overlap distributions all exhibit ±q EA peaks with q EA equal to unity for most disorder realizations.Variation within the interior demonstrates that the correlations between low-energy minima vary between J matrices.This is the non-self-averaging phenomenon inherent to SK spin glasses [5].The three J matrices are chosen to display a representative diversity of structure found in the overlap distributions.The J in Fig. 7(a) produces an overlap distribution that is dominated by a single lowenergy spin configuration.Other peaks (from different configurations) occur in only ∼10% of the trajectories.Figure 7(c) shows the other extreme in which many different spin configurations are found with multiple levels of clustering between states.This overlap matrix is indicative of a far glassier system.The character of most overlap matrices falls between these two for our system size of N = 15, such as the J in Fig. 7(b).
In the large-size limit, high peaks should be sparse in the overlap distribution because only a small set of thermodynamic states have significant weight.The peak positions do not average out into a smooth distribution between the goalposts.This is indicative of the lack of self-averaging manifest in these order parameter observables of the spin glass state.The order parameter that does average is the Parisi distribution [5], which we discuss in Sec.VIII below.

VII. EFFECTIVE TEMPERATURE
The overlap distributions do not appear to be consistent with an effective thermal equilibrium model.This is not surprising in this driven-dissipative quantum optical setting.Nevertheless, it is instructive to compare these distributions to those expected at equilibrium.Equilibrium overlap distributions can be constructed by assigning probabilities to spin states according to a Boltzmann factor exp(−E/k B T ), where E is the Ising energy of the spin state and T is an effective equilibrium temperature that serves as a fit parameter.The overlap between all pairs of states is then binned and weighted by their Boltz-mann factors.We perform a least-squares fit to each of the overlap distributions in Fig. 7 to extract T fit .The corresponding distributions are shown in red.The extracted temperatures are provided in units of T c , the largest eigenvalue of the J matrix averaged over all J matrices.This quantity corresponds to the critical temperature for the SK spin glass transition in the thermodynamic limit [4].In our finite system, T c corresponds to the average crossover temperature.
While the equilibrium model does capture the location of peaks in the overlap distribution, it is not able to quantitatively match their heights; instead, they seem to often be underestimated (overestimated) near the center (wings) of the distribution.The average T fit is 0.21(8)T c .Despite the lack of quantitative correspondence, the fitted temperatures are well-enough below T c to infer the presence of a low-energy ordered phase in this system, even when using realistic parameters.Indeed, the authors have observed such states in a related experimental system [13].

VIII. SPIN GLASS, FERROMAGNETIC, AND PARAMAGNETIC PHASES
Two order parameters are needed to distinguish between the different types of phases described by the Ising Hamiltonian.These are the J-averaged spin overlap distribution, also known as the Parisi order parameter, and the usual magnetization.The magnetization order parameter m = i ⟨σ x i ⟩/N is used to discriminate ferromagnetic from glass or paramagnetic ordering.It should approach ±1 at low temperature in a ferromagnetic state but be close to zero in the spin glass and paramagnetic phases.
The spin glass is distinguished from the paramagnet via the Parisi order parameter, the J-averaged overlap distribution.The average should form a smooth distribution for the SK spin glass.The paramagnet has a Parisi order-parameter distribution that is peaked around zero, while the spin glass and ferromagnet are peaked around ±q EA .Last, while the ferromagnet's Parisi distribution has no support between ±q EA , the spin glass has a "netwith-goalposts" structure of smooth interior support.
We average the overlap and magnetization distributions for 100 confocal J matrices to yield the aggregate distributions in Fig. 8(a).Distributions characteristic of a spin glass arise, consisting of extremal peaks at ±q EA bridged by a continuous interior for the overlap and a magnetization peaked near zero.This is consistent with the result expected from the Parisi solution to the SK spin glass [5].
A fit to the thermal equilibrium model yields a T fit = 0.21(2) T c , closely matching the average temperature found from the individual fits discussed above.The magnetization is well approximated by a centered binomial distribution, indicating that the local minima are uncorrelated with a ferromagnetic state.The standard devi- The same set of spin glass J matrices as in panel (a), but the system is rapidly quenched into the superradiant regime rather than ramped according to f (t).The overlap and magnetization both cluster around zero, indicative of a paramagnetic phase.The thermal fit is poorly constrained in this regime: Shown is the distribution for T fit = 5 T c, which bears similarity to the data.ation 0.31 of this distribution is close to that expected of the SK model at this system size, 1/ √ N ≈ 0.26.The difference from a Gaussian may indicate a small ferromagnetic remnant in the confocal J matrices that could be eliminated by placing spins further from the cavity midpoint [14].
By contrast, the ferromagnetic confocal J matrices reveal a very different behavior in Fig. 8(b).The ferromagnetic ensemble is constructed by using Gaussiandistributed spin positions with standard deviation of only 0.5w 0 in the cavity midplane.This leads to ferromagnetic J matrices with predominantly positive matrix elements and two global ground states corresponding to the two fully aligned spin states [14].Two-hundred quantum trajectories with the same ramp schedule as in Fig. 1c are used to construct the overlap and magnetization distribution per J matrix.The distributions are then averaged over the 100 matrices in the ensemble to produce the aggregate distributions in Fig. 8(b).The lack of support in the interior of the overlap and magnetization distributions indicates that only these two Z 2 -related spin states are found with high probability.This is consistent with a ferromagnetic phase.A least-squares fit to the thermal model yields a temperature T fit = 0.011(2) T c , where T c is twice the maximum eigenvalue for J matrices in the ferromagnetic regime [4].The T fit is lower than that found for the spin glass ensemble.This may be due to the larger energy gap to the ground state in the unfrustrated ferromagnetic J matrices.This makes it easier to maintain adiabaticity during the ramp, and fewer Landau-Zener transitions means a lower effective temperature.
A paramagnetic regime can be accessed by quenching the system into the superradiant regime rather than slowly ramping through the transition.In this case, adiabaticity is lost, and transitions into many excited states occur.Figure 8(c) shows the overlap and magnetization distributions that result from such a quench.The same spin glass J matrices as in Fig. 8(a) are considered, with 200 quantum trajectories per J, and all other parameters remain the same.Both the overlap and magnetization distributions are well approximated by centered binomial distributions of standard deviation 1/ √ N .This is indicative of a paramagnetic phase in which states are found at random.
Last, we present in Fig. 9 the dynamical evolution of the Parisi order parameter distribution for the spin glass.The distribution becomes, after around 2 ms, Fig. 8(a) in steady state.We also note that a finite-size scaling analysis of the Binder ratio 1 − ⟨q 4  αβ ⟩/(3⟨q 2 αβ ⟩ 2 ) is often used to pinpoint the exact location of the spin glass transition [49].However, the Binder ratio is ill-defined in this quantum system at early times because the overlap distribution begins as a delta function at t = 0, for which both the second and fourth order moments are zero.This happens because the spins begin aligned along σ z rather than σ x , a difficulty not encountered in typical equilibrium states of the classical Ising model.This makes the scaling analysis in this system more complicated, which we leave to future work.

IX. ULTRAMETRICITY
A prediction of Parisi's RSB ansatz for the SK spin glass solution is the formation of an ultrametric structure in the space of replicas [50].An ultrametric space is one satisfying the strong triangle inequality: Given any three points x, y, and z, the distances between those points should satisfy d(x, z) ≤ max[d(x, y), d(y, z)].It can be shown from this inequality that any triplet of points must form an isosceles triangle, either acute or equilateral.In the SK spin glass, replicas cluster into groups corresponding to low-energy local minima.Any triplet of replicas obeys this inequality in the thermodynamic limit, where the distance between replicas is the normalized Hamming distance, or equivalently d(α, β) = 1 − |q αβ |.
Numerical studies have verified that ultrametricity slowly emerges in system sizes up to 10 3 [51][52][53].The approach to ultrametricity is quantified by use of the metric K = (d max − d med )/σ(d), where d max is the largest distance in a given triplet of states and d med is the second largest (or the median).Their difference should be zero in an ultrametric space due to the isosceles condition.The difference is normalized by σ(d), the width of the distribution of distances between all states in the ultrametric space.
The overlap matrices and distributions in Figures 5-7 already exhibit the expected clustering of replicas into groups associated with local minima.Figure 10(a) demonstrates this even more clearly by plotting the associated dendrogram above the overlap matrix.The clustering of replicas into four primary groups is visible.
Figure 10(b) plots the J-averaged distribution of K as a function of system size.For each N , we generate 100 confocal J matrices in the spin glass regime and perform 200 quantum trajectories per J.For each J, the K distribution is computed between all triplets of tra-FIG.10.(a) The replica overlap matrix of 200 quantum trajectories for a confocal J matrix in the spin glass regime.The states cluster into one of four primary groups of states, which also appear as four primary blocks on the matrix diagonal.The above dendrogram shows the hierarchical clustering associated with the fracturing of the overlap matrix into four sectors, with various degrees of correlation between sectors.(b) The distribution of the K metric for ultrametricity, averaged over 100 realizations of J matrix for each system size.The distribution becomes increasingly peaked near zero, providing evidence of an ultrametric space emerging with increased system size.(c) The mean of the J-averaged K distribution as a function of system size, further showing the emergence of an ultrametric space.
jectories and the resulting distribution is averaged over J matrices.We begin the analysis at N = 8 because we find that low-energy local minima are not reliably present in smaller systems, with only a single cluster of states typically found.But even over the restricted range of N available, we find that the K distribution becomes increasingly narrow with increasing N .Moreover, Fig. 10(c) plots the mean of the K distribution with N , further showing a decrease in K with N .These constitute evidence for the emergence of ultrametricity in the system.Oscillations in ⟨K⟩ may arise from finite-size effects but appear to dampen with increasing N .We conclude that there is evidence for an approach to ultrametricity that is consistent with the significant finite-size effects found in SK spin glasses [52].

X. DISCUSSION
In summary, the transient formation of entangled states in a confocal cavity QED system allows RSB, concomitant with Ising symmetry breaking, to arise via the interplay of unitary Hamiltonian evolution and dissipative dynamics.The resulting Parisi distribution does not appear to be in equilibrium, which is expected given the driven-dissipative nature of the open quantum system.We note that previous work identified an effective temperature [14,54] associated with the multimode cavity QED dissipative dynamics.Such a temperature was found by considering detailed balance between energy-raising and energylowering processes in a system driven far above threshold.This effective temperature is much larger than T c in the spin-1/2 limit.(Although it is small in the semiclassical limit when M ≫ 1 [14].)Fortunately, however, this poses no obstacle to observing sufficiently lowenergy states and RSB because the ramp through threshold reaches a quasi-steady state much faster than the timescale ∝∆ 4 c /(g 2 ω 2 z κ) at which thermalization at T eff would occur [14].
Last, we note that the light-matter coupling strength required to reach the superradiant threshold for the spin-1/2 system is higher by a factor of √ M compared to the semiclassical limit.We estimate that this interaction strength could be achieved by Rydberg-dressing large atomic ensembles to yield a Rydberg blockade within each [31].This allows each spin ensemble to behave as if it were a single spin-1/2 degree of freedom while retaining the same collectively enhanced coupling strength ∝ √ M g 0 [55].Coupling to a Rydberg state could be realized through the addition of two pump lasers to the atomic level scheme in Fig. 1(d); Ref. [31] provides more details.Briefly, a laser at 780 nm would drive the |↑⟩ state to the atomic excited state 5 2 P 3/2 , while a blue beam at 479 nm would drive a transition from this excited state to the 100 2 S 1/2 Rydberg state.(Rydberg dressing inside optical cavities has been achieved [56].)The combined coupling terms produce a dark state that mixes the Rydberg and |↑⟩ states while avoiding the atomic excited state.Conservatively estimating an 8-µm average interatomic separation within a spin ensemble, a Rydberg-Rydberg interaction strength on the order of 200 MHz could be achieved.This should be sufficiently strong to push multiply excited Rydberg states far off resonance, resulting in an effective spin-1/2 degree of freedom for each spin ensemble.The spontaneous emission lifetime of the atoms is estimated to be greater than 10 ms-longer than the time scales shown above for RSB to emerge.This would allow RSB to be observed in a quantum optical context where implementations of, e.g., associative memory [14,17,[57][58][59][60][61]] might be realized.Doing so would provide experimental access to questions regarding how quantum effects might determine memory capacity and fidelity.
The research data supporting this publication can be accessed on the Harvard dataverse [62].
, where α is any complex number with real part greater than zero and w 0 is the waist of the fundamental mode.In a confocal cavity, a given resonance supports only the set of even modes or the set of odd modes.As such, it is useful to define the Green's function corresponding to either even modes (symmetrized) or odd modes (antisymmetrized).Choosing the even case, the symmetrized Green's function is defined as The confocal cavity-mediated interaction D(r, r ′ ) is expressed in terms of Green's functions [29,30]: The ideal interaction for a perfectly degenerate cavity with infinite mode support and delta-function-wide atomic ensembles is found by setting α = 0.The α = 0 limit corresponds to the J matrix in Eq. ( 2), where the term on the first line gives rise to delta function local and mirror interactions, while the other terms give rise to the nonlocal interaction.Allowing for α > 0 provides a good approximation for cavities with both mirror aberrations and finite-sized atomic ensembles [29].In this work, we use α = 0.02 to achieve a ratio of approximately ten between the local and nonlocal interactions, which roughly matches observations in recent confocal cavity experiments [10].This α yields local and mirror interactions with Gaussian waist much smaller than w 0 and a nonlocal interaction that has a large Gaussian envelope of waist much larger w 0 .While finite α does limit the maximum distance over which the nonlocal interaction can occur, it does not significantly affect the J matrices for this work, which considers atomic positions out to only ∼4w 0 .We note that the precise form of the confocal interaction depends on the details of both the atomic distribution and the nature of the cavity imperfections [29].However, these precise details matter little for the present work since they result in only small changes to the already disordered J matrices.
The degree of randomness in the J matrices produced by this cavity-mediated interaction was studied in depth in previous work [14].The elements of the J matrix become increasingly uncorrelated as w, the standard deviation of the atomic positions in the cavity midplane, becomes large compared with w 0 = 35 µm, the waist of the fundamental mode of the cavity.The correlation between randomly chosen J ij elements is less than one percent for w = 2w 0 , which is the value of w considered in the main text for generating glassy J matrices.This lack of correlation comes about because of the incommensurate periodic dependence of J ij on the positions r i , r j .The confocal J matrices produce eigenvalue spectra that approach a semicircle distribution, as expected for random matrices drawn from the Gaussian orthogonal ensemble (GOE), precisely like those of the SK model.At the system size N = 15 considered in the main text, the eigenvalue distribution is within 5% of the GOE semicircle distribution.

Appendix B: Derivation of semiclassical critical coupling strength
We now derive Eq. (3) for the critical coupling strength of the superradiant phase transition using linear stability analysis.The spin operators S α i are first mapped to bosonic operators b i through the Holstein-Primakoff transformation, where M is the number of atoms per ensemble.This transformation accurately models fluctuations around the normal phase when M is large.The original Hamiltonian in Eq. ( 1) is then transformed, up to a constant shift, to where N m is the total number of cavity modes and is the effective spin-photon coupling strength.The operator equations of motion, including cavity dissipation, are given by ȧµ The equations of motion are linear and thus directly solvable.To organize the set of operators, we introduce an operator-valued vector where the first 2N elements of u are the atomic operators followed by 2N m cavity operators.Using this notation, the equations of motion can be written in the concise form u = Au for a linear operator A.
The critical coupling strength g c can be found from the retarded Green's function, which describes the response of the system to an external drive.It takes , where θ(t) is the Heaviside step function.The Fourier transform is then defined by G R ij (ω) = dt e iωt G R ij (t).The Green's function is related to the linear operator A by [G R (ω)] −1 = S −1 (ω − iA) in the case of linear Heisenberg equations [40], where S ij = ⟨[v i (0), v † j (0)]⟩ are the equal time commutation relations.In this case, S = diag(+1, −1, +1, −1, • • • ) follows from canonical bosonic commutation relations.The full inverse Green's function, while large, has a simple 2×2 block form.The first few rows and columns of each block are shown below: We analyze the Green's function by first assigning blocks of the matrix: The diagonal matrix D spin is 2N ×2N with alternating elements ω z − ω, then ω z + ω.The matrix D cav is also diagonal of size 2N m ×2N m , with elements −∆ µ − ω − iκ, followed by −∆ µ + ω + iκ.In this expression, µ increases from 1 to N m .The matrix C describes coupling between the cavity modes and spin modes and is of size 2N ×2N m .It is most easily expressed in terms of 2×2 blocks given by We now determine when an instability in the normal phase occurs by considering the poles of the inverse Green's function.The poles dictate the characteristic response frequencies of the system, and thus the determination of when ω = 0 becomes a pole probes the global stability of the phase.This point can be found by considering when det [G R (ω = 0)] −1 crosses zero.This procedure relates to a normal mode analysis of the linear equations of motion, in which the system is stable only when all eigenvalues of A are greater than zero.The determinant can be written using the Schur complement [63] of ) is always positive.Thus, the instability condition simplifies to det D spin − CD −1 cav C T = 0.The matrix inside the determinant has a simple tensor product structure.At ω = 0 it is given by ) where ⊗ denotes the tensor product, I is the identity operator, and J is the cavity-mediated interaction connectivity matrix introduced in Eq. ( 2).
We must now determine when one of the eigenvalues of the above matrix crosses zero.The identity matrix simply shifts all eigenvalues by one.The eigenvectors of the total matrix now have the form v k ⊗ w, where v k is an eigenvector of J with eigenvalue λ k ≥ 0 and w is an eigenvector of the second matrix of ones.The second matrix has a zero eigenvalue and a nonzero eigenvalue 2 with eigenvector (1, 1)/ √ 2. We thus find that half of the eigenvalues of the matrix in Eq. (B9) are degenerate with value one, and the other half are given by The smallest value of g for which one of the above eigenvalues crossed zero occurs sets the critical coupling strength g c .The critical coupling thus depends on the largest eigenvalue λ max of the J matrix.Inserting λ max , setting the expression to zero, and solving for g yields the critical coupling strength of Eq. (3).

Appendix C: Derivation of the atom-only theory
We apply the method of Jäger et al. [43] to produce an atom-only theory for spins in a confocal cavity.The method accurately reproduces the low-energy spectrum of the single-mode driven-dissipative Dicke model, both below and above threshold.We extend the method to the multimode, multiple spin ensemble case described by the Hamiltonian in Eq. ( 1).
We must first find "effective fields" corresponding to operators in the spin Hilbert space that best approximate the effect of the cavity modes.There is one effective field Sµ for each cavity mode that is eliminated.The field may be time-dependent.The effective fields are chosen to satisfy the differential equations Solving for the effective fields begins with an ansatz of the form Inserting the ansatz into Eq.(C1) yields the following differential equations for the coefficients: This equation can be solved explicitly in the limit that the ramp function f (t) changes slowly over the cavity loss timescale 2π/κ ≈ 10 µs.This limit is well satisfied given that f (t) ramps over a period of 600 µs.In this limit, the differential equations are given by the solutions We now restrict ourselves to the degenerate case ∆ µ = ∆ C for all µ.The effective fields can then be written as with α ± given by Eq. ( 6).The effective fields can then be used to write the atomonly master equation.The Hamiltonian part is given by Inserting the effective fields from Eq. (C5) and recognizing the J matrix then leads to the atom-only Hamiltonian presented in Eq. ( 5).The full form of the master equation is now While correct, this expression is complicated to use because it involves an infinite sum of dissipation terms.This can be avoided by noting that there are only N linearly independent jump operators that can be created out of the N operators α + S x i + iα − S y i .As such, one can rewrite the dissipation term by expanding the effective fields.This yields where D[X, Y ] = 2XρY † − {Y † X, ρ} is the non-diagonal Lindblad superoperator and the Lindblad-Kossakowski matrix J is exactly the cavity-mediated interaction matrix defined in Eq. ( 2).Diagonalization of J brings the with collapse operators The square of the coefficients multiplying the collapse operators give the associated decoherence rates.A total decoherence rate per spin can be estimated by approximating the elements v k i of the normalized eigenvectors as uncorrelated random variables.Their variance should be 1/N to enforce unit normalization.We also approximate |α + | = 2 and α − = 0, which is valid well above threshold.The summed decoherence rate per spin can then be approximated as (κg Appendix D: Stochastic unraveling and reconstructing the spin measurement record The general formalism for homodyne unraveling has been well described by various authors [6][7][8].We provide a brief discussion of the specific measurement approach, and thus the associated unraveling scheme that we use in this work, which is based on invariance properties of the Lindbladian.
Spatial heterodyne detection is derived from the interference pattern of the LO and cavity light on a chargedcoupled device camera [11,33].While the measurement is a spatial heterodyne detection, meaning that the LO and cavity light have different propagation directions, the LO and cavity light possess the same optical frequency, as in a homodyne detection.We thus model the detection scheme as a balanced homodyne detection of the emitted cavity field.
To derive the unraveling corresponding to balanced homodyne measurements, we start by considering the atomonly master equation with Lindblad form where the collapse operators C k are given by Eq. ( 7).We can then manipulate the master equation to cast it in terms of the collapse operators corresponding to balanced homodyne detection.First, we write the Lindblad superoperators as a sum of two equal terms with rescaled collapse operators, We then use the shift-invariance property of the master equation, which admits shifts of a collapse operator C → C + aI if the Hamiltonian is also modified as H → H −i(a * C −aC † ).We perform shifts C k → C k ±iβI for each k and pair of collapse operators in Eq. (D2).
Here, β represents a real number proportional to the LO amplitude.Performing this shift yields the master equation (D3) The collapse operators (C k + iβ)/ √ 2 and (C k − iβ)/ √ 2 correspond to measurements of the two field quadratures, respectively.The Hamiltonian contribution from the two above collapse operators cancel out up to a global constant offset.The master equation is now unraveled into quantum trajectories using the standard quantum jump formalism [6] but with the above, shifted collapse operators.The extra field β means that the probability of a jump at each time step is higher, but the form of the collapse operator means that in such a jump, there is a smaller change to the state of the system.As one interpolates between β = 0 and β = ∞ this approach thus interpolates between quantum jumps (in terms of the original jump operators) and continuous quantum state diffusion.
The value of β used in our simulations is 0.1 √ κ.Recall that β corresponds to √ κ multiplied by the coherent state amplitude of the LO and the overlap between the LO and emitted cavity light.While this sounds small, the relevant quantity for determining how close the system is to the quantum state diffusion limit is the ratio of the number of detections per unit time to the rate associated with spin dynamics.The quantum state diffusion limit occurs when detections occur much more quickly than system dynamics.
The detection rate can be approximated for a given value of β [6] by The term in parentheses is the bare detection rate corresponding to the k'th collapse operator, where λ k is the k'th eigenvector of J.The term β 2 boosts the bare rate by the LO strength.The sum is taken over all collapse operators to approximate the total detection rate for each spin.This yields a detection rate of approximately 2 MHz at full ramp power when using typical J matrices in the spin glass regime; this translates to a timescale of about 0.5 µs.On the other hand, the spin dynamics typically occur no faster than ∼10 µs.Thus, we conclude that the dynamics are similar to those obtained in the diffusion limit.
The measurement records for each spin are constructed from the balanced homodyne signals where again v k i is the i'th element of the k'th eigenvector of the J matrix.The records s i (t) now change at a rate proportional to ⟨S x i ⟩, and thus, after sufficient integration time, their signs indicate the steady-state spin configuration.An example of the measurement records for a typical quantum trajectory is shown in Fig. 1(b).

Appendix E: Semiclassical equations of motion
We now derive the semiclassical equations of motion that describe a homodyne unraveling of the master equation.For these semiclassical calculations, we take the quantum state diffusion limit appropriate for large spins [6,8], corresponding to the limit β → ∞ in the previous section.
The stochastic differential equation describing the expectation value of an observable A is written in Itô form as where each dW k is the differential of an independent, Wiener process with ⟨dW k dW * k ⟩ = dt.The Wiener process can be either real, for homodyne detection, or complex for heterodyne detection.We consider a real Wiener process from here on for simplicity and without loss of generality.The semiclassical equations of motion are derived by first evaluating the above exact equation for S x/y/z i and simplifying the result through commutation relations.Any remaining product of spin operators A and B is then decomposed into a commutator and anticommutator, AB → [A, B]/2 + {A, B}/2.This step is necessary to achieve a semiclassical limit that retains stochastic terms from the homodyne detection.We then perform a mean-field decoupling of the anticommutator terms ⟨{A, B}⟩/2 → ⟨A⟩⟨B⟩ to arrive at semiclassical equations of motion: Above, α ± are as defined in Eq. ( 6).
Appendix F: The Parisi distribution in terms of quantum trajectories In this section, we first summarize the form of the Parisi distribution in a quantum spin glass, following Refs.[64,65].We then show how it can be understood in terms of an overlap operator in a doubled Hilbert space.This leads to the trajectory formulation of the overlap distribution presented in Eq. ( 11) of the main text.

Parisi order parameter for a quantum spin glass
For simplicity, we consider here the SK model in a transverse field.It is a close approximation of the more realistic model derived in Sec.II.The Hamiltonian is given by where σ x,y,z i are Pauli operators acting on one of the N total spins.The J ij couplings here are all-to-all with each element sampled independently from a Gaussian distribution with zero mean and variance σ 2 J /N .We consider here the equilibrium density matrix ρ = e −H/T /Z with partition function Z = Tr e −H/T .Application of the Suzuki-Trotter formula [45,46] allows for a reformulation of Z in terms of an equivalent classical model in a higher dimension.The classical energy is [47] F2) where s denotes the set of classical spin variables s i,τ = ±1.The term h c = T ln[coth(h q /LT )]/2 describes a nearest-neighbor coupling in the Trotter dimension indexed by τ , with periodic boundary conditions.The mapping becomes exact in the limit that the number sites L in the Trotter dimension tends to infinity.
Mapping the quantum partition function to an effective classical one unlocks the tools of replica theory.The free energy is evaluated via the "replica trick" log(Z) = lim n→0 (Z n − 1)/n, reducing the problem to computing the replicated partition function Z n for n replica spin systems s α .The free energy is then averaged over the quenched disorder matrices.This corresponds to the disorder-averaged, replicated partition function Z n given by E(s α , J) , (F3) where p(J ij ) = exp −N J 2 ij /2σ 2 J / √ 2πσ J and the sum is taken over all replica spin states.The disorder integrals can be computed exactly, which introduces a coupling between replicas.One finds that in the large-N limit the partition function can be expressed in terms of a local effective action with spins at different sites i decoupled, Z n = i e −Si .This gives a single-site action [64] The matrix Q αβ is the overlap order parameter that, after performing the disorder average, describes a coupling between replicas.The order parameter χ ∆τ is a replica-diagonal correlator describing a translationinvariant coupling in the Trotter dimension.Evaluation of Z n by the method of steepest descent gives selfconsistent equations for these order parameters, ⟨s α i,τ s α i,τ ′ ⟩ S , (F5) where ⟨A⟩ S = {s α } Ae −S /( {s α } e −S ) denotes an expectation value over the measure defined by S. While the overlap matrix may seem an abstract or purely theoretical object, Parisi realized [3] that Q αβ contains clear physical content: It describes the overlaps, or equivalently the distances, between the large number of distinct thermodynamic states in a single disorder realization.Finally, the Parisi order parameter P (q) is the distribution of elements of the Q αβ matrix, P (q) = lim (F6) The structure of P (q) [50,65] is modified in the quantum case by a transverse field [47,66,67].

Parisi distribution from the overlap operator
To establish a connection between Q αβ as predicted by replica theory and a physical observable of the quantum system, we consider the moments q (k) ≡ dqP (q)q k of the Parisi distribution.Using Eq. (F6) yields We now insert the self-consistency equation (F5) for Q αβ and make use of the fact that different sites decouple in the effective action, Eq. (F4).In the large-N limit, this yields The expectation value over the effective action S, in which the disorder has been integrated out, is now related back to expectation values at the level of individual disorder realizations.Following Ref. [68], each distinct replica index appearing in the expectation value under S corresponds to a distinct thermal average under the classical energy (F2).The result is then averaged to yield q (k) = [q (k) J ] J , where [•] J denotes the disorder average over J realizations and q (k) J is the moment associated with an individual disorder realization, given by q (k) (F8) Here, ⟨•⟩ E denotes a thermal average with respect to the classical energy in Eq. (F2).This can be related back to quantum expectation values as ⟨ j s ij ,τ ⟩ E = Tr ρ J j σ x ij , where ρ J is the equilibrium density matrix and the expression holds for any choice of τ [47].The key step is then to apply the simple relation Tr[X] 2 = Tr[X ⊗ X], yielding q (k) The term in parentheses to the power k is precisely the overlap operator O discussed in Eq. ( 9) of the main text, leading to relation q (k) = [⟨O k ⟩] J .The characteristic function φ(t) for the overlap distribution of each J is now straightforward to compute given the above simple form for the moments.Direct evaluation yields We obtain the overlap distribution through a Fourier transform of the characteristic function.Performing the disorder average then yields the Parisi distribution, where I is the identity matrix in the matrix exponential.
The exponential term can be recognized as the integral form of the Dirac delta function.To describe explicitly the action of the delta function, we expand the overlap operator in its eigenbasis as O = q qP q , where the sum runs over all N + 1 allowed values of the spin overlap q = −1, −1 + 2/N, • • • , 1.The operators P q are orthogonal projections onto the space of spin state pairs with overlap q.They mutually commute, [P q , P q ′ ] = 0, and have product relations P q P q ′ = P q δ qq ′ .Furthermore, the projectors span the full space of spin states and thus form a resolution of the identity, I = q P q .Inserting these forms into the integral expression of Eq. (F11) and evaluating the operator exponential, we find (q ′ − q) P q ′   (F12) = q ′ dt 2π exp [it (q ′ − q)] P q ′ = q ′ δ(q ′ − q)P q ′ .
Inserting this form back into Eq.(F11) produces the form of the overlap distribution presented in the main text, P (q) = q ′ δ(q − q ′ ) Tr (ρ J ⊗ ρ J )P q ′ J .
(F13) A variety of shapes are observed.Many distributions show interior support like peaks or a filling that is smooth.Interior support occurs in the most glassy J matrices where many nearly degenerate minima contribute to the low-energy manifold.Others show less pronounced interior structure.This occurs when the energy landscape is mostly dominated by a single ground state with energy much lower than any other local minimum.Even in these cases, there are often signs of smaller structures in the interior region arising from infrequently found local minima.This distinguishes the system from any ferromagnetic system, where only two goalpost peaks emerge from the single paramagnetic peak as the system cools.

FIG. 2 .
FIG. 2. (a-b) Example of two independent quantum trajectory simulations for the same system with a glassy J connectivity.The top panels show the dynamics of the N = 15 spin-1/2 system pumped through threshold.The spins begin to organize in the σ x i quadrature when the pump strength approaches the semiclassical normal-to-superradiant threshold at the time tc indicated by a vertical dashed line.Note, the quantities plotted differ from the simulated measurement records in Fig. 1b.(c) The same spin quantum trajectories from panel (a) are shown on (and within) the Bloch sphere.Note that here, the 15 traces are colored by radius rather than spin index.The red arrows show the general flow of the spin trajectories.(d) The averaged paths of the 15 spins, taken over 200 quantum trajectories.For each of the 15 spins, only trajectories with the same steady state are averaged together.(e-f) The bottom panels show entanglement entropy versus time of each spin for the trajectory simulations above; e.g., panel (e) pairs with panel (a).(g) The entanglement entropy per spin averaged over all 200 trajectories.

FIG. 3 .
FIG. 3. Example semiclassical trajectory simulation for the same frustrated J matrix and the same pump ramp schedule as in Fig. 2(a).(a) A sharper transition occurs with (b) dynamics restricted to the surface of the Bloch sphere.Panel (a) shows that the semiclassical transition, indicated by a vertical dashed line, is reached at time expected from Eq. (3).The red arrows in (b) show the general flow of spin trajectories.

FIG. 4 .
FIG. 4. (a) Plot of the energy of 20 independent semiclassical trajectories for the same J matrix as in Fig. 2. The energy is normalized to the instantaneous quantum ground-state energy E0.The N = 15 network has M = 10 5 spins per node.The shaded region is inaccessible to any unentangled state, and thus the semiclassical trajectories must climb over the barrier before reaching steady state.(b) By contrast, quantum trajectories can pass through the semiclassical energy barrier via entanglement between spins.This provides access to lower-energy steady states.Plotted are the normalized energies of 20 independent quantum trajectories for the same J matrix as above.Note the change in y-axis scale.The vertical dashed line marks the semiclassical threshold.The red (blue) colors in panel a (b) are chosen to be different shades to only help distinguish one trace from another.

FIG. 5 .
FIG.5.(a-e) Time evolution of q αβ for 200 quantum trajectories corresponding to the J matrix in Figs.2-4.Each panel shows the probability for each spin overlap value to occur at each given time during the ramp sequence.Error bars in this figure and in subsequent figures are estimated from bootstrap analysis; see App.G for details.(f-j) Time evolution of the full overlap matrix q αβ .The time that corresponds to each panel is the same as the panel above.The histograms in (a-e) can be recovered by binning the off-diagonal values in the overlap matrix.Each overlap matrix is ordered via an independent hierarchical clustering of spin states at each time.

FIG. 6 .
FIG. 6.(a) The energy of the six local minima of the Ising energy for the J matrix considered in Figs.2-5.Also plotted are the occurrence probabilities of those local minima as steady states for the 200 quantum trajectories.(b) The spin overlap matrix of these six local minima.(c) The same spin overlap distribution as in Fig. 5(j), but with annotated peaks.The notation X:Y indicates that the peak arises from the overlap of the local minima X and Y in the list of panel (a).Finite probability in the unlabeled bins is due to fluctuations of the overlap peaks around the labeled local minima.

FIG. 7 .
FIG. 7. Steady-state overlap distribution and overlap matrix for 200 quantum trajectories of three different confocal J matrices.The best-fit thermal distributions are shown in red.(a) Example of a J matrix that produces a dominant low-energy state; most trajectories find the same steady-state spin configuration.The temperature for the thermal fit is T fit = 0.17(1) T c.(b) A J matrix with multiple peaks and levels of structure in the overlap.The thermal fit yields T fit = 0.23(4) T c. (c) A more complex J matrix.Quantum trajectories find a large set of steady-state configurations.The thermal fit yields T fit = 0.17(3) T c.

FIG. 8 .
FIG. 8. Aggregate overlap and magnetization distributions in spin glass, ferromagnetic, and paramagnetic phases.The best-fitting thermal approximation is shown in red for all panels.(a) Confocal J matrices in the spin glass regime showing a smooth Parisi-like distribution and a magnetization peaked near zero.The thermal fit yields T fit = 0.21(2) T c.(b) Ferromagnetic regime showing the absence of interior support in the overlap and strong peaks in the magnetization.The thermal fit yields T fit = 0.011(2) T c. (c) The same set of spin glass J matrices as in panel (a), but the system is rapidly quenched into the superradiant regime rather than ramped according to f (t).The overlap and magnetization both cluster around zero, indicative of a paramagnetic phase.The thermal fit is poorly constrained in this regime: Shown is the distribution for T fit = 5 T c, which bears similarity to the data.

FIG. 9 .
FIG. 9. (a-e) Time evolution of the J-averaged Parisi distribution.The overlap distribution is averaged over all 100 J matrices in the confocal spin glass regime.The times are chosen to match those in Fig. 5.The steady-state distribution emerges after approximately 2 ms, demonstrating the distinctive goalpost peaks with a continuously filled interior.

FIG. 11 .
FIG.11.All 100 overlap distributions used in Fig.8(a).The y-axis is probability on a linear scale.The 100 J matrices used are all derived from that realizable in confocal cavities in the spin glass regime, as described in the main text.
The diagonalized atom-only collapse operators are non-Hermitian, in general, and Each record is initialized as h k (0) = 0. Every time that a jump occurs in (C k + iβ)/ √ 2, h k (t) is increased by one, while it is decreased by one every time a jump occurs in (C k − iβ)/ √ 2. The homodyne records do not yet reflect the spin measurements because the collapse operators C k are composed of a linear combination of spin operators from different sites.To construct the spin measurements records s i (t), the relation between homodyne records and spin states must be inverted by taking the linear combination