Reservoir-engineered spin squeezing: macroscopic even-odd effects and hybrid-systems implementations

We revisit the dissipative approach to producing and stabilizing spin-squeezed states of an ensemble of $N$ two-level systems, providing a detailed analysis of two surprising yet generic features of such protocols. The first is a macroscopic sensitivity of the steady state to whether $N$ is even or odd. We discuss how this effect can be avoided (if the goal is parity-insensitive squeezing), or could be exploited as a new kind of sensing modality to detect the addition or removal of a single spin. The second effect is an anomalous emergent long timescale and a"prethermalized"regime that occurs for even weak single-spin dephasing. This effect allows one to have strong spin squeezing over a long transient time even though the level of spin squeezing in the steady state is very small. We also discuss a general hybrid-systems approach for implementing dissipative spin squeezing that does not require squeezed input light or complex multi-level atoms, but instead makes use of bosonic reservoir-engineering ideas. Our protocol is compatible with a variety of platforms, including trapped ions, NV defect spins coupled to diamond optomechanical crystals, and spin ensembles coupled to superconducting microwave circuits.


I. INTRODUCTION
Among the most sought-after states in quantum metrology are spin-squeezed states, highly entangled states of spin-1/2 ensembles that enable parameter sensing with a sensitivity better than the standard quantum limit, even reaching fundamental Heisenberg-limit scaling [1,2]. The standard approach for producing these states is to unitarily evolve an initial product state under a collective spin-spin interaction Hamiltonian. While many interactions are possible, the most widely studied one is the one-axis twist (OAT) Hamiltonian [1], which has been realized in a number of ground-breaking experiments [3][4][5][6]. It unfortunately is not capable of achieving Heisenberg-limited squeezing even in the ideal case [2]. An alternate, more complex interaction Hamiltonian is the two-axis twist (TAT) Hamiltonian [1,[7][8][9][10], which, while more resource intensive, allows achieving Heisenberg-limited scaling.
While easy to understand, tailored unitaryevolution is not the only approach to spin squeezing. An alternative is to use the general strategy of reservoir engineering [11], where tailored dissipation is exploited to both produce and stabilize a nontrivial state of interest, i.e., a spin-squeezed state (see Fig. 1). The dissipative approach in principle has several advantages: the spin-squeezed state is stabilized in the steady state (as opposed to just prepared at a specific instant of time), the stabilization is largely insensitive to the initial state of the ensemble, and one can achieve Heisenberg-limited scaling.
FIG. 1. Schematic representation of a generic approach to generate dissipative spin squeezing by coupling spins to a bosonic mode that interacts with a squeezed reservoir. The squeezing rate experienced by the cavity is governed by the parameter κsqz, while g represents the spin-cavity coupling strength. Limiting factors to the protocol's performance are the intrinsic photon-loss rate κint, the local spin-relaxation rate γ rel , and the local spin-dephasing rate γ φ .
The dissipative stabilization of bosonic squeezed states has been studied extensively both theoretically [12][13][14] and experimentally [15][16][17][18][19][20]. Corresponding schemes for spin squeezing have also been studied theoretically. The earliest works analyzed schemes where atoms are directly illuminated with squeezed light. Both the cases of two-level atoms [21][22][23] and V -type multilevel atoms [24] were studied. More recently, it was shown theoretically that the same effective dissipative dynamics could be realized by using Raman processes in driven multilevel atoms coupled to a lossy cavity [9,25].
In this work, we revisit the dissipative approach to spin squeezing. Our work complements previous studies both by discussing a powerful, alternative method for implementing these schemes, as well as describing surprising phenomena that had not been fully analyzed in the past. In terms of implementation, we analyze a very general hybridsystems approach that harnesses bosonic dissipative squeezing. We consider a spin ensemble which is resonantly coupled to a cavity mode (via a standard Tavis-Cummings [26] interaction), which is in turn coupled to an effective squeezed reservoir (see Fig. 1). Previous proposals [22,24] suggested to implement this squeezed reservoir by driving the cavity with squeezed light, an approach which is limited by losses associated with the transport and injection of an externally prepared optical squeezed state. We show that there are also simpler methods to generate the effective squeezed reservoir, which can be implemented using only classical optical or microwaves drives by harnessing existing dissipative bosonic squeezing schemes. Specifically, we consider coupling the cavity to an auxiliary lossy degree of freedom (two level system or bosonic mode) which is driven simultaneously with imbalanced red-detuned and blue-detuned sideband drives. Such schemes produce an effective squeezed dissipator for the cavity and have been experimentally implemented in wide variety of platforms, including optomechanics [15], trapped ions [16] and superconducting circuits [20]. Since only classical radiation is required, this approach is insensitive to the aforementioned transport and insertion losses of squeezed radiation. We demonstrate that this hybrid-systems approach to dissipative spin squeezing can reach the Heisenberg limit, and also outperform OAT in the presence of single-spin T 1 decay. Note that unlike the Raman scheme of Ref. 25, which requires atoms with a specific four-level configuration, the approach here only requires standard two-level atoms, making it compatible with a wide variety of systems (including possibly solid-state systems such as ensembles of NV defect spins in diamond [27]).
Our work also analyzes surprising phenomena that were not fully discussed previously. Perhaps most striking is the extreme sensitivity of dissipative spin squeezing to the parity of the total number of spins N : the steady state is macroscopically different for N spins versus N + 1 spins. While this effect was implicitly contained in the results of Agarwal and Puri [21,22] (see Sec. VIII for a detailed discussion of the relation to previous works), we provide here a fully qualitative and quantitative analysis. We discuss how this effect can be avoided (if one wants strong spin squeezing independent of parity), and how it could also be exploited as a new kind of sensing modality. We also make a surprising connection to a non-dissipative many-body system, the antiferromagnetic Lipkin-Meshkov-Glick (LMG) model [28,29].
A second surprising and new phenomenon we describe is the interplay between collective dissipative-spin-squeezing dynamics and noncollective singlespin dephasing. As we show, this results in an extremely long relaxation timescale in the system (i.e., inverse dissipative gap) which grows with system size N . At a fundamental level, the effect has parallels to prethermalization behavior observed in weakly nonintegrable systems (see, e.g., [30,31]). At a practical level, we show that even infinitesimally weak single-spin dephasing dramatically impairs the steady-state spin squeezing to at most −3 dB. However, we also show that this need not be a limitation: large amounts of squeezing are possible in the prethermalized regime, i.e., at transient times parametrically shorter than the timescale required to reach the steady state, or by deliberately adding very small levels of single-spin relaxation. This effect may allow one to quickly generate strong squeezing even in parameter regimes that had previously been discarded based on the low level of steady-state squeezing.
Finally, we analyze the impact of imperfections in the reservoir-engineering process that lead to the engineered squeezed dissipation having a non-zero impurity and effective thermal occupancy. We reveal a striking sensitivity of spin squeezing to such imperfections if the squeezing strength is made too large.
The remainder of this paper is organized as follows: In Sec. II, we outline the key idea behind the standard approach to dissipative spin squeezing as well as summarize our generic protocol. In Sec. III, we explore the even-odd effect and briefly discuss connections to sensing. In Sec. IV, we carefully analyze the performance of our dissipative spin squeezing protocol in the presence of single-spin dissipation, showing that the steady-state squeezing it generates can outperform the transient squeezing produced by standard OAT. In Sec. V, we discuss the emergence of anomalously slow relaxation times and we introduce a dynamical-decoupling protocol to cancel the effect of inhomogeneous broadening. In Sec. VI, we analyze imperfections in the reservoir-engineering process, while in Sec. VII, we discuss in more detail how our protocol could be implemented in a variety of different physical systems. In Sec. VIII, we review previous works on dissipative spin squeezing and discuss their relation to our new findings. Conclusions and a summary are presented in Sec. IX.

II. MODEL AND THE BASIC DISSIPATIVE SQUEEZING PROTOCOL
The reservoir engineering approach to spin squeezing requires one to construct a nontrivial dissipative environment for the spins. In this section, we review the idealized spin-only quantum master equation that describes the needed dissipative dynamics [21,22]. We then present a more realistic model that corresponds to the generic, experimentally-friendly hybrid-systems setup sketched in Fig. 1, where a spin ensemble is coupled to a cavity (or other bosonic mode), which is in turn coupled to an engineered squeezed reservoir.
Throughout this paper, we quantify the amount of metrologically-useful spin squeezing (i.e., as relevant to a standard Ramsey measurement) using the Wineland parameter [2,32]. It is defined as where ∆Ŝ 2 ⊥ is the minimum variance in a direction perpendicular to the direction of the mean of the collective spin andŜ ≡ (Ŝ x ,Ŝ y ,Ŝ z ) is the vector of spin operators.

A. Idealized spin-only model
We consider the following quantum master equation acting on the Hilbert space of N spin-1/2 particles,ρ where we introduced the operator Σ[r] = cosh(r)Ŝ − − sinh(r)Ŝ + .
Here, Γ is the coupling rate to the engineered reservoir, r characterizes the squeezing strength, and D[ẑ]ρ =ẑρẑ † − {ẑ †ẑ ,ρ}/2 is the standard Lindblad dissipative superoperator. We also introduced the collective spin operatorsŜ ± =Ŝ x ± iŜ y witĥ denotes a standard Pauli matrix acting on the jth spin. Here,Σ[r] is analogous to a standard bosonic Bogoliubov annihilation operator, where bosonic raising and lowering operators have been replaced byŜ + andŜ − respectively. Similar to reservoir-engineered bosonic squeezing [13], the desired squeezed state will correspond to the vacuum of this operator and the squeezing parameter r characterizes the amount of squeezing, e −2r , of the vacuum fluctuations.
To be more explicit, Refs. 21 and 22 showed that, for even N , Eq. (2) has pure steady states that correspond to zero-eigenvalue eigenstates (i.e., "dark states") ofΣ[r],Σ Since Eq. (2) conserves the total angular momentum j, there is a dark state for each allowed value of j. Each |ψ dk [j; r] has a mean spin polarization in the z direction, and exhibits squeezing (anti-squeezing) ofŜ y (Ŝ x ). The choice of the squeezing axis is determined by the relative phase between theŜ + and S − terms in Eq. (3), which is chosen here to be −1.
If the system is initialized in an arbitrary state with a definite value of j, the dissipative dynamics will relax the system to a dark state in this subspace. For states in the maximum-angular-momentum subspace j = j max = N/2, the relaxation timescale (i.e., the inverse dissipative gap of the Liouvillian) is ∝ 1/N Γ, see Sec. III. Note that the dark states with j < j max are not unique, since the corresponding angular-momentum subspaces are degenerate [33]. However, if the initial state and the dynamics are invariant under permutation of spins, the system will only explore permutationally invariant states [34], and there is a unique dark state for each j subspace, see App. E 1.
As detailed in App. D, the dark states can be expressed in the form [22,23] where |j, m y denotes an eigenstate ofŜ 2 andŜ y , N (r) is a normalization constant, and we defined θ = ln tanh(r). In terms of the eigenstates |j, m ofŜ 2 andŜ z , these states read as follows.
where every second coefficient is nonzero, for k ∈ {0, ..., j}, and all other coefficients vanish, c −j (r) serves as a normalization constant. The parameter r controls the amount of squeezing in the steady state. If we initialize the system in an arbitrary state with j = N/2, the resulting pure steady state is squeezed, with ξ 2 R → 2/(N + 2) in the large-r limit. This corresponds to Heisenberglimited spin squeezing, and thus outperforms both the standard quantum limit (i.e., ξ 2 R = 1) as well as the maximum squeezing possible with an ideal OAT interaction (ξ 2 R ∝ 1/N 2/3 ). Note that a standard leading-order Holstein-Primakoff approximation could be used to map Eq. (2) to a bosonic squeezing dissipator; however, this would not let one understand the ultimate saturation of squeezing (with increasing r) to the Heisenberg-limited value.
B. Hybrid-systems approach to dissipative spin squeezing As noted in the introduction, previous studies have analyzed methods for realizing the dissipative dynamics in Eq. (2). These methods either required direct driving of spins with squeezed light [22,24] (which is experimentally challenging), or the use of Raman processes in structured four-level atoms [9,25] (which is not applicable to generic two-level systems). We present here an alternate, generic method that takes a hybrid-systems approach: a cavity (or other bosonic mode) is coupled both to an ensemble of two-level systems, as well as to an engineered, bosonic squeezed reservoir (see Fig. 1). As discussed, such a bosonic squeezed reservoir can be realized using only classical driving fields, and has been implemented in a variety of different experiments [15,16,20]. We discuss specific implementation strategies of this general approach in Sec. VII; here, we present the general structure of the overall quantum master equation.
To this end, we consider a spin ensemble that is resonantly coupled to a bosonic mode (with lowering operatorâ). In the rotating frame, the Hamiltonian isĤ where g is the spin-cavity coupling strength. We further assume that this mode is coupled both to an engineered squeezed reservoir (with coupling rate κ sqz and squeezing parameter r) as well as subject to unwanted zero-temperature loss (at rate κ int ). The quantum master equation is theṅ We have also included standard single-spin decay and dephasing dissipators (at rates γ rel and γ φ , respectively).
At a heuristic level, the cavity serves as a transducer that allows the spins to inherit the squeezed fluctuations produced by the bosonic squeezed reservoir. As the squeezed reservoir is engineered, we will treat r and κ sqz as tuneable parameters that can be optimized. In contrast, we will take the coupling g and the unwanted dissipation (i.e., κ int , γ φ , and γ rel ) to be fixed. This then motivates introducing singlespin cooperativities η k and collective cooperativities C k via: where k ∈ {φ, rel}. The goal will be to understand the optimal squeezing possible for a fixed value of C k . As we will show in Sec. IV, in the case where single-spin relaxation dominates over dephasing, the optimized dissipative scheme achieves steady-state squeezing scaling as ξ 2 R ∝ 1/ √ C rel . This is significantly better than the optimized transient OAT squeezing in this regime, which only scales as [35]. To connect our setup to the simpler quantum master equation (2), we consider the regime where the condition √ N g κ int +κ sqz holds, and we adiabatically eliminate the cavityâ. We obtain (see App. A) where we have defined and We see that the internal loss of the cavity results in a collective relaxation process for the spin ensemble; this is similar to OAT-based protocols that are derived using a strongly detuned cavity-spin ensemble system (in contrast to the resonant regime considered here).

A. Basic effect
A striking feature of the purely dissipative dynamics described by Eq. (2) is an extreme sensitivity to the parity of the number N of spins. As we will see, the steady state can be macroscopically different for N spins vs. N + 1 spins. While early work noted that the form of the steady state depends on parity [21,22], subsequent studies on achievable squeezing focused on the even-N case [23,25]. Our work reveals important new aspects of this parity effect. We show that by appropriate parameter tuning, one can avoid this effect, allowing steady-state squeezing that is near Heisenberg limited regardless of the parity of N . We also discuss a different regime where the even-odd effect could be used for a new sensing modality based on the macroscopic sensitivity to spin-number parity. Crucially, we show that there is Sketch of the steady state for (a) even N and (b) odd N . The size of the black circles represents the population of a level |j, m . For even N , a pure dark state exists for any squeezing parameter r because the jump operatorΣ leads to destructive interference between adjacent levels (blue arrows) such that every second level is unoccupied (dashed red lines). For odd N and large r, the interference condition cannot be satisfied for all levels (brown flashes) and the steady state is mixed. The two pure-state contributions with largest statistical weight are sketched here.
no long timescale associated with the emergence of this sensitivity to the addition or removal of a single spin. Note that the even-odd effect in dissipative spin squeezing has no counterpart in bosonic dissipative squeezing. We start with a simple intuitive picture that explains why the steady state of Eq. (2) is so sensitive to the parity of N . Recall that pure bosonic squeezed states are fully paired: they are superpositions of states having even photon numbers only [36]. A similar structure holds in our spin problem. We can think of the fully polarized state |j, m = −j as being the "vacuum", and a state |j, m = −j + q as having q excitations (i.e., q flipped spins). We thus see directly from Eq. (7) that, like bosonic squeezed states, the spin dark states |ψ dk [j; r] also only involve even numbers of excitations q.
Formally, in both the bosonic and spin problem, this paired structure leads to destructive interference that makes the state dark. WhenΣ[r] acts on a paired state, it creates a state having only odd number of excitations. For a given odd excitation number q odd , achieving a dark state requires destructive interference between the two pathways leading to q odd :Ŝ − could have acted on the state with (q odd + 1) excitations, orŜ + could have acted on the state with (q odd − 1) excitations. These destructive interference conditions can be directly used to derive the coefficients in Eq. (7) that determine |ψ dk [j; r] . This structure is shown schematically in Fig. 2

(a).
With this picture in mind, it is easy to see why we cannot have a pure dark state for odd N . In this case, the maximum number of excitations q max is odd. As such, the needed destructive interference is impossible to achieve. Starting with a fully paired state, we can create a state with q max excitations by acting withŜ + on |j, −j + (q max − 1) . However, there is no complementaryŜ − process, as there is no state with q max + 1 excitations. The best one can then do is to construct fully paired states that are only approximately dark due to this incomplete destructive interference [see Fig. 2 The net result of this "frustration" is dramatic: for odd N and large r, the dissipative steady state of Eq. (2) is impure and, moreover, exhibits no spin squeezing for large r. More specifically, for odd N , the steady-state squeezing diverges in the large-r limit, while the purity tends asymptotically to 1/3. This behaviour is shown explicitly in Fig. 3. One also sees that, for modest r, there is no appreciable even-odd effect: the odd-N steady state is almost pure and has the same squeezing as the even N case. This also follows from our heuristic picture: for small enough r, there is very little probability to have a large number of "excitations", and hence one is almost insensitive to the frustration resulting from the cut-off on maximum excitation number.
While our discussion has been focused on the ideal quantum master equation (2), the even-odd effect persists even in the presence of single-spin relaxation and dephasing [as described by Eq. (11) in the limit κ int → 0]. As discussed in App. D, observing the even-odd effect in the steady state requires the single-spin cooperativities η rel and η φ defined in Eq. (10) to be order unity or larger.
Finally, we note that the even-odd effect discussed here is distinct from the sensitivity to parity exhibited by unitary evolution under a OAT Hamiltonian H OAT = χŜ 2 x [37][38][39]. The unitary evolution generated byĤ OAT for a time π/2χ maps the initially fully polarized state |N/2, −N/2 to Greenberger Horne Zeilinger (GHZ) states oriented along orthogonal axes in phase space, depending on the parity of N . This coherent effect results in a strong sensitivity to parity at a particular instant in time; in contrast, in our system, we have a dissipative effect where the sensitivity manifests itself in the steadystate of the system. Moreover, in our case, the even vs. odd states are not equivalent up to a rotation, but differ both in their purity and the magnitude of their fluctuations. N and its purity (dash-dotted red line) approaches 1/3. Inset: Minimum Wineland parameter for even N (crosses) and odd N (dots) obtained at an optimal squeezing strength ropt. The squeezed spin component is alwaysŜy.

B. Parity-independent Heisenberg-limited squeezing
In most experimental situations, the even-odd effect will be a nuisance: one aims for strong steadystate squeezing without needing to control N at the single particle level. We therefore derive a quantitative estimate on the maximum squeezing parameter r that can be used without any parity sensitivity. For small r, the system and its steady state are well described by a Holstein-Primakoff approximation [40]; one recovers bosonic squeezing physics [41], which is independent of the parity of N . However, the correspondence between bosonic squeezing and spin squeezing will break down if the populations of the states |j, m ≈ j become nonzero. Using the steady-state occupation number of the Holstein-Primakoff bosons b †b ss = sinh 2 (r), one can estimate that this breakdown happens if the condition b †b ss ≈ N/2 holds. This yields the breakdown criterion which provides an estimate for the maximum squeezing parameter r possible with no even-odd effect. While working with e 2r N avoids parity effects, one might worry that this constraint precludes ever reaching Heisenberg-limited scaling of the steadystate spin squeezing. This is not the case. As shown in the inset of Fig. 3, the minimum Wineland parameter for odd N (obtained at a squeezing parameter r opt , which depends on N ) exhibits Heisenberg-like scaling, and the spin squeezing differs only by a constant prefactor ≈ 2.6 from the maximum achievable spin squeezing of ξ 2 R,HL = 2/(N + 2) for even N , which is obtained in the limit e 2ropt N .

C. Connections to the LMG Model
Despite first appearances, the extreme even-odd effect of our system is more than a nuisance. At a fundamental level, the effect has a surprising connection to a seemingly unrelated closed-system manybody model, the LMG model [28]. To see this, recall that in a quantum trajectories formulation of the master equation in Eq. (2), the evolution of a state vector in the absence of quantum jumps is governed by the non-Hermitian Hamiltonian (−i/2)Ĥ LMG , wherê x + e 2rŜ2 y +Ŝ z . (15) H LMG is precisely the Hamiltonian of the anisotropic antiferromagnetic LMG model [28], a generalized transverse field Ising model with all-to-all Ising couplings. For even N , we are thus dissipatively stabilizing the many-body ground states |ψ dk [j; r] of the antiferromagnetic LMG model [29,42] and converge to one of them depending on the total angular momentum j of the initial condition. Note that the physics here has crucial differences from the more studied ferromagnetic LMG model, which is also known to exhibit spin squeezing in its ground state [43,44] (but has no simple connection to a dissipative protocol). Focusing on the case where N is odd and r > 0, H LMG is positive and the steady state of Eq. (2) in a given total-angular-momentum subspace j can be written as (see App. D) where λ k and |ψ k are the ordered eigenvalues and eigenvectors ofĤ LMG . We can thus directly connect the properties of the odd-N steady state to the spectrum of the LMG Hamiltonian. Consider first the limit r → 0, whereĤ LMG →Ŝ 2 −Ŝ 2 z +Ŝ z . Then, the Hamiltonian has a unique ground state |ψ 0 → |ψ dk [j; 0] = |j, −j . Moreover, the groundstate energy is zero for any N and the gap to the double-degenerate first excited states is finite, i.e., lim r→0 λ 0 = 0 and lim r→0 λ 1,2 = 2j. As a result, the steady state is approximately pure even when N is odd, as |ψ 0 dominates the sum in Eq. (16).
In the opposite limit r → ∞, the LMG Hamiltonian is dominated by theŜ 2 y term,Ĥ LMG ≈ e 2rŜ2 y , and its eigenvalues are the eigenstates |j, m y ofŜ y with energy λ m ≈ m 2 e 2r . Now, there is no zero energy ground state for odd N (because m takes halfinteger values), the ground state is double degenerate, and the steady state converges to an incoherent mixture ofŜ y eigenstates, A direct computation shows that the purity converges to lim N →∞ Tr (ρ (N/2) ss ) 2 = 1/3. In the limit r → ∞, there is no mean spin polarization, but the variance ofŜ y remains finite, Ŝ 2 y ≥ 1/4. As a result, the Wineland parameter will diverge as shown in Fig. 3.
The connection to the LMG model thus provides useful intuition into the odd-N steady state. For even N , the dark state |ψ dk [j; r] remains an exact zero mode ofΣ † (r)Σ(r) for any value of r and interpolates smoothly between the limits |ψ dk [j; 0] = |j, −j and lim r→∞ |ψ dk [j; r] = |j, 0 y . In terms of the LMG model, this implies that for even N , the ground-state gap does not close as a function of r [45]. This feature of the antiferromagnetic LMG model has been discussed previously in the context of a closed system quantum phase transition [29,42,46].

D. Enhanced sensing
The dramatic even-odd sensitivity of the steady state, which has no counterpart in bosonic spin squeezing, could enable a new kind of sensing modality: it provides a means for detecting changes in N at the single-spin level. This kind of sensing has long been of interest for both fundamental studies and applications [47][48][49][50][51][52]; a recent experiment has even used dispersive sensing to measure real-time changes in atom number in an atomic ensemble dispersively coupled to a cavity [53]. Our dissipative setup could provide an alternative route for an analogous kind of sensing.
As discussed above, for a large squeezing parameter e 2r N , the squeezing of the collective steady state depends exponentially on the parity of N (see Fig. 3). A simpler quantity, the variance ofŜ y , also exhibits this strong sensitivity in the large r limit: for even N , it vanishes like N 2 e −4r /8 whereas, for odd N , it converges to the constant value N/π 2 if N 1. We thus see that measuringŜ 2 y provides a direct means for estimating the parity of N . Such collective spin fluctuation measurements have been implemented in variety of systems [54][55][56][57][58][59][60].
While the parity sensitivity is in principle a steady-state effect, the relatively fast relaxation timescale here means that it can be harnessed for real-time sensing. We stress that the strong sensitivity to parity does not come at the expense of a vanishingly small bandwidth: if a spin is suddenly lost, the relaxation time to the new opposite-parity steady state is (at worst) set by the inverse coupling rate 1/Γ. This timescale does not grow with system size [see inset of Fig. 4(a)]. The relaxation is even faster if one is in the maximum-j subspace; here, the relaxation rate is collectively enhanced by a factor of N .
We thus have a means for detecting spins leaving or decoupling from the cavity one by one, as each such event causes a large change inŜ 2 y [see Fig. 4(b)]. Note that the variance detection requires multiple repetitions of the measuement to distinguish even N from odd N : although the probability to obtain a measurement result with |m y | > 1/2 is negligible for even N in the large-r limit, it is only between 15 % and 19 % for odd N [cf. Eq. (16) and App. D]. One thus has to wait for a probabilistic measurement outcome with sufficiently large |m y | (the value depends on the detector resolution) to determine the spin-number parity unambiguously.
Imperfections of the squeezing process, e.g., an impure engineered reservoir, and local dissipation will reduce the visibility of the even-odd effect (see Sec. VI and App. D, respectively). In Sec. VII D, we discuss that even-odd effect can be observed in a state-of-the-art trapped-ion platform for N 10. This opens the exciting possiblity to experimentally verify the even-odd effect in spin squeezing, which has no counterpart in bosonic squeezing.

IV. ENHANCED PROTECTION AGAINST SINGLE-SPIN RELAXATION
The dissipative approach to spin squeezing also provides strong advantages when unwanted singlespin dissipation is included. In this section, we focus on the case where local relaxation is dominant, i.e., we study Eq. (11) in the limit γ rel = 0, γ φ → 0. For atomic systems, this can be viewed as a fundamental limit arising from spontaneous emission, whereas single-spin dephasing is a technical imperfection. As noted in Ref. 35, in this limit, standard OAT achieves an optimized squeezing that yields the scaling ξ 2 R ∼ C −1/3 rel for large N . This work also introduced an alternate Hamiltonian protocol involving two mutually-interacting spin ensembles, which could achieve a more favourable ξ 2 y is plotted for a system described by the ideal quantum master equation (2) with r = 2.5, starting from the state |N/2, −N/2 . The evolution is interrupted at randomly chosen times (black triangles), where a single (randomly chosen) spin is removed from the system. These spin-loss events cause the system to relax to a new steady state, leading to dramatic swings in the value of Ŝ 2 y after each loss event. Note the logarithmic scale used for the y axis.
at a specific optimized time. As we show below, our dissipative approach can achieve an identical scaling, but now for the steady state, and only using a single ensemble of standard two-level systems. We also show that this enhanced performance over OAT holds even for small-N ensembles. Note that singlespin dissipation was also studied in Ref. 25, but only for spontaneous emission in an ensemble of four-level atoms with a specific structure. This is distinct from the more generic model Eq. (11) we study.
Focusing on the limit of large N and a small singlespin cooperativity, we can approximate our system well using a standard mean-field theory based on linearizing the equations of motion for the system's covariance matrix. Solving these in the steady state and considering the limit of a sufficiently large r (see App. C), one finds that the steady-state squeezing is The numerator here describes unwanted heating by both single-spin relaxation and the collective decay γ coll associated with internal cavity loss. The only parameter left to optimize over is κ sqz , the coupling between the cavity and the squeezed reservoir, which enters Eq. (18) via Eqs. (12) and (13). There is a non-trivial minimum here. Suppression of unwanted collective heating requires a large κ sqz , as this reduces the ratio γ coll /Γ. In contrast, suppressing the effects of γ rel requires a large Γ and hence small κ sqz .
Minimizing with respect to κ sqz , we find where the optimal value of κ sqz satisfies We thus obtain an optimized squeezing that scales significantly better with collective cooperativity in this relaxation-dominated regime than the OAT result of ξ 2 rel . In App. C 2 we show numerical simulations of a more accurate non-linear mean-field theory that confirm these results. As we have stressed, the squeezing here is also achieved in the steady state (and not just at one optimal time). While we assumed a large value of r to derive these results, in practice one only needs exp(−2r) 1/ √ C rel for this scaling to hold.
The advantage over OAT in this relaxationdominated regime also persists for smaller-sized spin ensembles. To study this regime, we numerically solve Eq. (11) for the steady state. Figure 5 shows the obtained results for the steady-state squeezing (orange curve) as a function of N , where we have fixed g, κ int , and γ rel so that the single-spin cooperativity is η rel = 2. For each value of N , we optimize the parameters of the squeezed reservoir (κ sqz , r) to minimize the steady state ξ 2 R ; the optimized values are presented in App. F. For comparison, we also plot the optimized transient squeezing achievable using OAT (blue curve) in an identical cavity-spin system [27,35] (see App. G for details). For the OAT setup, there is no squeezed reservoir, κ sqz = 0, and there is a large detuning ∆ between the spins and cavity, which is optimized for each value of N . Figure 5 shows that, even for small N , the dissipative protocol yields an advantage over OAT. While for these small values of N and large η rel , the linearized mean-field theory scaling predictions are not expected to hold exactly, there is a qualitative agreement with the predicted power laws (as indicated by black dashed lines).
In App. H, we provide a brief performance analysis of a special case where κ int = 0. Mathematically, such a scenario is equivalent to a setup where one directly shines squeezed light onto the spin ensemble. We show that in the limit of large spin number, one can achieve the scaling of ξ 2 R ∝ (N Γ/γ rel ) −1 , although naturally, having either κ int = 0 or irradiating a spin ensemble directly, would likely be difficult to realize experimentally.

A. Pre-thermalization and emergent slow timescales
We now consider the effects of weak single-spin dephasing [i.e., the γ φ term in Eq. (11)] on our dissipative spin squeezing protocol. For very weak dephas-ing, such that γ φ < γ rel /N holds, the mean-field theory results of the previous section still provide a good description; one simply substitutes γ rel → γ rel + 2γ φ in Eqs. (19) and (20). The more interesting case is when dephasing is the dominant form of singlespin dissipation, but is still weak compared to the rate Γ associated with the collective spin squeezing dissipator (i.e., Γ γ φ γ rel ). In this case, the dynamics is surprisingly rich, exhibiting features reminiscent of prethermalization behavior observed in weakly nonintegrable systems [30,31]. Prethermalization is associated with approximately conserved quantities that can only be dynamically randomized on extremely long timescales; this results in an intermediate-time quasi-steady state whose form is contingent on the initial value of the conserved quantities. Here, a similar phenomenon arises, with total angular momentum playing the role of the approximately conserved quantity. We discuss this more in what follows.
Starting from an initial product state, we find that a seemingly tiny amount of single-spin dephasing is enough to completely destroy spin squeezing in the eventual steady state. Using a mean-field analysis, one can show that in the presence of arbitrarily weak but non-zero single-spin dephasing (and γ rel = γ coll = 0), the steady-state squeezing is bounded by −3 dB in the large N limit: where the optimal value is achieved with r = 1 8 ln N . Despite this, there exists an extremely long-lived intermediate-time regime (a quasi-steady state) where strong spin squeezing is observed. The system's dissipative dynamics is thus characterized by two vastly different timescales, as shown in Fig. 6. The system first evolves into a transient spinsqueezed state on a fast timescale ∝ 1/N Γ. In contrast, the eventual relaxation to the true steady state (which has minimal squeezing) occurs on a much slower timescale ∝ N/γ φ . For a large system size N , the ratio of these timescales can be dramatic. We also note that the slow relaxation time is parametrically slower than the single-spin dephasing time 1/γ φ .
The emergence of this surprisingly long timescale, and the corresponding fragility of the steady state to weak dephasing, are both surprising; we stress that single-spin relaxation (as discussed in the previous section) does not give rise to an analogous behavior. In App. E, we analyze this effect using Liouvillian perturbation theory [61] and develop an intuitive physical picture: Single-spin dephasing enables transitions between subspaces of different total angular momentum [34] such that an initial state in the  (11) for weak local dephasing, r = 1.0, and N = 50 spins (thick dashed orange line, γ φ /Γ = 0.005, γ rel /Γ = 0, and γ coll /Γ = 0). The final amount of steady-state spin squeezing is indicated by the thin dash-dotted orange line. Local dephasing deteriorates the amount of steady-state spin squeezing compared to an ideal system without local dissipation (solid blue line, γ φ /Γ = γ rel /Γ = γ coll /Γ = 0). Local relaxation counteracts this effect and partially restores the steady-state spin squeezing (dotted green line, γ φ /Γ = 0.005, γ rel /Γ = 0.001, and γ coll /Γ = 0). Note that the transient state is strongly spin squeezed even in the presence of local dissipation since the collective dissipator Σ induces spin squeezing on a short timescale ∝ 1/N Γ whereas the system approaches its steady state on a longer system-size-dependent timescale ∝ N/γ φ . (b) Wineland parameter calculated using the mean-field equations of motion detailed in App. B for N = 1000 spins and the same sets of dissipation rates as in (a). j = j max subspace evolves into a steady state populating subspaces with j < j max . The degeneracy of the j < j max subspaces gives rise to anomalously small matrix elements between the subspaces, which represent bottlenecks for the relaxation to the steady state.
We stress that the surprising impact of dephasing need not be problematic for experiments. The spin squeezing exhibited by the Wineland parameter ξ 2 R in the "prethermalized" intermediate-time regime is comparable to ξ 2 R of the steady state obtained in an ideal system without single-spin dissipation, as long as the conditions γ φ Γ, N γ φ are satisfied. Moreover, there is a simple but effective way to improve the spin squeezing of the steady state by deliberately adding a competing single-spin relaxation process γ rel . If this relaxation rate satisfies the condition γ rel γ φ /N , population will be pushed back to the large-angular-momentum subspaces, which decreases the steady-state Wineland parameter significantly and increases spin squeezing beyond the −3 dB limit, as shown in Fig. 6.

B. Inhomogeneous broadening and dynamical decoupling
In addition to the Markovian mechanism described in the previous section, in some platforms dephasing due to inhomogeneous broadening of the spin ensemble can also play a role. A major advan-tage of spin squeezing generated by OAT dynamics is that it is compatible with dynamical decoupling and thus allows for an effective cancellation of the impact of inhomogeneous broadening by a simple sequence of π pulses about the x axis [27]. At first glance, this does not seem to be the case for our dissipative scheme. However, we will show here that our dissipative scheme is in fact compatible with a slightly modified dynamical decoupling sequence.
Our starting point is a generalization of the Hamiltonian of Eq. (9), where ∆ j = ω j − ω 0 denotes the detuning of spin j from the resonance frequency of the bosonic modê a due to inhomogeneous broadening. Using average Hamiltonian theory [62], we will now derive an effective Hamiltonian for the two different decoupling sequences shown in Fig. 7, which are designed such that the effects of the ∆ j terms cancel out on average. Instead of transforming the state of the system at each decoupling pulse, it is more convenient to consider a Heisenberg picture where the Hamiltonian changes at each pulse, the so-called toggling frame [62].
A single π pulse about the x axis will flip the sign ofσ erators in the interaction term of Eq. (22), The first effect is desired and will cancel inhomogeneous broadening, but the second effect is unwanted because it will turn damping of the Bogoliubov modê Σ(r) into anti-damping. If one can control the coupling constant g as a function of time, one can switch off the undesired interaction after every second π pulse, as shown in Fig. 7(a), and obtains the average Hamiltonian where the inhomogeneous broadening has been canceled at the cost of a reduction of g by a factor of 2. Experimentally, the coupling g could be switched off by detuning the spins rapidly from the cavity. Instead of switching off the interaction between the spins and the bosonic mode for half of the period T , one could also use the dynamical decoupling sequence shown in Fig. 7(b). By applying a π rotation about the z axis, one can flip the sign of the second term without disturbing the first one, The sequence is terminated by a π rotation about the y axis which reverts all signs and restores the original Hamiltonian (22). If the waiting times between the pulses have a ratio of 2 : 1 : 1, both inhomogeneous broadening and the unwanted interaction terms will be canceled in the average Hamiltonian, Generating the additional π pulse about the z axis may seem challenging because it requires a controlled detuning of the spins from the cavity such that the accumulated phase is exactly π. Experimentally, this would likely be even more difficult than turning the coupling off for half a period, and one may conclude that this scheme is harder to implement than the first one. However, it is well-known from NMR that pulses about the z axis can also be realized using a so-called composite pulse which is a suitable combination of x and y rotations [63]. Specifically, a π pulse about the z axis can be decomposed into pulses about the x and y axes as follows: Dynamical decouping sequences to cancel inhomogeneous broadening in Eq. (22). Each sequence has a total duration time T , where κsqzT 1 and κintT 1, and is repeated multiple times. (a) A simple sequence of π pulses about the x axis will generate unwantedΣ † (r) antidamping terms [cf. Eq. (23)]. Therefore, the coupling to the spins, g(t), must be switched off after every second πx pulse. (b) By adding two additional π pulses about the y and z axes, the detrimental impact of the first πx pulse can be cancelled and the coupling g(t) can be kept constant.
which does not require any detuning between the spins and the cavity. Higher order contributions to the average HamiltonianĤ will become negligible if the period T of the decoupling sequence satisfies the conditions κ sqz T 1 and κ int T 1.

C. Non-uniform single-spin couplings
Another experimentally very relevant source of imperfections are inhomogeneities in the coupling strength g between the spins and the common bosonic mode, i.e., one has to replace the Hamiltonian in Eq. (9) witĥ This breaks the permutational symmetry of the spins reflected in the collective spin operatorsŜ ± . The standard strategy to analyze this effect is based on an expansion of the mean-field equations of motion around the average couplingḡ = N j=1 g j /N [64]. Variants of this approach have been applied to study superradiance [64,65], microwave quantum memories [66], and spin squeezing [67][68][69][70][71]. Nonuniform couplings then lead to an approximate collective model with a renormalized coupling parameter and a reduced effective length of the collective spin vector. Similar results will hold in our case if the squeezing parameter is not too large, e 2r N .
To analyze the case of strong squeezing, e 2r N , more powerful theoretical methods need to be de-veloped, which is an interesting and relevant subject for further study. We expect non-uniform single-spin couplings to reduce the visibility of the even-odd effect and to influence the prethermalization physics, since the suppression of transition rates crucially relies on the indistinguishability of the spins in the ensemble.

VI. IMPURE ENGINEERED RESERVOIR
In this section, we explore a different kind of imperfection that has not been studied in previous discussions of dissipative spin squeezing: the engineered reservoir may have an imperfect purity (or equivalently, mimic thermal squeezed light rather than vacuum squeezed light). At the most basic level, this corresponds to modifying our ideal quantum master equation (2) bẏ where n th ≥ 0 parameterizes the effective temperature of the squeezed reservoir. We will discuss below in Sec. VII and in App. I how this generic model can be related to more microscopic mechanisms (including collective decay, γ coll > 0). Even though the hybrid-systems reservoir-engineering approach we focus on is not limited by losses associated with the transport and injection of an externally prepared optical squeezed state, many reservoir engineering techniques will inevitably result in a n th = 0, hence it is important to understand the impact of this unwanted heating. Note that for n th = 0, the steady-state of the spin ensemble will necessarily be impure. This can have a deleterious impact on Ramsey-type sensing experiments if one is interested in signal phases that are not infinitesimally small (as has been discussed in the context of OAT [72,73]). For e 2r N , the spin squeezing described by Eq. (29) is well-approximated by a linearized bosonic master equation (via use of the Holstein-Primakoff approximation [40]). In this limit, one thus expects that a small n th will have only a small impact on the steady-state squeezing [13], i.e., The interesting question is whether n th also has innocuous effects for the larger values of r needed to approach the Heisenberg limit. Fig. 8(a) shows that this is not the case. The linearized bosonic theory breaks down if the squeezing parameter r is too large, with dramatic consequences: even small imperfections, n th 1, cause the steady-state Wineland parameter to strongly deviate from its ideal limit ξ 2 R,HL = 2/(N + 2). Further, one finds that the steady state Wineland parameter increases with increasing r (irrespective of the parity of N ), i.e., the steady-state squeezing is a non-monotonic function of r and exhibits a minimum at an optimal value r opt (which depends on N and n th ).
We thus have an important caveat: if the engineered squeezed reservoir has a non-zero effective temperature, increasing the reservoir squeezing parameter r does not result in ever-increasing steadystate spin squeezing. We stress that this is true even when N is even. Numerical simulations indicate that the minimum Wineland parameter at r opt almost follows a Heisenberg-like scaling with N , as shown in Fig. 8(b), but with a significantly larger prefactor than the ideal result ξ 2 R,HL = 2/(N + 2). Further numerical results exploring the parameter dependence of the optimal squeezing parameter r opt are given in App. J, while App. K briefly discusses scaling obtained using a mean-field theory approach.
Interestingly, Fig. 8(c) shows that, although the ratio ξ 2 R (r opt )/ξ 2 R,HL deviates from the bosonic expectation, the purity of the steady state at the optimal squeezing parameter r opt does closely follow the corresponding relation (2n th + 1) −1 valid for bosonic dissipative squeezing.
At a heuristic level, the much stronger sensitivity of dissipative spin squeezing to n th > 0 and the similar purity (in comparison to dissipative bosonic squeezing) can be understood by the following simplified picture (see also App. J). Similar to dissipative bosonic squeezing, the dominant contributions to the mixed steady state of Eq. (29) are the dark state ofΣ, |χ 0 = |ψ dk [r] , and the "first excited" state |χ 1 ∝Σ † |ψ dk [r] . Their statistical weights are in a thermal ratio p 1 /p 0 = n th /(n th + 1), which explains the similar n th -dependence of the purity for r ≤ r opt . However, the Ŝ 2 y variance of |χ 1 differs strongly from its counterpart in a bosonic squeezed state for r r opt : in this limit, the dark state converges to the m y = 0 eigenstate ofŜ y , |ψ dk [r] → |j, 0 y [cf. Eq. (5)], with a vanishing Ŝ 2 y variance. The operatorΣ † expressed in theŜ y basis contains spin raising and lowering operators, i.e., the state |χ 1 is an equal superposition of the |j, ±1 states and has a finite Ŝ 2 y variance. This leads to an increase of the Wineland parameter with increasing r as soon as the Ŝ 2 y variance approaches to its nonzero minimum value.
Thus, the exponential increase of the Wineland parameter for r ≥ r opt has a very similar origin as the corresponding effect in the the odd-N zeron th case discussed in Sec. III. These results are yet another demonstration of the fact that dissipative spin squeezing is more complicated than dissipative bosonic squeezing, due to the finite-dimensional Hilbert space and the intrinsic nonlinearity of spin systems.
A consequence of these findings is that, for large squeezing parameters r r opt , an impurity of the engineered reservoir reduces the even-odd effect. We will discuss the consequences for the even-odd effect below in Sec. VII D, where we will show that the even-odd effect can nevertheless be observed on state-of-the-art experimental platforms.

VII. HYBRID-SYSTEMS IMPLEMENTATION USING DISSIPATIVE BOSONIC SQUEEZING
As discussed in Sec. II B, the dissipative spinsqueezing setup described by the general quantum master equation (11) can be realized using standard two-level systems (unlike the more structured fourlevel atoms in Refs. [9,25]), and without requiring the use of non-classical squeezed input light. Instead, one harnesses a standard (resonant) Tavis-Cummings coupling between a spin ensemble and a bosonic mode, along with the dissipative squeezing of this bosonic mode which is engineered by coupling the bosonic mode to a lossy auxiliary mode that is driven by imbalanced red-detuned and blue-detuned sideband drives. The second element here has been experimentally realized in a variety of systems. In this section, we provide more details on the physical implementation of our hybrid-systems approach to dissipative spin squeezing in three promising platforms: trapped ions, solid-state optomechanical devices, and superconducting circuits.

A. Trapped ions
In trapped ions, the relevant spin degree of freedom usually corresponds to two metastable internal states (spin or orbital) of each individual ion. In contrast, the bosonic "cavity" mode corresponds to a collective motional mode of the ions [74] and the coupling parameter g now characterizes the spin-phonon coupling. Recent experiments have already utilized the spin-motion coupling for over 50 ions in a 2D Penning trap [57] and 1D linear Paul trap [75]. The desired Tavis-Cummings coupling is commonly realized by applying a laser field that is resonant with the red motional sideband of the spin-level transition (see, e.g., Ref. 76). Motional dissipation is, in turn, mediated by coupling the motional mode to a dipole-allowed transition of an ion.
To realize dissipative spin squeezing with N spins, we imagine a setup that consists of N + 1 ions. N of these ions make up the spin ensemble that we wish to squeeze; the remaining additional ion serves as a "cooler" ion that is used to dissipatively squeeze the collective motional mode, as shown in Fig. 9. A squeezed bath can be engineered by applying two laser fields that are resonant with the red and blue sideband transitions of the cooler ion [12], leading to an effective Hamiltonian where G ion ) is the red (blue) sideband coupling, andσ − is the lowering operator of the coolerion transition. The squeezing strength can be controlled by the ratio of the couplings, i.e., tanh(r) = |G ion |, and the squeezed axis is determined by their relative phase. We stress that the sideband transitions are implemented by classical drives, i.e., no squeezed radiation is required for the bath engineering. Reading out the spin-squeezed state can be performed by individually measuring each ion in the Pauliσ y basis. The collective spin variance can then be obtained from the statistics of the measurement outcomes collected in multiple runs.
To assess the practical requirement of our proposal, we adopt the state-of-the-art system parameters in the dissipative motional squeezing experiment by Kienzler et al. [16]. This experiment employs the bath engineering techniques described above to prepare the motion of a single trapped ion in a −12.6 dB squeezed ground state with up to 12% infidelity. Utilizing this scheme for spin squeezing requires two additions.
First, the reservoir engineering technique has to be extended to the collective motion of an ion chain.
Since the center-of-mass mode frequency is not altered by the number of ions in the trap, and the frequencies of other motional modes are well resolved [77], the same sideband drives used in the single-ion case can be applied to engineer the same squeezed bath for the collective center-of-mass mode. Note that, for incoherent electric field noise, the centerof-mass heating rate does not depend on the number of ions in the trap [78]. Second, spin-motion coupling has to be applied between the center-of-mass mode and the N system ions with a spin-phonon coupling strength g. To implement the collective spin dissipator in Eq. (11), a sufficient but not necessary condition is that the coupling strength g is sufficiently weak to guarantee the adiabatic elimination of the motional mode, while sufficiently strong such that qubit decoherence does not deteriorate the spin squeezing, i.e., κ sqz / √ N g √ κ sqz γ φ . From the time evolution of squeezed-state pumping in Ref. 16, we estimate the experimentally realized squeezed-reservoir coupling rate to be κ sqz ≈ 0.5 kHz. Moreover, γ φ ≈ 0.1 Hz has been realized in many experiments (e.g., Ref. 79). With these realistic parameters, a coupling g ≈ 30 Hz should fit in this regime for a modest chain with N 10 ions. Since spin-motion coupling as strong as kHz has been routinely implemented in trapped-ion experiments (e.g., Ref. 80), our desired range is thus well achievable by simply using a weaker drive.
The remaining issue is whether the thermal excitation due to motional heating will mask the desired even-odd effect. Motional heating leads to collective spin excitation and relaxation processes at equal rates γ heat . As discussed in App. I, this can be mapped onto a quantum master equation of the form (29) with an effective squeezing parameterr, i.e.,Σ(r) →Σ(r) wherẽ and with an effective thermal occupation number n th given by Kienzler et al. achieved a squeezed ground state with a fidelity larger than 88 % at −12.6, dB of squeezing [16], which corresponds to n th ≈ 0.14 or, equivalently, a motional heating rate γ heat /Γ ≈ 0.017. Note that this is a worst-case estimate of the fidelity since its reported value includes measurement errors. Consequently, the true value of γ heat in the experiment is likely lower. However, even with this pessimistic values the even-odd effect can be observed in experiments, as we will show below in Sec. VII D.

B. Solid-state spins in an optomechanical crystal
In solid-state platforms, the spin ensemble in our scheme could be realized using defect centers in a semiconductor, e.g., NV-center defect spins in diamond. These spins can be implanted in a structure which in turn realizes an optomechanical crystal: a patterned photonic crystal beam with a defect that localizes both a mechanical mode and an optical mode [81]. We note that high-Q diamond optomechanical crystals have been realized experimentally [82], with a recent experiment even integrating such a system with NV-center defect spins [83]. The localized mechanical mode plays the role of the bosonic "cavity" in Eq. (11). The spins and mechanical motion exhibit an intrinsic coupling due to the strain dependence of spin-level transitions, and the coupling could be further enhanced by incorporating the high strain sensitivity of excited states through phonon-assisted Raman transitions [83][84][85].
In this kind of setup, the optomechanical coupling between the localized mechanical and optical cavity mode provides a mechanism for the dissipative squeezing of the mechanical mode. If one is in the sideband-resolved regime (where the mechanical frequency is larger than the optical cavity decay rate), then this dissipative mechanical squeezing can be realized by driving the optical cavity by two control lasers that are resonant to the red and blue motional sidebands, respectively [41]. We stress that these drives are classical, coherent state drives. Ignoring the non-linear coupling that is usually negligibly weak in most platforms, the optomechanical coupling is well approximated bŷ where G OM |, which can be tuned by varying the amplitude of the driving tones. We note that this kind of dissipative squeezing of mechanical motion via optomechanics has been realized in several experiments [15,[17][18][19]. Our protocol thus provides a means of harnessing this capability to generate spin squeezing. Finally, in solid-state settings inhomogeneous broadening of the spin ensemble is almost always an issue; this is typically mitigated by using dynamical decoupling techniques. For spin-squeezing protocols based on OAT dynamics, a very simple decoupling sequence can be used, which repeatedly applies π pulses about the x axis [27]. This very simple strategy fails in our case because it transforms theΣ decay term in Eq. (11) into aΣ † anti-damping term. However, as discussed in Sec. V B, this unwanted excitation dynamics can be canceled by two additional π pulses about the z and y axis, which makes our protocol compatible with dynamical decoupling.

C. Superconducting microwave cavities
Superconducting microwave cavities and circuit QED are another promising class of systems for implementing our ideas. Our basic building block of a bosonic mode coupled to a spin ensemble could be realized by coupling a single microwave cavity mode to either a set of superconducting qubits [86][87][88][89], or to electronic spins in substrate (e.g., Bi donors implanted in Si [90,91]). The second ingredient, a mechanism for the dissipative generation of microwave squeezing, could also be implemented in different ways. One approach is to inject squeezed microwave radiation directly into the cavity using the output of a Josephson parametric amplifier [92][93][94]. This has already been achieved experimentally in Ref. 90, in a system where a cavity has been coupled to a spin ensemble. An alternate approach which has the advantage of not being limited by insertion losses (associated with transporting a squeezed state) is to mimic the same dissipative squeezing protocols used in optomechanics to squeeze a mechanical mode. This can be accomplished by coupling three microwave modes via a Josephson ring modulator [95], which generates a three-wave mixing term (p+p † )(â+â † )(b+b † ) between the modesâ,b, andp [96]. By driving the pump modep coherently at the sum and difference frequency of theâ andb modes, ω ± , one can engineer an interaction of the form of Eq. (34), where the prefactors G (±) OM depend on the strength of the drives at ω ± , respectively. Adiabatic elimination of the strongly-dampedb mode generates an effective squeezed bath for theâ mode as shown in Eq. (9). A recent experiment implementing this approach has demonstrated up to −8 dB of intracavity squeezing of theâ mode [20].

D. Experimental viability of the even-odd effect
Given the results of the Sec. VI, one may worry that a finite effective temperature of the engineered (a) Ratio between Ŝ 2 y for N = 8 vs. 9 as a function of the motional heating rate γ heat (red triangles). The dotted red line is a guide to the eye to highlight a ratio of unity, i.e., no even-odd effect. For each value of N and γ heat , we optimize the squeezing parameter over the range 0 ≤ e 2r ≤ 12 dB to minimize the variance Ŝ 2 y (blue data points). (b) Ratio of Ŝ 2 y for N vs. N + 1 spins as a function of N for an experimentally realistic parameter of γ heat /Γ = 0.017 (red triangles). Again, the dotted red line is a guide to the eye to highlight a ratio of unity, i.e., no even-odd effect. For each N , we optimize the squeezing parameter over the same range 0 ≤ e 2r ≤ 12 dB as in (a) to minimize the variance Ŝ 2 y (blue stars).
reservoir will substantially decrease spin squeezing and mask the even-odd effect introduced in Sec. III.
Here, we show that this is not always the case and, in particular, that the even-odd effect can be observed in experimentally accessible parameter regimes. For concreteness, we focus on the trapped-ion platform introduced in Sec. VII A, which is most mature. With a first experimental observation of the even-odd effect in mind, we consider a modest number of spins, N 10. In this platform, motional heating due to classical trap noise is the dominant source of imperfections of the engineered squeezed reservoir. As discussed in Sec. VII A and App. I, its impact can be modeled using Eq. (29) with an effective squeezing parameterr and an effective thermal occupation number n th given by Eqs. (32) and (33), respectively.
Compared to the bare squeezing parameter r (which is related to the amplitude of the red and blue sideband drives), the effective squeezing parameter strengthr is reduced by the motional heating rate γ heat . The effective thermal occupation number n th grows with increasing r, i.e., a large squeezing parameter reduces the purity of the squeezed state. For a given number N of ions, the variance Ŝ 2 y will thus take a minimum value at an optimal squeezing parameter r opt , whose value depends on N and γ heat /Γ, as shown in Fig. 10(a). Comparing these minimum variances for N = 8 vs. 9 ions, we find an even-odd difference of more than 20 % over a wide range of motional heating rates γ heat . The even-odd difference decreases with effective temperature n th and will disappear if the impurity of the squeezed reservoir is sufficiently large. Note that even for the upper bound of the motional heating rate derived in Sec. VII A, γ heat /Γ = 0.017, the even-odd difference is clearly visible.
An alternative way of probing the even-odd effect is shown in Fig. 10(b). There, the motional heating rate γ heat is kept fixed but the number of ions in the trap is varied. Again, we optimize r for each data point individually to minimize the variance Ŝ 2 y . The ratio between these variances shows pronounced oscillations for N ≤ 10 ions. Note that we restrict the optimal squeezing parameter e 2ropt in both scenarios to be smaller than the experimentally achievable −12.6 dB. Thus, our results suggest that the even-odd effect is realizable on state-of-the-art trapped ion platforms.
Instead of varying r for each N individually to minimize the variances of even and odd N independently, one could also use the same squeezing parameter r for a pair of N and N + 1 spins and maximize the ratio of their variances Ŝ 2 y (r) N / Ŝ 2 y (r) N +1 . In this case, the even-odd effect would be even more pronounced than the data shown in Fig. 10.

VIII. CONNECTION TO PREVIOUS WORKS
Here we review previous works on dissipative spin squeezing [9,[21][22][23][24][25], summarize their results and point out the differences to this work.
Agarwal and Puri discussed the idealized spinonly quantum master equation (2) with an additional collective decay term D Ŝ − ρ and a coherent drive ΩŜ + +Ω * Ŝ − [21,22]. They derived explicit expressions for the steady state and pointed out that the steady state of Eq. (2) is pure for even N (wherê Σ has a zero eigenvalue) and mixed for odd N (wherê Σ has only nonzero eigenvalues) [21]. Moreover, they discussed the pairwise excitation structure of the even-N steady state illustrated in Fig. 2(a) and showed numerical results for the odd-N population distribution in the regime e 2r N [21]. In a follow-up article [23], they calculated the Wineland spin-squeezing parameter ξ 2 R of a state of the form |ψ ∝ exp θŜ z exp −iπŜ y /2 |j, m , which contains the even-N steady state (5) as a special case, but not the odd-N steady state. They also discussed squeezing in the presence of a coherent drive [22].
Unlike our work, Agarwal and Puri did not evaluate the spin-squeezing properties of the odd-N steady state. Consequently, they did not find the dramatic difference in spin squeezing between the undriven even-N and odd-N steady states for e 2r N , which is one of the central results of this work. Note that the spin-squeezing parameter r required to see an even-odd effect will decrease if the number N of spins is lowered, which makes the regime e 2r N attainable on state-of-the-art experimental platforms, as discussed in Sec. VII. Therefore, another new and experimentally relevant result of our work is that the even-odd effect is not a mere mathematical complication, but a real and testable physical effect.
Agarwal and Puri suggested to generate spin squeezing by illuminating the spins with squeezed radiation, which is experimentally very challenging. In a later work, Kuzmich et al. proposed a related alternate scheme to generate spin squeezing using Vtype atoms illuminated by squeezed light [24]. They did not discuss dissipative spin squeezing and, consequently, did not comment on the even-odd effect at all. In contrast to the proposals by Agarwal et al. and Kuzmich et al., our work does not require any injection of nonclassical light, which lowers the experimental challenges significantly.
Dalla Torre et al. proposed another method to implement dissipative spin-squeezing dynamics in a specific multilevel atomic system using Raman transitions [25]. They use a pure dark state of the type of Eq. (6) to analyze the ideal situation when there is only collective loss of the form of Eq. (3). In addition, they discuss the impact of dissipation due to single-atom Raman scattering into free space. They do neither comment on an even-odd effect nor mention that their pure-state analysis is strictly speaking only valid in the case of even N .
Our work proves that, while this mathematical treatment is indeed admissible if squeezing is not too strong, e 2r N , substantial changes to the spinsqueezing physics will show up at e 2r N , which have not been discussed in the literature before. Regarding experimental versatility, Kuzmich et al. and Dalla Torre et al. discuss very specific implementations, which are not applicable to a generic ensemble of two-level systems. In contrast, we propose a more generic mechanism that has not been discussed in the literature before. Importantly, our approach is compatible with a wide range of systems including solid-state implementations (see Sec. VII) and is perhaps the most flexible and experimentally viable implementation. This versatility allowed us to discover the prethermalization physics discussed in Sec. V A, which is not present in the scheme of Ref. 25 because of their very specific decay mechanism that couples collective and local dissipation.
Finally, Borregaard et al. considered Λ-type atoms driven by multiple laser drives and identified a dissipative spin-squeezing scheme as its resonant limit [9]. They use this observation to interpret their results, but they do not discuss the dissipative spinsqueezing scheme in detail nor explore any of its consequences revealed in this work, like the even-odd physics and prethermalization.
Finally, we note that in contrast to our work, the previous works on dissipative spin squeezing reviewed above did not discuss or analyze the consequences of an imperfect engineered reservoir (see Sec. VI).

IX. CONCLUSIONS
In this work, we have revisited the reservoirengineering approach to preparing and stabilizing spin-squeezed states. We analyzed in detail a particular implementation strategy that had not previously been studied, but that is compatible with a number of experimental platforms: employ a hybrid-systems approach where one first uses bosonic reservoir-engineering techniques to stabilize a bosonic squeezed state, and then uses this state (via a standard Tavis-Cummings-type coupling) to dissipatively squeeze a spin ensemble. We also discussed how this approach compared favourably to the standard one-axis-twist method for spin squeezing in the presence of single-spin relaxation.
Our work also addressed fundamental aspects of dissipative spin squeezing, with a focus on two general but surprising phenomena. The first was an extreme, macroscopic sensitivity of the steady state to the parity of the number N of spins in the en-semble. We analyzed both how this effect could be avoided (if the goal was to generate spin squeezing without any parity sensitivity), and how it might be harnessed for a powerful new sensing modality. The second general effect we studied was the emergence of a surprisingly long slow timescale and "prethermalization" behavior when weak single-spin dephasing is added to our model. While the steady state in this regime exhibits at best limited squeezing, the intermediate time quasi-steady state can be highly squeezed. Moreover, the reduction of steady-state spin squeezing can be almost completely suppressed by deliberately introducing a small amount of singlespin relaxation.
Finally, we investigated the impact of an engineered reservoir stabilizing an impure steady state. We discovered a strong sensitivity of the Wineland parameter to impurity if the squeezing parameter r is large.
We hope our work will lay the groundwork for near-term experimental implementations of reservoir-engineered spin squeezing in a variety of systems. In future theoretical work, it will be interesting to explore extensions of the models analyzed here. For example, it is well known that collective Hamiltonian interactions that are not truly infinite range can still generate large amounts of spin squeezing [97,98]. Is the same true with dissipative spin-spin interactions, and if so, are the requirements more or less forgiving than in the coherent case? It would also be interesting to study in more detail the effects of disorder, e.g., due to inhomogeneous broadening, both on spin squeezing and on the parity-sensing scheme proposed here.

ACKNOWLEDGMENTS
This work was primary supported by the DARPA DRINQS program (agreement D18AC00014). We also acknowledge partial support by the University of Chicago Materials Research Science and Engineering Center, which is funded by the National Science  In this Appendix, we outline the derivation of the effective quantum master equation (11) of the main text. Our starting point is Eq. (9) describing a collection of spins interacting with a squeezed bosonic mode. For the moment, we ignore the terms in Eq. (9) describing local dissipation of the spins, Assuming that the cavity evolves on a much shorter timescale than the spins, we adiabatically eliminate the cavity by a projection operator technique [99] similar to the calculation outlined in Ref. 100. To this end, we split the quantum master equation into two superoperators, where L intρ is considered to be constant on the timescale defined by L cavρ . Using this approximation, we can formally solve Eq. (A3), Performing a Born approximation, we decompose the state asρ(t) ≈ρ sp (t) ⊗ρ ss cav , whereρ sp (t) is the reduced density matrix of the spin system andρ ss cav is the steady state of L cav . The equation of motion of the reduced spin density matrix is Inserting the explicit form of L int , we find that the integral on the right-hand side of Eq. (A5) depends on the cavity correlation functions Tr cav â † e Lcavtâρss cav = Tr cav âe Lcavtρss These correlation functions decay fast compared to the timescale on whichρ sp evolves, therefore, we can perform the Markov approximationρ sp (t ) ≈ρ sp (t) and rewrite Eq. (A5) as follows: Taking into account the remaining terms in Eq. (9) describing single-spin dissipation, we recover Eq. (11) of the main text.

Appendix B: Mean-field theory equations of motion
In this section, we provide the set of nonlinear equations of motion for symmetrized products of spin operators, for the effective spin-only model considered in the main text, namely Eq. (11). While such a system of equations is not closed, we neglect third-order cumulants (equivalently performing a 2 nd -order cumulant expansion) [101,102], which lets us approximate the third-order expectation values of various operators as We stress that the above approximation is applied to expectation values of symmetrized operators, defined according to the following convention: Given the initial state with spins completely polarized along the −z direction (i.e., Ŝ z = −N/2), the evolution is governed by the equations where Ĉ ZZ = Ŝ 2 z − Ŝ z Ŝ z . We stress that if we assume that Eq. (11) is a result of coupling the spin system to a cavity interacting with an engineered squeezed reservoir with photon loss κ int , then we have and as discussed in the main text and shown in detail in App. A.
Appendix C: Cooperativity scaling of the ξ 2

R parameter
In this Appendix, we provide a derivation of the cooperativity scaling of the Wineland parameter ξ 2 R . We concentrate our analysis on the weak dephasing limit, and start with the case where γ φ = 0 and where only the local decay γ rel as well as the collective cavity-induced decay γ coll are present. A scenario where local spin dephasing is dominant can lead to substantially altered behavior of the system, and is the subject of Sec. V A and App. E.

Analytical derivation
We begin by linearizing the mean-field-theory equations of motion shown in App. B by focusing on the limit where Ŝ z stays fixed at −N/2. This approximation closely reflects the true system dynamics when the spin number N is large and when the single cooperativities η φ or η rel are not much larger than unity, resulting in effective spin squeezing that is far from the Heisenberg limit. Hence, taking Ŝ z = −N/2 (i.e., spins keeping their polarization throughout the evolution and in the steady state), the Wineland parameter takes a simple form, which we can write using the results in App. B as Note that ξ 2 R gets smaller as r increases and hence, in what follows, we will take the limit r → ∞. It is worth pointing out, however, that choosing a finite r which satisfies exp (−2r) 1/ √ C rel is sufficient to reproduce the scaling of ξ 2 R derived below. In the large-r limit, we find Next, we use Eqs. (12) and (13) of the main text to rewrite the above expression as We consider a limit where N → ∞, while G = √ N g stays fixed. In such a case, the last term of the numerator can be dropped. Here it is crucial to point out that in an experimental setting, one will typically not have much control over κ int and γ rel , while κ sqz can be tuned at will through appropriate reservoir engineering (see Sec. VII). Hence, it is important to understand what value of κ sqz should be chosen to maximize the amount squeezing that this protocol can achieve. At first glance, one might think that choosing κ sqz as large as possible (i.e., κ sqz → ∞) is ideal as that maximizes the amount of bosonic squeezing that the spin-coupled cavity experiences. From Eq. (12), however, we see that such a choice will actually limit the value of Γ, which directly impacts the strength of squeezed-vacuum reservoir that the spins see [see Eq. (11) in the main text], resulting in the squeezing performance being strongly limited by the value of γ rel . Hence, as we shall see shortly, the right thing to do is to still choose κ sqz κ int , γ rel , but yet not too large, so that the Γ-controlled process is dominant over the local spin decay γ rel . To see this explicitly, we minimize Eq. (C4) with respect to κ sqz . Assuming N 1, this leads to where, in the second line, we used Eq. (10) of the main text to express the result in terms of the collective cooperativity C rel and then expanded in the limit of large C rel . Our above expression shows the C −1/2 rel cooperativity scaling for the dissipative protocol, which outperforms the C −1/3 rel behavior of the OAT method [35], in the case where spin decay is the dominant local noise process. The optimal value of κ sqz that results in Eq. (C5) reads which confirms the need for κ sqz κ int (given large C rel ), while also showing that it should not be infinitely large. Finally, we stress that ξ 2 R is ultimately limited by κ int /(κ int + κ sqz ), hence it is crucial that an appropriate κ sqz can be realized in an experimental setting.
While the above result has been calculated in the limit where γ φ = 0, a similar expression is valid when some local dephasing is present (i.e., γ φ = 0). In such a case, one can simply assume γ rel → γ rel + 2γ φ in Eq. (C4). As discussed in more detail in Sec. V A of the main text, however, this is only true when γ φ is not too large, namely when γ φ N γ rel . Otherwise, a local dephasing process can have a significant impact on the evolution and therefore dramatically limit the steady-state performance of the protocol.

Mean-field theory simulations
In this section, we present mean-field-theory simulations of the dissipative protocol obtained using R as a function of collective cooperativity C rel . The blue curve corresponds to ξ 2 R calculated by evolving the full (nonlinear) mean-field equations of motion for the dissipative system (see App. B). Here κint = 500g, γ rel = 0.04g, giving η rel = 0.2. The number N of spins is changed in order to vary C rel . At each blue point, both r and κsqz are optimized in order to minimize ξ 2 R . The orange dashed curve shows the corresponding fit (calculated over the three data points with the largest C rel ). The black dashed line describes the optimized squeezing of the engineered bosonic reservoir. The solid black line shows an ideal Heisenberg scaling 2/(N + 2). Finally, the black dotted curve shows the OAT scaling as calculated in [35].
the full (nonlinear) equations shown in App. B. We consider the case where local spin decay dominates over local spin dephasing, and in particular work in the limit of γ φ = 0.
The plot in Fig. 11, shows the scaling of the Wineland parameter as a function of the collective cooperativity C rel . The parameters are κ int = 500g, γ rel = 0.04g, giving η rel = 0.2, while the number N of spins is varied in order to modify C rel . At each blue point, both r and κ sqz are optimized in order to minimize ξ 2 R . The orange curve shows the corresponding fit, which is calculated using the three data points with largest C rel . We see good agreement with the cooperativity scaling discussed in the main text and derived in detail in the section above (where we have linearized the equations of motion). For comparison, the black dashed line describes the optimized squeezing of the engineered bosonic reservoir. The black solid line shows an ideal Heisenberg scaling 2/(N + 2). We also plot the dotted black curve, which corresponds to the squeezing one would get from the OAT protocol in the limit where γ rel dominates over γ φ (in the large C rel limit) -see Ref. 35. The simulations confirm that the dissipative protocol can indeed outperform the OAT approach.

Appendix D: Even-odd effect
In this Appendix, we briefly review previous results on the dissipative steady state of Eq. (2) in the main text and we derive Eq. (16) of the main text. We then comment on variance detection measurements required to use the even-odd effect as a sensor, and we discuss the impact of local dissipation.

Properties of the steady state
Agarwal and Puri derived that the steady state of Eq. (2) isρ if the Hermitian operatorΣ †Σ is invertible [22].
If N is even,Σ(r) has a zero eigenvalue in each subspace of angular momentum j and the associated eigenstates are the dark states |ψ dk [j; r] ∝ e θ(r)Ŝz |j, 0 y , 2. Using the even-odd effect for sensing As described in Sec. III D of the main text, the sensitivity of the steady state on the parity of the number N of spins can be used for sensing. Experimentally, sudden changes in the parity of N can be induced by various mechanisms. Trapped atoms can be physically lost from the trap by collisions with background gas, internal collisions, and photon-assisted processes [104]. If the spin-1/2 degree of freedom is a subspace of an atomic multi-level structure, undesired internal transitions can occur, which take the atom out of the spin-1/2 subspace and effectively remove it from the collective dynamics even though it may still be trapped [25]. Moreover, one could devise a system where the coupling strength of a single spin to the cavity and, thus, to the collective spin depends on an external parameter. A change of this single-spin coupling strength modifies the number of collective spins, which is collectively amplified and yields a large change of the steady state.
Note that such effective atom loss events do not change the collective expectation value Ŝ y = 0 of the distribution. However, the statistics of the fluctuations Ŝ 2 y depends on the parity, as shown in Fig. 12. The parity of N can thus be inferred by imposing a threshold condition on the variance Ŝ 2 y measured using spin-noise spectroscopy [54][55][56][57][58][59][60].

Impact of local dissipation
So far, our analysis of even-odd effects in the steady state has focused on the idealized case without any single-spin dissipation: γ rel = γ φ = γ coll = 0. We found that the Wineland parameters for even and odd N differ strongly in the regime e 2r N , as shown in Fig. 3 of the main text. Figure 13 shows that if local dissipation is taken into account, spin squeezing is reduced but the ratio between the Wineland parameters for even and odd N remains large. Moreover, for a fixed value of the squeezing parameter r, the ratio of the Wineland parameters in the presence of local dephasing can even be larger than the corresponding ratio obtained for γ φ = 0. For fixed local dissipation rates, the ratio is largest around the onset of the even-odd effect. At this optimum squeezing parameter r max , effective single-spin cooperativities much larger than unity, Γ/γ φ 1 or Γ/γ rel 1, are required to observe a ratio of the Wineland parameters greater than two.
Appendix E: Liouvillian perturbation theory of the slow timescale In this Appendix, we use Liouvillian perturbation theory [61] to analyze the emergence of the long relaxation timescale in the presence of local dephasing, which has been discussed in Sec. V A of the main text. We also provide a simple physical argument to understand this effect.
1. Hilbert space of N spin-1/2 systems and permutational invariance Addition of angular momenta of N spin-1/2 systems gives rise to N/2 + 1 subspaces of total angular momentum j, where j takes values between j max = N/2 and j min = 0 (1/2) if N is even (odd) [105]. For N > 2, all but the maximum-angularmomentum subspace are degenerate since there are multiple ways to combine N spin-1/2 systems to a total angular momentum j < N/2 [33] (for an illustration, see, e.g., Ref. 106). If local dissipative processes act identically on each spin-1/2 system, the equations of motion are invariant under permutation of the spins [34]. Consequently, if the system is initialized in a permutationally invariant state, e.g., any state in the subspace j = j max , the collective and local dissipative processes will preserve the permutational symmetry. Exploiting this symmetry, one can derive an effective quantum master equation which requires significantly less degrees of freedom to describe the system [34] and gives rise to efficient numerical simulation of large spin ensembles [107].

Analysis of the slow timescale
Our starting point, the quantum master equation (11) of the main text, belongs to the class of permutationally invariant systems described above. In the following, we focus on the case γ coll = γ rel = 0. Introducing the dimensionless time τ = Γt, the equation can be rewritten in the form dρ/dτ = L 0ρ + εL 1ρ , where we introduced the dimensionless superoperators and the dimensionless perturbation strength ε = 2γ φ /Γ. In the absence of local dephasing, ε = 0, the superoperator L 0 has N/2 + 1 different steady statesρ (j) 0 , each of them living in a different subspace of collective angular momentum j. Weak local dephasing, γ φ Γ, enables incoherent transitions between adjacent angular-momentum subspaces [34], which can be visualized as trajectory in a triangular (j, m) state space [108]. This perturbation lifts the degeneracy of the steady states and opens a new dissipative gap that determines the relaxation timescale towards the new, unique steady state.
The first-order corrections to the vanishing eigenvalues of L 0 are given by the eigenvalues of the tridiagonal matrix containing the transition rates j → j between collective angular momentum subspaces. Here,1 (j) is the identity operator in the angular-momentum subspace j. For even N , the transition rates are shown in Fig. 14(a). They depend on the structure of the dark state (6) given in the main text, where c (j) m are the expansion coefficients of the dark state and Γ (5,6) j;m,m are the transition rates derived in Ref. 34 using the notation introduced in Ref. 107. Note that, in our case, L 1 is dimensionless such that the transition rates (E4) and (E5) are dimensionless, too. The asymptotic decay rate, i.e., the absolute value of the smallest gap in the spectrum of M , is shown in Fig. 14(b).
For r = 0, the unperturbed steady states are the ground states of each angular-momentum subspace, ρ For r = 0, transition rates Γ j→j−1 are nonzero and dominate over the rate Γ N/2−1→N/2 if the condition r > 1/ √ N holds. As a consequence, an initial state in the maximum-angular-momentum subspace j = N/2 will undergo a directed hopping process towards lower angular momentum subspaces until it reaches a subspace j 0 where "downward" and "upward" rates are balanced, Γ j0→j0−1 ≈ Γ j0−1→j0 . Note that the downward rates Γ j→j−1 are almost constant as a function of j whereas the upward rates Γ j→j+1 depend strongly on j, as shown in Fig. 14(a). The asymptotic decay rate towards the steady state is proportional to 1/N if the transition rates in the vicinity of the equilibrium point j 0 scale proportional to 1/N . Inspection of the rates Γ (5,6) j;m,m listed in Ref. 107 shows that this is the case if j N/2 and m 0. For a given squeezing parameter r, these conditions can be fulfilled if N is sufficiently large, as shown in Fig. 14(b). Numerically, we find an exponent a ≈ 5, see inset of Fig. 14(b).
In the limit r → ∞, the asymptotic decay rate converges to the constant value 1/2.

Physical argument for the slow timescale
The existence of a bottleneck relaxation rate causing a 1/N scaling of the asymptotic decay rate for local dephasing can be understood by an intuitive argument. To explain it, we focus on the transition rate Γ N/2−1,−N/2+1→N/2,−N/2+1 , which is the bottleneck determining the asymptotic decay rate in the limit r < 1/ √ N . The states that are involved in this transition can be parametrized as [109,110] where p ∈ {0, . . . , N − 1}. Here, |j denotes the Nparticle state where the jth spin is in the excited state and all others are in the ground state. The p = 0 state has total angular momentum j = N/2, i.e., we can identify it with the state in the maximum-angular-momentum subspace. In contrast, the N − 1 states with p > 0 have total angular momentum j = N/2 − 1. Therefore, the index p > 0 allows us to label the N − 1 degenerate states in the j = N/2 − 1 subspace, |p ≡ |N/2 − 1, −N/2 + 1, p for p > 0 .
Local dephasing of spin n changes one sign in the superposition (E7), 1 2σ and thus creates an overlap between the orthogonal states |0 and |p > 0 that is proportional to 1/N , For identical dephasing processes on all N spins and for a collective initial state, i.e., a uniform statistical mixture of all N − 1 states |N/2 − 1, −N/2 + 1, p , the total upward transition rate between the two collective angular momentum subspaces is which is the bottleneck of the relaxation process and features the 1/N scaling with system size. Note that the corresponding downward rate is of the order of unity because we have to sum over all N − 1 possible target states |p > 0 , too, Also note that local relaxation does not lead to a similar emergence of the slow timescale because the overlap corresponding to Eq. (E11) will only be proportional to 1/ √ N and is thus canceled by the summation performed in Eq. (E12).

Appendix F: Optimal parameters in master equations
In this section, we show how the optimal protocol parameters vary as a function of increasing system size N in simulations from Fig. 5 of the main text. In the case of the dissipative protocols, the optimization included varying both r as well as κ sqz , whereas in the case of OAT, the spin-cavity detuning ∆ c (see App. G) is was varied. The results are shown in Fig. 15.
Appendix G: Effective One-Axis-Twist quantum master equation In this Appendix, we present the effective model that we consider when discussing the OAT protocol both in the main text and in App. C 2. In particular, following [27,35], we envision an ensemble of spins dispersively coupled to a bosonic cavity. After adiabatically eliminating the cavity, the spin-only quantum master equation can be approximated by [35] with and with ∆ c representing the cavity-spin detuning, g the cavity-spin coupling strength, κ int the decay rate of the cavity, and γ rel the local spin decay. We point out that we assume in the simulations that ∆ c is a tunable parameter, over which we optimize in order to maximize the amount of spin squeezing that the protocol can achieve.
Appendix H: Scaling of the Wineland parameter ξ 2 R in the limit η rel → ∞ When analyzing the performance of the dissipative spin-squeezing protocol in the main text, as one means of implementation, we envisioned engineering the required dissipator by coupling a spin ensemble to a lossy cavity that in turn interacts with an appropriately engineered squeezed bath. Furthermore, in our cooperativity scaling analysis (see Sec. IV and App. C) we investigated the limit of weak single spin cooperativity η rel ≤ 1. It is also interesting to consider a different asymptotic regime, where the internal cavity loss κ int is negligible, giving an extremely large η rel . We focus on the specific case where κ int = 0, and the only undesired dynamics is due to single-spin relaxation at a rate γ rel . Such a situation could be realized without any cavity, by directly irradiating an ensemble of two-level atoms with squeezed light. While this situation was analyzed in Refs. [22,24] the impacts of single spin relaxation were not studied.
The master equation in our chosen limit is thuṡ The key dimensionless parameter that describes the competition of the desired collective dissipative dynamics and the unwanted relaxation is Once again concentrating our attention on the large-N limit and fixing Ŝ z = −N/2, we can approximate the Wineland parameter using the mean-field equations of App. B as In the above expression we have already taken the limit r → ∞, which minimizes ξ 2 R . As one would expect, achievable squeezing increases asη gets larger, but more importantly we have that ξ 2 R ∝ 1/N . We can also define a quantity analogous to a collective cooperativity in this simplified system, which then lets us write assumingη ≤ 1 andC 1. In this Appendix, we show that the generalized model (29), describing an engineered reservoir that stabilizes an impure steady state, can capture the impact of collective excitation and relaxation if the squeezing parameter r and the effective thermal occupation number n th are adjusted properly. To proof this, we start with the fully general quantum master equationρ = ΓD Σ ρ + Γ D Σ † ρ which can model, e.g., interaction with an impure squeezing reservoir and a finite-temperature collective-decay reservoir if the ratios Γ/Γ and γ/γ are chosen properly. Equation (I1) is equivalent to a quantum master equation of the form given in Eq. (29),ρ =ΓD Σ ρ +Γ D Σ † ρ , where we defined a new squeezing parameterr and new decay ratesΓ,Γ as follows: 1 tanh(2r) = 1 tanh(2r) Γ = (Γ + Γ ) cosh(2r) + γ + γ 2 cosh(2r) Condition (I4) can be satisfied for arbitrary nonnegative rates Γ, Γ , γ, and γ . Collective excitation or decay, γ = 0 or γ = 0, respectively, will decrease the squeezing parameter, i.e., we always haver ≤ r. Sufficient (but not necessary) conditions to obtain nonnegative decay ratesΓ andΓ are Note that these conditions are satisfied if γ = γ , which is the case for the trapped-ion implementation dicussed in Sec. VII. There, we have Γ = 0 and γ = γ = γ heat Γ (see Sec. VII). Expanding the general results (I4) to (I6) in the small parameter γ heat /Γ, we find This corresponds to an impure squeezed reservoir with reduced squeezing parameterr < r, effective thermal occupation number n th = cosh(2r)γ heat /Γ, and unchanged decay rate Γ, as shown in Eqs. (32) and (33).
Similarly, for a perfect squeezing reservoir and small zero-temperature collective decay, i.e., Γ = i.e., the presence of collective decay can be understood in terms of an impure squeezing reservoir with decreased squeezing parameterr < r, effective temperature n th , and an enhanced decay rate Γ + γ coll . Ratio between the minimum steady-state Wineland parameter obtained by solving the mean-field equations of motion for the system governed by Eq. (K1), and the ideal Heisenberg-limited value ξ 2 R,HL = 2/(N + 2). Each dot corresponds a solution with optimized squeezing strength r, while the dashed curves are fits to a horizontal line (taken over the largest four values of N for each n th ).
of the main text, which we reproduce here for completeness: ρ = Γ(n th + 1)D Σ ρ + Γn th D Σ † ρ. (K1) Following App. B, we can once again write the corresponding mean-field equations setting third cumulants to zero, which can be readily solved numerically. Of particular interest is the scaling of the Wineland parameter ξ 2 R , with the spin number N . The results are presented in Fig. 17, where we show the ratio between the minimum steady-state Wineland parameter obtained by solving the meanfield equations of motion, and the ideal Heisenberglimited value ξ 2 R,HL = 2/(N + 2) (which full theory predicts at n th = 0). Each dot corresponds a solution with optimized squeezing strength r, while the dashed curves are fits to a horizontal line taken over the largest four values of N for each n th . From the plot it is clear that the mean-field theory predicts ∼ 1/N scaling at large N which is weakly skewed at smallest values of N that we consider. These results are consistent with our full master equation simulations presented in Sec. VI of the main text. Similarly to a simple bosonic theory prediction, we also find that the minimum steady state value of Wineland parameter scales linearly with growing n th . In particular, to a good approximation it satisfies the phenomenological equation ξ 2 R ≈ 2.8(1 + 2n th ) exp(−2r opt ), (K2) with r opt separately optimized for each n th . Finally, we point out that, as can be seen from the blue curve of Fig. 17, the mean-field theory solution does not correctly predict the prefactor of the 1/N in the case of n th = 0 (i.e., where the corresponding ratio ξ 2 R (r opt )/ξ 2 R,HL should be equal to 1). This is a somewhat expected behavior, as in parameter regimes where the ξ 2 R may be either at, or near its Heisenberg limited value, the state of the spin ensemble is not Gaussian, which may substantially lower the ability of mean-field theory to quantitatively describe its behavior.