Phonon-induced spin-spin interactions in diamond nanostructures:application to spin squeezing

# Phonon-induced spin-spin interactions in diamond nanostructures: application to spin squeezing

S. D. Bennett    N. Y. Yao    J. Otterbach    P. Zoller    P. Rabl    M. D. Lukin Physics Department, Harvard University, Cambridge, Massachusetts 02138, USA Institute for Quantum Optics and Quantum Information, Austrian Academy of Sciences, 6020 Innsbruck, Austria Institute for Theoretical Physics, University of Innsbruck, 6020 Innsbruck, Austria Institute of Atomic and Subatomic Physics, TU Wien, Stadionallee 2, 1020 Wien, Austria
July 17, 2019July 17, 2019
July 17, 2019July 17, 2019
###### Abstract

We propose and analyze a novel mechanism for long-range spin-spin interactions in diamond nanostructures. The interactions between electronic spins, associated with nitrogen-vacancy centers in diamond, are mediated by their coupling via strain to the vibrational mode of a diamond mechanical nanoresonator. This coupling results in phonon-mediated effective spin-spin interactions that can be used to generate squeezed states of a spin ensemble. We show that spin dephasing and relaxation can be largely suppressed, allowing for substantial spin squeezing under realistic experimental conditions. Our approach has implications for spin-ensemble magnetometry, as well as phonon-mediated quantum information processing with spin qubits.

###### pacs:
07.10.Cm, 71.55.-i, 42.50.Dv

Electronic spins associated with nitrogen-vacancy (NV) centers in diamond exhibit long coherence times and optical addressability, motivating extensive research on NV-based quantum information and sensing applications. Recent experiments have demonstrated coupling of NV electronic spins to nuclear spins Jelezko2004 (); Childress2006 (), entanglement with photons Togan2010 (), as well as single spin Maze2008 (); Balasubramanian2008 () and ensemble Acosta2009 (); Pham2011 () magnetometry. An outstanding challenge is the realization of controlled interactions between several NV centers, required for quantum gates or to generate entangled spin states for quantum-enhanced sensing. One approach toward this goal is to couple NV centers to a resonant optical Englund2010 (); Faraon2011 () or mechanical Rabl2010 (); Arcizet2011 (); Kolkowitz2012 () mode; this is particularly appealing in light of rapid progress in the fabrication of diamond nanostructures with improved optical and mechanical properties Zalalutdinov2011 (); Burek2012 (); Hausmann2012 (); Ovartchaiyapong2012 (); Tao2012 ().

In this Letter, we describe a new approach for effective spin-spin interactions between NV centers based on strain-induced coupling to a vibrational mode of a diamond resonator. We consider an ensemble of NV centers embedded in a single crystal diamond nanobeam, as depicted in Fig. 1a. When the beam flexes, it strains the diamond lattice which in turn couples directly to the spin triplet states in the NV electronic ground state Maze2011 (); Doherty2012 (). For a thin beam of length m, this strain-induced spin-phonon coupling can allow for coherent effective spin-spin interactions mediated by virtual phonons. Based on these effective interactions, we explore the possibility to generate spin squeezing of an NV ensemble embedded in the nanobeam. We account for spin dephasing and mechanical dissipation, and describe how spin echo techniques and mechanical driving can be used to suppress the dominant decoherence processes while preserving the coherent spin-spin interactions. Using these techniques we find that significant spin squeezing can be achieved with realistic experimental parameters. Our results have implications for NV ensemble magnetometry, and provide a new route toward controlled long-range spin-spin interactions.

Model.—The electronic ground state of the negatively charged NV center is a spin triplet with spin states labeled by as shown in Fig. 1b. In the presence of external electric and magnetic fields and , the Hamiltonian for a single NV is Doherty2012 ()

 HNV=(D0+d∥Ez)S2z+μBgs→S⋅→B−d⊥[Ex(SxSy+SySx)+Ey(S2x−S2y)],Ê (1)

where GHz is the zero field splitting, , is the Bohr magneton, and () is the ground state electric dipole moment in the direction parallel (perpendicular) to the NV axis Vanoort1990 (); Dolde2011 ().

