We theoretically analyze an optically spin-polarized collection of atoms, which serves as a basis for atomic sensors. Assuming that the intrinsic atomic spin projection noise is negligible, we provide the closed-form autocorrelation function and the power spectral density (PSD) of the solution to a noisy version of an optically pumped Bloch equation, wherein each component of the external magnetic field is subjected to white noise. We conclude that noise in the bias B-field direction does not affect the autocorrelation function, up to first order in white noise covariance amplitudes. Moreover, the noise terms for the remaining two axes make different contributions to the magnetic noise-driven spin PSD; in particular, the contribution corresponding to noises perpendicular to the probing direction dominates at high frequencies. Some results concerning the second (and higher)-order terms are given. In particular, we anticipate a decrease in the effective Larmor frequency despite an increase in the magnetic field magnitude in the case of anisotropic transversal B-field noises. The analytic results are supported by Monte Carlo simulations employing the Euler–Maruyama method. The analytic methodology is applied to the case of a Bell–Bloom magnetometer, which reveals a non-linearity in the PSD of the magnetometer output and also a broadening effect due to magnetic field noise.

Magnetic resonances, such as nuclear magnetic resonance (NMR) and electron paramagnetic resonance (EPR), serve as basic principles for modern atomic sensors, such as atom spin gyroscopes (ASGs)1,2 and optically pumped atomic magnetometers (OPAMs),3–5 respectively. The sensitivity of an atomic sensor is determined by various noise sources, such as the photon shot noise, the atomic spin projection noise (intrinsic spin noise), and the noise induced by undesired magnetic fields from the surrounding electronics (e.g., electric heaters) and the environment. Among these, photon shot noise and atomic spin projection noise set the fundamental sensitivity limit of the atomic sensor. Atomic spin projection noise can be reduced by utilizing spin-squeezed states, which are produced by quantum nondemolition measurements.6,7 Photon shot noise can be reduced by utilizing squeezed light.8,9 However, the sensitivities tend to be degraded from the fundamental limit due to magnetic field noises. Such noises from various electronics can be suppressed by employing sensor components generating low levels of magnetic field noise. For example, electric heaters for heating up atom cells can be designed to minimize the resulting magnetic interference. Such designs can be realized by double-layered polyimide films10 or by employing the SOS-CMOS technology.11 On the other hand, OPAMs operated in the Earth’s magnetic field cannot avoid the white-like noise from their environments, including the surrounding electronics. In ASGs and OPAMs, as the spin collections are polarized by optical pumping, the magnetic field noise-induced spin noise can dominate the intrinsic atomic spin projection noise. For a fully polarized spin system, the atomic spin projection noise is given as σin=F/(2Na), where Na is the number of atoms in the collection and F does not carry units of .12 In contrast, the magnetic field noise-induced spin noise depends on the degree of the spin polarization. For example, we show that the magnetic field noise-induced spin noise can be written as σexγPΓxx+Γyy/(4Γ2), given that the Larmor frequency is much larger than the transversal spin relaxation rate Γ2, where P = 2Sz is the degree of polarization, γ is the gyromagnetic ratio, and Γxx and Γyy are the magnetic field noise amplitudes along the x and y directions as formally defined in Eq. (6): see Eq. (11c). In this paper, we focus on the regime where σexσin. The precise mathematical formulation for such a regime will become clearer later after we compare the noise powers for these two types of noises in Sec. II: see Eq. (12).

Although the field-noise-driven spin precession is complex, we can still analytically investigate the spin noise in the presence of white magnetic field noise and provide some new insights into atomic sensors. In this regard, we investigate the response of a magnetic resonance system subjected to external B-field white noise by analyzing the corresponding stochastic differential equation.

In this paper, we focus on the Bloch equation, which is one of the simplest phenomenological equations describing magnetic resonance,2,13–16
dSdt=γS×B+δBRSs2ΓreS,
(1)
where S is the average spin (e.g., the average on a hot 87Rb vapor cell), γ is the gyromagnetic ratio, R is the optical pumping rate, s is the photon spin vector,17 and Γre is the spin relaxation matrix. The arrow notation emphasizes that the corresponding vectors are three-dimensional. In this paper, we assume that δB is white noise, B is along the z-direction, and the pumping beam is also z-directional with a σ+ polarization, i.e., s=ẑ. This equation includes an optical pumping term so that it describes a motion of polarized spins in the presence of the white magnetic field noise. We note that the stochastic Bloch equation for unpolarized spins was previously analyzed in Ref. 16 by means of the stochastic path integrals. As it turns out in this paper, in the spin-polarized case Sz = P/2, even with a small external magnetic field noise, the spin noise is mainly determined by the fluctuations of the polarized spin driven by the magnetic field noise. A similar remark applies when we compare this paper to a related experiment on an activated crystal.18 We also note that as Eq. (1) focuses on spin noise driven by a noisy magnetic field, it contrasts to studies such as Ref. 19, which incorporated spontaneous fluctuations due to, e.g., collisions between atoms, or Ref. 20, which focused on noise due to atomic diffusion. We note that studies tend to describe the power spectral density (PSD) of the spin noise Sx, arising from various noise sources, as a Lorentzian function or a sum of nearly Lorentzian functions; see Refs. 12 and 2123.

To the best of our knowledge, no one has explicitly noted the fact that the PSD should depend on the noise direction (for cases in which such a direction is fixed): x-directional noises and y-directional noises lead to different PSDs of Sx. A close inspection of the magnetic resonance system reveals that such anisotropy would depend on the frequency range of interest. For the high-frequency range of the PSD, which will correspond to high-frequency magnetic field noises, the Larmor precession axis B+δB is hastily altered along the field noise direction before a spin component along that direction is sufficiently created. Hence, the spin motion is mainly distributed on the axis perpendicular to the field noise axis [see Fig. 1(a)]. In other words, spin noises parallel to the magnetic field noise are suppressed. On the other hand, when a low-frequency magnetic field noise comes in, the precession axis does not change hastily; therefore, sufficient spin precession can occur, enabling spin noise parallel to the magnetic field noise to come into the picture [see Fig. 1(b)].

FIG. 1.

(a) Intuition of the system’s response to high-frequency components of an x-directional white noise. The precession axis (red) is blurred, leading to transitory motions of the spin vector (blue). Being transitory, these motions are mainly perpendicular to the noise direction. (b) Intuition of the response to low-frequency components. In this case, the precession axis is somewhat fixed; therefore, the spin vector can precess around a well-defined axis for times much shorter than the fastest timescale of the low-frequency noise. Therefore, spin motions parallel to the noise direction occur. (c) A sample transversal spin trajectory corresponding to a high-frequency noise and a low-frequency noise, drawn by a computer simulation. The horizontal axis S of the graph denotes the spin component parallel to the field noise direction, whereas the vertical axis S corresponds to the component perpendicular to the noise. Note that ⟨SS⟩ appears to be non-zero, as can be verified from Eq. (10b).

FIG. 1.

