Engineering Graph States of Atomic Ensembles by Photon-Mediated Entanglement

Graph states are versatile resources for quantum computation and quantum-enhanced measurement. Their generation illustrates a high level of control over entanglement. We report on the generation of continuous-variable graph states of atomic spin ensembles, which form the nodes of the graph. The edges represent the entanglement structure, which we program by combining global photon-mediated interactions in an optical cavity with local spin rotations. By tuning the entanglement between two subsystems, we either localize correlations within each subsystem or enable Einstein-Podolsky-Rosen steering. We further engineer a four-mode square graph state, highlighting the flexibility of our approach. Our method is scalable to larger and more complex graphs, laying groundwork for measurement-based quantum computation and advanced protocols in quantum metrology.

Graph states are versatile resources for quantum computation and quantum-enhanced measurement.Their generation illustrates a high level of control over entanglement.We report on the generation of continuous-variable graph states of atomic spin ensembles, which form the nodes of the graph.The edges represent the entanglement structure, which we program by combining global photon-mediated interactions in an optical cavity with local spin rotations.By tuning the entanglement between two subsystems, we either localize correlations within each subsystem or enable Einstein-Podolsky-Rosen steering.We further engineer a four-mode square graph state, highlighting the flexibility of our approach.Our method is scalable to larger and more complex graphs, laying groundwork for measurement-based quantum computation and advanced protocols in quantum metrology.
Entanglement is a key resource for enabling quantum computation and advancing precision measurements towards fundamental limits.Crucial to these applications is the ability to controllably and scalably generate quantum correlations among many particles.A leading platform for achieving these ends are systems of cold atoms.Here, entangled states of over 20 atoms, such as cluster states with applications in quantum computation, have been generated by bottom-up approaches using local interactions [1].Conversely, global interactions among 10 2 to 10 5 atoms have been applied to prepare collective entangled states, including squeezed states [2][3][4][5][6][7] that enable enhanced precision in clocks [5,6,8,9] and interferometers [7,10].Such states, featuring symmetric correlations between all atom pairs, have been generated by collisions in Bose-Einstein condensates [2,3] and by photon-mediated interactions in optical cavities [5][6][7].
Atoms in cavities offer a particularly versatile platform for scalable generation of entanglement [5-8, 11, 12], with a single mode of light serving as an interface for correlating the atoms across millimeter-scale distances.In this setting, entanglement between spatial modes of an atomic gas has been achieved by splitting a global squeezed state into distinct subensembles [13], building on past work with optically dense ensembles in free space [14] and with spinor condensates [15][16][17][18][19]. Combining such top-down generation of entanglement with advances in local control and detection [20][21][22] provides the opportunity to engineer and probe richer spatial structures of entanglement, with applications in multimode quantum sensing [23], multiparameter estimation [24], and quantum computation [25].
A paradigmatic class of multimode entangled states are graph states [26], universal resources for quantum computation [25] with broader applications in quantum metrology [23] and in simulations of condensed-matter physics [27].These states, also known as cluster states, derive their name from a graph that defines the entangle- * Authors contributed equally ment structure, with edges representing correlations between nodes that may represent either individual qubits or continuous-variable degrees of freedom.Discretevariable graph states have been generated with superconducting qubits [28], trapped ions [29], and Rydberg atoms [1], while continuous-variable graph states have been prepared in photonic systems [30,31].Hitherto unexplored are opportunities for combining the benefits of light and matter to engineer graph states with flexible connectivity and long-lived information storage in atomic states.
Here, we report on the generation of programmable multimode entanglement in an array of four atomic ensembles coupled to an optical cavity.To control the structure of entanglement, we intersperse global interactions with local spin rotations.These two ingredients provide control over the strength of entanglement between subsystems and thereby enable a general protocol for preparing graph states.As a minimal instance, we prepare and characterize a two-mode graph state that exhibits Einstein-Podosky-Rosen (EPR) steering, a strong form of entanglement which is a resource for quantum teleportation and which has previously been demonstrated using collisional interactions in Bose-Einstein condensates [15][16][17]19].To illustrate the versatility of our protocol, we further construct a four-mode square graph state.Our work offers a blueprint for scalable generation of resource states for continuous-variable quantum computation and multimode quantum metrology.
As the mechanism for generating global entanglement, we implement cavity-mediated spin-nematic squeezing of spin-1 atoms [32].When a drive field is applied to the cavity (Fig. 1a), photons mediate spin-exchange interactions [33], and the system is governed by the Hamiltonian Here, F denotes the collective spin of all N atoms in the cavity, with spin length F ≤ N , and χ quantifies the collective interaction strength.In the second term, q parameterizes the quadratic Zeeman energy, proportional to the difference Q 0 = N 1 + N −1 − N 0 between the pop-  b, The resulting spin-nematic squeezing is visualized on a spherical phase space spanned by the collective spin-1 observables {F x , Q yz , Q 0 }.For short interaction times, the dynamics can be described on an effective two-dimensional phase-space spanned by the conjugate observables {x, p}.c, Combining the global interactions with local spin rotations allows for engineering a variety of entanglement structures, such as entanglement localized to selected subsystems and graph states with up to four nodes.