Motion of the diamond nanoresonator changes the local strain at the position of the NV center, which results in an effective, strain-induced electric field Doherty2012 (). We are interested in the near-resonant coupling of a single resonant mode of the nanobeam to the transition of the NV, with Zeeman splitting , as shown in Fig. 1b,c. The perpendicular component of strain mixes the states. For small beam displacements, the strain is linear in its position and we write , where is the destruction operator of the resonant mechanical mode of frequency , and is the perpendicular strain resulting from the zero point motion of the beam. We note that the parallel component of strain shifts both states relative to Acosta2010 (); however, with near-resonant coupling and preparation in the subspace, the state remains unpopulated and parallel strain plays no role in what follows. Within this two-level subspace, the interaction of each NV is , where is the Pauli operator of the th NV center and is the single phonon coupling strength. For many NV centers we introduce collective spin operators, and , which satisfy the usual angular momentum commutation relations. The total system Hamiltonian can then be written as

 H=ωma†a+ΔBJz+g(a†J−+aJ+), (2)

which describes a Tavis-Cummings type interaction between an ensemble of spins and a single mechanical mode. In Eq. (2) we have assumed uniform coupling of each spin to the mechanical mode for simplicity. In general the coupling may be nonuniform and we discuss this further below.

To estimate the coupling strength , we calculate the strain for a given mechanical mode and use the experimentally obtained stress coupling of Hz Pa in the NV ground state Togan2011 (); si (). We take a doubly clamped diamond beam (see Fig. 1a) with dimensions such that Euler-Bernoulli thin beam elasticity theory is valid LL (). For NV centers located near the surface of the beam we obtain si ()

 g2π≈180(ℏL3w√ρE)1/2GHz, (3)

where is the mass density and is the Young’s modulus of diamond. For a beam of dimensions m we obtain a vibrational frequency GHz and coupling kHz. While this is smaller than the strain coupling MHz expected for electronic excited states of defect centers Habraken2012 (); Soykal2011 () or quantum dots Wilson-Rae2004 (), we benefit from the much longer spin coherence time in the ground state. An important figure of merit is the single spin cooperativity , where is the mechanical damping rate and is the equilibrium phonon occupation number at temperature . Assuming , ms and K, we obtain a single spin cooperativity of . This can be further increased by reducing the dimensions of the nanobeam and operating at lower temperatures.

Spin squeezing.—In the dispersive regime, , virtual excitations of the mechanical mode result in effective interactions between the otherwise decoupled spins. In this limit, can be approximately diagonalized by the transformation with . To order this yields an effective Hamiltonian,

 (4)

where is the phonon-mediated spin-spin coupling strength. Rewriting , and provided the total angular momentum is conserved, we obtain a term corresponding to the one-axis twisting Hamiltonian Kitagawa1993 ().

To generate a spin squeezed state, we initialize the ensemble in a coherent spin state (CSS) along the axis of the collective Bloch sphere. The CSS satisfies and has equal transverse variances, . This can be achieved using optical pumping and global rotations of the spins with microwave fields Taylor2008 (). The squeezing term describes a precession of the collective spin about the axis at a rate proportional to , resulting in a shearing of the uncertainty distribution and a reduced spin variance in one direction as shown in Fig. 1d. This is quantified by the squeezing parameter Wineland1992 (); Ma2011 (),

 ξ2=2J⟨ΔJ2min⟩⟨Jx⟩2, (5)

where is the minimum spin uncertainty with and . The preparation of a spin squeezed state, characterized by , has direct implications for NV ensemble magnetometry applications, since it would enable magnetic field sensing with a precision below the projection noise limit Wineland1992 ().

We now consider spin squeezing in the presence of realistic decoherence. In addition to the coherent dynamics described by , we account for mechanical dissipation and spin dephasing using a master equation si ()

 ˙ρ= −i[−λ2J2z+(ΔB+λa†a)Jz,ρ]+12T2∑iD[σiz]ρ +Γγ(¯nth+1)D[J−]+Γγ¯nthD[J+], (6)

where and the single spin dephasing is assumed to be Markovian for simplicity (see below). Note that we absorbed a shift of into , and ignored single spin relaxation as can be several minutes at low temperatures Jarmola2012 (). The second line describes collective spin relaxation induced by mechanical dissipation, with . Finally, the phonon number shifts the spin frequency, acting as an effective fluctuating magnetic field which leads to additional dephasing.