(a) Intuition of the system’s response to high-frequency components of an x-directional white noise. The precession axis (red) is blurred, leading to transitory motions of the spin vector (blue). Being transitory, these motions are mainly perpendicular to the noise direction. (b) Intuition of the response to low-frequency components. In this case, the precession axis is somewhat fixed; therefore, the spin vector can precess around a well-defined axis for times much shorter than the fastest timescale of the low-frequency noise. Therefore, spin motions parallel to the noise direction occur. (c) A sample transversal spin trajectory corresponding to a high-frequency noise and a low-frequency noise, drawn by a computer simulation. The horizontal axis S of the graph denotes the spin component parallel to the field noise direction, whereas the vertical axis S corresponds to the component perpendicular to the noise. Note that ⟨SS⟩ appears to be non-zero, as can be verified from Eq. (10b).

Close modal

We note that such a frequency-dependence can be confirmed by computer simulations. Figure 1(c) shows a sample transversal spin noise trajectory corresponding to two types of transversal magnetic field noises: a high-frequency noise (blue) and a low-frequency noise (purple). Details of the simulation are found in  Appendix A.

In this paper, we validate the above-mentioned intuitions by analytically solving Eq. (1) up to first order in B-field noise covariance amplitudes (these amplitudes will be denoted as Γαβ, where α, β ∈ {x, y, z}). We give corrections to the PSD, which are second order in the covariance amplitudes. We believe that such higher-order information cannot be obtained if we approximate Eq. (1) by a Langevin equation, which is arguably the simplest approach for modeling noise: assuming δBz(t) = 0 for all t and approximating Sz to be a constant transform the equation into a Langevin form. Such higher-order terms exhibit some qualitatively different features from their first-order counterparts: the higher-order line shapes are much more complicated, and they exhibit longitudinal–transversal coupling, which is absent in the first-order terms. Finally, the analytic methodology is applied to the case of a Bell–Bloom magnetometer, which is an instance of an atomic sensor utilizing a collection of atomic spins.

Summary: We consider the model given by Eq. (2). The autocorrelation function of Sx is provided by Eq. (5), for which the initial conditions can be given arbitrarily. If the system is at its steady-state, we have the steady-state conditions [Eq. (10)] (a suggestive notation t0 → ∞ for steady-state values is employed in this equation). Plugging in Eq. (10) into Eq. (5), we obtain the main formulas in Eq. (11). In particular, Eq. (11b) formally proves the intuition suggested in Fig. 1. Equation (11c) leads to Eq. (12), which quantifies the regime where σexσin as alluded in Sec. I.

Above-mentioned equations are valid up to first order in the field noise amplitudes Γxx, Γxy, … , which are defined in Eq. (6). Higher-order corrections are presented in Subsection II B. In particular, second-order corrections to the power spectral density (PSD) of Sx are presented by Eq. (20), for which we assume Eq. (16). In addition, we conclude that Gaussian white noise leads to a decrease in the Larmor frequency [Eq. (23)]; this happens when an anisotropic B-field noise is injected perpendicularly to the nominal B-field axis. An intuitive explanation on this phenomenon is provided in Fig. 3.

We consider a situation in which B=Bẑ+δBx(t)δBy(t)δBz(t)T, where δBα(t) are mean-zero white noise. We note that the average spin S of the collection experiences a spatially homogeneous magnetic field fluctuation as we are implicitly assuming that the size of the atomic vapor cell is much smaller than the spatial scale of B-field fluctuations, in contrast to the inhomogeneous fluctuations, which lead to additional spin relaxations.24 We assume a diagonal relaxation matrix with Γ2, Γ2, and Γ1 as its three consecutive diagonal entries, where Γ2 and Γ1 are the transversal and longitudinal spin relaxation rates excluding the contributions of optical pumping. Note that the relaxation matrix is transversally (x, y) isotropic, although it is not fully isotropic as Γ1 may differ from Γ2.13 In this case, the Bloch equation can be written in the following matrix form:
ddtSxSySz=Γ2RγB0γBΓ2R000Γ1R+γ0δBzδByδBz0δBxδByδBx0SxSySz+00R/2.
(2)
From now on, we denote Γ1+R and Γ2+R as Γ1 and Γ2, respectively. We assume γB ≠ 0 and Γ1, Γ2 to be positive. Denoting the right-hand side of Eq. (2) as (A1+A2(t))S+c and by defining a coordinate change as S̃(t)exptA1c00S1, the equation is transformed into the following homogeneous form (boldface zeros denote zero matrices with varying dimensions, whereas non-bold zeros denote 1 × 1-dimensional zero blocks):
dS̃dt=exptA1c00A2000exptA1c00S̃(t),
(3)
where A2000 is meant to be a 4 × 4 matrix. Hereafter, we denote the right-hand side of the above equation as Ã(t)S̃(t).
We have the following second-order approximation for the autocorrelation function of S̃, where 1 denotes the 4 × 4 identity matrix and h.o.t. means higher-order terms:
S̃(t0+τ)S̃(t0)T=1+t0t0+τt0tÃ(t)Ã(t)dtdt+(h.o.t.)×S̃(t0)S̃(t0)T.
(4)

From now on, the analysis proceeds in two steps. The first step is to calculate the time-evolution operator in the parentheses in Eq. (4), which is achieved in Eq. (5). The second step is to calculate the initial value S̃(t0)S̃(t0)T when the system reaches the steady state (denoted by t0 → ∞); the result of this step is summarized in Eq. (10). In addition, we note that an additional analysis involving the higher-order terms in Eq. (4) will be discussed later.

Calculating the time-evolution operator and returning to the original coordinates, we have the following expression for the autocorrelation function of Sx:
Sx(t0+τ)Sx(t0)=eΓ2τcos(τω)+γ2(ΓxxΓyy)4ωsin(τω)γ2(Γxx+Γyy+2Γzz)4τcos(τω)Sx(t0)2+eΓ2τ1+γ2Γxy2ωsin(τω)γ2(Γxx+Γyy+2Γzz)4τsin(τω)×Sx(t0)Sy(t0)+CSx(t0)Sz(t0)+DSx(t0),
(5)
where we denoted the Larmor frequency γB as ω. Here, we assumed that the noises are white with their covariances given by delta functions as
δBα(t)δBβ(t)=Γαβδ(tt)
(6)
for symmetric coefficients Γαβ; that is, the correlation time of the B-field noises is much shorter than any relevant time scales of the system. C and D are terms that vanish when the coefficients Γαβ for αβ are zero (the full expressions are given in  Appendix C).

Approximation (5) for the autocorrelation function of Sx is a priori valid only for short evolution times τ with τÃ1 (or more precisely, t0t0+τdtÃ(t)1). Meanwhile, we argue that Eq. (5) also holds for τ large enough so that Γ1τ, Γ2τ ≫ 1. First, the right-hand side of the equation clearly vanishes whenever Γ1τ, Γ2τ ≫ 1, at least up to first order in Γαβ. In addition, the left-hand side can be approximated by ⟨Sx(t0 + τ)⟩⟨Sx(t0)⟩ for such large τ because physical intuition tells us that Sx(t0) and Sx(t0 + τ) will be statistically uncorrelated; any deviation from the expectation at time t0 will be exponentially damped out in time, thereby not affecting the spin at time t0 + τ. Therefore, both sides of the equation vanish up to first order in Γαβ, assuming that the system is at its steady-state. We emphasize that this argument hinges on the change of coordinates into S̃, which resulted in the damping exponential eΓ2τ terms shown in the right-hand side of Eq. (5).