Local Spin Rotations
ulations N m of atoms in the m = ±1 and m = 0 Zeeman states.
We visualize the collective spin dynamics in a spherical phase space, analogous to the Bloch sphere, for spin-1 observables (Fig. 1b).We focus on a system initialized with all atoms in m = 0, i.e., polarized along the Q 0 axis.The effect of the cavity-mediated interactions is to twist the quasiprobability distribution of this initial state about the F x axis, inducing squeezing [34].Simultaneously, the quadratic Zeeman effect generates socalled spinor rotations about the Q 0 axis, mapping states along F x to polarized states of the quadrupole operator Q yz after a rotation of 90 • .The early-time dynamics explored in our experiments are well described by approximating a patch of the sphere as a two-dimensional phase space spanned by the conjugate observables x = F x / √ CN and p = Q yz / √ CN , which are normalized such that the Heisenberg uncertainty relation for x and p is Var (x) Var (p) ≥ 1.The contrast C, set by the commutator |⟨[F x , Q yz ]⟩| = 2CN , accounts for imperfect polarization along the Q 0 axis.
We engineer entanglement in an array of four atomic ensembles (Fig. 1a), each containing up to 5 × 10 3 Rubidium-87 atoms in the f = 1 hyperfine manifold.The ensembles are placed near the center of a near-concentric optical cavity with a Rayleigh range of 0.9 mm and are spaced by 250 µm.Applying a drive field to the cavity for 50 µs generates spin-nematic squeezing in the symmetric mode that directly couples to the cavity.To read out each ensemble i in a specified quadrature x i cos ϕ − p i sin ϕ, we map this quadrature onto the spin component F x via a spinor rotation by an angle ϕ.A subsequent spin rotation converts this signal into a population difference between Zeeman states, which we detect by fluorescence imaging.
To verify the generation of spin-nematic squeezing, we measure the variance ζ 2 = Var (x cos ϕ − p sin ϕ) for the symmetric mode x + = i x i /2 of all four ensembles.As shown in Fig. 2a, we measure a minimum value ζ 2 = 0.52 ± 0.07, limited primarily by technical noise [35].We confirm the presence of entanglement by evaluating the Wineland squeezing parameter ξ 2 = ζ 2 /C = 0.63 ± 0.08.Values below the standard quantum limit (SQL) ξ 2 = 1, shown by the dashed line at ζ 2 = C, indicate enhanced metrological sensitivity compared to any unentangled state of N atoms [10,35,36].We calibrate N from measurements of the atomic projection noise (Extended Data Fig. 1) and determine C from measured populations in the three Zeeman states (Methods Sec.F).
To demonstrate that only the symmetric mode couples to the cavity, we also evaluate the variance ζ 2 for the mode x − = (x L − x R )/ √ 2 which is anti-symmetric under the exchange of the left two ensembles x L and the right two ensembles x R .As expected, the variance for the antisymmetric mode shows no statistically significant dependence on ϕ and has an average value ζ 2 = 1.14 ± 0.04 near the quantum projection noise level.
We confirm the long-range character of the entanglement by evaluating a witness for entanglement [37] between the left and right subsystems, Here, x ′ + denotes the squeezed quadrature in the symmetric mode and p ′ − is the corresponding conjugate observable in the anti-symmetric mode.Generically W can take on any value since x ′ + and p ′ − commute.However, in the absence of correlations between the left and right subsystems, their independent Heisenberg uncertainty relations impose the constraint W ≥ 1, such that values W < 1 imply entanglement.The uncertainty product from the data in Fig. 2a is W = 0.55 ± 0.10, witnessing entanglement between the left and right subsystems.Consistent with the entanglement between subsystems, we observe a degradation in squeezing when measuring each subsystem individually, as shown in Fig. 2b.To further highlight that the left and right subsystems are in locally mixed states, we quantify the increase in phase space area due to the mutual information between them.For Gaussian states, the phase space area A m = ζ min ζ max for a mode m is the product of the standard deviations of the squeezed and anti-squeezed quadratures.Local measurements that discard correlations between the left and right subsystems yield a total phase space volume A L A R = 3.7 ± 0.4, larger than the total phase space volume A + A − = 2.2 ± 0.3 for global measurements of the symmetric and anti-symmetric modes.This emphasizes the loss of information when ignoring correlations between the local subsystems.
To optimize squeezing within each subsystem, e.g., for applications in spatially resolved sensing, the correlations between subsystems should be removed while maintaining the entanglement internal to each subsystem.Combining the global spin-nematic squeezing with local rotations provides the requisite control of the entanglement structure.To disentangle the left and right subsystems, we perform a sequence akin to spin echo, as shown in Fig. 3a.Between two pulses of interactions, we rotate the spins of the right subsystem by 180 • by optically imprinting a local vector ac Stark shift (Methods Sec.D).The effect is to cancel out interactions between the two subsystems, leaving only local squeezing (Fig. 3c).The scheme can equivalently be viewed as squeezing both the symmetric and anti-symmetric modes in the same quadrature (Fig. 3b).
More broadly, applying a sequence of squeezing operations in the basis of collective modes enables control over the spatial structure of entanglement via the relative orientations of the squeezed quadratures.Whereas a relative phase Φ = 0 between the squeezed quadratures of the symmetric and antisymmetric modes disentangles the left and right subsystems, the entanglement between subsystems can alternatively be maximized by introducing a relative phase Φ = 90 • via a spinor rotation in the sequence shown in Fig. 3a.The 90 • phase improves the entanglement witness W in Eq. ( 2) by producing simultaneous squeezing of both x ′ + and p ′ − .The resulting variances, shown in Fig. 3d, yield an entanglement witness W = 0.23 ± 0.05.The presence of squeezing in both orthogonal quadratures is indicative of entanglement of the paradigmatic Einstein-Podolsky-Rosen (EPR) type.
A notable feature of the EPR entangled state is its capacity for steering, in which measurements of one subsystem can predict measurements of both quadratures of the other subsystem to better than the local Heisenberg uncertainty product.Steering is a stricter condition than entanglement and enables teleportation of quantum information [38].To witness the left subsystem steering the right, we use measurements of the left subsystem to estimate x ′ R and p ′ R and calculate the error of the inference after subtracting a small detection noise contribution (see Methods Sec.G).The product of conditional variances Var 18 is less than one, the local Heisenberg uncertainty bound.The comparable witness for the right subsystem steering the left is 0.66 ± 0.18.We thus establish bidirectional steering at the 92% confidence level, which justifies identifying the state as a continuous-variable EPR state.
Our preparation of the EPR state constitutes a minimal instance of a scalable protocol for preparing graph states, in which the edges of the graph denote quantum correlations between conjugate observables on connected sites.Mathematically, this defining property of an ideal graph state can be expressed as where the adjacency matrix A encodes the connectivity of the graph and we implicitly sum over the repeated index j.As a general recipe for preparing a specified graph state, we diagonalize the adjacency matrix A to obtain a set of eigenvectors representing collective modes that should be squeezed.For each eigenmode m, the corresponding eigenvalue λ m specifies the orientation ϕ m = arccot λ m of the squeezed quadrature.q Fig. 3. Tunable entanglement: from local squeezing to EPR correlations.a, Scheme for controlling the strength of entanglement between left (L) and right (R) subsystems of the four-site array.After squeezing the symmetric mode (red), we transfer the squeezing into the anti-symmetric mode (blue) by applying a local 180 • spin rotation (green) to the right subsystem.Next, a global spinor rotation (purple) adjusts the angle of the squeezed quadrature.Finally, a second interaction pulse produces squeezing in the symmetric mode.The relative angle Φ between the squeezed axes of the collective modes determines the form of entanglement.b, To disentangle the left and right subsystems, we choose a relative phase Φ = 0 between the squeezed axes of the symmetric (red circles) and anti-symmetric (blue squares) modes.c, Entanglement internal to each subsystem manifests in variances ζ 2 = 0.41 ± 0.06 and ζ 2 = 0.38 ± 0.07 for the left and right subsystems (yellow squares and purple circles), respectively.d, To generate EPR entanglement between the left and right subsystems, we choose a relative angle Φ = 90 • between squeezed quadratures of the collective modes.The variances ζ 2 = 0.50 ± 0.07 and ζ 2 = 0.46 ± 0.08 for orthogonal quadratures of the symmetric and anti-symmetric modes yield an entanglement witness W = 0.23 ± 0.05 < 1. e, Representation of the resulting EPR entangled state as a graph state, corroborated by the reconstructed correlation matrix Corr (xi, pj).
The graph representing the two-mode EPR state is shown in Fig. 3e and corresponds to an adjacency matrix Diagonalizing A yields a state preparation protocol that matches the scheme of Fig. 3a: the eigenmodes of A are the symmetric and antisymmetric modes, while the eigenvalues λ ± = ±1 indicate that the squeezed quadratures should be oriented at ϕ ± = ±45 • , consistent up to a global rotation with the squeezing curves in Fig. 3d.Henceforth we work in a globally rotated basis chosen to orient the squeezed quadratures at the angles ϕ m .To visualize the equivalence of squeezing the collective modes with engineering the graph of entanglement, we use the data from Fig. 3d to reconstruct the correlations between conjugate variables in the two subsystems where Cov (•, •) denotes the covariance (Methods Sec.I).These correlations, shown in Fig. 3e, agree with the adjacency matrix A.
We additionally directly probe the graph of the EPR state by measuring the variances of the nullifers n i = p i − A ij x j in Eq. (3).As the ideal limit of zero variance requires infinitely strong squeezing, a practical definition of a graph state is that the variances of the nullifiers should approach zero in the limit of perfect squeezing.Defining normalized variances such that v i = 1 for a coherent state, our state preparation protocol theoretically produces variances assuming equal squeezing of all eigenmodes.Experimentally, we access each nullifier n i by performing a local 90 • spinor rotation on subsystem i.For the two-mode EPR state, with n L = p L − x R and n R = p R − x L , we measure variances v L = 0.53 ± 0.11 and v R = 0.36 ± 0.09 (Extended Data Fig. 3), directly confirming the entanglement structure specified by the graph.
To illustrate the scaling to more complex graphs, we produce the square graph state shown in Fig. 4a, with adjacency matrix The eigenbasis of A is shown in Fig. 4b.The eigenvalues λ m = (2, 0, 0, −2) specify squeezing angles ϕ m = (27 • , 90 • , 90 • , 153 • ) for the four eigenmodes.We sequentially couple each eigenmode to the cavity with the aid of local spin rotations, analogously to the scheme in Fig. 3a, squeezing the desired quadrature of each mode via global cavity-mediated interactions followed by a global spinor rotation (Extended Data Fig. 2).The result is shown in Fig. 4b, where the orientation of the squeezed quadrature for each eigenmode is within 5 • of the target squeezing angle ϕ m .Reconstructing the correlations Corr (x i , p j ) between sites from these measurements of the collective modes yields the matrix shown in Fig. 4c, which is consistent with the target adjacency matrix.We additionally directly measure the nullifiers n i for the square graph state.Their normalized variances v i , listed in Fig. 4e, have an average value 0.63 ± 0.07 consistent with the squeezing ζ 2 of the collective modes.Each nullifier further satisfies a condition v i < 0.94 ruling out separability into the independent nodes of the graph (Methods Sec.K), highlighting the presence of spatial entanglement between the four ensembles.
Our scheme for preparing graph states generalizes to any method of generating global entanglement that can be combined with local rotations.The approach is scalable to larger arrays, requiring only M squeezing operations to prepare arbitrary M -node graph states.For atoms in a cavity, the rate of each squeezing operation is collectively enhanced by the number of modes, such that the total interaction time required is independent of array size [35].Similarly, the degree of squeezing per mode is fundamentally limited only by the collective cooperativity per ensemble.Combining our approach with cavity-mediated generation of non-Gaussian states [11,12,39] or atom counting [40][41][42][43] opens prospects in continuous-variable quantum computation.Proposals for fault-tolerant measurement based quantum computation with continuous-variable cluster states require initial squeezing of 20.5 dB [44], which is comparable to the strongest demonstrated cavity-based spin squeezing [8].Continuous-variable cluster states additionally enable novel forms of quantum-enhanced measurement [23,24], including simultaneous sensing of displacements in conjugate variables [45] with applications including vector magnetometry.
Our protocol can be extended to a variety of platforms where either bosonic modes or qubits form the nodes of the graph and a central ancilla mediates collective interactions.Opportunities include generating continuousvariable graph states in multimode optomechanical systems [46], or in superconducting circuits featuring multiple microwave or acoustic modes coupled to a single qubit [47,48]; and discrete-variable graph states of individual atoms, superconducting qubits [49], or ions [50] with photon-or phonon-mediated interactions.Our approach offers the benefit of programmable connectivity and prospects for leveraging the central ancilla to perform quantum non-demolition measurements with applications in computation, error correction, and continuous quantum sensing.