Let us for the moment ignore fluctuations of the phonon number ; we address these in detail below. Starting from the CSS , we plot the squeezing parameter in Fig. 2a for an ensemble of spins and several values of , in the presence of dephasing and collective relaxation . Here we calculated by solving Eq. (Phonon-induced spin-spin interactions in diamond nanostructures: application to spin squeezing) using an approximate numerical approach treating and separately, and verified that the approximation agrees with exact results for small si (). To estimate the minimum squeezing, we linearize the equations of motion for the averages and variances of the collective spin operators (see dashed lines in Fig. 2a). From these linearized equations, in the limits of interest, , and to leading order in both sources of decoherence, we obtain approximately

 ξ2≃4Γγ¯nthJλ2t+tT2. (7)

Optimizing and the detuning , we obtain the optimal squeezing parameter,

 ξ2opt≃2√Jη, (8)

at time , similar to results for atomic systems SchleierSmith2010 (); Leroux2010 (); Leroux2012 (). Note that for non-Markovian dephasing, the scaling is even more favorable Marcos2010 (). In Fig. 2b we plot the scaling of the squeezing parameter with for small but finite decoherence, and find agreement with Eq. (8). For comparison we also plot the unitary result in the absence of decoherence, scaling as and limited by the Bloch sphere curvature Kitagawa1993 ().

Phonon number fluctuations.—In Eq. (4) we see that the phonon number couples to , leading to additional dephasing due to thermal number fluctuations. On the other hand, this same coupling can also lead to additional spin squeezing from cavity feedback, by driving the mechanical mode SchleierSmith2010 (); Leroux2010 (); Leroux2012 (). In the following, we consider a twofold approach to mitigate thermal spin dephasing while preserving the optimal squeezing. First, we apply a sequence of global spin echo control pulses to suppress dephasing from low-frequency thermal fluctuations. This also extends the effective coherence time of single NV spins Taylor2008 (). Second, we consider driving the mechanical mode, and identify conditions when this results in a net improvement of the squeezing.

To simultaneously account for thermal dephasing, driven feedback squeezing, and spin control pulse sequences, we write the interaction term in Eq. (4) in the so-called “toggling frame” Haeberlen1968 (),

 Hint(t)=λJzf(t)δn(t). (9)

The function periodically inverts the sign of the interaction as shown in the inset of Fig. 3a, describing the inversion of the collective spin with each pulse of the spin echo sequence. Phonon number fluctuations are described by , where is the mean phonon number and we have omitted an average frequency shift proportional to in Eq. (9). The number fluctuation spectrum is plotted in Fig. 3a for a driven oscillator coupled to a thermal bath si ().

We calculate the required spin moments within the Gaussian approximation for phonon number fluctuations, and obtain si ()

 ⟨J+(t)⟩=e−χ⟨e−iμ(Jz−1/2)J+(0)⟩, (10)

and similar results for and . In Eq. (10) the dephasing parameter and effective squeezing via are given by

 χ =λ2∫dω2πF(ωτ)ω2¯Sn(ω), (11) μ =λ2∫dω2πK(ωτ)ω2An(ω), (12)

where and . The filter function describes the effect of the spin echo pulse sequence with time between pulses Martinis2003 (); Uhrig2007 (); Cywinski2008 (). The function plays the analogous role for the effective squeezing described by , and is related to by a Kramers-Kronig relation si (). We plot and for a sequence of pulses in Fig. 3a.

Discussion.—We now consider the impact of thermal fluctuations on the achievable squeezing. The noise spectrum is symmetric around . Without spin echo control pulses, this low frequency noise results in nonexponential decay of the spin coherence, (with ), familiar from single qubit decoherence Taylor2008 (); deSousa2009 (). The inhomogeneous thermal dephasing time is , severely limiting the possibility of spin squeezing. In particular, at time we find that squeezing is prohibited when si (). However, one can overcome this low frequency thermal noise using spin echo. By applying a sequence of equally spaced global -pulses to the spins during precession of total time , we obtain , suggesting that thermal dephasing can be made negligible relative to both and . For a sufficiently large number of pulses, , we recover the optimal squeezing in Eqs. (7) and (8).

Adding a mechanical drive can further enhance squeezing via feedback; however, it also increases phonon number fluctuations, contributing to additional dephasing. We consider a detuned external drive of frequency , leading to two additional peaks in at , as shown in Fig. 3a. The area under the left [right] peak scales as [], where is the mean phonon number due to the drive at zero temperature. The symmetric and antisymmetric parts of this noise contribute to dephasing and squeezing as described by Eqs. (11) and (12). Choosing the interval between pulses, we obtain additional dephasing and effective squeezing with . In the limit , the effects of the drive dominate over and and we recover the ideal scaling given in Eq. (8), even with a small number of echo pulses. This is shown in Fig. 3b,c where we see that the optimal squeezing improves with increasing for a fixed number of pulses .

