A hybrid-systems approach to spin squeezing using a highly dissipative ancillary system

Squeezed states of spin systems are an important entangled resource for quantum technologies, particularly quantum metrology and sensing. Here we consider the generation of spin squeezed states by interacting the spins with a dissipative ancillary system. We show that spin squeezing can be generated in this model by two different mechanisms: one-axis twisting and driven collective relaxation. We can interpolate between the two mechanisms by simply adjusting the detuning between the dissipative ancillary system and the spin system. Interestingly, we find that for both mechanisms, ancillary system dissipation need not be considered an imperfection in our model, but plays a positive role in spin squeezing. To assess the feasibility of spin squeezing we consider two different implementations with superconducting circuits. We conclude that it is experimentally feasible to generate a squeezed state of hundreds of spins either by one-axis twisting or by driven collective relaxation.


I. INTRODUCTION
The generation of non-classical states of large quantum systems has attracted significant attention due to the potential of such states in emerging quantum technologies [1,2], such as quantum metrology and sensing [3][4][5][6][7][8]. For instance, it is well known that highly entangled states of N spin-1/2 particles, such as spin squeezed states, can -in principle -be exploited to increase the precision of some measurements by a factor that scales with N 1/2 compared to the best precision that is achievable with a separable state [4-6, 9, 10]. Interestingly, an improvement in the precision is even possible in the presence of certain types of realistic decoherence, although the scaling of the improvement is reduced to N 1/4 [11][12][13]. The motivation of the work described here is the generation of such spin squeezed states, starting from an easily prepared separable state of the spin system. Solid state spin defects, such as nitrogen vacancy centres or electron donor spins in silicon, are particularly promising candidate spin systems due to their long coherence times [14][15][16]. However, to generate entanglement it is clear that we require some sort of interaction between the spins. Although it has been proposed that this can be achieved using the natural magnetic dipoledipole interaction between the spins [17], in practice this is difficult because any spin will interact very weakly with a distant spin (the strength of the dipole-dipole coupling between two spins scales as r −3 where r is the distance between the two spins). Since this interaction is weak it will be challenging to generate highly entangled states within the spin coherence time. Instead, we adopt a hybrid-systems approach where the spins are allowed to interact with an auxilliary system. This interaction with the auxilliary system can induce coupling between the * dooleysh@nii.ac.jp spins (including long range interactions between distant spins) which can then be exploited to generate the necessary entanglement. This approach has been used to experimentally generate few-qubit entangled states of many different systems, for example trapped ions [18,19], Rydberg atoms [20] and superconducting qubits [21,22]. In this context, the auxilliary system is sometimes called a "quantum bus" [23]. However, these experiments are typically limited to few qubit systems. Also, some schemes [7,24] need significant entanglement with the ancillary system throughout the interaction. This means that they are limited by the requirement that the ancillary system must have a coherence time that is longer than the duration of the entanglement. In this paper we consider the interaction of a short-lived ancillary system with a longlived spin system, and we show that this hybrid-systems approach can be used to generate relatively large spin squeezed states. This is a typical feature of the hybridsystems approach: the strengths of both the auxilliary system and the system of interest are exploited to generate dynamics that would be difficult to generate with either system individially [25].
We structure our paper as follows. In section II we give our measure of spin squeezing and we introduce the two spin squeezing mechanisms that are relevant to this paper: (i) spin squeezing by one-axis twisting (OAT) [26,27], and (ii) spin squeezing by driven collective relaxation (DCR) [28,29]. In section III we describe our model and we adiabatically eliminate the ancillary system to obtain an effective master equation for the spin system. We show that both of the spin squeezing mechanisms, OAT and DCR, emerge from these effective dynamics. For concreteness, we focus on two different implementations of the model, one with a superconducting flux qubit playing the role of the ancillary system, and the other with a superconducting microwave resonator. In IV we consider spin squeezing by OAT, including the effect of realistic imperfections in the model, such as dissipation of the ancillary system, inhomogeneity in the spin energies due to fluctuations in their local magnetic fields, and inhomogeneity in the couplings between the ancillary system and the spins. Interestingly, we find that ancillary system dissipation need not be considered an "imperfection" in the model. The spin squeezing is very robust to such decoherence and, perhaps counter-intuitively, moderate dissipation can even improve the spin squeezing by OAT. Such "dissipation-assisted" spin squeezing is an interesting effect since it is unusual for spin squeezing by OAT to be improved by adding dissipation to a part of the system. We also find that the inhomogeneity in the couplings can be reduced to a negligible level by a judicious experimental setup. Inhomogeneity of the spin energies can be compensated by dynamical decoupling. We find that a pulse sequence known as concatenated-XY8 effectively preserves spin squeezing. However, a drawback of pulsed dynamical decoupling is that, in practice, each pulse in a sequence introduces errors that damage the spin squeezing. At the end of section IV we show that driving the spin system enables spin squeezing by OAT without the need for a dynamical decoupling pulse sequence. In section V we consider spin squeezing by DCR, including the effect of realistic imperfections in the model. We show that for squeezing by DCR, standard pulse sequences are not effective in overcoming inhomogeneity in the spin energies, but we present a novel pulse sequence that preserves the spin squeezing. We conclude that it is experimentally feasible to generate squeezed states of hundreds of spins, either by OAT or by DCR. Finally -for our chosen model parameters -we estimate the improvement in precision that this can give in magnetic field sensing.