A. Definition of spin and quadrupole operators
While for spin-1/2 particles all single-particle spin operators can be written as a linear combination of the dipole moments f x , f y , and f z , the space of spin-1 operators additionally includes quadrupole operators defined as 3 δ αβ where α, β ∈ {x, y, z} and δ αβ is the Kronecker delta function [3].For plotting the state on the generalized Bloch sphere, we use the operator q 0 = q zz + 1 3 , which quantifies the population difference between the m = 0 state and the m = ±1 states.We additionally construct collective observables corresponding to each spin-1 operator in a system of N atoms.

B. State Preparation
To prepare the array of four atomic ensembles in an optical cavity, we initially load 87 Rb atoms from a magnetooptical trap into an array of optical dipole traps, each with a waist of 6 µm.After optically pumping the atoms into the |f = 2, m = −2⟩ state, the ensembles are transferred into a 1560 nm intracavity optical lattice.Further details of the trapping procedue are described in Ref. [20].The atoms are then evaporatively cooled by decreasing the lattice depth from U 0 = h × 14 MHz to U 0 = h × 175 kHz in 200 ms.A series of composite microwave pulses [51] is used to transfer the atoms from |2, −2⟩ to |1, 0⟩.Any remaining population in the |1, ±1⟩ states is removed by first transferring this population into the |f = 2⟩ manifold using microwave pulses, and then applying resonant light to push and heat the |f = 2⟩ population out of the lattice.The lattice is then ramped up to a depth of U 0 = h × 25 MHz to minimize atom loss and increase confinement during the interaction phase of the sequence, yielding a final temperature in the lattice of 80 µK.During the interaction phase of the experiment, the ratio of the lattice depth to atomic temperature is U 0 /(k B T ) = 15 for an ensemble at the center of the cavity.

C. Interactions and Cavity Parameters
The spin-exchange interactions between atoms are mediated by a near-concentric Fabry-Perot cavity with length 2R − d, where R = 2.5 cm is the radius of curvature of the mirrors and d = 70 µm.The drive field is detuned from the 5S 1/2 , f = 1 → 5P 3/2 transition by ∆ = −2π × 9.5 GHz, after accounting for the ac Stark shift of the excited state due to the 1560 nm lattice.At the drive wavelength of 780 nm, the cavity mode has a Rayleigh range z R = 0.93 mm and waist w 0 = 15 µm, resulting in a vacuum Rabi frequency 2g = 2π ×3.0 MHz.
We parameterize the dispersive atom-light coupling by the vector ac Stark shift per intracavity photon, which for a maximally coupled atom is Ω 0 = − g 2 6∆ = 2π×41 Hz.As the array of atomic ensembles spans a length of 750 µm along the cavity axis, centered at the focus of the cavity mode, the maximally coupled ensembles experience a 30% larger Stark shift than the two minimally coupled ensembles.In addition, thermal motion of the atoms in the lattice means that the average atom experiences a reduced single-photon Stark shift compared with an on-axis atom at an antinode, resulting in a thermally averaged single-photon Stark shift Ω = 2π × 27 Hz.
Our method of generating cavity-mediated interactions is described in Refs.[33,52].The interactions are controlled by a drive field detuned from cavity resonance by an amount δ c .This corresponds to detunings δ ± = δ c ∓ ω z from two virtual Raman processes in which a collective spin flip is accompanied by emission of a photon into a cavity, where ω z is the Zeeman splitting.Rescattering of this photon into the drive mode is accompanied by a second collective spin flip, producing resonant spin-exchange processes of collective interaction strength where N is the total number of atoms and n is the intracavity photon number [52].We operate in a magnetic field of 4.1 G perpendicular to the cavity axis, corresponding to a Larmor frequency ω z = 2π ×2.9 MHz.The drive light is typically detuned by 2π × 4.2 MHz from the shifted cavity resonance, so that δ − = −2π × 1.3 MHz and δ + = −2π × 7.1 MHz.We define a total interaction strength χ = χ − + χ + .The drive light produces a typical intracavity photon number n = 800.A representative atom number N = 1.5 × 10 4 yields a collective interaction strength χ = −2π × 4 kHz.Exact parameters for each data set are detailed in Extended Data Table 1.
The parameters were selected to optimize squeezing, as discussed in Sec.IV of the supplement.

D. Global and Local Control over Spin Orientation
To access different quadratures of the squeezed states generated in our experiments and to adjust the relative squeezing angles of the collective modes, we apply global rotations about the Q 0 axis by two different methods.In the first method we let the system evolve under the quadratic Zeeman shift q = 2π × 1.2 kHz.Alternatively, we apply a detuned 2π microwave pulse on the hyperfine clock transtion |f = 1, m = 0⟩ ↔ |f = 2, m = 0⟩.For a suitable choice of detuning δ mw and microwave Rabi frequency Ω mw , the imparted phase is ϕ = π(1 − δ mw / Ω 2 mw + δ 2 mw ).This latter technique reduces the time required to rotate the orientation of the squeezed state before the final readout, since the Rabi frequency Ω mw = 2π × 7.5 kHz is much larger than the quadratic Zeeman shift.However, inhomogeneities in the microwave Rabi frequency on different ensembles can lead to unwanted population transfer from |1, 0⟩ to |2, 0⟩, which shifts the cavity resonance for subsequent interaction periods.Therefore, in sequences employing multiple drive field pulses to squeeze different collective modes, we use only the rotation under quadratic Zeeman shift to adjust the squeezing angle.
We apply local spin rotations around F y to read out the observables x and p, and rotations about F z to transform between collective modes.For these rotations, we use circularly polarized light that is blue-detuned from the 5S 1/2 , f = 1 → 5P 3/2 transition by 120 GHz.The laser beam is perpendicular to the cavity axes and is focused to individually address a single atomic ensemble, which we select by controlling the position of the beam via an acousto-optical deflector (AOD).The angle between the magnetic field, which defines our quantization axis, and the propagation direction of the laser is chosen to be 70 • .The circular component parallel to the magnetic field induces a vector ac Stark shift that acts as an artificial magnetic field, generating local rotations about F z .Rotations by 180 • about F z flip the sign of both F x and Q yz on selected ensembles.We thus utilize these rotations to transfer squeezing between orthogonal collective modes, as shown in Fig. 3a of the main text.For this transfer we simultaneously address two ensembles and induce the required spin rotation in approximately 18 µs.
The same laser allows for driving Raman transitions within the f = 1 hyperfine manifold, as the circular polarization component orthogonal to the magnetic field acts as an effective transverse field.Specifically, we use an arbitrary waveform generator to modulate the drive amplitude of an acousto-optical modulator, and thus the power of the laser, at the Larmor frequency.This rf modulation induces spin rotations about an axis in the F x −F y plane.Since there is no prior phase reference, we define the rotation to be around F y so that a π/2 pulse maps F x onto a population difference between Zeeman states.
To avoid differential evolution of the spinor phase ϕ, we typically perform global Raman rotations by simultaneously addressing all four ensembles (except for the direct measurement of the nullifiers described in Sec.J).In this setting, we achieve a global Rabi frequency of Ω Raman = 2π × 12.5 kHz.