To sum up, Eq. (5) is valid for values of τ satisfying τÃ1 or Γ1τ, Γ2τ ≫ 1. As a result, it will be valid for all τ ≥ 0 as long as these two regions overlap: see Fig. 2. Such an overlap happens when the noise is small enough (i.e., Ã/Γ11 and Ã/Γ21) since in this case the region τÃ1 depicted in Fig. 2 will be shifted toward τ → ∞.

FIG. 2.

Equation (5) is valid for all τ ≥ 0 for sufficiently small noises.

FIG. 2.

Equation (5) is valid for all τ ≥ 0 for sufficiently small noises.

Close modal

In any case, Eq. (4) is exact if we consider all the higher-order terms.25 Later, in this section, the effect of the first non-zero higher-order term is considered; if the resulting correction terms are much smaller than the first-order terms in Γαβ, we can, in retrospect, conclude that the noise was small enough for the approximation in (4) to be valid.

So far, we did not assume that the system is in its steady state. Now, we determine the relevant steady state values Sx(t0)2, ⟨Sx(t0)Sy(t0)⟩, ⟨Sx(t0)Sz(t0)⟩, and ⟨Sx(t0)⟩ after all the irrelevant transients are damped out (i.e., t0 → ∞). This can be achieved by enlarging the main vector variable further to include the second-order terms S(t)Sx(t)Sx(t)Sx(t)Sy(t)Sz(t)Sz(t)T as in the following equation:
ddtSS1=A11+1A1c1+1c00A1c000+A21+1A2000A20000SS1.
(7)
Here, 1 denotes the 3 × 3 identity matrix and ⊗ denotes the Kronecker product. Denoting the terms in the parentheses as A1+A2 and defining a coordinate change S̃(t)etA1S(t)S(t)1T, we conclude that the second-order approximation for S̃(t) similar to Eq. (4) coincides with the first-order approximate solution to the following differential equation [Eq. (8b)]. Although such an equivalence is a priori valid only for sufficiently short time evolutions, this limitation is superficial because (8b) is independent of the initial time t0,
ddtS̃(t)=etA1Γ̃etA1S̃(t),
(8a)
Γ̃t0tA2(t)A2(t)dt.
(8b)
Here, Γ̃ is a 13 × 13 constant matrix that is independent of t > t0. It is a matrix that has linear combinations of γ2Γαβ as its components (the full expression is given in  Appendix C).
In conclusion, up to first order in Γαβ (i.e., up to second order in the noise amplitudes), we have the following equation in the original coordinates:
ddtSS1=(A1+Γ̃)SS1.
(9)

Incidentally, every eigenvalue of the 12 × 12 submatrix of A1+Γ̃ at the top left corner has a negative real part when Γαβ = 0; this mathematically proves that the expectations S and S will indeed converge to an equilibrium, at least for sufficiently small Γαβ values.

These steady state values can be obtained by solving a 12 × 12 linear equation; some are provided in the following equations, written up to first order in Γαβ. We used the suggestive symbol ∞ in Eq. (10) to emphasize that the equation represents the steady state values,
Sx()Sx()=γ2R216Γ12Γ2ω2+Γ22×ω2Γxx2Γ2ωΓxy+(ω2+2Γ22)Γyy,
(10a)
Sx()Sy()=γ2R2ωΓxx2Γ2ΓxyωΓyy16Γ12ω2+Γ22,
(10b)
Sx()Sz()=γ2R2ωΓyz+Γ2Γzx8Γ12ω2+Γ22,
(10c)
Sx()=γ2RωΓyz+Γ2Γzx4Γ1ω2+Γ22,
(10d)
Sz()=R2Γ1γ2R4Γ12Γxx+Γyy.
(10e)
Note that the xy-product given in Eq. (10b) vanishes when Γxx = Γyy and Γxy = 0, as expected due to symmetry. In addition, since ω2Γxx2Γ2ωΓxy+Γ22Γyy0, Eq. (10a) gives non-negative values. Moreover, the steady state value for Sz has a first-order correction term, which cannot be easily justified by the Langevin approach; the Langevin approach, in order to be justified, requires Sz to be approximated by its a priori equilibrium value R/2Γ1. Finally, we note that (10c) and (10d) vanish whenever the correlation coefficients Γyz and Γzx are zero.
Plugging the relevant steady state values into (5), we finally obtain the autocorrelation function at steady state. The autocorrelation function, the one-sided PSD, and the total power of Sx(t) are given in the following equations, which are written up to first order in Γαβ. It turns out that there is no first-order effect of non-zero Γzz, Γyz, or Γzx. That is, the system under consideration is insusceptible to longitudinal noises,
Sx(t0+τ)Sx(t0)=Cω2Γxx2Γ2ωΓxy+(2Γ22+ω2)Γyycos(ωτ)+CΓ2(ωΓxx2Γ2ΓxyωΓyy)×sin(ωτ)for τ0, at steady state, 
(11a)
PSDSx(f)=40Sx(t0+τ)Sx(t0)cos(2πfτ)dτ=γ2R2ω2Γxx2ωΓ2Γxy+(Γ22+Ω2)Γyy2Γ12(Ωω)2+Γ22(Ω+ω)2+Γ22,
(11b)
Sx2=0PSDSx(f)df=γ2R2ω2Γxx2Γ2ωΓxy+(ω2+2Γ22)Γyy16Γ12Γ2ω2+Γ22.
(11c)
In Eq. (11), f denotes the argument of the PSD function, Cγ2R2eΓ2τ16Γ12Γ2Γ22+ω2, and Ω is a shorthand for 2πf. Note that (11c) coincides with (10a) as expected. Equation (11b) follows from the Wiener–Khinchin theorem. From Eq. (11c), we can derive the following condition for the intrinsic spin noise to be dominated by the field-noise-driven spin noise:
γ2R2Γxx or yyΓ12Γ21Na.
(12)
Here, Na denotes the number of atoms in the collection.12 In the case of a cubic 87-rubidium cell of side length 5 mm and temperature 80 °C and assuming that the pumping and the relaxation parameters are given as in Sec. III, this will hold for Γxx or yy ≥ 10−5 nT2 ms, which corresponds to δBα100fT/Hz, where δBα2Γαα. Note that by applying P = 2SzR1 and δBα2=Γαα to Eq. (11c) and then by taking square roots, we arrive at the expression σexγPδBx2+δBy2/(4Γ2), which was previously noted in Sec. I. Such a slight change in notation explicitly reveals the role of the degree P of the spin polarization.

In Eq. (11c), the contribution of Γyy to the total power is larger than that of Γxx. Our result agrees with physical intuition; for example, an x-directional noise tilts the Larmor precession axis slightly along the x axis so that the precessions projected on the xy-plane will be somewhat compressed along the x axis, leading to smaller Sx magnitude compared to the case of a y-directional noise input. Moreover, as the Larmor frequency ω grows larger, discrimination between the x and y directions will be gradually blurred so that the discrepancy becomes negligible for ω ≫Γ2.

