Scale-invariant spin dynamics and the quantum limits of field sensing

The quantum limits of measurement is a rich topic of both fundamental and practical interest. The theory of these limits informs many other topics, including the statistics of parameter estimation [1], the geometry of quantum states [2], entanglement in many-body systems [3], quantum information processing [4, 5], and quantum non-locality [6, 7]. Understanding quantum measurement effects has led to improved sensitivity in gravitational wave detectors [8–10] and progress toward similar improvements in measurements of time [11–15], dc magnetic fields [16, 17] and radio-frequency fields [18, 19].

The vast majority of prior work on quantum sensitivity limits concerns the problem of linear interferometric parameter estimation. For example, the standard quantum limit (SQL) ⟨δϕ²⟩ ⩾ 1/N and the Heisenberg limit (HL) ⟨δϕ²⟩ ⩾ 1/N² constrain linear estimation of a phase ϕ given the resource of N non-interacting two-level systems. These dimensionless limits acquire units, e.g. length or time, through implementation-dependent scale factors, e.g. a wavelength or a transition frequency. Because these scale factors, as well as the available N, can vary greatly from one implementation to another, such dimensionless limits do not by themselves provide benchmarks by which to compare different sensor implementations.

Here we study a qualitatively different kind of quantum sensitivity limit, one that contains no implementation-specific scale factors, and no makes reference to available resources, only to the quantity to be measured and to the method of measurement. To see what form such a limit could take, consider a sensor that measures the field B in a volume V over an observation time T, and gives a reading B_obs = B_true + δB, where B_true is the true value of the field and δB is the measurement error a zero-mean random variable if the sensor is properly calibrated. The mean apparent magnetostatic energy in the sensor volume is ${E}_{\mathrm{o}\mathrm{b}\mathrm{s}}=\langle {B}_{\mathrm{o}\mathrm{b}\mathrm{s}}^{2}\rangle V/\left(2{\mu }_{0}\right)={B}_{\mathrm{t}\mathrm{r}\mathrm{u}\mathrm{e}}^{2}V/\left(2{\mu }_{0}\right)+\langle \delta {B}^{2}\rangle V/\left(2{\mu }_{0}\right)$, where the second term expresses the sensor's so-called 'energy resolution,' which more properly can be identified as the bias of the apparent energy.

Allowing for averaging of independent measurements of duration T, ⟨δB²⟩T is a figure of merit (smaller is better) that gauges the remaining error after unit total acquisition time. Allowing also for averaging over spatial region, ⟨δB²⟩VT is the relevant figure of merit. It is traditional to express this via the 'energy resolution per bandwidth'

$$\begin{equation}{E}_{\mathrm{R}}\equiv \frac{\langle \delta {B}^{2}\rangle VT}{2{\mu }_{0}}=\frac{{S}_{B}\left(0\right)V}{2{\mu }_{0}}\end{equation} \tag{ 1 }$$

where S_x(0) is the low-frequency limit of the power spectral density of variable x. We seek a limit of the form ${E}_{\mathrm{R}}{\geqslant}\mathcal{S}$, where $\mathcal{S}$ is a constant with units of action. As a purely empirical observation, several technologies to sense low-frequency magnetic fields approach E_R = ℏ, while to date none surpasses it [20]. Such a limit makes no reference to available resources, and the scale factors are the fundamental constants μ₀ and ℏ. This describes a quantum limit resembling in some ways the relation ΔE_xΔB_y ⩾ 2πℏcV^−4/3 derived by W Heisenberg [[21], pp 48–54] where c is the speed of light and ΔE_x and ΔB_y are uncertainties of orthogonal components of the electric and magnetic field, respectively. Whereas Heisenberg's result bounds the joint uncertainty of a pair of observables and allows ΔB = 0 if the electric field E is completely uncertain, our results bound the magnetic field uncertainty while making no reference to the electric field.