E. Readout and Fluorescence Imaging
We characterize the multimode entangled states in our experiment by state-sensitive fluorescence imaging.To read out a specified quadrature in the x − p plane (where x ∝ F x and p ∝ Q yz ), we first perform a global spinor rotation by a variable angle ϕ and subsequently perform a 90 • spin rotation about F y to convert F x to F z .The implementations of these rotations are described in Sec.D. To ensure that the rotation angle stays close to 90 • during the whole duration of our experiments, we calibrate the frequency of the Rabi oscillation every hour.
For the data shown in Figs. 2 and 3 of the main text, where each subsystem (left and right) consists of two atomic ensembles, we modify the readout to minimize the impact of global technical fluctuations.Specifically, we apply a local 180 • rotation about one of the ensembles in each subsystem prior to the final spin rotation, thereby mapping the symmetric mode onto one that involves a differential measurement of F z between ensembles.Similarly, the anti-symmetric mode is mapped onto a mode that remains robust against technical noise.
To measure the atomic state populations, we collect a sequence of four images, with one detecting any population in the f = 2 hyperfine manifold and the remaining three images detecting the populations in the three magnetic substates m = 0, ±1 within the f = 1 manifold.For this portion of the experimental sequence, we lower the power of the 1560 nm trapping laser to reduce the ac Stark shift of the electronically excited 5P 3/2 state and reconfine the atoms in the microtraps.We apply two counter-propagating laser beams resonant with the f = 2 → f ′ = 3 transition of the D 2 line and collect the resulting fluorescence signal on an EMCCD camera.To avoid interference of the two imaging beams, we switch them on one at a time for 3 µs each and alternate between the two beams for 126 µs per image.After this time, most of the atoms in f = 2 have escaped the trapping potential due to heating, and we switch on one of the imaging beams for 150 µs to remove any residual atoms in f = 2.To measure the atoms in the remaining states, we transfer the population in the desired state to the f = 2 manifold via microwave pulses and repeat the imaging sequence above.To reduce the sensitivity of this transfer to magnetic field noise and microwave power fluctuations, we use a composite pulse that involves a sequence of four microwave pulses with different relative phases [51].
To calibrate the conversion from fluorescence signal to atom number, we employ a measurement of the atomic projection noise.We prepare N atoms in a superposition of m = ±1 by initializing all atoms in m = 0 and then rotating by 90 • about F y .To isolate the projection noise, we vary the atom number N = N +1 + N −1 and measure the variance of the population difference N +1 − N −1 .Extended Data Figure 1 shows these data in units of camera counts for each of the four collective modes.For each mode, we perform a polynomial fit, where c m denotes the signal from atoms in state m.The linear coefficient a 1 = r + g includes the count-to-atom conversion factor r and an additional contribution g ≪ r from photon shot noise, augmented by the excess noise of the EMCCD.From the fit value a 1 = 415 ± 6 and an independent calibration of g = 20, we obtain the conversion factor r = 395 ± 6 counts/atom.This calibration is consistent with an independent measurement of the dispersive cavity shift δ N = 4N Ω due to N atoms.
The fit offset a 0 and quadratic component a 2 provide information, respectively, about the imaging noise floor and technical noise in the fluorescence readout.The quadratic component of the fits in Extended Data Fig. 1 determines the atom number N ∼ 1/a 2 at which technical fluctuations become comparable to the projection noise.For the mode with the highest technical fluctuations, we find a quadratic component a 2 = 5 × 10 −5 .We therefore limit the maximal atom number in the experiment to N ≲ 2×10 4 to ensure that projection noise dominates over technical fluctuations.For our typical atom numbers, the background noise a 0 is small compared with the photon shot noise, the latter being equivalent to a fraction g/r ≈ 0.05 of the atomic projection noise.For the direct measurement of the nullifiers in Fig. 4 and the EPR steering, we subtract the photon shot noise contribution from the measured variances.

F. Measurement of Contrast C
To compute the normalized variance ζ 2 = Var (F x ) /(CN ), we extract the contrast C from the same data set as the variance of F x .Specifically, in terms of the Zeeman state populations N ′ m after the readout spin rotation, we measure both the spin component The quadrupole moment Q xx is directly proportional to the contrast C in our Larmor-invariant system.
In the most general case, the contrast C may be expressed exactly in terms of the collective quadrupole moments as Due to the Larmor invariance of the states generated in our experiment [35], ⟨Q xx ⟩ = ⟨Q yy ⟩.Further, the three quadrupole moments sum to zero, as can be seen from the definitions in Sec. A. We can thus reexpress the contrast as We use this expression to normalize all variances reported in the main text.

G. Steering Criterion
To confirm Einstein-Podolsky-Rosen (EPR) steering, we show that a measurement on the right subsystem can be used to infer the measurement results in the left subsystem with a higher precision than permitted by the local Heisenberg uncertainty relation.To calculate the error of the inference of an observable O of the left subsystem conditioned on measurements of the right subsystem, we find weights g i that minimize the conditional variance where i indexes ensembles within the right subsystem and the weights g i capture inhomogeneities in coupling for different ensembles.For the EPR-steered state, these variances are minimized for the x ′ and p ′ observables.We measure EPR steering in both directions, requiring inferences in two directions and two quadratures.The values of all of the conditional variances, after subtracting a small photon shot noise contribution as calibrated in Sec.E, are summarized in Extended Data Table 2.
We also report the optimal values of g i for each inference.For most of the inferences, higher weight is given to the ensemble closest to the center of the cavity, which we attribute to the difference in atom-light coupling for different ensembles.