II. BACKGROUND
Our primary system of interest thoughout this paper is an ensemble of N spin-1/2 particles. The collective spin operators for this system areĴ µ = 1 2 N i=1σ are the Pauli operators for the i'th spin with µ ∈ {x, y, z}. The mean spin vector for an arbitrary state ρ s of the spin system is the expectation value Ĵ = Tr(ρ s Ĵ ) where Ĵ = (Ĵ x ,Ĵ y ,Ĵ z ) is a vector of operators. We denote the unit vector in the direction of the mean spin as n = Ĵ /| Ĵ |. We quantify spin squeezing of a state ρ s with the Wineland squeezing parameter [10], where the minimisation is over all unit vectors n ⊥ that are perpendicular to the mean spin direction and Var( n ⊥ · Ĵ ) is the variance of the operator n ⊥ · Ĵ . For a spin coherent state of the form [30] where |↑ i and |↓ i are the eigenstates ofσ (i) z , we have | Ĵ | = N/2 and min n ⊥ Var( n ⊥ · Ĵ ) = N/4. Hence, by the definition above we have ξ 2 = 1 for a spin coherent state. A state is spin squeezed if ξ 2 < 1, implying that the variance of the operator n ⊥ · Ĵ (for some choice of n ⊥ ) is less than that of a spin coherent state. This is illustrated in Fig. 1 where we plot the Q-function for a spin coherent state and for two example spin squeezed states. It is also known that the squeezing parameter ξ 2 is an entanglement witness, meaning that ξ 2 < 1 implies that the (possibly mixed) state ρ s is entangled. Moreover, the parameter ξ has a specific operational meaning in that it is the ratio of the phase sensitivity of the state ρ s to that of a spin coherent state in a Ramsey inteferometric measurement [10].
There are various possible mechanisms for the generation of spin squeezing starting from a spin coherent state. In this paper we are particularly interested in two of these: (i) spin squeezing by one-axis twisting (OAT) [26,27] and (ii) spin squeezing by driven collective relaxation (DCR) [28,29].
(i) OAT is generated by the evolutioṅ where the Hamiltonian is quadratic in one of the collective spin operators, for example,Ĥ oat = χĴ 2 z . It is well known that for large N this leads to optimum spin squeezing value ξ 2 ∼ N −2/3 after evolution time [26,27]: For the one-axis twisting HamiltonianĤ oat = χĴ 2 z the most spin squeezing is achieved for an initial spin coherent state that is on the equator of the Bloch sphere of each spin, i.e., |θ, φ = π 2 , φ . In Fig. 1(b) we plot the Q-function for a one-axis-twisted state for N = 40.
(ii) Spin squeezing by DCR is induced by the Lindblad master equatioṅ ± are the collective spin raising and lowering operators and the superoperator D is defined as The parameter Ω = (Ω x , Ω y , 0) is a transverse magnetic field applied to the spin system, and γ is the collective spin relaxation rate. Although the model described by Eq. 5 has been well-studied [31][32][33], spin squeezing by this mechanism has been explored only recently [28,29]. It has been shown that any initial spin coherent state |θ, φ relaxes to a steady state [31,32].
For an appropriate value of the transverse field | Ω|, this steady state is squeezed [28]. In contrast, the analagous driven, dissipative dynamics for a bosonic mode leads to steady states that are always coherent states rather than squeezed states [34]. In Eq. 5, the steady state with the most squeezing is achieved for | Ω| of the order N γ [28,29]. In the following we see that both the OAT and the DCR spin squeezing mechanisms emerge in the interaction of an ancillary system with a spin ensemble.

III. MODEL
We model the interaction of a spin ensemble with a dissipative ancillary system by the master equatioṅ whereĤ =Ĥ s (t) +Ĥ anc + H int and: x cos(ωt + η), (7) The operatorÂ = d−2 n=0 √ n + 1 |n n + 1| is a lowering operator for the d-dimensional ancillary system, where |n is a basis for the ancilla state space. If, for example, the ancillary system is a qubit (d = 2), the opera-torÂ is the qubit lowering operatorσ − = (σ x − iσ y )/2, which has the commutation relation [σ + ,σ − ] =σ z . If the ancillary system is a bosonic mode (d → ∞), the operatorÂ is the annihilation operatorâ, with the commutation relation [â,â † ] = 1. The ancillary system frequency is ω anc and its relaxation rate is γ with corresponding relaxation time γ −1 . In Eq. 7 for the spin ensemble, each spin may have a different frequency ω (i) with an averageω = 1 [ω (i) −ω] 2 /N . Also, the spins are driven at the average spin frequencyω by a classical field of amplitude Ω and phase η. The interaction Hamiltonian Eq. 9 describes the coupling between the spins and the ancillary system, where each spin may have a different coupling λ (i) , with an average couplingλ and standard deviation δλ.
Rotating to an interaction frame defined by the unitary we make a rotating wave approximation which gives the master Eq. 6 but with the new Hamiltonian: where ∆ = ω anc −ω is the detuning between the ancillary system frequency and the average of the spin frequencies, n η = (cos η, sin η, 0) is a unit vector in the equatorial plane, and we have definedĤ . In Eq. 10 the Hamiltonian is separated into an "inhomogeneous" part represented by the Hamiltonian termsĤ IB andĤ IC and a "homogeneous" part represented by the remaining terms of Eq. 10. The subscript "IB" onĤ IB stands for "inhomogeneous broadening". If each spin has the same frequency ω (i) =ω (equal to the average value), the inhomogeneous broadening termĤ IB vanishes. Similarly, the subscript "IC" onĤ IC stands for "inhomogeneous couplings", and if each spin is equally coupled to the flux qubit we have λ (i) =λ and the inhomogeneous coupling termĤ IC vanishes.
We note that spin relaxation and spin dephasing have been neglected in the model described above. If the spins are, for example, an ensemble of donor spins in silicon then this is a reasonable assumption since the spin dephasing time is of the order of seconds and the relaxation time is of the order of tens of minutes at low temperatures [16]. Also, long coherence times of ∼ 30 ms have been achieved for ensembles of nitrogen-vacancy centres in diamond by the use of dynamical decoupling [15].

A. Effective dynamics
To see how spin squeezing is generated in this model, we first define the parameter Γ = ∆ 2 + γ 2 /4. If the ancillary system is initially in its ground state |0 and if Γ satisfies the conditions we can adiabatically eliminate the ancillary system. The result is the effective master equation (see Appendix A for details) [35,36]: where ρ s = Tr anc (ρ) is the reduced state of the spin system and γ eff =λ 2 γ/Γ 2 is the collective-spin relaxation rate. The effective Hamiltonian is [37] H eff = Ω n η · Ĵ + χ effĴ 2 z − χ effĴz − χ eff Ĵ · Ĵ , (13) where χ eff = ∆λ 2 /Γ 2 . For clarity we have neglected the inhomogeneous couplings, i.e., λ (i) =λ in Eq. 12 and Eq. 13, although the effect of both inhomogeneous couplings and inhomogeneous broadening will be assessed in later numerics. These effective dynamics have features of both spin squeezing mechanisms that were discussed in section II, that is, squeezing by OAT and by DCR. The term χ effĴ 2 z in the effective Hamiltonian is an OAT term, while the collective relaxation in Eq. 12, in combination with spin drive Ω n η · Ĵ , are the necessary ingredients for squeezing by DCR (for comparison, see Eq. 5). The two spin squeezing mechanisms appear independently in two different regimes of the effective master Eq. 12. The OAT regime emerges when the OAT coefficient is much larger than the collective relaxation rate, χ eff γ eff [35]. By comparing the expressions for χ eff and γ eff it is easy to see that this reduces to the condition that the detuning should be much larger than the ancillary system relaxation rate, ∆ γ. On the other hand, the DCR regime emerges when the collective relaxation dominates the OAT, χ eff γ eff , which corresponds to the condition ∆ γ. We note that the DCR regime includes, for example, the case where the ancillary system and the spins are resonant (∆ = 0), in which case the effective master Eq. 12 is in the same form as Eq. 5.