Finally, we discuss our assumption of uniform coupling strength in Eq. (2). This is an important practical issue, as we expect the coupling to individual spins to be inhomogeneous in experiment due to the spatial variation of strain in the beam. Nonetheless, even with nonuniform coupling, we still obtain squeezing of a collective spin with a reduced effective total spin , provided . First, we note that inhomogeneous magnetic fields resulting in nonuniform detuning are compensated by spin echo. Second, for a distribution of coupling strengths , the effective length of the collective spin is for the direct squeezing term, and for feedback squeezing with a mechanical drive. Similar conclusions were reached in atomic and nuclear systems SchleierSmith2010 (); Leroux2010 (); Leroux2012 (); Rudner2011 (). In the case of direct squeezing, it is important that the sign of the ’s is the same to avoid cancellation; this is automatically achieved by using NV centers implanted on the top of the beam. For beam dimensions m analyzed above, we estimate that NV centers can be embedded without being perturbed by direct magnetic dipole-dipole interactions. A reduction of the effective spin length by factor still leaves , sufficient to observe spin squeezing.

Conclusions.—We have shown that direct spin-phonon coupling in diamond can be used to prepare spin squeezed states of an NV ensemble embedded in a nanoresonator, even in the presence of dephasing and mechanical dissipation. With further reductions in temperature, beam dimensions, and spin decoherence rates, the regime of large single spin cooperativity could be achieved. This would allow for coherent phonon-mediated interactions and quantum gates between two spins embedded in the same resonator via , and coupling over larger distances by phononic channels Habraken2012 ().

Acknowledgments.—The authors gratefully acknowledge discussions with Shimon Kolkowitz and Quirin Unterreithmeier. This work was supported by NSF, CUA, DARPA, NSERC, HQOC, DOE, the Packard Foundation, the EU project AQUTE and the Austrian Science Fund (FWF) through SFB FOQUS and the START grant Y 591-N16.

## References

• (1) F. Jelezko et al., Phys. Rev. Lett. 93, 130501 (2004).
• (2) L. Childress et al., Science 314, 281 (2006).
• (3) E. Togan et al., Nature 466, 730 (2010).
• (4) J. R. Maze et al., Nature 455, 644 (2008).
• (5) G. Balasubramanian et al., Nature 455, 648 (2008).
• (6) V. M. Acosta et al., Phys. Rev. B 80, 115202 (2009).
• (7) L. M. Pham et al., New J. Phys. 13, 045021 (2011).
• (8) D. Englund et al., Nano Lett. 10, 3922 (2010).
• (9) A. Faraon et al., Nature Photon. 5, 301 (2011).
• (10) P. Rabl et al., Nature Phys. 6, 602 (2010).
• (11) O. Arcizet et al., Nature Phys. 7, 879 (2011).
• (12) S. Kolkowitz et al., Science 335, 1603 (2012).
• (13) B. J. M. Hausmann et al., Phys. Status Solidi A 209, 1619 (2012).
• (14) P. Ovartchaiyapong et al., Appl. Phys. Lett. 101, 163505 (2012).
• (15) M. K. Zalalutdinov et al., Nano Lett. 11, 4304 (2011).
• (16) M. J. Burek et al., Nano Lett. 12, 6084 (2012).
• (17) Y. Tao et al., arXiv:1212.1347.
• (18) J. R. Maze et al., New J. Phys. 13, 025025 (2011).
• (19) M. Doherty et al., Phys. Rev. B 85, 205203 (2012).
• (20) E. Vanoort and M. Glasbeek, Chem. Phys. Lett. 168, 529 (1990).
• (21) F. Dolde et al., Nature Phys. 7, 459 (2011).
• (22) V. M. Acosta et al., Phys. Rev. Lett. 104 (2010).
• (23) See supplementary information.
• (24) E. Togan et al., Nature 478, 497 (2011).
• (25) L. D. Landau and E. M. Lifshitz, Theory of Elasticity (Butterworth-Heinemann, Oxford, 1986).
• (26) S. J. M. Habraken et al., New J. Phys. 14, 115004 (2012).
• (27) O. O. Soykal, R. Ruskov, and C. Tahan, Phys. Rev. Lett. 107 (2011).
• (28) I. Wilson-Rae, P. Zoller, and A. Imamoglu, Phys. Rev. Lett. 92, 075507 (2004).
• (29) M. Kitagawa and M. Ueda, Phys. Rev. A 47, 5138 (1993).
• (30) J. M. Taylor et al., Nature Phys. 4, 810 (2008).
• (31) D. J. Wineland et al., Phys. Rev. A 46, R6797 (1992).
• (32) J. Ma et al., Phys. Rep. 509, 89 (2011).
• (33) A. Jarmola et al., Phys. Rev. Lett. 108, 197601 (2012).
• (34) M. H. Schleier-Smith, I. D. Leroux, and V. Vuletic, Phys. Rev. A 81, 021804(R) (2010).
• (35) I. D. Leroux, M. H. Schleier-Smith, and V. Vuletic, Phys. Rev. Lett., 104, 073602 (2010).
• (36) I. D. Leroux et al., Phys. Rev. A 85, 013803 (2012).
• (37) D. Marcos et al., Phys. Rev. Lett. 105, 210501 (2010).
• (38) U. Haeberlen and J. Waugh, Phys. Rev. 175, 453 (1968).
• (39) L. Cywiński et al., Phys. Rev. B 77, 174509 (2008).
• (40) J. M. Martinis et al., Phys. Rev. B 67, 094510 (2003).
• (41) G. S. Uhrig, Phys. Rev. Lett. 98, 100504 (2007).
• (42) R. de Sousa, Top. Appl. Phys. 115, 183 (2009).
• (43) M. Rudner et al., Phys. Rev. Lett. 107, 206806 (2011).