H. Graph State Generation
Our prescription for preparing graph states rests on diagonalizing the adjacency matrix A, with the resulting eigenvectors specifying collective modes to squeeze and the eigenvalues specifying the squeezed quadratures.Formally, the diagonal matrix Λ of eigenvalues λ m is related to the adjacency matrix A by where the columns of V represent the eigenmodes.In terms of the quadrupole operators x = (x 1 , . . ., x M ) and p = (p 1 , . . ., p M ) on the individual sites, each eigenmode is parameterized by collective quadrature operators x = V x and p = V p. Re-expressing Eq. ( 3) in terms of these collective modes, shows that the antisqueezed axis for each mode lies along the line pm = λ m xm .Thus, the squeezed quadrature is oriented at a spinor phase ϕ m = arccot λ m .The experimental sequence for preparing the square graph state is presented in Extended Data Fig. 2. For this graph, where the columns have eigenvalues λ m = (2, 0, 0, −2) corresponding to the angles ϕ m = (27 We choose angles Φ 1,2,3 for the global spinor rotations so that each mode is squeezed along the appropriate axis at the end of the sequence.Schematically, we incorporate the global spinor rotations that occur during the pair creation dynamics and Larmor rotations into these angles.Our approach of squeezing the eigenmodes of the adjacency matrix allows for generating arbitrary graph states.In the most general case, the eigenmodes may have weighted couplings to the cavity, which could be controlled via the positions or populations of the array sites, or by driving the cavity from the side with a spatially patterned field.However, even with equally weighted couplings to the cavity, a wide variety of graphs are accessible by squeezing eigenmodes of the form V jm = exp(iφ jm )/ √ M .The phases φ jm can be imprinted by local spin rotations, generalizing the 180 • rotations applied in this work.For translation-invariant graphs, the eigenmodes are spin waves with φ jm = j(2πm/M ), and a magnetic field gradient suffices to transform between them.

I. Correlation Matrix Reconstruction
The definition of a graph state in Eq. (3) considers an ideal limit of infinite squeezing.In the following, we elaborate on the definition of the adjacency matrix for realistic states with finite squeezing and show that the square graph state generated in our experiment is consistent with this definition.The adjacency matrix that best describes a given state is the one that minimizes Var (p − Ax), which is given by Since A is necessarily symmetric we also have A ij ∝ Cov (x i , p j ).Thus, the correlations between sites in the x and p bases directly reveal the adjacency matrix.
To reconstruct the correlation matrices in Figures 3e  and 4c from measurements of the collective modes, we use a transformation of basis to express the covariance matrix in Eq. (5) as where we sum over the repeated indices m and m ′ .The variances of x and p transform analogously.
In the case of equal couplings to the cavity and equal atom number in each ensemble, the eigenmodes are independent, and Cov (x, x), Cov (p, p), and Cov (x, p) are all diagonal.We use this assumption to extract all relevant information about the state by measuring the covariance matrix for each individual eigenmode.From the variances ζ 2 min,m and ζ 2 max,m in each collective mode and the orientation ϕ m of the squeezed quadrature, we calculate the covariance matrix as where R is a 2 × 2 rotation matrix.

J. Direct Nullifier Measurements
To confirm the efficacy of our graph state generation method, we directly measure the nullifiers n = p − Ax and their variances.This measurement requires simultaneously measuring some sites in the p basis and others in the x basis.To perform this readout for the two-mode graph state, we first apply a variable spinor rotation ϕ to set the measurement basis globally and then apply a 90 • readout rotation about F y only to ensembles 1 and 2. This sequence maps the observable Q L = −x L cos ϕ + p L sin ϕ onto a population difference between Zeeman states.Subsequent evolution under the quadratic Zeeman shift thus only affects the measurement basis in ensembles 3 and 4.After a 90 • rotation about the Q 0 axis, we apply a second Raman rotation to the remaining ensembles to enable readout of the observable P R = x R sin ϕ + p R cos ϕ.The results are shown in Extended Data Fig. 3a, where the red/blue data points represent normalized variances of the sum/difference N ± = Q L ± P R .The nullifiers are given by n To extract the nullifiers for the square graph state, the procedure is the same except that we measure sites 1 and 3 in the x basis (at ϕ = 0) and apply the additional 90 degree rotation about Q 0 to sites 2 and 4. Thus on sites 1 and 3 we read out Q 1,3 = −x 1,3 cos ϕ + p 1,3 sin ϕ and on sites 2 and 4 we read out P 2,4 = x 2,4 sin ϕ + p 2,4 cos ϕ.
To extract the nullifiers, we construct the following four observables,   Extended Data Fig. 1.Imaging calibration.To calibrate the count-to-atom conversion factor r, we measure the fluctuations of the difference in counts cm from atoms in states m = ±1 as a function of the average total counts ⟨c+1 + c−1⟩ from atoms in these two states.The blue dashed line is the polynomial fit of Eq. 9, where the linear coefficient a1 = r + g accounts for the atomic projection noise and a small contribution g ≪ r from photon shot noise.The black dashed line represents the atomic projection noise for r = 395 counts/atom.

Engineering Graph States of Atomic Ensembles by Photon-Mediated Entanglement: Supplementary Information
This supplement provides supporting derivations pertaining to the spin-nematic squeezing dynamics and the verification of entanglement, as well as supporting measurements.Section I presents matrix representations of spin-1 operators and motivates the generalized Bloch sphere for visualization of the spin-1 dynamics.Section II presents a derivation of the metrological squeezing parameter for spin-f atoms.Section III provides an alternative determination of the contrast used in calculating normalized variances and metrological squeezing.In Sec.IV we analytically derive the initial dynamics of spin-nematic squeezing and discuss both technical and fundamental limits on the achievable multimode squeezing.

I. SPIN-1 ALGEBRA AND GENERALIZED BLOCH SPHERE
The matrix forms of the spin-1 dipole operators are From these matrices and the definitions in Methods Sec.A, we can also readily construct the matrix representations of the quadrupole operators q αβ = f α f β + f β f α − 4 3 I 3 δ αβ and q 0 = q zz + 1 3 I 3 , where I 3 is the identity matrix and δ αβ is the Kronecker delta function.For example, We visualize the spin-1 dynamics in an approximate SU(2) subspace spanned by F x , Q yz , and Q 0 , where denote collective observables for the system of N atoms.This visualization is motivated by the commutation relations and Q 0 .For a Larmor-invariant system with only a small side-mode population, ⟨D⟩ = (N +1 + N −1 )/2 ≪ Q 0 .This justifies the use of Q 0 , which is directly accessible from the atomic populations N m , as an approximation for the exact SU(2) subspace {F x , Q yz , Q 0 − D}.Identical dynamics occur in the subspace spanned by F y , Q xz , and Q 0 -where the analogous approximation applies-as our system is Larmor invariant.

II. ENTANGLEMENT DETECTION VIA METROLOGICAL GAIN
Squeezing within a single mode can be used to directly witness entanglement when the generated states circumvent limits on the metrological sensitivity of unentangled states.The appropriate measure of squeezing for this context is the Wineland parameter ξ 2 defined in the main text, which quantifies enhanced sensitivity to rotations [10,36] in our case on the quadrupole spin sphere with axes {F x , Q yz , Q 0 }.Entanglement is necessary for reaching values of this metrological squeezing parameter that are below the standard quantum limit (SQL) ξ 2 = 1.
As the seminal work defining the Wineland squeezing parameter focused on ensembles of spin-1/2 particles [36], we show here how the connection between metrological gain and entanglement generalizes to systems with larger internal spin f , including the spin-1 atoms in this work [54].Generically, the derivation of a Wineland criterion for entanglement follows from considering the metrological task of estimating a parameter λ in a unitary transformation of the form U = e iλH .To derive a criterion suitable for detecting spin-nematic squeezing, we assume without loss of generality that the squeezed quadrature is F x and consider the resulting enhancement in sensitivity to rotations generated by H = Q yz .
We first derive the standard quantum limit on the sensitivity attainable to rotations about Q yz using an unentangled state of N atoms.For any initial state specified by a density matrix ρ, the quantum Cramer-Rao bound limits the precision of estimating the angle λ in a single trial to ∆λ ≥ 1 where is the quantum Fisher information [55,56], and the inequality in Eq. ( S5) is saturated by pure states.For product states, the quantum Fisher information for detecting rotations about Q yz is limited to where in the last inequality the single-atom spin f = 1 dictates that Var (q yz i ) ≤ f 2 .Thus, the minimum imprecision of a measurement of λ using an unentangled state is ∆λ SQL = 1/(2f √ N ).The Wineland parameter compares the sensitivity of an arbitrary state to the standard quantum limit.Assuming that the state is used to estimate λ via the dependence of ⟨F x ⟩ on the rotation about Q yz , the uncertainty ∆λ is related to the spin variance Var (F x ) by where the rate of change d ⟨F x ⟩ /dλ of ⟨F x ⟩ under the Hamiltonian H = Q yz is given by The ratio of the measurement variance (∆λ) 2 to the standard quantum limit (∆λ SQL ) 2 is thus For our atoms with internal spin f = 1, Eq. (S9) simplifies to ξ 2 = Var (F x ) /(C 2 N ).The Wineland parameter is related to the normalized variance ζ 2 plotted in the main text as ξ 2 = ζ 2 /C, so that the SQL is represented by a line at ζ 2 = C.

III. CALIBRATING CONTRAST AND CONFIRMING LARMOR INVARIANCE
To calculate the normalized variances presented in the main text, we determine the commutator CN from the average populations of the Zeeman states after the readout spin rotation, as described in the Methods (Sec.F).The ability to determine both the noise Var (F x ) and the commutator CN = |⟨[F x , Q yz ]⟩| /2 from the same data set is a feature of the Larmor-invariant spin-1 system and is advantageous for minimizing sensitivity to any drifts in experimental parameters.However, an alternative approach is to determine the contrast C by a separate interferometric measurement, as suggested by the role of the metrological squeezing parameter in quantifying enhanced interferometric sensitivity.Here, we present such a measurement of the contrast C Rabi of a Rabi oscillation on the {F x , Q yz , Q 0 } sphere, which corroborates the method used in the main text and confirms the assumption of Larmor invariance.
We determine the contrast C Rabi from the amplitude of a Rabi oscillation after spin-nematic squeezing.For this measurement, we prepare the atoms in m F = 0 before driving the cavity for 50 µs to induce pair-creation dynamics.Up to this point, the sequence is the same as that used for the data in Fig. 2 of the main text.After the dynamics, we apply a spin rotation by an angle θ = Ω Raman t.Assuming Larmor invariance of the state, this operation is equivalent to a rotation around F x by an angle 2θ on the {F x , Q yz , Q 0 } sphere (Fig. S1a).We then evaluate the normalized population imbalance ))/N as function of θ.The full expression for the resulting Rabi oscillation in terms of the initial expectation values Q 0 , ⟨Q yz ⟩, and