B. Realistic parameters
To assess the feasibility of spin squeezing we choose some reasonable parameters for our model. We consider two different implementations: (i) a superconducting flux qubit (FQ) coupled to an ensemble of nitrogen-vacancy (NV) centres in diamond, and (ii) a superconducting microwave resonator (MR) coupled to an ensemble of electron donor spins in silicon.

Superconducting flux qubit and nitrogen-vacancy centres
The spin parametersω and δω depend on the type of spin system. Here we take the spins to be NV centres in diamond. Although the NV centre ground state is spin-1 [38], a magnetic field can be applied to detune one of the spin sub-levels so that -to a good approximation -the NV centre can be considered spin-1/2. Due to the NV zero-field splitting, we haveω ≈ 2π × 3 GHz [39], and based on recent experimental results we estimate δω ≈ 2π × 3 kHz [40]. The spin drive parameters Ω and η are experimentally tunable.
For a superconducting FQ, the ancillary system oper-atorÂ in Eqs. 8 and 9 is the qubit lowering operatorσ − . We assume that the FQ is tuned so that its two persistent current states are degenerate, with tunnel splitting ω anc = ω FQ between them. The detuning ∆ = ω FQ −ω can be varied experimentally by changing the flux qubit tunnel splitting ω FQ [41,42].
The coupling strength is determined by the magnetic field B(y i , z i ) that is generated by the FQ at the position (0, y i , z i ) of the i'th NV centre (with axes as shown in Fig. 2). We assume a square FQ of length 3 µm, wire thickness 0.1 µm and wire height 0.2 µm (see figure Fig. 2(a)) and we assume a uniform critical current of I = 1.4 µA. Based on these values, the coupling strength λ(y, z) in the interior of the FQ can be estimated by the Biot-Savart law [43]. This is shown in the contour plot in Fig. 2(a). For NV centres positioned near the middle of the FQ the coupling is relatively homogeneous across a broad area (the blue region in Fig. 2(a)). Assuming that the NV centres are contained in a diamond sample of volume 1.58 × 1.58 × 0.2 µm 3 with NV density 10 15 cm −3 gives a total of N = 500 nitrogen-vacancy centres randomly placed throughout the diamond sample. We find numerically that in this case the average coupling isλ ≈ 2π × 12 kHz with standard deviation δλ ≈ 2π × 1 kHz. We note that coherent coupling between a FQ and an ensemble of NV centres has been demonstrated experimentally with a similar coupling strength [44].
Finally, we assume that the FQ relaxation rate is γ = 1 MHz, corresponding to the relaxation time γ −1 = 1 µs. This is a reasonable estimate, since relaxation times an order of magnitude longer than this have been measured in recent experiments with flux qubits [45,46]. We will find it useful to write both the adjustable detuning ∆ and the relaxation rate γ as a proportion of the collective couplingλN . For the relaxation rate this gives γ = 0.0265 ×λN , whereλ = 2π × 12 kHz and N = 500, as determined above. We use these expressions for ∆ and γ when our numerical simulations are limited to small numbers of spins, N . This is useful because with these expressions the condition Γ = ∆ 2 + γ 2 /4 λ N in Eq. 11 is satisfied by the same proportion for any value of N , and we can extrapolate from our small-N numerical results to our larger estimated value N = 500.

Superconducting microwave resonator and donor spins in silicon
The coupling of a MR to a spin system is much weaker than the coupling of a FQ to a spin system. For this reason, we take the spin system in this case to be an ensemble of electron spins of phosphorus atoms doped in a silicon crystal, since longer coherence times have been measured for these spins than for NV centres [16]. The donor electron spin frequencyω is determined by the electron Zeeman splitting. The donor electron spins interact with the donor nuclear spins via a hyperfine in- teraction of 2π × 118 MHz. However, by polarizing the nuclear spins this interaction can be regarded as a contribution to the electron Zeeman splitting [47]. With an additional externally applied magnetic field ∼ 100 mT we haveω ≈ 2π × 3 GHz [48] (similar to the NV zerofield splitting). Based on experimental results [16], we assume δω = 2π × 15 Hz, which is much smaller than for NV centres. The spin drive parameters Ω and η are experimentally tunable.
For a superconducting MR, the operatorÂ in Eqs. 8 and 9 is the bosonic lowering operatorâ. The detuning ∆ = ω MW −ω can be varied experimentally by adjusting the MR frequency ω MW . The coupling strength λ (i) = g e µ B |B(y i , z i )|/(2 ) is determined by the magnetic field B(y i , z i ) that is generated by the MR at the position (0, y i , z i ) of the i'th donor spin. We assume that the wire of the MR has length 1.5 mm, width 2 µm, height 50 nm (see Fig. 2(b)), a penetration depth of 90 nm, and an inductance L = 1.5 nH. Based on these values, the coupling strength is shown in the contour plot in Fig. 2(b). There is a region of relatively homogeneous coupling directly above the wire (the green area between the two ridges in the contour plot). We suppose that a silicon sample of dimensions 1 mm × 2 µm × 50 nm and donor spin density 1.2×10 14 cm −3 is placed in this region at a distance of 100 nm from the resonator. With these values we estimate that there are N = 1.2 × 10 4 spins placed randomly throughout the silicon sample, with an average couplingλ ≈ 2π × 56 Hz and a standard deviation δλ ≈ 2π × 4 Hz [49]. This is a considerably weaker coupling than for the FQ and NV implementation in the previous section. However, measured coherence times for donor spins in silicon are much longer than those for NV centres [16]. The relaxation rate of the resonator is γ = ω MR /Q ≈ 2π × 0.34 MHz, assuming a resonator quality factor Q = 4.5 × 10 4 and the frequency ω MR ≈ 2π × 3 GHz. Forλ = 2π × 56 Hz and N = 1.2 × 10 4 this relaxation rate can be expressed as γ = 0.1λN . We use this value when our numerical simulations are restricted to small values of N .
In the following sections we numerically investigate the spin squeezing, using realistic parameters as far as possible. For all of our simulations we use the master Eq. 6 with the Hamiltonian Eq. 10, that is, the master equation before the approximations that leads to the effective master Eq. 12. This gives us meaningful results even when the approximation conditions in Eq. 11 are not well satisfied.