## Appendix A Supplemental information

### a.1 Coupling strength

We assume that the NV axis is aligned with both the magnetic field and with the direction of beam deflection, so that the longitudinal strain due to deflection is entirely perpendicular to the NV axis. From experiment Togan2010supp (), the splitting of the states with stress is Hz/Pa. We convert this into the deformation potential coupling frequency, GHz, using the Young’s modulus of diamond, GPa. Next we calculate the strain at the NV center using elasticity theory. The equation for the bending mode of a thin beam is

 ρA∂2tϕ(z,t)=−EI∂4zϕ(z,t) (13)

where is the transverse displacement in the direction and is along the beam. Here is the mass density, is the cross sectional area, and is the moment of inertia. The solutions are of the form where

 u(z)=1√N[coskz−coshkz−(coskL−coshkL)(sinkL−sinhkL)(sinkz−sinhkz)], (14)

which satisfies the boundary conditions for a doubly clamped beam. The allowed wavenumbers are given by the solutions of , and the corresponding eigenfrequencies are

 ωn=k2n√EIρA. (15)

We normalize the modefunction of the fundamental mode by setting the free energy stored in the beam to the zero point energy,

 W=12EI∫L0dz(∂2u0∂z2)2=ℏω02. (16)

Integrating by parts we obtain the normalization condition,

 ∫L0dzu20=ℏρAω0. (17)

If the NV center lies at the midpoint along the beam, , and at a distance from the neutral axis of the beam, the strain due to the zero point motion of the fundamental mode is

 ϵ0=−r0∂2zu0(L/2)∼r0√ℏρAω027L5/2∼52r0t√ℏL3w√Eρ, (18)

where we Eq. (15) and for the fundamental mode. The coupling strength is given by the deformation potential and the strain due to zero point motion,

 g2π=Ξ2πℏϵ0∼180 GHz⋅2r0t√ℏL3w√Eρ. (19)

For an NV near the surface of the beam, and we obtain Eq. (4) of the main text.

### a.2 Effective squeezing Hamiltonian from spin-phonon coupling

In this section we provide additional details on deriving in Eq. (4) from the original in Eq. (1). Assuming that the magnetic field is aligned along the NV axis, , and defining and , we can rewrite as ()

 HNV=(D0+d∥Ez)S2z+gsμBBSz−d⊥2(E+S2++E−S2−), (20)