Cavity Drive Squeezing
Raman Rotation .Direct measurement of contrast.a, After applying the cavity drive field for 50 µs to generate spin-nematic squeezing, we perform a rotation about an arbitrary axis with Larmor phase φL in F x − F y plane.When φL = 0 this rotation corresponds to Rabi oscillations that are visualized on the {F x , Q yz , Q 0 } sphere.b, Following the rotation, we measure the normalized imbalance ))/N as a function of Rabi oscillation time t.The dashed line is a fit of the model in Eq. (S12).From this fit, we extract the contrast The final term arises from the fact that {F x , Q yz , Q 0 } form an approximate rather than exact SU(2) subspace (see Sec. I) and is small, ⟨D⟩ ≈ (N +1 +N −1 )/2 ≪ Q 0 , for the early-time squeezing dynamics in our system.The contrast of the Rabi oscillation can be expressed in terms of the operators prior to the spin rotation as which provides an alternative measurement of the contrast C as defined in the main text.
To extract the contrast C Rabi , we fit measurements of the population imbalance Q 0 (t) after a Raman pulse of duration t with a function of the form where Ω Raman is the Rabi frequency and γ is the Rabi dephasing rate.The amplitude A and offset d together account for both the final term in Eq. (S10) and incoherent collisional spin-exchange interactions between the readout rotation and imaging, where the latter is dominant in our system.We use the fit to extract the contrast The primary source of imperfect contrast is spin-changing collisions that occur in the time (20 ms) between the Raman rotation and the fluorescence readout.These collisions permit states polarized with either all atoms in m = 0 or all atoms in m = ±1 to decay into a mixture of all three states.This process reduces contrast and introduces an offset to the Rabi oscillation, but does not directly affect the measurement of Var (F x ), since the collisions preserve the population difference N ′ +1 − N ′ −1 from which we infer the polarization F x in the experiment.Additionally, dephasing from inhomogeneity in the Raman rotation and detuning due to a finite ratio Ω Raman /q reduce the measured value of Q 0 θ=90 • .The measurement of contrast C used for evaluating normalized variances depends only on Q 0 θ=90 • and not on Q 0 θ=0 • , and is thus more strongly impacted by these effects, leading to a 3% difference between the two measurements as shown in Fig. S1c.This difference indicates that the normalized variances and squeezing parameters reported in the main text are conservative estimates.
In calibrating the contrast C via Rabi oscillations, we perform the Raman rotation about an arbitrary axis in the F x − F y plane, as there is no prior phase reference in the experiment.We thus implicitly assume Larmor invariance of the generated state.We confirm this assumption by performing repeated measurements of the population imbalance at fixed rotation angle θ while varying the Larmor phase φ L .We perform these measurements with a readout rotation of θ = 90 • where the measured imbalance is most sensitive to changes in the Larmor phase, and we observe no dependence on phase over a full Larmor oscillation period (see Fig. S1d).

IV. SQUEEZING DYNAMICS
We present a theoretical formulation of the spin-nematic squeezing dynamics.This model provides a basis to estimate the contributions of technical noise in our measured squeezing, as well as the effects of cavity dissipation.Finally, we discuss fundamental limits set by the cavity cooperativity on the degree of multimode squeezing attainable in our protocol for preparing graph states.

A. Equations of Motion
We analyze the dynamics of spin-nematic squeezing for a system initialized with all atoms in m = 0, focusing on the experimentally relevant regime of early times where the population in states m = ±1 remains small (N +1 +N −1 ≪ N ).To allow for effects of finite atomic temperature, we begin by writing down a Hamiltonian that incorporates nonuniform coupling to the cavity mode, Here the collective spin F defined in the main text is replaced with a weighted collective spin F = i w i f i , which includes a correction for inhomogeneous cavity couplings w i ∝ Ω i , where Ω i denotes ac Stark shift per intracavity photon experienced by the i th atom.The weights w i are normalized such that ⟨w i ⟩ = 1.We describe the early-time dynamics in the two dimensional subspace spanned by the weighted spin operators F x and Q yz .During these early times, commutators relevant to the dynamics are Identical dynamics occur in the subspace spanned by F y and Q xz so that state remains invariant under global spin rotations about F z .The linear equations of motion for F x and Q yz can be solved exactly.This system has eigenvalues ±λ where λ = −q(q + 2χ).The corresponding solutions are In general, these two operators are not orthogonal unless χ = −q.The expectation value and variance of any observable F ϕ = cos ϕ F x − sin ϕ Q yz can be calculated from these operators by noting that F ϕ (t) = aF −λ (t) + bF +λ (t) where a and b are real coefficients that are independent of time and may be solved for from the expressions at time t = 0.The variance for a particular spinor angle ϕ at time t is then given by where the cross term proportional to ab again highlights that F −λ (t) and F +λ (t) are in general not orthogonal.At t = 0 the system is in a coherent state with variance F ϕ (0) 2 = N at the level of projection-noise for all values of ϕ.This condition yields the constraint which reduces to a 2 + b 2 = 1 when χ = −q.

B. Technical Limitations on Squeezing
The model of the dynamics in Sec.IV A provides a foundation for estimating the effect of experimental noise on the amount of observed squeezing.Two primary limits are fluctuations in the collective interaction strength χ and inhomogeneous coupling of the thermal distribution of atoms to the cavity, where the latter effect introduces a slight discrepancy between the collective mode squeezed by the cavity and the observable detected in fluorescence imaging.

Fluctuations in Interaction Strength
The collective interaction strength χ varies the angle of squeezing as ϕ min = arctan(−q/λ), and so fluctuations in χ act to wash out the squeezing.We can write the variance as a function of spinor phase as To find the effect of fluctuations ∆χ in interaction strength, we expand Eq. (S18) around the angle ϕ min of optimum squeezing and write our expression in terms of χ.To leading order in ∆χ, we are left with the additional noise ∆ζ 2 interaction from interaction strength fluctuations, which is given by (S19) This is a lower bound on the added noise, since the angle ϕ min at finite times has larger fluctuations than in the t > λ −1 limit.In principle, added noise from interaction strength fluctuations can be suppressed by working in the regime of χ ≫ q.However, we do not operate in this limit because setting χ ∼ q maximizes the squeezing rate λ relative to the fundamental cavity dissipation, as we shall see in Sec.IV C. The fluctuations in interaction strength χ in our experiment arise from variations in the number of intracavity photons n, the number of coupled atoms N , or the detunings from the two virtual Raman processes δ ± (see Eq. ( 8)).These sources of noise are correlated, as the number of intracavity photons not only depends on the input drive strength n i , but also depends on the detuning from cavity resonance.The detunings δ c and δ ± in turn depend on the atom number N due to the dispersive shift δ N = 4ΩN of the cavity resonance induced by the atoms.The two direct sources of fluctuations in χ are then the number of input photons n i and the atom number N , which lead to total fluctuations ∆χ χ where evaluates to α ≈ 2 for our parameters.We stabilize the drive input power to ensure ∆n i /n i < 5%, and we reduce fluctuations in atom number by post-selecting so that ∆N/N < 5% within each data set.At ∆χ/χ < 10%, and a typical value of ζ 2 max = 10, using the values from Methods Sec.C, the additional noise from interaction strength fluctuations is ∆ζ 2 interaction ≈ 0.02.

Inhomogeneous Atom-Cavity Coupling
Our fluorescence imaging measures a uniformly weighted collective spin F x , while the cavity couples to the inhomogeneously weighted collective spin F x defined in Sec.IV A. Any width in the distribution of coupling weights w i manifests itself in reduced squeezing.Without loss of generality, we assume F x is the squeezed observable and compute the projection of the measured observable on the squeezed observable: The excess noise is given by the magnitude of the remaining component, The variance in couplings comes primarily from the thermal distribution of the atomic states.We parameterize the temperature by the ratio β = U 0 /(k B T ) of the lattice depth to the atomic temperature.Assuming a harmonic trap, the excess noise Eq. (S23) for a single lattice site is given by ∆ζ In the low temperature limit β → ∞, to leading order the added noise is