We can obtain a first implementation-independent limit from a well-known analysis of dc superconducting quantum interference devices (dc SQUIDs) by Tesche and Clarke (TC) [22]. Considering a lumped-circuit model for dc SQUID magnetometers with resistively-shunted Josephson junctions, TC computed the sensitivity, i.e. power spectral density of the equivalent noise, for an optimized device. At zero temperature, the sensitivity is limited by zero-point current fluctuations in the shunt resistances, to give S_Φ(0)/(2L) ⩾ ℏ where Φ is the flux through the SQUID loop, and L is the loop inductance [23, 24]. With careful construction, small dc SQUID devices have reported S_Φ(0)/(2L) as low as 2ℏ [25–27].

The implementation-dependent factor L can be eliminated by noting that a wire loop has Φ = BA and $L=\sqrt{A}{\mu }_{0}/\alpha$, where A is the loop area and α is a wire-geometry factor of order unity¹. This gives a limit that concerns only field geometry and time: when measuring the field on a patch of area A in a time T with a dc-SQUID, the limiting sensitivity is equivalent to minimum energy per bandwidth of ⟨δB²⟩A^3/2T/(2μ₀) ⩾ αℏ.

In what follows we illustrate and derive an analogous limit for a second important field sensing technology, the spin-precession sensor. As the name suggests, such devices detect magnetic fields by the precession induced in an ensemble of spins. Notable examples include alkali vapors [28, 29], nitrogen-vacancy centers in diamond [30, 31], and spinor Bose–Einstein condensates [32, 33]. One might expect such sensors to be described by the 'quantum metrology' analysis of linear interferometry [5], which gives rise to the SQL and HL given above. If this were the case, there could be no energy resolution limit, because the spin number density ρ = N/V could be taken to infinity, such that 1/N and thus ⟨δB²⟩ approach zero for fixed V, leaving vanishing E_R. At high densities interactions cannot be neglected, however, and the linear interferometry results are not directly applicable [34].

Prior works on quantum sensing with interacting particles have considered scaling [35, 36] and optimization [34, 37] scenarios. Here we show how the mutual magnetic interaction of spins in a sensing ensemble can produce a self-similar spin dynamics resulting in implementation-independent energy resolution. Numerically, we find also $\mathcal{S}=\alpha \hslash$, in agreement with the TC analysis. Prior modeling of high-density NV-center ensembles [38] and high-density alkali vapors [39, 40] have noted invariance of sensitivity with respect to spin density.

We consider a generic prepare-evolve-project sensing protocol, in which an ensemble of N spins {s_i}, each with spin quantum number s, is initialized in a product state $\vert {{\Psi}}_{0}\rangle ={\vert {\psi }_{0}\rangle }^{\otimes N}$, allowed to evolve under a Hamiltonian H containing as parameters the unknown positions of the spins {x_i} and a field component $\mathcal{B}$ to be estimated, e.g. the dc or rf field amplitude. The spins are detected at time T by projection of the total spin S ≡ ∑_is_i onto a direction n.

This model describes an idealized spin-precession sensor, without real-world complications such as the four possible spatial orientations of NV centers in a diamond lattice. We focus on magnetic dipole–dipole interactions among the sensing spins as a source of spin depolarization. In contrast to the much studied effects of a nuclear spin bath [41], dipolar coupling among the sensing spins is intrinsic to the sensor rather than a feature of the environment [42], and will dominate at high densities [43]. We also note that some proposed spin-precession sensors employ other many-body effects including ferromagnetism [44] and quantum degeneracy [33] that are not included here.

Taking V as constant, $T\langle \delta {\mathcal{B}}^{2}\rangle$ determines the energy resolution. For large N, S ⋅ n is nearly Gaussian, and optimal estimation [45] of $\mathcal{B}$ can be understood using propagation of error: if Γ_S is the covariance matrix of S, with elements ⟨S_iS_j + S_jS_i⟩/2 − ⟨S_i⟩⟨S_j⟩, the sensitivity limit is (see appendix A)

$$\begin{equation}T\langle \delta {\mathcal{B}}^{2}\rangle {\geqslant}{\mathrm{min}}_{\mathbf{n},T}\;T\frac{n\cdot {{\Gamma}}_{S}\left(T\right)\cdot n}{{\left\vert {\partial }_{\mathcal{B}}\langle \mathbf{S}\cdot \mathbf{n}\rangle \right\vert }^{2}}.\end{equation} \tag{ 2 }$$