Note that we have an additional “high-pass” term Ω2Γyy in the numerator of Eq. (11b). This term dominates the Γxx term when ω ≪Ω. Therefore, y-directional noise is somewhat more “high-passed” compared to x-directional noise. This agrees with the intuitive predictions presented in Sec. I. In case ω ≪Ω, we imagine S as being along the z-direction, wherein some transitory displacements δS appear due to a “blurring” of the precession axis [Fig. 1(a)]. On the other hand, in case ω ≫Ω, the precession axis B+δB can now be imagined to be permanently tilted [Fig. 1(b)]. In this case, there is enough time for the precession to evolve so that x-directional magnetic noise can now contribute to the PSD of Sx, as long as Γ2 is not too large to wipe out the precession. (Such a wipe-out could happen when the B-field is sufficiently small, in the order of a few tens of nanoteslas in the case of 87Rb.) Therefore, in this case, the relative strength between Γ2 and ω will be an issue; indeed, it is the case as can be seen in Eq. (11b).

In addition, we note that there is no first-order contribution of Γzz to the PSD. This is clear since without any x or y-directional perturbation to the bias field, the spin will be fixed along the z-direction. However, if some transversal noise enters into the system, small precessions will occur, therefore enabling the z-directional noise to affect the whole picture. However, this will be a second-order effect involving the products ΓxxΓzz or ΓyyΓzz.

We note that a directional transversal noise δB(t)=δB(t)cosθx̂+sinθŷ with ⟨δB(t)δB(t′)⟩ = Γδ(tt′) results in the following noise power, which is clearly non-negative:
Sx2=γ2R2ΓΓ22+ω2Γ22cos(2θ)Γ2ωsin(2θ)16Γ12Γ2ω2+Γ22.
(13)
From the equation, it is clear that the noise-induced signal power is minimized at θ=π/4argω+iΓ2/2.
Finally, we note a non-Lorentzian behavior of the first-order PSD given in Eq. (11b). First, note that the two-sided PSD, i.e., PSDSx(f)/2, equals the following expression for all Ω:
γ2R216Γ12ω2+Γ222ω2Γxx+2Γ22Γyyω(ΓxxΓyy)Ω(Ωω)2+Γ22+2ω2Γxx+2Γ22Γyy+ω(ΓxxΓyy)Ω(Ω+ω)2+Γ22,
(14)
where we assumed Γxy = 0 for the sake of clarity. Each fraction in the parentheses has a Ω-proportional term in its numerator, deviating it from a Lorentzian function. This additional term, at Ω ≃±ω, is always comparable to the Lorentzian term in the case of an x-directional noise; in the case of a y-directional noise, it may even dominate the Lorentzian term, provided that Γ2ω, leading to a non-Lorentzian behavior of the PSD. However, if we focus on the near-resonance region so that we could replace Ω’s in the numerators by ω, the equation reduces to a sum of Lorentzians as
γ2R216Γ12Γxx+Γyy(Ωω)2+Γ22+3ΓxxΓyy(Ω+ω)2+Γ22
(15)
in cases where Γ2ω and Γxx ≃Γyy.
Now, we point out some alterations to the theory required for the second-order analysis. For this purpose, we will assume that third-order moments of the external noise vanish, and fourth-order moments satisfy the following equation; these assumptions were motivated by Wick’s theorem,
δBx(t1)δBx(t2)δBx(t3)δBx(t4)=Γxx2δ(t1t2)δ(t3t4)+Γxx2δ(t1t3)δ(t2t4)+Γxx2δ(t1t4)δ(t2t3).
(16)
Although the above equation is written for δBx, we assume the analogous equations for any fourfold products of the noise terms. Such assumptions are required to compute some higher-order terms mentioned in Eq. (4). In particular, our second-order analysis will require a quadruple integral term in Eq. (4),
t0t0+τt0t1t0t2t0t3Ã(t1)Ã(t2)Ã(t3)Ã(t4)dt4dt3dt2dt1,
(17)
whereas its triple-integral counterpart vanishes.
First, we point out that
RF(t1,t2,t3,t4)δ(t1t2)δ(t3t4)d4t=14t0t0+τt0t1F(t1,t1,t3,t3)dt3dt1,
(18)
where R denotes the region defined by t0t4 ≤ ⋯ ≤ t1t0 + τ, while similar expressions with δ(t1t3)δ(t2t4) or δ(t1t4)δ(t2t3) terms vanish; this follows by approximating δ by a sequence of smooth functions and estimating the volume of the support of the integrand. (Number 4 in the right-hand side will be replaced by powers of two for the analogous higher-order equations.) As a result, the first non-zero higher-order term in Eq. (4) is simplified as the following expression:
γ44t0t0+τt0t1et1KDe(t1t3)KDet3Kdt3dt1,
(19)
where K=A1c00 and D reduces, for the case in which Γxy = Γyz = Γzx = 0, to the 4 × 4 diagonal matrix with its consecutive diagonal entries given by −Γyy − Γzz, −Γzz − Γxx, −Γxx − Γyy, and 0 (the general form is given in  Appendix C). A similar inspection shows that Eq. (8b) is already correct up to second order. After this, the analysis proceeds verbatim; however, one should check that Eq. (5) is satisfied up to second order whenever Γ1τ, Γ2τ ≫ 1. This can be verified by utilizing formulas (10d) and (C2).
We can observe that the second-order coefficients for the PSD corresponding to ΓxxΓzz, ΓyyΓzz, and ΓxxΓyy are indeed non-zero (details in  Appendix C), whereas the coefficient corresponding to Γzz2 is zero. Moreover, the coefficients for Γxx2 and Γyy2 are both negative at resonance. For completeness, we inform the reader that the correction terms corresponding to Γxx2 and Γyy2 are given as ω2Γxx2+(Γ22+Ω2)Γyy2F, where F is given as follows:
F=γ4R2Ω4Ω22ω2Γ1Γ22Γ22+ω2+Γ22ω2+Γ1Γ2+Γ22/2Γ13(Ωω)2+Γ222×(Ω+ω)2+Γ222.
(20)

Since we gave the second-order terms corresponding to Γxx2 and Γyy2 and the first-order terms corresponding to Γxx and Γyy, we can now quantify the smallness of noise required for the first-order approximation to be valid: it is when the second-order terms are much smaller than the corresponding first-order terms for all Ω in the range of interest.

Now, we investigate on the higher-order terms of the spin autocorrelation function in the presence of Gaussian white noises, albeit the corresponding assumptions for the higher-order moments will be increasingly difficult to verify, in practice. Although it is computationally difficult to obtain the exact analytic formulas for higher-order terms, we can proceed to obtain the eigenvalues of an ODE representation for the spin autocorrelation function, which was previously approximated by a Dyson series in Eq. (4). This enables one to extract various effective relaxations and frequency shift for the spin variable.