C. Cavity Dissipation
The unitary dynamics described in Methods Sec.IV A are modified by decay channels inherent to any real cavity system.As described in Sec.C and Refs.[20,33], spin-exchange interactions in the cavity are mediated by a virtual process in which atoms collectively scatter photons from a vertically polarized drive field into a horizontally polarized cavity mode.For coherent interactions, these horizontally polarized photons are subsequently scattered back into the vertical drive mode, allowing for unitary transfer of information among the atoms.However, in practice, photons may also be lost before completing the unitary dynamics and thereby carry away quantum information.
A photon may be lost either due to the finite cavity lifetime or by atomic scattering into free space.In the case of cavity decay, the loss of a photon is accompanied by creation or annihilation of a collective spin excitation.This decay channel is described by the Lindblad operators L ± = √ γ ± F ± and has a characteristic strength Γ coll = 2N (γ + + γ − ) that, in analogy to the collective interaction strength χ, is enhanced by the number of atoms.The ratio of the collective decay to the collective interaction strength is determined by the detuning δ − from the dominant Raman process in our experiment.To quantify the impact of the collective decay process on the squeezing, we write down the Lindblad equation of motion for the squeezed quadrature with the two loss operators L ± = √ γ ± F ± .Under the simplifying assumptions of uniform couplings (F ± = F ± ), χ ≈ −q and sufficiently early-time dynamics, the equation of motion for squeezing is The steady state of this equation leads to a bound on the variance due to collective decay of The variance exponentially decays to this bound at a rate proportional to exp(−2λ).At finite times, the effect of the bound is mathematically equivalent to mixing the ideal squeezed state achieved under unitary dynamics with vacuum fluctuations on a beam splitter with transmission 1 − ∆ζ 2 coll .Additionally, photons may be lost due to free-space scattering at a rate Γ sc per atom.Free-space scattering is not a collective process, as the scattered photons carry away information about individual atoms.On cavity resonance, free-space scattering is thus suppressed with respect to interactions by a factor N η/k, where N η is the collective cooperativity for a cycling transition and the numerical factor k = 96 includes the strengths of the atomic transitions in our level scheme [20].Overall, the rate of free-space scattering in the limit δ − > κ is then The effect of a scattering event is to erase correlations between the atom that scatters a photon and the remaining atoms.This erasure of correlations adds noise to the squeezed quadrature at a rate where α = 2 for the worst case where an atom is projected into the m = 0 state.However, while collective decay only impacts the mode coupled to the cavity, free-space scattering continues to impact all modes that have already been  squeezed.On average each mode is impacted by scattering for a total duration τ × M/2, where τ is the duration of each squeezing pulse, yielding a noise contribution While free-space scattering additionally reduces the contrast by a factor of approximately e −M τ Γsc/2 , in the regime of strong squeezing the added noise is the dominant limitation.The impact of the two noise contributions ∆ζ 2 coll and ∆ζ 2 sc can be minimized by choosing optimal values of the interaction strength χ and detuning δ − .An interaction strength of χ = −q optimizes the speed of the coherent dynamics λ/χ.Optimizing the two limits for the mode with maximal scattering, with τ λ = 1, yields a detuning of δ − = κ N η/(192M ), balancing the impact of the two loss channels to minimize their combined effect.In practice, we experimentally optimize squeezing, which takes into account additional noise sources, resulting in a slight deviation from the theoretical optimum.
Having derived expressions for both collective decay and free-space scattering, we summarize the impact on the experiments in the current work in Sec.IV D and derive fundamental limits on the scaling of the squeezing in Sec.IV E.

D. Summary of Noise Contributions
We summarize the impact of all noise processes limiting our squeezing in Table S1.The effects of cavity decay, freespace scattering, and coupling variation are all mathematically equivalent to mixing the squeezed quantum state with zero point fluctuations, as if on a beam splitter.Starting from the minimal possible variance given unitary dynamics, 1/ζ 2 max , each process results in a factor of (1 − ∆ζ 2 i ) reduction in the amount by which the state is squeezed.We calculate the combined effect of these noise sources ∆ζ 2  i , assuming that they are all independent, as ∆ζ Photon shot noise and interaction strength noise behave differently, since they directly add noise rather than degrading the state toward the standard quantum limit.These terms are added at the end using standard propagation of uncertainty.Finally, we divide by the Rabi oscillation contrast C Rabi = 0.89 measured in Sec.III to obtain the expected normalized variance ζ 2 .In principle, working at larger atom number increases the collective cooperativity, decreasing the relative effects of cavity dissipation.However, in addition to the noise sources in Table S1, different collective modes are sensitive to technical noise in the readout procedure, as discussed in Sec.E. The example data presented in Table S1 are measured in the ↑↓↑↓ mode, which has negligible technical noise (see Extended Data Fig. 1).For Figs. 2 and 3 of the main text we squeeze up to two collective modes, and these modes can always be mapped to the two collective modes with minimal technical noise using local Larmor rotations.In this case we choose an atom number of N = 1.5 × 10 4 , which is limited by the density of atoms in the trap, as any collisional spin exchange interactions are incoherent with the photon-mediated interactions, reducing the effective spin length.For the square graph state all 4 collective modes need to be squeezed, and the technical noise in our projection noise calibration becomes a relevant parameter, so we reduce the atom number in Fig. 4 of the main text to N = 8 × 10 3 .

E. Fundamental Scaling
Fundamental limits to the degree of squeezing attainable by global cavity-mediated interactions among N atoms are governed by the collective cooperativity N η [57], where η = 4g 2 /(κΓ) is the single-atom cooperativity.In this section, we derive limits on squeezing multiple collective modes which demonstrate the scalability of our graph-state preparation protocol.We focus first on the specific case of the spin-nematic squeezing employed in this work, and additionally comment on generalizations to other methods of cavity-mediated spin squeezing.
The fundamental limit on squeezing set by the cavity cooperativity arises from two dissipation processes: loss of photons from the cavity mode that mediates interactions; and scattering of photons into free space.These loss processes are parameterized by the rates Γ coll and Γ sc in Eq. (S30).The effect of the free-space scattering is proportional to the number of modes M because scattering at any point during the protocol can reduce squeezing.Conversely, collective decay does not depend on M because it acts only on the collective mode coupled to the cavity.To place free-space scattering and cavity loss on an even footing, we imagine dividing the squeezing of each collective mode into multiple segments interleaved with squeezing of the other collective modes, so that the scattering is interspersed with the coherent dynamics.Any scattering loss while addressing each of the M modes should then be included as a component of Eq. (S28).The full equation of motion for the variance of each collective mode is thus The squeezing is optimized by choosing the drive-cavity detuning that minimizes the total contribution to Eq. (S34) from scattering and cavity decay: This optimum detuning is set by the collective cooperativity per mode N η/M and leads to an overall variance While these derivations focus on variance ζ 2 , the same scaling applies to the Wineland squeezing parameter ξ 2 .Equation (S36) shows that our protocol for generating arbitrary M -node graph states by squeezing M collective modes yields fixed squeezing at constant atom number per ensemble (N/M ), independent of the number of graph nodes.The protocol requires order M squeezed modes and order M sets of rotations.However, at fixed atom number per ensemble, with increasing M the increased collective interaction strength causes each mode to be squeezed faster, so that the total interaction time remains fixed.The protocol can thus be scaled to larger arrays of ensembles, limited only by the spatial extent of the cavity mode and the fidelity of local addressing.
The limit on squeezing ξ 2 ∝ 1/ (N/M )η set by the collective cooperativity per mode generalizes to a wide variety of methods of cavity spin squeezing, including approaches employing either photon-mediated interactions [57,58] or quantum non-demolition measurements [8].Improvements to both the numerical prefactor and the overall scaling with cooperativity are possible, however, by a suitable choice of atomic level scheme.Notably, for squeezing on a cycling transition, the scaling for the single-mode case improves to ξ 2 ∝ 1/(N η) [7,59].Optimizing the scheme for the collective entangling operations may facilitate future work seeking the error-correction threshold of −10 log ξ 2 = 20.5 dB [44], or generating discrete-variable graph states in arrays of single atoms.

