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.

## I. INTRODUCTION

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 films^{10} 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 $\sigma in=F/(2Na)$, where *N*_{a} 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 $\sigma ex\u2243\gamma P\Gamma xx+\Gamma yy/(4\Gamma 2)$, given that the Larmor frequency is much larger than the transversal spin relaxation rate Γ_{2}, where *P* = 2*S*_{z} 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.

^{2,13–16}

^{87}Rb vapor cell),

*γ*is the gyromagnetic ratio,

*R*is the optical pumping rate, $s\u20d7$ 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 $\delta B\u20d7$ is white noise, $B\u20d7$ is along the z-direction, and the pumping beam is also

*z*-directional with a

*σ*

^{+}polarization, i.e., $s\u20d7=z\u0302$. 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

*S*

_{z}=

*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

*S*

_{x}, arising from various noise sources, as a Lorentzian function or a sum of nearly Lorentzian functions; see Refs. 12 and 21–23.

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 *S*_{x}. 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\u20d7+\delta B\u20d7$ 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)].

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 *δB*_{z}(*t*) = 0 for all *t* and approximating *S*_{z} 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.

## II. THEORY

Summary: We consider the model given by Eq. (2). The autocorrelation function of *S*_{x} 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 *t*_{0} → ∞ 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 *S*_{x} 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.

### A. First-order analysis

*δB*

_{α}(

*t*) are mean-zero white noise. We note that the average spin $S\u20d7$ 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 $\Gamma 2\u2032$, $\Gamma 2\u2032$, and $\Gamma 1\u2032$ as its three consecutive diagonal entries, where $\Gamma 2\u2032$ and $\Gamma 1\u2032$ 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 $\Gamma 1\u2032$ may differ from $\Gamma 2\u2032$.

^{13}In this case, the Bloch equation can be written in the following matrix form:

_{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\u20d7+c\u20d7$ and by defining a coordinate change as $S\u0303(t)\u2254exp\u2212tA1c\u20d700S\u20d71$, 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):

**1**denotes the 4 × 4 identity matrix and h.o.t. means higher-order terms:

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 $\u27e8S\u0303(t0)S\u0303(t0)T\u27e9$ when the system reaches the steady state (denoted by *t*_{0} → ∞); 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.

*S*

_{x}:

*γB*as

*ω*. Here, we assumed that the noises are white with their covariances given by delta functions as

_{αβ}; 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 *S*_{x} is *a priori* valid only for short evolution times *τ* with $\tau A\u0303\u226a1$ (or more precisely, $\u222bt0t0+\tau dtA\u0303(t)\u226a1$). 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 ⟨*S*_{x}(*t*_{0} + *τ*)⟩⟨*S*_{x}(*t*_{0})⟩ for such large *τ* because physical intuition tells us that *S*_{x}(*t*_{0}) and *S*_{x}(*t*_{0} + *τ*) will be statistically uncorrelated; any deviation from the expectation at time *t*_{0} will be exponentially damped out in time, thereby not affecting the spin at time *t*_{0} + *τ*. 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\u0303$, which resulted in the damping exponential $e\u2212\Gamma 2\tau $ terms shown in the right-hand side of Eq. (5).

To sum up, Eq. (5) is valid for values of *τ* satisfying $\tau A\u0303\u226a1$ 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., $A\u0303/\Gamma 1\u226a1$ and $A\u0303/\Gamma 2\u226a1$) since in this case the region $\tau A\u0303\u226a1$ depicted in Fig. 2 will be shifted toward *τ* → ∞.

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.

*S*

_{x}(

*t*

_{0})

*S*

_{y}(

*t*

_{0})⟩, ⟨

*S*

_{x}(

*t*

_{0})

*S*

_{z}(

*t*

_{0})⟩, and ⟨

*S*

_{x}(

*t*

_{0})⟩ after all the irrelevant transients are damped out (i.e.,

*t*

_{0}→ ∞). This can be achieved by enlarging the main vector variable further to include the second-order terms $S(t)\u2254Sx(t)Sx(t)Sx(t)Sy(t)\cdots Sz(t)Sz(t)T$ as in the following equation:

**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\u0303(t)\u2254e\u2212tA1$$S(t)S\u20d7(t)1T$, we conclude that the second-order approximation for $\u27e8S\u0303(t)\u27e9$ 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

*t*

_{0},

*t*>

*t*

_{0}. It is a matrix that has linear combinations of

*γ*

^{2}Γ

_{αβ}as its components (the full expression is given in Appendix C).

_{αβ}(i.e., up to second order in the noise amplitudes), we have the following equation in the original coordinates:

Incidentally, every eigenvalue of the 12 × 12 submatrix of $A1+\Gamma \u0303$ at the top left corner has a negative real part when Γ_{αβ} = 0; this mathematically proves that the expectations $\u27e8S\u27e9$ and $\u27e8S\u20d7\u27e9$ will indeed converge to an equilibrium, at least for sufficiently small Γ_{αβ} values.

_{αβ}. We used the suggestive symbol ∞ in Eq. (10) to emphasize that the equation represents the steady state values,

*xy*-product given in Eq. (10b) vanishes when Γ

_{xx}= Γ

_{yy}and Γ

_{xy}= 0, as expected due to symmetry. In addition, since $\omega 2\Gamma xx\u22122\Gamma 2\omega \Gamma xy+\Gamma 22\Gamma yy\u22650$, Eq. (10a) gives non-negative values. Moreover, the steady state value for

*S*

_{z}has a first-order correction term, which cannot be easily justified by the Langevin approach; the Langevin approach, in order to be justified, requires

*S*

_{z}to be approximated by its

*a priori*equilibrium value $R/2\Gamma 1$. Finally, we note that (10c) and (10d) vanish whenever the correlation coefficients Γ

_{yz}and Γ

_{zx}are zero.

*S*

_{x}(

*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,

*f*denotes the argument of the PSD function, $C\u2254\gamma 2R2e\u2212\Gamma 2\tau 16\Gamma 12\Gamma 2\Gamma 22+\omega 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:

*N*

_{a}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}nT

^{2}ms, which corresponds to $\delta B\alpha \u2265100fT/Hz$, where $\delta B\alpha 2\u2254\Gamma \alpha \alpha $. Note that by applying

*P*= 2

*S*

_{z}≈

*R*/Γ

_{1}and $\delta B\alpha 2=\Gamma \alpha \alpha $ to Eq. (11c) and then by taking square roots, we arrive at the expression $\sigma ex\u2243\gamma P\delta Bx2+\delta By2/(4\Gamma 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 *S*_{x} 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\u20d7$ as being along the *z*-direction, wherein some transitory displacements $\delta S\u20d7$ appear due to a “blurring” of the precession axis [Fig. 1(a)]. On the other hand, in case *ω* ≫Ω, the precession axis $B\u20d7+\delta B\u20d7$ 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 *S*_{x}, 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 ^{87}Rb.) 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}.

*δB*(

*t*)

*δB*(

*t*′)⟩ = Γ

_{⊥}

*δ*(

*t*−

*t*′) results in the following noise power, which is clearly non-negative:

_{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

_{2}≪

*ω*and Γ

_{xx}≃Γ

_{yy}.

### B. Higher-order analysis

*δB*

_{x}, 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),

*t*

_{0}≤

*t*

_{4}≤ ⋯ ≤

*t*

_{1}≤

*t*

_{0}+

*τ*, while similar expressions with

*δ*(

*t*

_{1}−

*t*

_{3})

*δ*(

*t*

_{2}−

*t*

_{4}) or

*δ*(

*t*

_{1}−

*t*

_{4})

*δ*(

*t*

_{2}−

*t*

_{3}) 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:

_{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).

_{xx}Γ

_{zz}, Γ

_{yy}Γ

_{zz}, and Γ

_{xx}Γ

_{yy}are indeed non-zero (details in Appendix C), whereas the coefficient corresponding to $\Gamma zz2$ is zero. Moreover, the coefficients for $\Gamma xx2$ and $\Gamma yy2$ are both negative at resonance. For completeness, we inform the reader that the correction terms corresponding to $\Gamma xx2$ and $\Gamma yy2$ are given as $\u2212\omega 2\Gamma xx2+(\Gamma 22+\Omega 2)\Gamma yy2F$, where $F$ is given as follows:

Since we gave the second-order terms corresponding to $\Gamma xx2$ and $\Gamma 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.

_{yz}= Γ

_{zx}= 0, we have the following eigenvalues:

_{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\u20d7<B\u20d7+\delta B\u20d7(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\u20d7(t)$ with respect to two different precession axes $B\u20d7tot=B\u20d7+\delta B\u20d7(t)$.

The intuition presented as in Fig. 3 is implicitly relying on the assumption that $\delta B\u20d7(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 $(\Gamma xx\u2212\Gamma 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 $\gamma S\u20d7\xd7\delta B\u20d7(t)$ is always orthogonal to $S\u20d7$; it does not lead to an additional relaxation of the squared norm $S\u20d7\u22c5S\u20d7$.

### C. Implications on optically pumped atomic magnetometers

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 $\gamma S\u20d7\xd7B\u20d7$ and a damping term $\u2212\Gamma reS\u20d7$. 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.

*ω*=

*γ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:

*R*(

*t*) =

*R*

_{0}+

*R*

_{1}cos(

*ωt*), $B\u20d7=[0,By,0]T$, $\delta B\u20d7=[0,\delta By(t),0]T$, and $\Gamma re=diag[\Gamma 2\u2032,\Gamma 1\u2032,\Gamma 2\u2032]$. Here,

*ω*=

*γB*

_{y}, and henceforth, we will denote $\Gamma i\u2032+R0$ as Γ

_{i}.

*S*

_{+}=

*S*

_{z}+

*iS*

_{x}and performs a coordinate change

*S*

_{+}=

*e*

^{−iωt}

*K*

_{+}for some

*K*

_{+}=

*K*

_{1}(

*t*) +

*iK*

_{2}(

*t*), the transversal components of Eq. (24) become

_{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[\u2212\Gamma ,\u2212\Gamma ,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

*πf*and

*δ*(Ω) denotes the delta function. From the $2\Gamma 2+\gamma 2\Gamma 2+4\Omega 2$ term in the denominator, it is clear that there is a broadening $\Gamma 2\u2192\Gamma 2+\gamma 22\Gamma $. Aside from this broadening, we can see a curious non-linearity of the form $\Gamma /\Gamma 2+\gamma 2\Gamma $, which becomes significant when

*γ*

^{2}Γ is comparable to Γ

_{2}. Apparently, the damping Γ

_{2}hinders the impact of this non-linearity.

*S*

_{x}cos(

*ωt*), which is equal to

*K*

_{2}is ∼Γ

_{2}by Eq. (27) and the same estimate holds for

*K*

_{1}, the PSD of the magnetometer output is approximately given by the PSD of

*K*

_{2}divided by 4 (assuming Γ

_{2}≪

*ω*and an appropriate choice of the LPF).

## III. COMPARISON BETWEEN THE THEORY AND STOCHASTIC SIMULATIONS

*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 $\delta B\u20d7=\delta Bn\u0302$ with $\delta B(t)\delta B(t\u2032)=\Gamma \delta (t\u2212t\u2032)$, the Bloch equation [Eq. (1)] can be rewritten in the following form, which is convenient for computer simulations:

*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 ^{87}Rb). 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 nT^{2} 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\Gamma 1z\u0302$. 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 2^{17} for each sample trajectory corresponding to a discretized Brownian motion (details are provided in Appendix B).

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.

*B*both

*in silico*and in theory:

As a further illustration, we investigated the situation in which the noise direction is given by the unit vector $n\u0302=cos\u2061\theta x\u0302+sin\u2061\theta y\u0302$ 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 nT^{2} ⋅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.

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.

## IV. CONCLUSIONS

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 Γ.

## ACKNOWLEDGMENTS

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.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

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

## DATA AVAILABILITY

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

### APPENDIX A: SIMULATION PARAMETERS FOR Fig. 1(c)

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 $\delta B\u20d7=\delta Bn\u0302$, where $n\u0302\u22a5z\u0302$, satisfied ⟨*δB*(*t*)*δB*(*t*′)⟩ = (1 nT^{2} ms) *δ*(*t* − *t*′). Indeed, noise-parallel spin magnitudes were smaller than the noise-perpendicular spin magnitudes when the input noise consisted of high-frequency components.

### APPENDIX B: DETERMINATION OF TIME STEP SIZE FOR EULER–MARUYAMA METHOD

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 2^{25} 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/2^{16}. With such a time step size, each Brownian increment Δ*W* for a single time step of the method was obtained by summing 2^{25−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/2^{6} 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/2^{25}, 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.

### APPENDIX C: OMITTED FORMULAS

_{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

*S*

_{x}, 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

*S*

_{x}:

*a*≔ Γ

_{xx},

*b*≔ Γ

_{yy},

*c*≔ Γ

_{zz},

*d*≔ Γ

_{xy},

*e*≔ Γ

_{yz}, and

*f*≔ Γ

_{zx}to save space. In addition, in Eq. (19),

*S*

_{x}corresponding to the products Γ

_{xx}Γ

_{zz}, Γ

_{yy}Γ

_{zz}, and Γ

_{xx}Γ

_{yy}. They are given as

## REFERENCES

^{87}Rb–

^{129}Xe/

^{131}Xe atom spin gyroscopes

^{87}Rb

*Optical Magnetometry*

^{3}He and

^{129}Xe

*Stochastic Modelling and Applied Probability*

*An Introduction to Stochastic Differential Equations*