IV. SPIN SQUEEZING BY OAT
A. Ideal case, with ancilla relaxation First, we consider spin squeezing in the OAT regime Γ ≈ ∆ γ, assuming that the initial state is the spin coherent state |θ, φ = π 2 , 0 , since this state is in the class of states π 2 , φ that leads to the most squeezing by the OAT term χ effĴ 2 z . To prepare the state π 2 , 0 we make use of the fact that a general spin coherent state |θ, φ =R(θ, φ) |0, 0 can be prepared by rotating each spin from the state [30]. This rotation can be implemented by applying an electomagnetic pulse to the spin system. The state |θ, φ = |0, 0 is itself easily prepared, e.g., by cooling (it is the ground state of the spin Hamiltonian Eq. 7 when Ω = 0) or, for NV centres, by optical pumping [39,50]. After the state π 2 , 0 =R( π 2 , 0) |0, 0 has been prepared, unitary evolution byĤ eff leads to spin squeezing.
As mentioned in the previous section, the detuning ∆ and the spin drive Ω are experimentally tunable parameters in both of the considered implementations. To get a comprehensive picture of which values of these parameters lead to spin squeezing we plot min t ξ 2 , the minimum spin squeezing that is achieved across all evolution times t, as a function of ∆ and Ω. This is shown in Fig. 3(a) for the FQ and NV model assuming (for the moment) the ideal case where there is no flux qubit relaxation (γ = 0), no inhomogeneous broadening (δω = 0) and homogeneous coupling of the spins to the flux qubit (δλ = 0). Fig. 3(a) shows that there is significant spin squeezing (the dark red region) in the lower right portion of the plot where Γ ≈ ∆ Ω and Γ ≈ ∆ λ N , corresponding the the regime of validity of the effective Hamiltonian Eq. 13. In Fig. 3(b) we include a realistic amount of flux qubit relaxation (γ = 0.0265×λN ) and we see that the spin squeezing is very robust to this kind of decoherence. The corresponding plots for the resonator model are not shown as they are qualitatively similar to Figs. 3(a,b).
In the OAT regime (∆ γ) the evolution time required to reach the optimal spin squeezing is given by Eq. 4. From this equation we estimate where we have substituted the expression for the OAT coefficient χ eff and we have used Γ ≈ ∆. It may appear from this expression for t opt that the optimum squeezing time decreases with the number of spins N , but this is not so, since we require Γ ≈ ∆ λ N for the effective Hamiltonian Eq. 13 to be valid. For a detuning ∆ = kλN λ N for some k 1, this translates to an optimal squeezing time t opt = 3 1/6 kN 1/3 /λ. We see that the optimum squeezing time actually increases as the number of spins N increases. However, the scaling is N 1/3 so that t opt is not too large for moderate values of N . It was determined above that N = 500 was a realistic number of NV centres that could be coupled to the FQ for our chosen FQ dimensions. In this case, substituting ∆ = 20λN andλ = 2π × 12 kHz, we estimate t opt = 2.5 ms, which is within the spin coherence times achieved in recent experiments with ensembles of nitrogen-vacancy centres in diamond [15]. In Fig. 3(c) the solid black line shows the time evolution of the spin squeezing parameter, assuming the large detuning ∆ = 20λN , zero ancillary system relaxation γ = 0, and zero spin drive Ω = 0. For comparison, the dashed black line shows the spin squeezing for relaxation γ = 0.0265 ×λN . The spin squeezing is almost indistinguishable from the γ = 0 case, confirming that this spin squeezing mechanism is very robust to realistic levels of ancillary system relaxation. Also, the horizontal dotted black line shows the level of optimum spin squeezing ξ 2 ∼ N −2/3 for perfect OAT, i.e., the optimum squeezing by the effective Hamiltonian Eq. 13 (or, alternatively, by Eq. 3). The minimum of the solid black line reaches this optimum since the dynamics are well approximated by the effective Hamiltonian for the large detuning ∆ = 20λN . For the microwave resonator and donor spins in silicon implementation, we estimated N = 1.2 × 10 4 andλ = 2π × 56 Hz. Substituting these values, along with ∆ = 20λN , gives t opt = 1.6 s, which is within the spin coherence time ∼ 10 s of phosphorus donor spins in silicon at low temperatures with spin echo [16].
Although these optimum squeezing times are within the achievable spin coherence times in both implementations, it is desirable to decrease t opt . Interestingly, we observe that for a given number of spins we can significantly reduce the optimum squeezing time as follows. Since the optimum squeezing time t opt scales linearly with the detuning ∆ (see Eq. 14), the spins can be squeezed more quickly by decreasing the detuning. Decreasing the detuning comes at a cost, however: we need ∆ to be large enough to satisfy our approximation condition Γ ≈ ∆ λ N . This leads to a tradeoff: if we decrease the detuning we can squeeze more quickly at the expense of a worse approximation to the effective Hamiltonian, and conversely, if we increase the detuning we have a better approximation but at the expense of a longer wait for the optimal squeezing. This can be seen by comparing the solid black line and the solid red line in Fig. 3(c) for the FQ and NV model. The detuning ∆ = 2λN for the solid red line, is an order of magnitude smaller than ∆ = 20λN for the solid black line, so that the optimum spin squeezing is achieved an order of magnitude faster. However, the minimum spin squeezing min t ξ 2 is degraded compared to the optimal squeezing (the horizontal black dotted line). This difference between the optimum and the minumum of the solid red curve shows the importance of using the full master Eq. 6 instead of the effective master Eq. 12 for smaller values of ∆. Interestingly, if ∆ is not quite large enough to satisfy ∆ λ N , flux qubit relaxation can significantly improve the approximation so long as the effect of collective relaxation on the OAT is still negligible, i.e, provided that ∆ γ. This is shown by the dotted red line in Fig. 3(c) for the FQ and NV model with γ = 0.0265 ×λN . With this realistic amount of flux qubit relaxation the squeezing can be significantly improved compared to γ = 0 (the solid red line). This improvement of the spin squeezing by ancillary system relaxation may be surprising on first sight, since in most models any form of decoherence is an unwanted influence on the dynamics. However, the effect of the relaxation is to suppress excitation of the ancillary system. This, in turn, inhibits entanglement between the ancillary system and the spin system, which would be damaging to the spin squeezing. If we choose ∆ = 2λN the optimum squeezing time for N = 500 spins with the FQ and NV implementation is reduced to t opt = 250 µs. Similarly, the optimum squeezing time for the MR and donor spins implementation is decreased to t opt = 160 ms.