Fig. 1 .
Fig. 1.Programmable entanglement in an array of four atomic ensembles within an optical cavity.a, Initializing all atoms in the m = 0 state and driving the cavity with light induces creation of correlated atom pairs in states m = ±1.b,The resulting spin-nematic squeezing is visualized on a spherical phase space spanned by the collective spin-1 observables {F x , Q yz , Q 0 }.For short interaction times, the dynamics can be described on an effective two-dimensional phase-space spanned by the conjugate observables {x, p}.c, Combining the global interactions with local spin rotations allows for engineering a variety of entanglement structures, such as entanglement localized to selected subsystems and graph states with up to four nodes.

Fig. 2 .
Fig.2.Global squeezing and entanglement between subsystems.a, Cavity-mediated interactions lead to squeezing of the symmetric mode (red circles) below the standard quantum limit (SQL, dashed line).The anti-symmetric mode (blue squares) does not couple to the cavity and remains approximately in a coherent state.Multiplying values of the variance ζ 2 for the squeezed quadrature x ′ + of the symmetric mode and the orthogonal quadrature p ′ − of the antisymmetric mode yields the entanglement witness W . Inset: green ellipse shows area √ W , smaller than dashed circular region representing minimum-uncertainty unentangled state.b, Analyzing the left and right subsystems separately (yellow squares and purple circles) yields a degradation in squeezing, consistent with neglecting information contained in correlations between the subsystems.Error bars show 1 s.d.confidence intervals extracted via jackknife resampling.Shaded curves show the 1 s.d.confidence intervals of sinusoidal fits to the data.

Fig. 4 .
Fig. 4. Generation of a square graph state.a, Diagram of four-mode square graph state and theoretical correlation matrix Corr (xi, pj) ∝ Aij.b, Left: schematic illustration of eigenmodes of the adjacency matrix A and the corresponding squeezing ellipses, with orientations specified by eigenvalues λm = cot ϕm.Right: measured variances ζ 2 in the four eigenmodes, showing squeezing at the specified spinor phases ϕm (black dashed lines).Error bars show 1 s.d.confidence interval.c, Top: directly measured variances vi of the nullifiers, with schematics showing central node i (dark gray circle) and neighbors (black circles) contributing to each nullifier.Bottom: correlation matrix reconstructed from the measurement results in b.

such that n 1 , 3 =
N 1,3 (90 • ) and n 2,4 = N 2,4 (0 • ).The measured normalized variances as a function of the initial spinor rotation are shown in Extended Data Fig. 3b.The nullifier variances reported in the main text are obtained from the data in Figs. 3 by subtracting a small detection noise contribution, as described in Sec.E.

. 2 .
Sample sequence for generating the 4-mode square graph state by squeezing collective modes.Bottom four rows show the state of each eigenmode throughout the entire pulse sequence.The spinor angles Φ1,2,3 = (0, 117 • , −54 • ) are chosen such that, at the end of the sequence, each eigenmode is squeezed along the axis specified by the corresponding eigenvalue of the adjacency matrix.

. 3 .
Direct measurements of the nullifiers, as described in Methods Sec.J. a, Nullifiers for the two-mode EPR state.The spinor phase ϕ gives the basis in which the left ensemble pair is measured, while the right ensemble pair experiences an additional 90 • of spinor evolution.The nullifier for the left subsystem nL = pL − xR is extracted from N− (blue) at ϕ = 90 • .The nullifier for the right subsystem nR = pR − xL is extracted from N+ (red) at ϕ = 0 • .b, Nullifiers for square graph state.We extract the nullifier variances shown in Fig.4cof the main text from N1,3 (blue, purple) at ϕ = 90 • and N2,4 (red, yellow) at ϕ = 0 • .

)
Fig.S1.Direct measurement of contrast.a, After applying the cavity drive field for 50 µs to generate spin-nematic squeezing, we perform a rotation about an arbitrary axis with Larmor phase φL in F x − F y plane.When φL = 0 this rotation corresponds to Rabi oscillations that are visualized on the {F x , Q yz , Q 0 } sphere.b, Following the rotation, we measure the normalized imbalance−Q 0 θ /N = (N ′ 0 − (N ′ +1 + N ′ −1))/N as a function of Rabi oscillation time t.The dashed line is a fit of the model in Eq. (S12).From this fit, we extract the contrastC Rabi = Q 0 θ=0 • − Q 0 θ=90 • /(2N ) =0.89 ± 0.01.c, We compare the contrast C Rabi to the value C = (3Q 0 θ=90 • − 1)/(2N ) computed as derived in Methods Sec.F. The percent-level difference indicated by the average ratio C Rabi /C = 1.03 ± 0.01 (black dashed line) is attributable to the dephasing of the Rabi oscillation at rate γ, which results in conservative estimates of normalized variances and squeezing parameters in the main text.d, Measurements of population imbalance after a rotation by θ = 90 • confirm that the calibration of contrast is independent of the choice of Larmor phase φL.Data shown are the average over typically 50 iterations.
Fig.S1.Direct measurement of contrast.a, After applying the cavity drive field for 50 µs to generate spin-nematic squeezing, we perform a rotation about an arbitrary axis with Larmor phase φL in F x − F y plane.When φL = 0 this rotation corresponds to Rabi oscillations that are visualized on the {F x , Q yz , Q 0 } sphere.b, Following the rotation, we measure the normalized imbalance−Q 0 θ /N = (N ′ 0 − (N ′ +1 + N ′ −1))/N as a function of Rabi oscillation time t.The dashed line is a fit of the model in Eq. (S12).From this fit, we extract the contrastC Rabi = Q 0 θ=0 • − Q 0 θ=90 • /(2N ) =0.89 ± 0.01.c, We compare the contrast C Rabi to the value C = (3Q 0 θ=90 • − 1)/(2N ) computed as derived in Methods Sec.F. The percent-level difference indicated by the average ratio C Rabi /C = 1.03 ± 0.01 (black dashed line) is attributable to the dephasing of the Rabi oscillation at rate γ, which results in conservative estimates of normalized variances and squeezing parameters in the main text.d, Measurements of population imbalance after a rotation by θ = 90 • confirm that the calibration of contrast is independent of the choice of Larmor phase φL.Data shown are the average over typically 50 iterations.
and [Q 0 , F x ] = 2iQ yz .The Heisenberg equations dO/dt = i[H, O] for both spin observables in this space are d dt near cavity center with an inverse temperature of β = 15, Eq.(S24) limits squeezing to ζ 2 > 0.08.We directly measure the distribution of couplings w i via microwave spectroscopy, probing the ac Stark shift induced by the drive field on the hyperfine clock transition |f = 1, m f = 0⟩ → |2, 0⟩.The drive light, detuned from atomic resonance by ∆ = −2π × 9.5 GHz, induces a differential ac Stark shift that is directly proportional to the weight w i for each atom.We measure the distribution of Stark shifts at different drive intensities, as shown in Fig. S2.The measured spectra are well fit by a model of a thermal distribution with inverse temperature β = 15, exhibiting variances and means that directly corroborate the bound established in Eqs.(S23)-(S24).
Fig. S2.Measurement of the drive-induced ac Stark shift on the |f = 1, m f = 0⟩ → |2, 0⟩ transition.Red, blue, and gold indicate increasing intracavity photon number.The measured distributions of Stark shifts (circles) are well described by a numerical fit based on a thermal distribution of atomic positions with a ratio β = 15 of trap depth to temperature.

Table 1 .
Summary of experimental parameters for cavity mediated interactions: detuning δ− of cavity drive field from two-photon resonance, total atom number N , collective interaction strength χ and interaction time τ .

Table 2 .
Summary of EPR steering values.To measure EPR steering between different subsystems, we need to infer the value of the left subsystem in the x ′ and p ′ quadratures from measurements of the right subsystem, and vice versa.Variances representing the error of each inference, and the resulting steering witnesses, are presented.

Table S1 .
Summary of noise sources contributing to the variances in Figs.2 and 3of the main text, along with the relevant experimental parameters from Extended Data Table1.The method of calculating the expected variance from the individual noise contributions is described in Sec.IV D.