where is the effective electric field due to strain. We quantize the perpendicular strain field, and , where is the strain due to the zero point motion of the resonant mode. Next we focus on the two-level subspace only, and assume that transitions to state are not allowed due to the large zero field splitting . For the th spin we write Pauli operators and , and within this two-level subspace the interaction for a single NV is

 Hi=ΔB2σzi+g(σ+ia+a†σ−i)+ωma†a, (21)

where , is the energy between , and we included the mechanical oscillator of frequency . Summing Eq. (21) for many NVs coupled to the same mode with uniform coupling strength we obtain

 H=ΔBJz+g(a†J−+aJ+)+ωma†a, (22)

which is Eq. (2) of the main text. To obtain Eq. (4), we first rewrite in the rotating frame at the mechanical frequency ,

 H=ΔJz+g(a†J−+aJ+), (23)

where . Next we apply the transformation , with , and to order we obtain

 Heff≃ΔJz+g2Δ(J+J−+2a†aJz). (24)

Transforming back to the nonrotating frame yields Eq. (4) of the main text.

### a.3 Individual spin dephasing and phonon-induced relaxation

#### a.3.1 Individual spin dephasing from intrinsic T2

Each NV spin experiences intrinsic decoherence in the absence of the mechanical mode. Individual relaxation () processes are due to lattice phonons; at low temperature can be s and we ignore it Jarmola2012supp (). However, we include intrinsic single spin dephasing, which arises from magnetic noise of C nuclear spins in the diamond lattice. In practice, single spin dephasing may be nonexponential Taylor2008supp (); deSousa2009supp (), but for simplicity we approximate the effect of single spin dephasing by an effective Markovian master equation with dephasing rate ,

 ˙ρ=12T2∑i[σziρσzi−ρ]. (25)

#### a.3.2 Collective phonon-induced spin relaxation

The transformation used to obtain the effective squeezing Hamiltonian also introduces a relaxation channel for the collective spin by admixing phonon and spin degrees of freedom. Mechanical dissipation is described by the master equation for the system density matrix ,

 (26)

Transforming and using the transformation in the main text, we obtain effective spin relaxation terms in the master equation,

 ˙ρ=Γγ(¯nth+1)[J−ρJ+−12(J+J−ρ+ρJ+J−)]+Γγ¯nth[J+ρJ−−12(J−J+ρ+ρJ−J+)] (27)

where . From Eqs. (25) and (27) we calculate the equations for spin averages and variances accounting for individual dephasing and collective relaxation using .

#### a.3.3 Squeezing estimate from linearized equations for spin averages and variances

Here we sketch the derivation of the estimated optimal squeezing given in Eq. (8) in the main text. In order to treat the squeezing Hamiltonian, collective relaxation and spin dephasing on equal footing, we linearize the equations for the spin averages and variances. This corresponds to expanding in the small error from decoherence at short times and ignoring the curvature of the Bloch sphere for sufficiently short times, when the spin uncertainty distribution remains on a locally flat region of the Bloch sphere. To linearize the equations we assume that all (connected) correlations of order higher than two vanish. The linearized equations are valid for short times, so we also make use of the initial conditions in the spin coherent state at , which are and . Here we define the covariance operator , while in the main text we refer only to its average, . Within these approximations, and using Eqs. (25) and(27), the linearized equations for the spin averages required to calculate the squeezing parameter are

 ∂t⟨Jx⟩ =−Γ2⟨Jx⟩ (28) ∂t⟨Jy⟩ =λJ⟨Jz⟩−Γ2⟨Jy⟩−Γγ(¯nth+12)⟨Jy⟩+Γγ⟨Cyz⟩ (29) ∂t⟨Jz⟩ (30) ∂t⟨J2y⟩ (31) ∂t⟨J2z⟩ (32) ∂t⟨Cyz⟩ (33)

This linear set of equations can be directly solved. The full analytic solutions are lengthy so we simply plot the numerical solution for the squeezing parameter (blue dash-dotted in Fig. 4), which agrees with exact numerics at short times.

We use these linearized equations to estimate the scaling of the optimal squeezing parameter (see Eqs. (7) and (8) in the main text). First, we solve Eqs. (28-33) to second order in . Second, we calculate (see Eq. (5) in the main text) from the resulting spin averages, and simplify the result in the limit of interest, , and assuming as required for significant squeezing. Third, we assume that all sources of decoherence are small, and expand in the errors and . Within these approximations we obtain

 ξ2≃1+4JΓγ¯ntht(Jλt)2+(5Γγ¯nth+T−12)t, (34)