Here the expectation is defined as

$$\begin{equation}\langle \mathcal{A}\rangle =\int \langle {{\Psi}}_{0}\vert {U}^{{\dagger}}\left(t\right)\mathcal{A}U\left(t\right)\vert {{\Psi}}_{0}\rangle {P}_{\rho }\left(\left\{{\mathbf{\text{x}}}_{i}\right\}\right)\;\mathrm{d}\left\{{\mathbf{\text{x}}}_{i}\right\},\end{equation} \tag{ 3 }$$

and includes classical averaging over configurations {x_i} with probability density P_ρ({x}), parametrized by ρ. For example, {x_i} could be distributed as a Poisson point process (PPP), in which an infinitesimal volume dV contains a spin with probability ρ dV. U is the solution to the Schrödinger equation and thus depends on {x_i} and $\mathcal{B}$. Due to coherent signal accumulation at short times and decoherence at long times, the minimum in equation (2) occurs at a finite time T = T_opt. Examples (explained in detail below) are shown in figure 1.

Zoom In Zoom Out Reset image size

Considering a uniform magnetic field B with a constant component B_dc along the z-axis, and including the MDDI interaction (indicated with the subscript dd), the Hamiltonian is [[46], p 103]

$$\begin{equation}H={H}_{L}+{H}_{\mathrm{d}\mathrm{d}}\end{equation} \tag{ 4 }$$

$$\begin{equation}{H}_{L}\equiv -\gamma \hslash \sum _{i}{\mathbf{s}}_{i}\cdot \mathbf{B}\end{equation} \tag{ 5 }$$

$$\begin{align}\hfill {H}_{\mathrm{d}\mathrm{d}}& \equiv \sum _{i\ne j}\frac{{\gamma }^{2}{\hslash }^{2}{\mu }_{0}}{4\pi {r}_{ij}^{3}}\left[{\mathbf{s}}_{i}\cdot {\mathbf{s}}_{j}-3\left({\mathbf{s}}_{i}\cdot {\mathbb{R}}_{ij}\right)\left({\mathbf{s}}_{j}\cdot {\mathbb{R}}_{ij}\right)\right]\hfill \\ \hfill & =\hslash \sum _{i\ne j}\frac{{{\Upsilon}}_{dd}}{{s}^{2}\rho {r}_{ij}^{3}}\left[{\mathbf{s}}_{i}\cdot {\mathbf{s}}_{j}-3\left({\mathbf{s}}_{i}\cdot {\mathbb{R}}_{ij}\right)\left({\mathbf{s}}_{j}\cdot {\mathbb{R}}_{ij}\right)\right],\hfill \end{align} \tag{ 6 }$$

where γ is the gyromagnetic ratio, ${\mathbb{R}}_{ij}={\mathbf{r}}_{ij}/{r}_{ij}$, r_ij ≡ x_i − x_j and r_ij ≡ |r_ij|. ϒ_dd ≡ s²γ²ℏμ₀ρ/(4π) is the strength of the dipole–dipole coupling expressed as an angular frequency.

We now explore the general features of this problem through numerical simulation. Due to the rapid fall-off of ${H}_{\mathrm{d}\mathrm{d}}^{\left(ij\right)}$ with r_ij, and the fact that ${H}_{\mathrm{d}\mathrm{d}}^{\left(ij\right)}$ vanishes when averaged over a sphere of constant r_ij, the dynamics of any given spin s_i will be determined mostly by its nearest neighbors and by B. This motivates the following approximation: we consider the full system of N spins as an ensemble of Q = N/M clusters of M spins each, with each cluster evolving independently under H. The collective spin is then $\mathbf{S}={\sum }_{q=1}^{Q}{\mathbf{S}}^{\left(q\right)}$, where ${\mathbf{S}}^{\left(q\right)}={\sum }_{{\mathbf{s}}_{i}\in {c}^{\left(q\right)}}{\mathbf{s}}_{i}$ is the total spin of cluster c^(q). Within each cluster, positions $\left\{{\mathrm{x}}_{i}^{\left(q\right)}\right\}$ are assigned by finding the M − 1 closest, PPP-distributed points to ${\mathrm{x}}_{1}^{\left(q\right)}$, which is taken as the origin. For a product-state initial condition, the S^(q) are independent, so that ${{\Gamma}}_{S}={\sum }_{q}{{\Gamma}}_{{S}^{\left(q\right)}}$. We compute U(T) using equation (4) and matrix exponentiation for Q ∼ 10⁴ clusters, to find ${{\Gamma}}_{{S}^{\left(q\right)}}\left(T\right)$, ⟨S^(q)(T)⟩, and its derivatives.