B. Realistic case, with dynamical decoupling
We now consider the effect of inhomogeneous broadening and inhomogeneous couplings on the OAT spin squeezing. In this case, since the state space dimension increases exponentially in the number of spins, N , our numerics are restricted to a small number of spins, N = 6. In Fig. 4(a) the solid red line shows the spin squeezing including the effect of inhomogeneous cou-plingsĤ IC and inhomogeneous broadeningĤ IB in the FQ and NV model, using the value δλ = 2π × 1 kHz that was estimated for the standard deviation of the couplings λ (i) , and spin frequencies ω i chosen at random from a Gaussian distribution with standard deviation δω = 2π × 3 kHz. The dynamics are averaged over 100 evolutions to remove fluctuations due to the randomness of the ω i . Similarly, the solid red line in Fig. 4(b) shows the spin squeezing including inhomogeneities with the MR model parameters. For both implementations, we see that the spin squeezing is degraded, although there is a small amount of spin squeezing at very short times for the MR model parameters in Fig. 4(b). Further numerics have shown that the decay of spin squeezing is primarily due to the inhomogeneous broadening term H IB rather than inhomogeneous coupling termĤ IC . In fact, for the relatively homogeneous couplings achieved in the setups illustrated in Fig. 2, the inhomogeneous coupling termĤ IC can be safely neglected for timescales of interest. To understand the damaging effect of the inhomogeneous broadening we note that the Hamiltonian H IB causes each spin to evolve around its Bloch sphere at a different rate determined by the frequency ω i −ω. For a Gaussian distibution of ω i , this dephasing leads to Gaussian decay of | Ĵ |, the magnitude of the mean spin vector, with a decay time (δω) −1 . Since the Wineland squeezing parameter ξ 2 is inversely proportional to | Ĵ | 2 , decay of the mean spin vector leads to an increase in the squeezing parameter.
It is well known that for evolution byĤ IB , a single πpulse at some time τ leads to a spin echo at the time 2τ , where the π-pulse is an instantaneous rotationR(π, φ) = exp[−iπ(Ĵ x sin φ −Ĵ y cos φ)] by an angle π about an axis on the equator of the Bloch sphere of each spin. Crucially, the π-pulse also commutes with the OAT operatorĴ 2 z so that a π-pulse at some time τ has the effect of undoing the inhomogeneous broadening at time 2τ without affecting the spin squeezing by OAT [51]. To see this we consider the effective master Eq. 12 in the OAT regime ∆ γ. Assuming Ω = 0, δλ = 0 and neglecting collective relaxation, the spin system evolves bẏ If at time τ we apply the π-pulse operatorR(π, φ), the state is transformed to ρ s (τ ) =R(π, φ)ρ s (τ )R † (π, φ) and the evolution equation for the following period of time t > τ is: Operating on Eq. 16 on the left byR † (π, φ) and on the right byR(π, φ) gives, for t > τ , the evolution equation: where we have usedR † (π, φ)Ĥ IBR (π, φ) = −Ĥ IB ,