The time evolution of the autocorrelation function can be summarized by the ODE,
ddtS(t)1S(t0)T1=A1c00+γ22D×S(t)1S(t0)T1,
(21)
although we still need to calculate the steady state values as before [as in Eq. (10)]. Such an effective renormalization of the Bloch operator by a noise term D was previously noted in Ref. 16, wherein the effects of the longitudinal noise component were investigated. Here, we may extend the results in Ref. 16 by computing the eigenvalues of the renormalized Bloch operator. In the case of Γyz = Γzx = 0, we have the following eigenvalues:
Γ1γ22(Γxx+Γyy)
and
Γ2γ24(Γxx+Γyy+2Γzz)±iωd,
(22)
where
ωd=ω2γ416ΓxxΓyy2+4Γxy2.
(23)
The first (or latter) formula in Eq. (22) describes the longitudinal (or transversal) relaxation. They show that the longitudinal relaxation rate is broadened by Γxx and Γyy, whereas the transversal relaxation rate is broadened by any Γαα. For the unpolarized spin case, a similar formula was derived in Ref. 16.

Moreover, the transversal noise coefficients lead to a frequency shift toward the DC (0 Hz) direction, which is not obvious at first glance as B<B+δB(t) in the case of a transversal noise. This becomes plausible when we recall that the transversal noises change the B-field direction so that the trajectory of the spin vector wobbles (in the longitudinal direction) during the precession. If this trajectory is projected onto the transversal plane, the resulting precession angle may observed to be smaller. This intuition is captured by an example in Fig. 3, which describes a repeated rotation of S(t) with respect to two different precession axes Btot=B+δB(t).

FIG. 3.

The initial spin vector S0=(1/2,0,0) is consecutively rotated around two different B-field axes Btot=B+δB(t). The rotations are assumed to be counterclockwise. Here, B=(0,0,10), and δB(t) suddenly changes from (−10, 0, 0) to (10, 0, 0) (arb. unit). The first rotation (red solid) points upward, whereas the second rotation (black solid) points downward. The blue dashed curve corresponds to the xy-planar circle of radius 1/2 centered at the origin. The dashed-dotted segments represent the axes of rotations. Blue dots correspond to (±1/2, 0, 0). Each rotation corresponds to a rotation angle of 90°, corrected for the instantaneous increase of Larmor frequency. As the blue dots indicate, the over precession angle is less than 180°. That is, the longitudinal wobbling of the spin vector trajectory rendered the overall precession angle to appear smaller.

FIG. 3.

The initial spin vector S0=(1/2,0,0) is consecutively rotated around two different B-field axes Btot=B+δB(t). The rotations are assumed to be counterclockwise. Here, B=(0,0,10), and δB(t) suddenly changes from (−10, 0, 0) to (10, 0, 0) (arb. unit). The first rotation (red solid) points upward, whereas the second rotation (black solid) points downward. The blue dashed curve corresponds to the xy-planar circle of radius 1/2 centered at the origin. The dashed-dotted segments represent the axes of rotations. Blue dots correspond to (±1/2, 0, 0). Each rotation corresponds to a rotation angle of 90°, corrected for the instantaneous increase of Larmor frequency. As the blue dots indicate, the over precession angle is less than 180°. That is, the longitudinal wobbling of the spin vector trajectory rendered the overall precession angle to appear smaller.

Close modal

The intuition presented as in Fig. 3 is implicitly relying on the assumption that δB(t) has a definite axis (here, the x axis). As a result, one might suspect that such an anisotropy of noise might be important for the predicted Larmor frequency shift. Indeed, this is the case as indicated by the difference term (ΓxxΓyy)2 appearing in Eq. (23). Moreover, as can be seen from Eq. (23), if the transversally anisotropic noise is very large, ωd will be an imaginary number, and in this case, the autocorrelation function is overdamped; there will be no oscillations in the autocorrelation function. The longitudinal–transversal correlations Γzx and Γyz do not seem to admit straightforward effective descriptions (however, see  Appendix C).

We note that although the relaxation rates of the autocorrelation function (effective relaxation rates) are altered, the spin relaxation rates themselves are not altered by B-field noises; after all, once the B-field noise δB(t) is given, the cross product term γS×δB(t) is always orthogonal to S; it does not lead to an additional relaxation of the squared norm SS.

It is plausible that the above-described methodology can be used to extract noise response information from other variants of the Bloch equation, which contain a multiplicative Bloch term γS×B and a damping term ΓreS. We present an example of such equations, which is also directly relevant to real world applications. From the representative equation [Eq. (1)], extensions to various atomic sensor dynamics will depend on the precise physical setup and signal processing details; as an example, we present an analysis for the case of a Bell–Bloom magnetometer.

Recall that as the Larmor frequency ω = γB is directly proportional to B, a collection of atomic spins serves as a basis for an optically pumped atomic magnetometer (OPAM). In particular, we consider a Bell–Bloom magnetometer,3 which is described by the following Bloch equation:
dSdt=γS×B+δBR(t)Sẑ2ΓreS,
(24)
where R(t) = R0 + R1 cos(ωt), B=[0,By,0]T, δB=[0,δBy(t),0]T, and Γre=diag[Γ2,Γ1,Γ2]. Here, ω = γBy, and henceforth, we will denote Γi+R0 as Γi.
If one defines a transversal variable S+ = Sz + iSx and performs a coordinate change S+ = eiωtK+ for some K+ = K1(t) + iK2(t), the transversal components of Eq. (24) become
dK+dt=iγδBy(t)K+Γ2K++R14R12eiωt+eiωtK++R0eiωt2+R1ei2ωt4.
(25)
Assuming a rotating wave approximation, we have
dK1dt=Γ2K1+γδByK2+R14,
(26a)
dK2dt=Γ2K2γδBK1,
(26b)
which has the same structure compared to Eq. (2). Assuming Γ2 > 0 and Gaussian white noise, we can analyze Eq. (26) in a parallel fashion to Eq. (2). For example, the D matrix analogous to that appearing in Eq. (19) will be given as D=diag[Γ,Γ,0], where Γ = Γyy is defined similarly as in the previous analyses. The exactly same enlargement scheme as in Eq. (7) is employed. The final results (including all orders) are given as
PSDK2(f)=γ2R12Γ2Γ2Γ2+γ2Γ2Γ2+γ2Γ2+4Ω2
(27)
and
PSDK1(f)=γ2ΓPSDK2(f)2Γ2+γ2Γ+πR12δ(Ω)2Γ2+γ2Γ2,
(28)
where Ω ≔ 2πf and δ(Ω) denotes the delta function. From the 2Γ2+γ2Γ2+4Ω2 term in the denominator, it is clear that there is a broadening Γ2Γ2+γ22Γ. Aside from this broadening, we can see a curious non-linearity of the form Γ/Γ2+γ2Γ, which becomes significant when γ2Γ is comparable to Γ2. Apparently, the damping Γ2 hinders the impact of this non-linearity.
In practice, the signal Sx cos(ωt), which is equal to
12+e2iωt+e2iωt4K2(t)+ie2iωte2iωt4K1(t),
(29)
is fed into a low-pass filter (LPF), the output of which can be considered as the final magnetometer output. Since the bandwidth of K2 is ∼Γ2 by Eq. (27) and the same estimate holds for K1, the PSD of the magnetometer output is approximately given by the PSD of K2 divided by 4 (assuming Γ2ω and an appropriate choice of the LPF).
In this section, we numerically deal with cases in which the external noise is given as Gaussian white noise along the x direction or the y direction. We employ the Euler–Maruyama method, a direct generalization of the Euler method to the stochastic case, to numerically solve the corresponding stochastic differential equations. We refer to Ref. 26 for a less technical review of the Euler–Maruyama method, Sec. 10.2 of Ref. 27 for a convergence proof of the method, and Refs. 28 and 29 for an introduction to stochastic differential equations. In case δB=δBn̂ with δB(t)δB(t)=Γδ(tt), the Bloch equation [Eq. (1)] can be rewritten in the following form, which is convenient for computer simulations:
dS=Γ2ω0ωΓ2000Γ1S+R2ẑdt+γΓ0nznynz0nxnynx0SdW,
(30)
where W denotes the standard one-dimensional Wiener process (Brownian motion).