where the first term is the result from linearized unitary squeezing, and the remaining terms are the lowest order corrections in both sources of decoherence. We further approximate , valid for sufficiently large detuning , and , valid self-consistently at the optimal squeezing time and in the relevant limit . Within these approximations we obtain Eq. (7) in the main text. Finally, we optimize with respect to , obtaining Eq. (8) and given in the main text.

#### a.3.4 Combining individual dephasing and collective relaxation: numerics

As discussed in the main text, in the absence of a mechanical drive we can neglect phonon number fluctuations for a sufficiently large number of pulses . In this case the remaining sources of decoherence are intrinsic single spin dephasing and collective relaxation induced by mechanical dissipation. These sources of decoherence are simple to treat separately but difficult to treat simultaneously for a large number of spins. To calculate the solid blue curves in Fig. 2 of the main text, we treat the combination of both sources of decoherence as approximately independent, valid provided both are small enough to still allow spin squeezing. To calculate the spin averages needed for the squeezing parameter, we first account for collective relaxation using the Dicke state basis in which total is conserved. We then account for individual dephasing by multiplying the resulting averages by dephasing factors such as where is the result of the Dicke state calculation. Each step would be numerically exact in the absence of the other source of decoherence; thus we expect that this procedure provides a good approximation if all errors are small. To verify the accuracy of the approach, we compare the result with exact numerics calculated by numerically integrating the full master equation for small in Fig. 4.

### a.4 Phonon number fluctutions

In this section we consider fluctuations of the phonon number, . We start by rewriting the effective Hamiltonian [see Eq. (4) in main text] for the collective spin coupled to a driven oscillator in the frame rotating at the mechanical drive frequency,

 Heff=ΔJz+λJza†a+λ2J2z−δa†a+Ω(a+a†) (35)

where is the detuning of the magnetic transition frequency from the drive, and is the drive detuning from the mechanical frequency. The amplitude of the drive is and we have made the rotating wave approximation. Our aim is to find the effect of the oscillator on the spin to second order in (within the Gaussian approximation). For this we require the number fluctuation spectrum of a damped, driven, thermal oscillator in the absence of coupling to the spin.

#### a.4.1 Number fluctuations of a driven thermal mode

To calculate the effective dephasing from number fluctuations, we first need the power spectral density of phonon number fluctuations,

 Sn(ω)=∫dteiωt⟨δn(t)δn(0)⟩, (36)

where , and the average is taken with respect to the oscillator in thermal equilibrium with its environment. In the absence of coupling, , the Langevin equation for the driven thermal mode in the frame of the classical drive frequency and within the rotating wave approximation is

 ˙a(t)=(iδ−γ2)a(t)+Ω+√γξ(t) (37)

The solution is , where is the coherent amplitude due to the drive, and

 d(t)=√γ∫t−∞dt′eiδ−γ/2)(t−t′)ξ(t′) (38)

describes thermal and quantum fluctuations. The mean phonon number is the sum of driven and thermal parts, , where

 ¯ndr=∣∣α2∣∣=Ω2δ2+γ2/4 (39)

and the thermal occupation is . Using Eq. (38) we find the two-time correlations,

 ⟨d†(t)d(0)⟩ =¯nthe−iδte−γ|t|/2, (40) ⟨d(t)d†(0)⟩ =(¯nth+1)eiδte−γ|t|/2, (41)

and from these we can calculate the full spectrum of driven thermal number fluctuations. The correlation using Wick’s theorem and is

 ⟨⟨δn(t)δn(0)⟩⟩= ⟨n(t)n(0)⟩−¯n2 = (42) +⟨d†(t)d(0)⟩⟨d(t)d†(0)⟩. (43)

Using Eq. (38) and taking the Fourier transform, we find that the number fluctuation spectrum for a driven thermal oscillator is given by

 Sn(ω)=γ¯ndr[¯nth(ω−δ)2+γ2/4+¯nth+1(ω+δ)2+γ2/4]+2γ¯nth(¯nth+1)ω2+γ2. (44)

#### a.4.2 Effects of number fluctuations on the spin in the Gaussian approximation

From the spectrum of number fluctuations we can calculate the effect of number fluctutions on the spin dephasing and squeezing. We write the full Hamiltonian as , where describes the driven damped oscillator, and

 H0=ΔJz+λ2J2z (45)