We first consider the case of dc magnetometry, in which B(t) = (0, 0, B_dc), the initial state is |ϕ₀⟩ = |+x⟩, and $\mathcal{B}={B}_{\mathrm{d}\mathrm{c}}$. Representative results are shown in figure 1. The mean amplitude of oscillation shows a steady and seemingly exponential decline, while the elements of the covariance matrix saturate, with the result that the imprecision reaches a minimum at which $\langle \delta {B}_{\mathrm{d}\mathrm{c}}^{2}\rangle VT/\left(2{\mu }_{0}\hslash \right)\approx 0.7$ at a finite time T_optϒ_dd ≈ 0.5. As shown in figure 2, the limiting sensitivity improves with increasing M, but saturates at about M = 6. Similarly, s larger than 1/2 provides an advantage that appears to saturate about s = 1. Simulations (not shown) with other conditions, including different spin quantum number, gyromagnetic ratio, field strength and density, find very similar limiting sensitivities, strongly suggesting an implementation-independent limit for fixed spin-precession sensors.

Zoom In Zoom Out Reset image size

To understand this limit, we return to analytical methods. We shift to a coordinate system rotating about the field at the Larmor frequency, in effect working in the interaction picture [47] with H_dd as interaction Hamiltonian. Rotating-frame spin operators ${\tilde {\boldsymbol{s}}}_{i}\left(t\right)\equiv {U}_{L}\left(t\right){\mathbf{s}}_{i}{U}_{L}^{{\dagger}}\left(t\right)$ are defined via the unitary U_L(t) ≡ exp[iH_Lt/ℏ], which is the inverse of time evolution under H_L alone. U_L(t) produces Larmor precession ${\tilde {\boldsymbol{s}}}_{i}\left(t\right)={R}_{z}\left({\omega }_{L}t\right){\boldsymbol{s}}_{i}\left(t\right)$, where the matrix R_z(θ) describes rotation about axis z by angle θ. With this transformation, time evolution is governed by the rotating-frame Hamiltonian

$$\begin{align}\hfill \tilde {H}& ={U}_{L}\left(t\right)H{U}_{L}^{{\dagger}}\left(t\right)+i\dot {{U}_{L}}\left(t\right){U}_{L}^{{\dagger}}\left(t\right)\hfill \\ \hfill & ={U}_{L}\left(t\right)H{U}_{L}^{{\dagger}}\left(t\right)-{H}_{L}\hfill \\ \hfill & =\hslash \sum _{i\ne j}\frac{{{\Upsilon}}_{dd}}{{s}^{2}\rho {r}_{ij}^{3}}\left[{\tilde {\mathbf{s}}}_{i}\cdot {\tilde {\mathbf{s}}}_{j}-3\left({\tilde {\mathbf{s}}}_{i}\cdot {\mathbb{R}}_{ij}\right)\left({\tilde {\mathbf{s}}}_{j}\cdot {\mathbb{R}}_{ij}\right)\right].\hfill \end{align} \tag{ 7 }$$

In the rotating frame, the spins evolve only under their mutual coupling, which consists of two terms: the spin exchange term $\propto {\tilde {\boldsymbol{s}}}_{i}\cdot {\tilde {\boldsymbol{s}}}_{j}$ has no explicit time dependence and conserves total angular momentum, whereas the term $\propto \left({\tilde {\boldsymbol{s}}}_{i}\cdot {\mathbb{R}}_{ij}\right)\left({\tilde {\boldsymbol{s}}}_{j}\cdot {\mathbb{R}}_{ij}\right)$ allows angular momentum to escape 'to the lattice' and has coefficients that oscillate at the Larmor frequency. We note that the sensitivity limit given in equation (2) is unaffected by the frame shift, due to optimization over the read-out direction n.