Numerical analyses were performed with R = 8 ms−1, Γ1 = 25/3 ms−1, Γ2 = 10 ms−1, and γ = 2π × 0.007 rad/(ms nT) (this γ value corresponds to the case of 87Rb). The z-directional magnetic field was set to 100 nT [Fig. 4(a)] or 1000 nT [Fig. 4(b)], and the noises were characterized by Γαα = 1 nT2 ms, where α = x or y; this corresponds to ∼32 pT/Hz. The stochastic Bloch equation was integrated for a time period of 12 ms, starting from its (zeroth-order) macroscopic equilibrium R2Γ1ẑ. The last 6 ms of the approximate solutions were used for further investigations, such as PSD calculations by Welch’s method. The number of forward time steps was at least 217 for each sample trajectory corresponding to a discretized Brownian motion (details are provided in  Appendix B).

FIG. 4.

(a) Theory vs Monte Carlo where B = 100 nT. Black plots correspond to the x-directional noise input, whereas blue plots correspond to the y-directional noise input. Monte Carlo simulations were repeated 100 times, and the averaged results are shown. The red dashed line corresponds to the Larmor frequency. The theory and the Monte Carlo simulations gave almost the same results. (b) This is for the case in which B = 1000 nT, with all other things being equal. We note that it tended to be much more difficult to deal with cases with stronger B-fields; the Larmor oscillations became quite rapid so that shorter time step sizes seemed to be required.

FIG. 4.

(a) Theory vs Monte Carlo where B = 100 nT. Black plots correspond to the x-directional noise input, whereas blue plots correspond to the y-directional noise input. Monte Carlo simulations were repeated 100 times, and the averaged results are shown. The red dashed line corresponds to the Larmor frequency. The theory and the Monte Carlo simulations gave almost the same results. (b) This is for the case in which B = 1000 nT, with all other things being equal. We note that it tended to be much more difficult to deal with cases with stronger B-fields; the Larmor oscillations became quite rapid so that shorter time step sizes seemed to be required.

Close modal

The results (Fig. 4) indicate that our analytic results are indeed correct. Note the relatively “high-passing” nature of the y axis compared to the x axis; this provides evidence to our intuition-based predictions.

Moreover, we plot (Fig. 5) the following ratio as a function of B both in silico and in theory:
PSD at resonance when noise is along y axisPSD at resonance when noise is along x axis.
(31)
Figure 5 shows that the theory and the numerical experiment are consistent.
FIG. 5.

Theory vs Monte Carlo for the ratio described in Eq. (31). Each point comes from an average of 50 independent Monte Carlo experiments.

FIG. 5.

Theory vs Monte Carlo for the ratio described in Eq. (31). Each point comes from an average of 50 independent Monte Carlo experiments.

Close modal

As a further illustration, we investigated the situation in which the noise direction is given by the unit vector n̂=cosθx̂+sinθŷ by numerically verifying the proposed equation [Eq. (13)]. Now, the longitudinal B-field is set to B = 227.36 nT (so that ω = Γ2), whereas the noise is characterized by Γ = 1 nT2 ⋅ms. The Monte Carlo results were obtained by directly calculating the time-average value of Sx2. The results are given in Fig. 6, which show agreement between the theory and the experiment.

FIG. 6.

Theory vs Monte Carlo for noise power described in Eq. (13). Here, the B-field is set to make ω = Γ2. Each point comes from an average of 50 independent Monte Carlo experiments. The red dashed-dotted line indicates for the theoretical minimum.

FIG. 6.

Theory vs Monte Carlo for noise power described in Eq. (13). Here, the B-field is set to make ω = Γ2. Each point comes from an average of 50 independent Monte Carlo experiments. The red dashed-dotted line indicates for the theoretical minimum.

Close modal

Finally, the second-order results are supported by numerical calculations. With the same relaxation and pumping parameters as before and with B = 500 nT, we plot the PSD value at resonance as a function of Γxx or Γyy (Fig. 7). For these calculations, we employed the first-order-corrected steady state value given in Eq. (10e) as the initial condition. Note that the result shows that the second-order analysis is indeed more accurate. We finally note that because both second-order coefficients at resonance are always negative (for non-zero R), the higher-order terms should contribute positively to the PSD at resonance to prevent the value from being negative.

FIG. 7.

Theory (up to first and second order) vs Monte Carlo for PSD values at resonance as a function of Γαα. Here, B = 500 nT (a) X-directional noises; (b) y-directional noises. Each point comes from an average of 300 independent Monte Carlo experiments.

FIG. 7.

Theory (up to first and second order) vs Monte Carlo for PSD values at resonance as a function of Γαα. Here, B = 500 nT (a) X-directional noises; (b) y-directional noises. Each point comes from an average of 300 independent Monte Carlo experiments.

Close modal

We have solved, by a time-dependent perturbative approach, the optically pumped stochastic Bloch equation describing a macroscopic spin collection subjected to small magnetic field white noises in the presence of an external bias magnetic field and a σ+ pumping beam both along the z direction. We have focused on the regime where the magnetic field-noise-driven spin precession dominates the effect of the intrinsic spin noise so that the latter can be considered negligible. The analytic forms for the autocorrelation function and the corresponding power spectral density (PSD) of the macroscopic spin along the probing axis (x axis in our case) are given up to first order. The analytic formula of the first-order PSD is of non-Lorentzian nature. Moreover, it reveals anisotropy in the effects of the magnetic field noises; the effect of x-directional magnetic field noises dominates at the low frequency range of the PSD (given ω ≫ Γ2), whereas the y-directional magnetic field noises prevail at high frequencies Ω ≫ ω. The z-directional magnetic field noises, represented by Γzz, Γzx, and Γyz, contribute to the PSD via second-order corrections. A closer inspection shows that anisotropic transversal B-field noises will lead to an effective Larmor frequency shift toward the DC direction when one observes the PSD; this is somewhat surprising because the magnetic field magnitude only increases in the case of transversal B-field noises. Our analytic results are supported by Monte Carlo simulations approximately solving the stochastic Bloch equation. Our analysis will be useful for understanding the behavior of a spin collection in a noisy magnetic field, which is a situation naturally arising in studies of atomic magnetometers; to emphasize this point, we indeed applied the analytic methodology to the case of a Bell–Bloom magnetometer. The results for the magnetometer indicate a noise-induced broadening effect and a nonlinear dependency of the magnetometer output spectrum [Eq. (27)] with respect to Γ.