describes the spin including the constant effective squeezing term. The coupling in the interaction picture and in the toggling frame is

 Hint(t)=λf(t)δn(t)Jz, (46)

where is time-independent as it commutes with the full Hamiltonian. We have included the function to describe spin echo, which effectively inverts the sign of the interaction with each pulse.

The equation of motion for the operator in the interaction picture is

 ˙J+(t)=iλf(t)[Jzδn(t),J+(t)]. (47)

We integrate this formally, insert the solution, and take the average with respect to the oscillator to get

 ˙J+(t)=−λ2f(t)∫t0f(t′)⟨[Jzδn(t)[Jzδn(t′),J+(t′)]]⟩osc (48)

where denotes averaging over the oscillator degrees of freedom. Note that we neglected additional noise terms; these play no role as we will only be interested in taking the average at the end. Next, we neglect the time dependence of under the integral, as it is higher order in , . Expanding the commutators we obtain

 ˙J+ =−λ2f(t)∫t0f(t′)∫dω2πe−iω(t−t′)[JzJ+Sn(ω)−J+JzSn(−ω)], (49)

using Eq. (36). Defining the symmetric and antisymmetric parts of the number fluctuation spectrum,

 ¯Sn(ω)=12[Sn(ω)+Sn(−ω)],An(ω)=Sn(ω)−Sn(−ω), (50)

we can rewrite Eq. (49) as

 ˙J+=−λ2f(t)∫t0f(t′)∫dω2πe−iω(t−t′)×[¯Sn(ω)+An(ω)(Jz−12)]J+. (51)

Solving Eq. (51) and finally taking the average with respect to spin degrees of freedom, we obtain

 ⟨J+(t)⟩=e−χ⟨eiμ(Jz−12)J+(0)⟩ (52)

where is the average over all degrees of freedom, and

 χ =λ2∫dω2π¯Sn(ω)∫dt1∫dt2e−iω(t1−t2)θ(t1−t2)f(t1)f(t2), (53) μ =iλ2∫dω2πAn(ω)∫dt1∫dt2e−iω(t1−t2)θ(t1−t2)f(t1)f(t2). (54)

Here all integration limits are from to , and the time integration limits are accounted for in and the step function . Similarly, we obtain the other averages needed to calculate the squeezing,

 ⟨J2+(t)⟩ =e−4χ⟨e2iμ(Jz−1)J2+(0)⟩, (55) ⟨J+(t)Jz(t)⟩ =e−χ⟨eiμ(Jz−12)J+(0)Jz(0)⟩. (56)

By comparing with the spin evolution under unitary one-axis twisting, we see that describes spin squeezing with an effective squeezing coefficient . The parameter describes collective dephasing.

To evaluate and for a given pulse sequence, we next define to rewrite the double time integral as

 It= ∫dt1∫dt2e−iω(t1−t2)θ(t1−t2)f(t1)f(t2) (57) =∫dω12π|f(ω1)|2[πδ(ω+ω1)−iω+ω1]. (58)

We define the filter function for pulse sequence with time between pulses,

 F(ωτ)=ω22|f(ω)|2. (59)

The dephasing term is

 χ =2λ2∫dω2π¯Sn(ω)∫dω12πF(ω1τ)ω21[πδ(ω+ω1)−iω+ω1] (60)

Since and are both even in , the imaginary part of the integrand is odd and integrates to zero. As a result is real and we obtain

 χ=λ2∫dω2πF(ωτ)ω2¯Sn(ω). (61)

The coherent term is

 μ =iλ2∫dω2πAn(ω)∫dω12πF(ω1τ)ω21[πδ(ω+ω1)−iω+ω1]. (62)

Since is odd, in this case the imaginary part (involving the -function) is zero. The real part is

 μ=λ2∫dω2πK(ωτ)ω2An(ω), (63)

where we defined

 K(ωτ)=2ω2∫dω12πF(ω1τ)ω21(ω+ω1). (64)

and satisfy a Kramers-Kronig relation (with a factor of from the definitions).

#### a.4.3 Dephasing from purely thermal oscillator

From Eq. (44), the number fluctuation spectrum of a purely thermal oscillator in the frame of the mechanical drive is

 Sthn(ω)=2γ¯nth(¯nth+1)ω2+γ2. (65)

The spectrum is symmetric in frequency and thus . We obtain the dephasing from thermal flucuations by inserting Eq. (65) in Eq. (