The periodicity of the rotating-frame Hamiltonian motivates a Kapitza approach in which we divide the dynamics of $\left\{{\tilde {\boldsymbol{s}}}_{i}\right\}$ into a slowly-varying 'secular' part and a 'micro-motion' part oscillating at ω_L. The micromotion is smaller than the secular part by a factor ∼ϒ_dd/ω_L, and for sufficient ω_L becomes negligible (see appendix B). The secular motion is governed by the Larmor-cycle-averaged Hamiltonian

$$\begin{align}\hfill \overline{H}\left(t\right)& \equiv \frac{1}{{T}_{L}}{\int }_{t}^{t+{T}_{L}}\mathrm{d}{t}^{\prime }\;\tilde {H}\left({t}^{\prime }\right)\hfill \\ \hfill & =\hslash \sum _{i\ne j}\frac{{{\Upsilon}}_{dd}}{{s}^{2}\rho {r}_{ij}^{3}}\frac{1-3{\mathbb{R}}_{ij,z}^{2}}{2}\left(3{\tilde {s}}_{i,z}{\tilde {s}}_{j,z}-{\tilde {\mathbf{s}}}_{i}\cdot {\tilde {\mathbf{s}}}_{j}\right),\hfill \end{align} \tag{ 8 }$$

where the subscript z indicates the component along , i.e., along the dc field. The x and y components are lost in the cycle average. We note that ω_L no longer appears.

We can now understand the effect of density in the secular regime, using a strategy from renormalization group (RG) theory. We imagine dividing the sensor volume into λ equal sub-volumes, while also increasing the density by a factor λ. We indicate post-transformation quantities with primes, e.g. ρ' = λρ. If P_ρ is self-similar, in the sense that the statistical distribution of {λ^1/3(x'_i − x'_j)} within a sub-volume the same as that of {x_i − x_j} in the full volume, and again assuming edge effects are negligible, the sub-volumes now represent λ independent, reduced-scale realizations of the original sensor. A PPP for example has such self-similarity.