This work was supported by the Agency for Defense Development Grant funded by the Korean Government (No. 915001102). The authors thank Ji Hoon Yoon for having a fruitful discussion with the authors about the topic.

The authors have no conflicts to disclose.

Minwoo M. Kim: Formal analysis (lead); Methodology (lead); Writing – original draft (lead). Sangkyung Lee: Supervision (lead); Writing – original draft (supporting).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

For Fig. 1(c), the simulation was again done using the Euler–Maruyama method. B was set to 500 nT; other parameters were the same as in Sec. III. The input noises were simulated by applying high-pass and low-pass filters to a prescribed transversal Gaussian white noise. These filtering processes were done by MATLAB R2022a highpass and lowpass functions with their cutoff frequencies being 10 times the Larmor frequency and 0.1 times the Larmor frequency, respectively.

The white noise δB=δBn̂, where n̂ẑ, satisfied ⟨δB(t)δB(t′)⟩ = (1 nT2 ms) δ(tt′). Indeed, noise-parallel spin magnitudes were smaller than the noise-perpendicular spin magnitudes when the input noise consisted of high-frequency components.

The time step sizes for the numerical simulations were determined by a relative error criterion. First, the whole time interval [0,12] (ms) was divided into 225 subintervals of equal lengths and a discretized Brownian motion was generated accordingly. Then, we executed the Euler–Maruyama method with time step size Δt = 12/216. With such a time step size, each Brownian increment ΔW for a single time step of the method was obtained by summing 225−16 = 512 consecutive entries of the discretized Brownian motion. Then, we reduced the time step size by a factor of 1/2, and the whole Euler–Maruyama method was repeated; such half-reductions of the time step were performed until the relative errors between the previous calculation and the current calculation were always less than 0.01 at time instances t = 12 − i × 12/26 for i = 0, 1, …, 29. This condition was checked separately for all three components of the spin vector (so that 30 × 3 = 90 relative error conditions were checked in total). If the condition failed until Δt reached 12/225, the corresponding discretized Brownian motion was discarded since it led to failure. The above-described procedure was also used for determining the time step size for Fig. 1(c); after the time step size was determined, the corresponding ΔW sequence was calculated and then fed into the high-pass and low-pass filters described in  Appendix A.