FIG. 4. For both (a) the FQ and NV implementation, and (b) the MR and donor spins implementation, the spin squeezing is badly degraded by inhomogeneous broadening
(the solid red lines) but can be completely recovered by dynamical decoupling with a concatenated-XY8 pulse sequence (dashed red lines). The XY8 pulse sequence is C =πx-πy-πx-πy-πy-πx-πy-πx, whereπx represents the π-pulseR(π, π/2) = exp[−iπĴx],πy represents the π-pulsê R(π, 0) = exp[iπĴy], and each dash represents free evolution for a time τ . The concatenated-XY8 pulse sequence is [15] C-πx-C-πy-C-πx-C-πy-C-πy-C-πx-C-πy-C-πx, a sequence of 72 π-pulses in total, ending atλt ≈ 50 in both (a) and (b). Each line in the figures is averaged over 100 evolutions to remove random fluctuations. Horizontal dotted lines show the optimum spin squeezing by OAT. (Fig. (a)  and the important propertyR † (π, φ)Ĵ 2 zR (π, φ) =Ĵ 2 z . Comparing Eq. 15 and Eq. 17 shows that the effect of the π-pulse is to reverse the sign of the inhomogeneous broadening HamiltonianĤ IB in the following period of evolution, without changing the OAT operatorĴ 2 z . Eqs. 15 and 17 are easily solved to give the combined unitary evolution operator: for times t > τ . The operators in the two exponents above commute, so that at t = 2τ we havê We see that at this time there is a spin-echo (the inhomogeneous broadening HamiltonianĤ IB has been cancelled), but that the OAT squeezing is unaffected. However, in a real system the higher order terms in the effective Hamiltonian will also contribute to the dynamics.
We have performed numerics that show that for evolution by the Hamiltonian Eq. 10, which includes higher-order terms, a single π-pulse does not completely recover the spin squeezing even if the approximation conditions Eq. 11 are satisfied. This is because the single π-pulse does not refocus the spin dephasing due to higher order inhomogeneous terms in the effective dynamics. To fully preserve the spin squeezing a more complicated dynamical decoupling pulse sequence is required. We have tried various pulse sequences numerically and found that spin squeezing can be completely recovered for a sequence of alternatingR(π, 0) andR(π, π/2) pulses, as shown in Fig.  4 (the dashed red line). The pulse sequence, known as concatenated-XY8, has recently been implemented experimentally to increase the coherence time of an ensemble of nitrogen-vacancy centres to ∼ 30 ms [15]. Finally, we note that this model for OAT has been discussed by previous authors for an ancillary bosonic mode [35,51]. However, compared to previous work, we highlight the robustness of spin squeezing to ancillary system dissipation, and we demonstrate that ancillary system dissipation can even play a positive role in the generation of spin squeezed states. We have also demonstrated the feasibility of spin squeezing in the two implementations considered, and we have shown that inhomogeneous broadening can be overcome by the concatenated-XY8 pulse sequence.

C. Realistic case, without dynamical decoupling pulses
A practical challenge in the spin squeezing method outlined above is the application of accurate dynamical decoupling pulses to the spin system. For example, we have assumed that each π-pulse has no errors and that it can be implemented instantaneously. In reality, however, a π-pulse cannot be implemented instantaneously, and if there are errors in each of the pulses in a sequence, these errors may accumulate, leading to degradation of the spin squeezing. Moreover, the preparation of the spin coherent state |θ, φ = π 2 , 0 requires a pulse that rotates each of the spins equally. If each spin is rotated by a slightly different angle, this introduces inhomogeneous broadening to the system and damages the spin squeezing. In this section we suggest an alternative approach that generates spin squeezing by OAT without the need for dynamical decoupling pulses and starting from the spin coherent state |θ, φ = |0, 0 that can be prepared without applying a pulse to rotate the spins.
To show how this works, we derive another effective Hamiltonian in the OAT regime (Γ ≈ ∆ γ) starting fromĤ eff (Eq. 13). First, we rotateĤ eff to an interaction frame defined by the unitary transformation In the rotating frame the Hamiltonian iŝ If the parameter Ω is large enough to satisfy the condition Ω N χ eff we can make a rotating wave approximation by neglecting quickly oscillating terms in Eq. 20. The resulting effective Hamiltonian is: Intuitively, this is the effective Hamiltonian that results from averagingĤ eff over rapid oscillations around the n ηaxis due to the large drive Ω. Compared to the effective HamiltonianĤ eff (Eq. 13), the key feature ofĤ eff (Eq. 21) is that the OAT term is ( n η · Ĵ ) 2 rather thanĴ 2 z . For instance, if η = 0 we have ( n 0 · Ĵ ) 2 =Ĵ 2 x , or if η = π/2 we have ( n π/2 · Ĵ ) 2 =Ĵ 2 y . Regardless of the value of the phase η of the driving field, preparation of the spin coherent state |θ, φ = |0, 0 will lead to the most spin squeezing with this OAT term. This is convenient because the state |θ, φ = |0, 0 can typically be prepared by cooling or by optical pumping, without the need for electromagnetic pulses to rotate the spins. Interestingly, the ability to easily change the axis of the spin squeezing term ( n η · Ĵ ) 2 by changing η can be exploited to increase the spin squeezing to the Heisenberg limit using optimal control techniques [52]. We note that in our scheme, this does not require control pulses as in [52,53], but can be achieved by simply shifting the phase η of the spin driving field.
By substituting the expression for χ eff into the approximation condition Ω N χ eff that leads to Eq. 21, the condition becomes where on the right hand side we have used the approximation Γ ≈ ∆, which is valid in the OAT regime. Since (from Eq. 11) we also have Γ ≈ ∆ Nλ, the condition Eq. 22 is easily satisfied for spin drive Ω comparable to (or greater than) the coupling strengthλ. We note, however, that although Ω should be large we must maintain Γ ≈ ∆ Ω to ensure consistency with our previous approximation conditions in Eq. 11.
We now present some numerical results. Neglecting inhomogeneous broadening and inhomogeneous couplings, we plot the minimum spin squeezing min t ξ 2 as a function of the experimentally adjustable parameters ∆ and Ω for the initial state |θ, φ = |0, 0 and for a realistic amount of ancillary system relaxation γ. This is shown for the FQ and NV model in Fig. 5(a). The correspoding plot for the MR model is not shown as it is qualitatively similar. We see in Fig. 5(a) that there are a wide range of values of ∆ and Ω that give significant spin squeezing. In Fig. 5(b) we plot the time evolution of the spin squeezing parameter ξ 2 for several choices of the detuning ∆ and the flux qubit drive Ω, again for the FQ and NV model parameters. Interestingly, even if the condition ∆ λ N is not well-satisfied, e.g. for ∆ = 2Nλ, it is still possible to achieve a level of spin squeezing that is comparable to the optimal squeezing (the horizontal dotted black line in Fig. 5(b)). This is somewhat surprising since we expect such a decrease in ∆ to damage the spin squeezing, as in Fig. 3(c). However, since the spin squeezing is not degraded, it is preferable to choose the smaller value ∆ = 2λN as the squeezing dynamics are faster in this case. With this in mind, we estimate the optimum squeezing time for the FQ and NV implementation, with ∆ = 2λN , N = 500 andλ = 2π × 12 kHz. We find that t opt = 500 µs. This is a factor of two longer than the corresponding time for the OAT dynamics in section IV A because the OAT coefficient χ eff /2 inĤ eff (Eq. 21) is a factor of two smaller than the OAT coefficient inĤ eff (Eq. 13). Similarly, in the microwave resonator model, the optimum squeezing time for N = 1.2 × 10 4 spins and detuning ∆ = 2λN is estimated to be t opt = 320 ms.
Finally, we consider inhomogeneous broadening and inhomogeneous couplings. If, in addition to Eq. 22, we have Ω δω, the inhomogeneous broadening Hamiltonian in the interaction frame,Û † (t)Ĥ IBÛ (t), is quickly oscillating and is suppressed in the rotating wave approximation. This is plotted for the FQ and NV device in Fig.  6(a) where the dashed lines show the spin squeezing for N = 6 spins interacting with a dissipative flux qubit with δλ = 2π × 1 kHz and δω = 2π × 3 kHz. For ∆ = 20λN (the dashed black line), the spin squeezing is slightly degraded compared to the ideal spin squeezing (the solid black line). This is due to higher order inhomogeneous terms that are not suppressed by the spin drive. However, for smaller detuning ∆ = 4λN the spin squeezing is achieved more quickly and inhomogeneous broadening is almost completely suppressed (the dashed red line). Fig.  6(b) shows the corresponding plot for the MR and donor spins implementation.

V. SPIN SQUEEZING BY DCR
In this section we consider spin squeezing in the DCR regime of our effective master Eq. 12, that is, when ∆ γ ≈ Γ. To easily access this parameter regime we assume that the ancillary system is highly dissipative with γ = 20λN . We note that this value of ancillary system relaxation is almost three orders of magnitude stronger than the value γ = 0.0265 ×λN that we used in the previous section for the FQ and NV model, and 200 times larger than γ = 0.1λN for the MR and donor spins model.

A. Ideal case
Again, we begin by neglecting inhomogeneous broadening and inhomogeneous couplings. In Fig. 7(a), for the FQ and NV model, we plot min t ξ 2 as a function of the detuning ∆ and the spin drive Ω for the easily prepared initial state |θ, φ = |0, 0 . We see that for a range of values of ∆ and Ω there is significant spin squeezing during the evolution. In Fig. 7(b) we plot the time evolution of the spin squeezing for various choices of ∆ and Ω. As expected [28,29], we see steady state spin squeezing for a carefully chosen value of the spin drive Ω (the black line, Fig. 7(b)). We have verified numerically that for these parameters the steady state of the master Eq. 6 with the Hamiltonian Eq. 10 is indeed squeezed.

B. Realistic case, dynamical decoupling
We now consider the effect of inhomogeneous broadening and inhomogeneous couplings on DCR spin squeezing. In this case our numerics are limited to a small number of spins N = 6. For simplicity we also assume that ∆ = 0, i.e., the ancillary system and the spin system are resonant, so that the effective Hamiltonian Eq. 13 only includes the spin drive termĤ eff = Ω n η · Ĵ , and the effective master equation is of the form Eq. 5. For the FQ and NV model, the red line in Fig. 8(a) shows that for δω = 2π × 3 kHz and δλ = 2π × 1 kHz the spin squeezing is quickly degraded. The red line in Fig. 8(b) shows that the inhomogeneities have a similar effect on spin squeezing for the parameters in the MR model. Further numerics have shown that, as with OAT, this is primarily due to the inhomogeneous broadeningĤ IB . Unfortunately, the dynamical decoupling approach that was taken to protect spin squeezing against inhomogenous broadening in the previous section for OAT will not work for DCR. This is because for spin squeezing by OAT, the π-pulse operator R(π, φ) = exp[−iπ(Ĵ x sin φ−Ĵ y cos φ)] has the convenient property that it commutes with the OAT operatorĴ 2 z so thatR † (π, φ)Ĵ 2 zR (π, φ) =Ĵ 2 z and the spin squeezing is not disrupted by the π-pulse. For DCR, however, the π-pulse operatorR(π, φ) will disrupt the DCR squeezing mechanism. To see this we start from the effective master Eq. 12 in the DCR regime ∆ γ (assuming ∆ = 0 and δλ = 0): If at time τ we apply the π-pulse operatorR(π, φ), the state is transformed to ρ s (τ ) =R(π, φ)ρ s (τ )R † (π, φ) and the master equation for the following period of time t > τ is:ρ Operating on Eq. 24 on the left byR † (π, φ) and on the right byR(π, φ) gives (for t > τ ) the evolution equation: where we have usedR † (π, φ)Ĥ IBR (π, φ) = −Ĥ IB and R † (π, φ)Ĵ −R (π, φ) = −e −2iφĴ + . Comparing Eq. 23 and Eq. 25 shows that the effect of the π-pulse is to reverse the sign of the inhomogeneous broadening Hamiltonian H IB in the following period of evolution. However, unlike for OAT, the DCR spin squeezing mechanism is also disrupted by the π-pulse, since the collective relaxation operator is transformed fromĴ − toĴ + . For example, if the system was in the steady state of the DCR master equation before the π-pulse, then after the π-pulse it will be far from the steady state. However, the desired effect can be achieved by a reflection of each spin at time τ through a plane in the Bloch sphere that contains the z-axis and the vector n η . Such a reflection is implemented by the complex conjugation operatorV , with the assumption that for each spin the matrix elements ofσ (i) z , n η · σ (i) and n η+π/2 · σ (i) are: This gives:Vσ Applying the complex conjugation operator to the spin state at time τ transforms the state to ρ s (τ ) = V ρ s (τ )V −1 . The master equation for the following period of evolution is: Operating on Eq. 30 on the left byV −1 and on the right byV gives, for t > τ , the evolution equation: where we have Comparing Eq. 23 and Eq. 31 we see that the sign of the inhomogeneous broadening HamiltonianĤ IB is reversed and the Lindblad term is unchanged, as desired. However, the spin drive term is also transformed from Ω n η · Ĵ to −Ω n η · Ĵ . This can easily be corrected by shifting the phase of the spin drive η → η + π so that −Ω n η · Ĵ → Ω n η · Ĵ , which finally gives: Eq. 32 is identical to Eq. 23, apart from a reversal of the inhomogeneous broadening. We note that this operation must be repeated many times, with a short free evolution time τ between each operation, since the dynamics before the reflection and phase-shift does not commute with the dynamics afterwards. Finally, we note that the reflection is an unphysical transformation, since the complex conjugation operator is anti-unitary. However, for some states, we can apply rotations that have the same effect as a reflection. For example, reflecting each spin of the spin coherent state |θ, π/2 in the xz-plane of its Bloch sphere gives the state |θ, −π/2 . Clearly, this transformation can also be implemented by a rotationR(−2θ, π/2) = exp[i2θĴ x ] of each spin around the x-axis. The angle of rotation 2θ depends on the spin coherent state parameter θ. For simplicity, we assume that the spin system is prepared in a spin coherent state that is "close" to the steady state. We take this to be the state |θ ss , φ ss where the angle θ ss = cos −1 (−2 Ĵ z ss /N ) is determined by the expectation value Ĵ z ss in the steady state and φ ss = tan −1 ( Ĵ y ss / Ĵ x ss ) = η + π 2 . This simplifies the procedure because the rotation that mimics the reflection of the state does not have to be changed as the system evolves. The rotation isR(−2θ ss , η), which, for example, transforms the state θ ss , η + π 2 to θ ss , η − π 2 , the reflection of the state through the plane that contains the z-axis and the vector n η . The result of repeating this operation many times is plotted in the dashed red lines, Fig. 8. We see that the inhomogeneous broadening can be significantly suppressed by this procedure and that some spin squeezing can be recovered. We note that it may be possible to get a further improvement by alternating reflections of the state in the plane containing the z-axis and n η with those in the orthogonal plane containing the z-axis and n η × z. This would be analagous to the alternating π-pulses around two orthogonal axes in the concatenated-XY8 pulse sequence in section IV B.

VI. CONCLUSION
We have shown that a single model -the interaction of a spin system with a dissipative ancillary system -can lead to spin squeezing by two distinct mechanisms: oneaxis twisting and driven collective relaxation. In either case, spin squeezing is generated even though the ancillary system coherence time is much smaller than the duration of the squeezing process. This is possible because we have adiabatically eliminated the ancillary system, which stays close to its ground state throughout the squeezing process. We focus on two possible implementations, with either a superconducting flux qubit or a superconducting microwave resonator playing the role of the ancillary system. With dynamical decoupling we have shown numerically that both spin squeezing mechanims are robust to inhomogeneities in the model. In practice, the dynamical decoupling pulses may introduce errors that reduce the spin squeezing. However, we have also shown that by driving the spin system it is possible to generate robust OAT spin squeezing without the need for electromagnetic pulses. This flexibility -squeezing can be generated in disparate parameter regimes and for a variety of practical requirements -is, we believe, a strength of this model. We conclude that spin squeezing of hundreds of solid state electron spins should be experimentally feasible in this model with current or near future technology. Finally, we estimate the sensitivity of a magnetic field measurement that can be achieved using a spin system prepared in a squeezed state. We assume that the squeezed state is prepared by OAT, since this leads to more spin squeezing than DCR. We also assume that the spin system is undergoing non-Markovian dephasing during the field sensing period. In this case, using the recent results of [13], we estimate that for our FQ and NV implementation a magnetic field B can be measured with sensitivity δB √ T = 1.4 pT/ √ Hz where T is the total sensing time (see Appendix B for details). This is a factor of ∼ 2.7 times improvement over the best sensitivity that can be achieved using a separable state of the spins. Such a magnetic field sensor would also have a very good spatial resolution ∼ 2 µm, as determined by the size of the diamond sample. For the MR and donor spin implementation, we estimate δB √ T = 10 fT/ √ Hz, a factor of ∼ 4.1 improvement over the best sensitivity that can be achieved using a separable state of the donor spins (see Appendix B for details), and a spatial resolution ∼ 1 mm. There is better sensitivity for the MR implementation than for the FQ implementation since it employs a higher number of spins, and because these spins are donor electrons in silicon, which have much longer coherence times than NV centres [16]. However, the spatial resolution is worse since the sample coupled to the MR is larger than the sample coupled to the FQ (see Fig. 2). By using a flux qubit ancillary system with N = 500 donor spins in silicon we could combine the best features of both implementations, giving a sensitivity δB √ T = 75 fT/ √ Hz and a spatial resolution ∼ 2 µm. In Fig. 9 we show how such a magnetometer would compare with the reported sensitivies of some existing state-of-the-art magnetometers. Relative to these existing devices, our proposed magnetometer would occupy an unexplored region of high sensitivity and high spatial resolution. The sensitivity, and the advantage of using a squeezed state instead of a separable state, can be improved by generating a squeezed state with a larger number of spins. To do this, the challenge is to couple a larger number of spins to the ancillary system, while maintaining long coherence times and strong, relatively homogeneous coupling. We note that this should be possible as experimental techniques advance. For example, if a superconducting circuit can be arranged in a Helmholz coil configuration, that is, with two superconducting solenoids separated by the solenoid radius, then the coupling to a spin ensemble placed between the solenoids will be stronger, and more homogeneous over a larger spatial region than for the setups in Fig. 2. This would enable the generation of much larger spin squeezed states.
To derive the effective master Eq. 12 we follow the procedure developed by Reiter and Sorensen [36] for adiabatic elimination in an open quantum system. The starting point is the master Eq. 6,ρ = − i [Ĥ, ρ] + γD[Â](ρ), where the Lindblad operatorÂ = d−2 n=0 √ n + 1 |n n + 1| represents dissipation of energy into the ground state of the ancillary system. The projector onto the ancillary system ground state is denoted P g = |0 0|, whileP e = I −P g = d−1 n=1 |n n| projects onto the excited subspace. This allows the Hamiltonian to be written asĤ = (P g +P e )Ĥ(P g +P e ) =V g +V e +V + +V − , wherê V g =P gĤPg = ( Ω n η · Ĵ +Ĥ IB ) ⊗P g , V e =P eĤPe = ( Ω n η · Ĵ +Ĥ IB ) ⊗P e + ∆ d−1 n=1 n |n n| , describe the dynamics within the ground and excited subspaces respectively, whilê give the dynamics that connects the two subspaces.
The Reiter-Sorensen procedure gives a prescription for the derivation of an effective master equation under the assumptions that the dynamics due toV g andV ± are much slower than either the dissipative dynamics or the dynamics due toV e . In our model this requirement is satisfied if Γ λ N , Γ Ω and Γ δω, where Γ = ∆ 2 + γ 2 /4. According to the Reiter-Sorensen procedure, the resulting effective master equation isρ = − i [V eff , ρ] + D[L eff ](ρ) with effective Hamiltonian and Lindblad operatorŝ We refer the reader to Ref. [36] for a derivation of the effective operators Eqs. A1, A2 and A3. (The essence of the approximation is that the dynamics due toV g and V ± are perturbatively small compared to the dynamics due to the (non-Hermitian) HamiltonianV NH , allowing adiabatic elimination of the excited subspace.) Since we have already assumed that Γ max{Ω, δω}, we can approximatê Substituting Eq. A5 into Eqs. A1 and A2 gives: We may ignore the projectorP g in Eqs. A6 and A7 since the ancillary system remains in its ground state