For a given configuration {x_i'} = {λ^−1/3x_i}, the Hamiltonian is ${\overline{H}}^{\prime }=\lambda \overline{H}$, implying a speed-up of the rotating-frame dynamics by a factor λ. When averaged over {x_i}, this produces faster evolution of spin means $\langle {\tilde {s}}_{i}^{\prime }\left(t\right)\rangle =\langle {\tilde {s}}_{i}\left(\lambda t\right)\rangle$ and correlators $\langle {\tilde {s}}_{i}^{\prime }\left(t\right){\tilde {s}}_{j}^{\prime }\left(t\right)\rangle =\langle {\tilde {s}}_{i}\left(\lambda t\right){\tilde {s}}_{j}\left(\lambda t\right)\rangle$. Considering then the collective spin S, which sums the λ sub-volumes, we find $\langle {\tilde {S}}^{\prime }\left(t\right)\rangle =\lambda \langle \tilde {S}\left(\lambda t\right)\rangle$ and Γ_S'(t) = λΓ_S(λt). Inserting into equation (2), we find T^'_opt = λ⁻¹T_opt and thus the same limiting sensitivity, independent of λ and thus of ρ.

The specific value of γ is similarly irrelevant: a change to γ' = λγ, through dynamics and quantum noise scaling, gives $\langle {\tilde {S}}^{\prime }\left(t\right)\rangle =\lambda \langle \tilde {S}\left(\lambda t\right)\rangle$, Γ_S'(t) = λΓ_S(λt), T^'_opt = λ⁻¹T_opt, and again the same limiting sensitivity. We see now that, within the secular regime ω_L ≫ ϒ_dd, the only factors that can influence E_R are s and fundamental constants. This analytic result confirms the existence of an implementation-independent limit suggested by the numerical calculations.

Returning to numerical methods, we can show that the secular regime limit is in fact the global optimum. Again we use the cluster simulation. We study a range of ω_L/ϒ_dd, see figure 3, and find that outside of the secular regime $T\langle \delta {B}_{\mathrm{d}\mathrm{c}}^{2}\rangle$ is about a factor of two larger than in the secular regime.

Zoom In Zoom Out Reset image size

Using $\overline{H}$ we can also consider RF magnetometry [19], which requires only minor modification to the above discussion. Now B(t) = (B_RF cos ωt, B_RF sin ωt, B_dc), the initial state is |ϕ₀⟩ = |+z⟩, and the unknown is $\mathcal{B}={B}_{\mathrm{R}\mathrm{F}}$. As before, we take {x_i} to be PPP-distributed. In the rotating frame, the RF field appears fixed, and contributes a term $-\gamma \hslash {B}_{\mathrm{R}\mathrm{F}}{\sum }_{i}{\tilde {\boldsymbol{s}}}_{i}\cdot \hat{x}$ to both equations (7) and (8). For γ|B_RF| ≪ ϒ_dd, i.e. for weak-field detection, the contribution of this term to the dynamics is small, and the scale-invariance arguments proceed as above. Here also, this shows analytically that there is an implementation-independent limit on the sensitivity. Using equation (8) we numerically evaluate ${{\Gamma}}_{\tilde {S}}$ and ${\partial }_{{B}_{\mathrm{R}\mathrm{F}}}\langle \tilde {\boldsymbol{S}}\cdot \boldsymbol{n}\rangle$ by the cluster expansion as above. Results, shown in figure 4, indicate a lower-bound of E_R/ℏ ≈ 1/4.

Zoom In Zoom Out Reset image size

Thus far, the discussion has concerned only spins with fixed positions. For mobile spins, we consider a set of trajectories {x_i(t)}, and the RG argument proceeds as above if the trajectories' distribution is self-similar in the sense that {λ^1/3[x'_i(t) − x'_j(t)]} has the same distribution as {x_i(λt) − x_j(λt)}. This describes sub-diffusive (Δx)³ ∝ Δt transport, as opposed to diffusive (Δx)² ∝ Δt or ballistic Δx ∝ Δt transport. As such, we cannot directly use the RG argument to establish an implementation-independent limit for vapor- or gas-phase spin-precession sensors. Nonetheless, we recover scale-invariance in two scenarios.

First, precisely because the decoherence rate ϒ_dd ∝ ρ grows faster with density than does the t ∝ ρ^2/3 (diffusive) or t ∝ ρ^1/3 (ballistic) transport time across the inter-particle spacing, for sufficiently large densities (or sufficiently slow transport) the spins will appear effectively immobile and can be treated as fixed, so that the conclusions given above for fixed spins apply. Second, in the opposite extreme of highly-mobile, weakly-coupled spins, we can expect short-range collisional processes to cause spin depolarization faster than do long-range spin–spin interactions. Such collisions produce a decoherence rate Γ ∝ ρ in both diffusive and ballistic regimes [48]. As a consequence, E_R becomes independent of density ρ. The collisional processes that produce decoherence do not, however, appear to scale with γ as would be required for a species-independent sensitivity limit [49]. For this regime, we thus expect a density-independent limit for any given sensing species, but we do not expect these limits to be the same for different species.

Conclusions

Acknowledgments

Appendix A.: Optimal readout

To efficiently evaluate equation (2), it is convenient to note that by parametrizing n = (cos θ, sin θ, 0), the variance of the estimate of B can be written

$$\begin{align}\hfill \langle \delta {B}^{2}\rangle & ={\mathrm{min}}_{\mathbf{n}}\frac{n{{\Gamma}}_{S}{n}^{T}}{{\left\vert {\partial }_{\mathcal{B}}\langle \mathbf{S}\cdot \mathbf{n}\rangle \right\vert }^{2}}\hfill \\ \hfill & ={\mathrm{min}}_{\theta }\frac{\left(\mathrm{cos} \theta ,\mathrm{sin} \theta ,0\right){{\Gamma}}_{S}{\left(\mathrm{cos} \theta ,\mathrm{sin} \theta ,0\right)}^{T}}{{\left[{\partial }_{\mathcal{B}}\left(\langle {S}_{x}\rangle \mathrm{cos} \theta +\langle {S}_{y}\rangle \mathrm{sin} \theta \right)\right]}^{2}}\hfill \\ \hfill & =\frac{\vert {\Gamma}{\vert }_{2}}{Z}\hfill \end{align} \tag{ A.1 }$$

where

$$\begin{equation}Z\equiv {\left({\partial }_{\mathcal{B}}\langle {S}_{x}\rangle \right)}^{2}\mathrm{v}\mathrm{a}\mathrm{r}\;{S}_{y}-2\;\mathrm{c}\mathrm{o}\mathrm{v}\left({S}_{x},{S}_{y}\right)\left({\partial }_{\mathcal{B}}\langle {S}_{x}\rangle \right)\left({\partial }_{\mathcal{B}}\langle {S}_{y}\rangle \right)+{\left({\partial }_{\mathcal{B}}\langle {S}_{y}\rangle \right)}^{2}\mathrm{v}\mathrm{a}\mathrm{r}\;{S}_{x}\end{equation} \tag{ A.2 }$$

and |⋅|₂ indicates the determinant of the upper left 2 × 2 sub-matrix.

Appendix B.: Kapitza-theory dynamics

To understand the conditions under which spin-interaction dynamics will simplify due to Larmor precession, we adapt the classical Kapitza method [55, 56] to the spin problem at hand. We begin by writing the equations of motion for the spins, in the frame rotating at ω_L and thus governed by Hamiltonian of equation (6). We have

$$\begin{equation}\frac{\mathrm{d}}{\mathrm{d}t}\tilde {\mathbf{s}}\left(i\right)=\sum _{k\ne i}\frac{{\gamma }^{2}{\hslash }^{2}{\mu }_{0}}{4\pi {r}_{ik}^{3}}\left[3\left({\tilde {\mathbf{s}}}^{\left(i\right)}{\times}{\tilde {\mathbb{R}}}^{\left(\mathrm{ik}\right)}\right)\left({\tilde {\mathbf{s}}}^{\left(k\right)}\cdot {\tilde {\mathbb{R}}}^{\left(\mathrm{ik}\right)}\right)+\left({\tilde {\mathbf{s}}}^{\left(i\right)}{\times}{\tilde {\mathbf{s}}}^{\left(k\right)}\right)\right]\end{equation} \tag{ B.1 }$$

where $\tilde {\boldsymbol{r}}\equiv {R}_{z}\left({\omega }_{L}t\right)\boldsymbol{r}$ is the rotated vector joining the spins. The precession period is T ≡ 2π/ω_L.

Collecting all spin components into a single vector $\boldsymbol{z}={\bigoplus}_{i}{\tilde {\boldsymbol{s}}}^{\left(i\right)}$, we note that equation (B.1) has the form

$$\begin{equation}\frac{\mathrm{d}}{\mathrm{d}t}{z}_{i}={\mathcal{R}}_{ijk}\left(t\right){z}_{j}{z}_{k}\end{equation} \tag{ B.2 }$$

where the tensor of coefficients $\mathcal{R}$ is periodic: $\mathcal{R}\left(t+T\right)=\mathcal{R}\left(t\right)$. It is convenient to identify a time t_α as the start of a cycle, and divide $\mathcal{R}\left(t\right)$ as $\mathcal{R}\left(t\right)=\overline{\mathcal{R}}+\circ{\mathcal{R}}\left(t\right)$ where the cycle-averaged part is $\overline{\mathcal{R}}\equiv {T}^{-1}{\int }_{{t}_{\alpha }}^{{t}_{\alpha }+T}\mathcal{R}\left(t\right)\;\mathrm{d}t$.

We write z = p + q, where p is the slowly-varying 'secular motion' and q is the small and rapidly-varying 'micro-motion', defined as the solution of

$$\begin{equation}\frac{\mathrm{d}}{\mathrm{d}t}{q}_{i}={\circ{\mathcal{R}}}_{ijk}\left(t\right)\end{equation} \tag{ B.3 }$$

with zero cycle-average: ${\int }_{{t}_{\alpha }}^{{t}_{\alpha }+T}\boldsymbol{q}\left(t\right)\;\mathrm{d}t=\mathbf{0}$. We can write the formal solution

$$\begin{equation}{q}_{i}\left(t\right)={\int }_{{t}_{\alpha }}^{t}{\circ{\mathcal{R}}}_{ijk}\left({t}^{\prime }\right)\mathrm{d}{t}^{\prime }\;\;{p}_{j}\left({t}_{\alpha }\right){p}_{k}\left({t}_{\alpha }\right).\end{equation} \tag{ B.4 }$$

We note that q ∼ T, and thus q becomes small for large ω_L.

Using the smallness of q we expand the rhs of equation (B.2), as applies to the time period t ∈ [t_α, t_α + T) to find

$$\begin{equation}\frac{\mathrm{d}}{\mathrm{d}t}{z}_{i}={\mathcal{R}}_{ijk}\left(t\right){z}_{j}\left({t}_{\alpha }\right){z}_{k}\left({t}_{\alpha }\right)+{q}_{l}\left(t\right){\left[{\partial }_{{z}_{l}}{\mathcal{R}}_{ijk}\left(t\right){z}_{j}{z}_{k}\right]}_{\mathbf{z}=\mathbf{p}\left({t}_{\alpha }\right)}+O{\left(q\right)}^{2}.\end{equation} \tag{ B.5 }$$

We drop the doubly-small O(q)² term and integrate over one cycle to find the cycle-averaged rate of change

$$\begin{equation}\frac{\mathrm{d}{p}_{i}}{\mathrm{d}t}\approx \frac{{\Delta}{p}_{i}}{T}={\overline{\mathcal{R}}}_{ijk}\left(t\right){p}_{j}\left({t}_{\alpha }\right){p}_{k}\left({t}_{\alpha }\right)+\frac{1}{T}{\int }_{{t}_{\alpha }}^{{t}_{\alpha }+T}\mathrm{d}t{\int }_{{t}_{\alpha }}^{t}\mathrm{d}{t}^{\prime }{\circ{\mathcal{R}}}_{ijk}\left({t}^{\prime }\right){p}_{j}\left({t}_{\alpha }\right){p}_{k}\left({t}_{\alpha }\right){\left[{\partial }_{{z}_{l}}{\mathcal{R}}_{ijk}\left(t\right){z}_{j}{z}_{k}\right]}_{\mathbf{z}=\mathbf{p}\left({t}_{\alpha }\right)},\end{equation} \tag{ B.6 }$$

which now refers only to p. The second term describes the effect of micromotion on the secular dynamics.

We note that the first term in equation (B.6) scales as ϒ_dd, and contains both the ${\tilde {\boldsymbol{s}}}^{\left(i\right)}{\times}{\tilde {\boldsymbol{s}}}^{\left(k\right)}$ factor that produces spin-exchange and $\left({\tilde {\boldsymbol{s}}}^{\left(i\right)}{\times}{\hat{z}\tilde {\mathbb{R}}}_{z}^{\left(ik\right)}\right)\left({\tilde {\boldsymbol{s}}}^{\left(k\right)}\cdot {\hat{z}\tilde {\mathbb{R}}}_{z}^{\left(ik\right)}\right)$ obtained by cycle-averaging $\left({\tilde {\boldsymbol{s}}}^{\left(i\right)}{\times}{\tilde {\mathbb{R}}}^{\left(ik\right)}\right)\left({\tilde {\boldsymbol{s}}}^{\left(k\right)}\cdot {\tilde {\mathbb{R}}}^{\left(ik\right)}\right)$. This latter factor is responsible for loss of angular momentum to the center of mass degrees of freedom. The second term in equation (B.6) scales as ${{{\Upsilon}}_{dd}}^{2}/{\omega }_{L}$, i.e., smaller than the first by a factor ϒ_dd/ω_L. It is this smallness that justifies using the cycle-averaged Hamiltonian of equation (8) for large ω_L.