C and D, appeared in Eq. (5), are given as
C=γ2eΓ1τeΓ2τcos(ωτ)ωΓyzΓ1Γzx+Γ2Γzx+eΓ2τΓ1ΓyzΓ2Γyz+ωΓzxsin(ωτ)2ω2+Γ1Γ22
(C1)
and
D=γ2R4Γ1ω2+Γ22ω2+Γ1Γ22ωΓyz+Γ2Γzxω2+Γ1Γ22eΓ1τω2+Γ22ωΓyzΓ1Γzx+Γ2ΓzxeΓ2τΓ1Γ1Γ2ΓzxΓ22Γzx+Γ1Γyzω2Γ2ωΓyz+ω2Γzxcos(ωτ)eΓ2τΓ1Γ1Γ2ΓyzΓ22ΓyzΓ1Γzxω+2Γ2ωΓzx+ω2Γyzsin(ωτ),
(C2)
respectively. Somewhat interestingly, Dγ2RωΓyz+Γ2Γzx4Γ1ω2+Γ22 for Γ1τ, Γ2τ ≫ 1, and this expression equals the right-hand side of Eq. (10d). In addition, note that Eqs. (C1) and (C2) exhibit some longitudinal–transversal coupling terms involving Γzx or Γyz, which are not multiplied by a sinusoidal term. Therefore, if these coupling terms are Fourier transformed to give the PSD of Sx, this will result in a DC noise component centered at 0 Hz, which is somewhat qualitatively similar to the double peak phenomenon noted in Ref. 30. Indeed, if one performs the relevant second-order calculations, the following second-order DC term occurs in the autocorrelation function of Sx:
γ4R216Γ12Γ22+ω22Γ2Γzx+ωΓyz2.
(C3)
Γ̃, appeared in Eq. (8b), is given as
Γ̃=γ2A000B0000,
(C4)
where
A=bcd2f2d2cef2ebd2a2b2ce2cd2fef2df2e2a2bc2efd2bdf2d2cea2b2cd2f2e2fdcd2fd2ace2fe2aefd2f2e2ab2c2dae2f2ebe2fda2bc2d2f2ef2dfe2ad2ab2c2e2bdf2dae2f2e2ab
(C5)
and
B=b2c2d2f2d2a2c2e2f2e2a2b2,
(C6)
wherein we defined a ≔ Γxx, b ≔ Γyy, c ≔ Γzz, d ≔ Γxy, e ≔ Γyz, and f ≔ Γzx to save space. In addition, in Eq. (19),
D=ΓyyΓzzΓxyΓzx0ΓxyΓzzΓxxΓyz0ΓzxΓyzΓxxΓyy00000.
(C7)
Finally, we provide the second-order corrections for the PSD of Sx corresponding to the products ΓxxΓzz, ΓyyΓzz, and ΓxxΓyy. They are given as
γ4R2C1ΓxxΓzz+C2ΓyyΓzz8Γ12Γ2ω2+Γ22(Ωω)2+Γ222(Ω+ω)2+Γ222
(C8)
and
γ4R2C3ΓxxΓyy2Γ13(Ωω)2+Γ222(Ω+ω)2+Γ222,
(C9)
where Ci are given by the following equations:
C1=ω2+2Γ22Ω6ω4+3Γ22ω26Γ24Ω4ω6+9Γ24ω26Γ26Ω2+ω87Γ22ω615Γ24ω45Γ26ω2+2Γ28,
(C10a)
C2=ω2Ω6ω43Γ22ω2+4Γ24Ω4ω6+22Γ22ω4+13Γ24ω2+8Γ26Ω2+ω8+11Γ22ω6+15Γ24ω4+Γ26ω24Γ28,
(C10b)
C3=Ω6+3Γ22ω2Ω4ω42Γ2(2Γ1+Γ2)ω23Γ24Ω2+ω2+Γ223.
(C10c)
1
S.
Lee
,
S. H.
Yim
,
T. H.
Kim
,
Z.
Kim
, and
K.
Shim
, “
Lock-in-detection in 87Rb–129Xe/131Xe atom spin gyroscopes
,”
J. Phys. B: At. Mol. Opt. Phys.
53
,
035502
(
2020
).
2
T. G.
Walker
and
M. S.
Larsen
, “
Spin-exchange-pumped NMR gyros
,”
Adv. At., Mol., Opt. Phys.
65
,
373
401
(
2016
).
3
W. E.
Bell
and
A. L.
Bloom
, “
Optically driven spin precession
,”
Phys. Rev. Lett.
6
(
6
),
280
281
(
1961
).
4
T. M.
Tierney
,
N.
Holmes
,
S.
Mellor
,
J. D.
López
,
G.
Roberts
,
R. M.
Hill
,
E.
Boto
,
J.
Leggett
,
V.
Shah
,
M. J.
Brookes
,
R.
Bowtell
, and
G. R.
Barnes
, “
Optically pumped magnetometers: From quantum origins to multi-channel magnetoencephalography
,”
NeuroImage
199
,
598
608
(
2019
).
5
R.
Jiménez-Martínez
,
W. C.
Griffith
,
S.
Knappe
,
J.
Kitching
, and
M.
Prouty
, “
High-bandwidth optical magnetometer
,”
J. Opt. Soc. Am. B
29
,
3398
3403
(
2012
).
6
R. J.
Sewell
,
M.
Koschorreck
,
M.
Napolitano
,
B.
Dubost
,
N.
Behbood
, and
M. W.
Mitchell
, “
Magnetic sensitivity beyond the projection noise limit by spin squeezing
,”
Phys. Rev. Lett.
109
(
25
),
253605
(
2012
).
7
J. M.
Geremia
,
J. K.
Stockton
, and
H.
Mabuchi
, “
Suppression of spin projection noise in broadband atomic magnetometry
,”
Phys. Rev. Lett.
94
(
20
),
203002
(
2005
).
8
F.
Wolfgramm
,
A.
Cerè
,
F. A.
Beduini
,
A.
Predojević
,
M.
Koschorreck
, and
M. W.
Mitchell
, “
Squeezed-light optical magnetometry
,”
Phys. Rev. Lett.
105
,
053601
(
2010
).
9
C.
Troullinou
,
R.
Jiménez-Martínez
,
J.
Kong
,
V. G.
Lucivero
, and
M. W.
Mitchell
, “
Squeezed-light enhancement and backaction evasion in a high sensitivity optically pumped magnetometer
,”
Phys. Rev. Lett.
127
,
193601
(
2021
).
10
S. H.
Yim
,
Z.
Kim
,
S.
Lee
,
T. H.
Kim
, and
K. M.
Shim
, “
Note: Double-layered polyimide film heater with low magnetic field generation
,”
Rev. Sci. Instrum.
89
,
116102
(
2018
).
11
H.
Korth
,
K.
Strohbehn
,
F.
Tejada
,
A. G.
Andreou
,
J.
Kitching
,
S.
Knappe
,
S. J.
Lehtonen
,
S. M.
London
, and
M.
Kafel
, “
Miniature atomic scalar magnetometer for space based on the rubidium isotope 87Rb
,”
J. Geophys. Res.: Space Phys.
121
(
8
),
7870
7880
, (
2016
).
12
M. V.
Romalis
, “
Quantum noise in atomic magnetometers
,” in
Optical Magnetometry
(
Cambridge University Press
,
2013
).
13
F.
Bloch
, “
Nuclear induction
,”
Phys. Rev.
70
(
7-8
),
460
474
(
1946
).
14
S. J.
Seltzer
and
M. V.
Romalis
, “
Unshielded three-axis vector operation of a spin-exchange-relaxation-free atomic magnetometer
,”
Appl. Phys. Lett.
85
,
4804
4806
(
2004
).
15
S. J.
Seltzer
, “
Developments in alkali-metal atomic magnetometry
,” Ph. D. thesis (
Princeton University
,
2008
).
16
F.
Li
,
S. A.
Crooker
, and
N. A.
Sinitsyn
, “
Higher-order spin-noise spectroscopy of atomic spins in fluctuating external fields
,”
Phys. Rev. A
93
,
033814
(
2016
).
17
S.
Appelt
,
A. B. A.
Baranga
,
C. J.
Erickson
,
M. V.
Romalis
,
A. R.
Young
, and
W.
Happer
, “
Theory of spin-exchange optical pumping of 3He and 129Xe
,”
Phys. Rev. A
58
,
1412
(
1998
).
18
M. M.
Sharipova
,
A. N.
Kamenskii
,
I. I.
Ryzhov
,
M. Y.
Petrov
,
G. G.
Kozlov
,
A.
Greilich
,
M.
Bayer
, and
V. S.
Zapasskii
, “
Stimulated spin noise in an activated crystal
,”
J. Appl. Phys.
126
,
143901
(
2019
).
19
K.
Mouloudakis
,
G.
Vasilakis
,
V. G.
Lucivero
,
J.
Kong
,
I. K.
Kominis
, and
M. W.
Mitchell
, “
Effects of spin-exchange collisions on the fluctuation spectra of hot alkali-metal vapors
,”
Phys. Rev. A
106
,
023112
(
2022
).
20
V. G.
Lucivero
,
N. D.
McDonough
,
N.
Dural
, and
M. V.
Romalis
, “
Correlation function of spin noise due to atomic diffusion
,”
Phys. Rev. A
96
,
062702
(
2017
).
21
E. L.
Ivchenko
, “
Fluctuations of spin polarization of free carriers in semiconductors
,”
Fiz. Tekh. Poluprovodn.
7
,
1489
(
1973
) [Sov. Phys.-Semicond. 7, 998 (1973)].
22
N. A.
Sinitsyn
and
Y. V.
Pershin
, “
The theory of spin noise spectroscopy: A review
,”
Rep. Prog. Phys.
79
,
106501
(
2016
).
23
V.
Shah
,
G.
Vasilakis
, and
M. V.
Romalis
, “
High bandwidth atomic magnetometery with continuous quantum nondemolition measurements
,”
Phys. Rev. Lett.
104
,
013601
(
2010
).
24
G. D.
Cates
,
S. R.
Schaefer
, and
W.
Happer
, “
Relaxation of spins due to field inhomogeneities in gaseous samples at low magnetic fields and low pressures
,”
Phys. Rev. A
37
,
2877
(
1988
).
25
W. J.
Rugh
,
Linear System Theory
(
Pearson
,
1995
).
26
D. J.
Higham
, “
An algorithmic introduction to numerical simulation of stochastic differential equations
,”
SIAM Rev.
43
(
3
),
525
546
(
2001
).
27
P. E.
Kloeden
and
E.
Platen
, “
Numerical solution of stochastic differential equations
,” in
Stochastic Modelling and Applied Probability
(
Springer
,
1995
).
28
L. C.
Evans
,
An Introduction to Stochastic Differential Equations
(
American Mathematical Society
,
2014
).
29
N. G.
van Kampen
,
Stochastic Processes in Physics and Chemistry
(
Elsevier
,
2007
).
30
S. V.
Poltavtsev
,
I. I.
Ryzhov
,
M. M.
Glazov
,
G. G.
Kozlov
,
V. S.
Zapasskii
,
A. V.
Kavokin
,
P. G.
Lagoudakis
,
D. S.
Smirnov
, and
E. L.
Ivchenko
, “
Spin noise spectroscopy of a single quantum well microcavity
,”
Phys. Rev. B
89
,
081304(R)
(
2014
).