Spectral analysis of spin noise in an optically spin-polarized stochastic Bloch equation driven by noisy magnetic fields

We provide a closed-form autocorrelation function and power spectral density (PSD) of the solution, along a prescribed probing direction, to a noisy version of an optically pumped Bloch equation wherein each component of the external magnetic field is subject to (possibly correlated) white noise. We conclude that, up to first order in the white noise covariance amplitudes, noise in the bias B-field direction does not affect the autocorrelation function. 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, and an effective Larmor frequency shift caused by anisotropic transversal B-field noises, towards the DC direction, is revealed. The analytic results are supported by Monte Carlo simulations employing the Euler-Maruyama method.


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][4][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 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/(2N a ) 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 σ ex ≃ |γP | δB 2 x + δB 2 y /(4 √ Γ 2 ) given that the Larmor frequency is much larger than the transversal spin relaxation rate Γ 2 , where P = 2S z is the degree of polarization, γ is the gyromagnetic ratio, and δB x and δB y are the magnetic field noise amplitudes along the x and y directions, respectively.In this paper, we focus on the regime where σ ex ≫ σ in .The precise mathematical formulation for such a regime, including the interpretation for δB x and δB y , will become clearer later after we compare the noise powers for these two types of noises in Section 2.
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 insight into atomic sensors.In this regard, we investigate the response of a magnetic resonance system subject 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][14][15][16], where ⃗ S is the average spin (e.g., the average on a hot 87 Rb 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.
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 [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 [19], which incorporated spontaneous fluctuations due to, e.g., collisions between atoms, or [20], which focused on noise due to atomic diffusion.
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: [12,[21][22][23].Moreover, 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); we can intuitively speculate that x-directional noises and y-directional noises lead to different PSDs of S x , because we are observing the x-direction.
Moreover, 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. 1a).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. 1b).
We note that such a frequency-dependence can be confirmed by computer simulations.Fig. 1c 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 the B-field noise covariance amplitudes (these amplitudes will be denoted as Γ αβ ).Moreover, we show that the PSD of S x cannot be written as a sum of nearly Lorentzian functions.Finally, we give some second-order corrections to the PSD.We believe that such higher order information cannot be obtained if we approximate Eq. ( 1) by a Langevin equation: assuming δB z (t) = 0 for all t and approximating S z to be a constant transforms the equation into a Langevin form.Such higher order terms exhibit some qualitatively different features from their first-order counterparts: the higher order lineshapes are much more complicated and they exhibit longitudinaltransversal coupling which is absent in the firstorder terms.

Theory
We consider a situation in which ⃗ B = Bẑ + δB x (t) δB y (t) δB z (t) T , where δB α (t) are meanzero noise terms with an approximately white power spectrum.We note that the average spin ⃗ S of the collection experiences a spatially homogeneous magnetic field fluctuation as we are implicitly assuming the size of the atomic vapor cell is much smaller than the spatial scale of the 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  In this case the precession axis is somewhat fixed; therefore, the spin vector can precess around a well-defined axis for a while.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 ⟨S ∥ S ⊥ ⟩ appears to be non-zero, as can be verified from Eq. (9b).
(x, y) isotropic, although it is not fully isotropic as Γ ′ 1 may differ from Γ ′ 2 [13].We moreover assume ⃗ s = ẑ, i.e., a σ + pumping beam.In this case, the Bloch equation can be written in the following matrix form.
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 righthand side of Eq. ( 2) as (A 1 + A 2 (t)) ⃗ S + ⃗ c, and by defining a coordinate change as S(t , the equation is transformed into the following homogeneous form. Hereafter, we denote the right-hand-side of the above equation as A(t) S(t).
We have the following second-order approximation for the autocorrelation function of S, where 1 denotes the 4-by-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 timeevolution operator in the parenthesis in Eq. ( 4), which is achieved in Eq. (5).The second step is to calculate the initial value ⟨ S(t 0 ) S(t 0 ) T ⟩ when the system reaches steady state (i.e., t 0 → ∞); the result of this step is summarized in Eq. ( 9).Also, 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 S x : 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 ′ )⟩ = Γ αβ δ(t−t ′ ) for symmetric coefficients Γ αβ ; that is, the correlation time of the B-field noises are much shorter than any relevant time scales of the system.C and D are complicated correlation 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 τ A ≪ 1 (or more precisely, t0+τ t0 dt A(t) ≪ 1).Meanwhile, we argue that equation ( 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 Γ αβ .Also, 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 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 τ A ≪ 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, since in this case the region τ A ≪ 1 depicted in Fig. 2 will be shifted towards τ → ∞.
In any case, equation ( 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.5) is valid for all τ ≥ 0 for sufficiently small noises.
So far, we did not assume that the system is in its steady state.Now, we determine the relevant steady state values ⟨S x (t 0 ) 2 ⟩, ⟨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 T as in the following equation.
Here, 1 denotes the 3-by-3 identity matrix and ⊗ denotes the Kronecker product.Denoting the terms in the parenthesis as A 1 + A 2 , and defining a coordinate change S(t) := e −tA1 S(t) ⃗ S(t) 1 T , we conclude that the second-order approximation for ⟨ S(t)⟩ similar to equation ( 4) coincides with the first-order approximate solution to the following differential equation (7).Although such an equivalence is a priori valid only for sufficiently short time evolutions, this limitation is superficial because ( 7) is independent of the initial time t 0 .
Here, Γ is a 13-by-13 constant matrix that is independent of t > t 0 .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.
Incidentally, every eigenvalue of the 12-by-12 submatrix of A 1 + Γ 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-by-12 linear equation; some are provided in the following equations, written up to first order in Γ αβ .We denoted the value of the time variable as ∞ to emphasize that Eq. ( 9) represents the steady state values.
Note that the xy-product given in Eq. (9b) vanishes when Γ xx = Γ yy and Γ xy = 0, as expected due to symmetry.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Γ 1 .Finally, we note that (9c) and (9d) 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 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 nonzero Γ zz , Γ yz , or Γ zx .That is, the system under consideration is insusceptible to longitudinal noises.
In Eq. (10), and Ω := 2πf .Note that (10c) coincides with (9a) as expected.Eq. (10b) follows from the Wiener-Khinchin theorem.From Eq. (10c), we can derive the following condition for the intrinsic spin noise to be dominated by the field noise-driven spin noise.
Here, N a denotes the number of atoms in the collection [12].In 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 Section 3, this will hold for Γ xx or yy ≥ 10 −5 nT 2 • ms which corresponds to δB α ≥ 100 fT/ √ Hz where δB 2 α := Γ αα .Note that by applying P = 2S z ≈ R/Γ 1 and δB 2 α = Γ αα to Eq. (10c) and then by taking square roots, we arrive at the expression σ ex ≃ |γP | δB 2 x + δB 2 y /(4 √ Γ 2 ) which was previously noted in the Introduction 1.Such a slight change in notation explicitly reveals the role of the degree P of the spin polarization.
In equation (10c), 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 equation (10b).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 the Introduction.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. 1a).On the other hand, in case ω ≫ Ω, the precession axis ⃗ B + δ ⃗ B can now be imagined to be permanently tilted (Fig. 1b).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 equation (10b).
Also, 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 .
From the equation, it is clear that the noiseinduced signal power is minimized at θ = π/4 − arg (ω + iΓ 2 )/2.Finally, we note a non-Lorentzian behavior of the first-order PSD given in Eq. (10b).First, note that the two-sided PSD, i.e.PSD Sx (f )/2, equals the following expression for all Ω: 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 case of an x-directional noise; in 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 the Ωs in the numerators by ω, the equation reduces to a sum of Lorentzians as 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: Although the above equation is written for δ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 secondorder analysis will require a quadruple integral term in equation ( 4 where R denotes the region defined by t 0 ≤ t 4 ≤ t 3 ≤ t 2 ≤ 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.(The number 4 in the right-handside will be replaced by powers of two for the analogous higher-order equations.)As a result, the first non-zero higher order term in equation ( 4) is simplified as the following expression: A similar inspection shows that Eq. ( 7) 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 the formulae (9d) and (26).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 Γ 2 zz is zero.Moreover, the coefficients for Γ 2  xx and Γ 2 yy are both negative at resonance.For completeness, we inform the reader that the correction terms corresponding to Γ Since we gave the second-order terms corresponding to Γ 2 xx and Γ 2 yy , 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 formulae 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 although we still need to calculate the steady state values as before (as in Eq. ( 9)).Such an effective renormalization of the Bloch operator by a noise term D was previously noted in [16], wherein the effects of longitudinal noise component were investigated.Here, we may extend the results in [16] by computing the eigenvalues of the renormalized Bloch operator.In the case of Γ yz = Γ zx = 0, we have the following eigenvalues: where The first (or latter) formula in Eq. ( 21) 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 towards the DC (0 Hz) direction, which is not obvious at first glance as in 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 ⃗ B tot = ⃗ B + δ ⃗ B(t).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. (22).Moreover, as can be seen from Eq. ( 22), 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).
xa x i s 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 (solid red) points upward, whereas the second rotation (solid black) points downward.The dashed blue curve corresponds to the xy-planar circle of radius 1/2 centered at the origin.The dash-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.
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 ⃗ S • ⃗ S.

Comparison between the theory and stochastic simulations
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 [26] for a less technical review of the Euler-Maruyama method, section 10.2 of [27] for a convergence proof of the method, and [28,29] for an introduction to stochastic differential equations.In case δ ⃗ B = δBn with ⟨δB(t)δB(t ′ )⟩ = Γδ(t − t ′ ), the Bloch equation ( 1) can be re-written in the following form, which is convenient for computer simulations: 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 87 Rb).The z-directional magnetic field was set to 100 nT (Fig. 4a) or 1000 nT (Fig. 4b), and the noises were characterized by Γ αα = 1 nT 2 • ms, where α = x or y; this corresponds to approximately 32 pT/ √ Hz.The stochastic Bloch equation was integrated for a time period of 12 ms, starting from its (zerothorder) macroscopic equilibrium R 2Γ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 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 intuitionbased predictions.(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.
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-axis PSD at resonance when noise is along x-axis .
(24) Fig. 5 shows that the theory and the numerical experiment are consistent.Fig. 5: Theory vs. Monte Carlo for the ratio described in Eq. (24).Each point comes from an average of 50 independent Monte Carlo experiments.
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 (12).Now, the longitudinal B-field is set to B = 227.36nT (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 S 2 x .The results are given in Fig. 6, which show an agreement between the theory and the experiment.6: Theory vs. Monte Carlo, for noise power described in Eq. (12).Here, the B-field is set to make ω = Γ 2 Each point comes from an average of 50 independent Monte Carlo experiments.The dash-dotted red line indicates for the theoretical minimum.
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 equation (9e) 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.

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 a 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 ydirectional 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 towards the DC direction when one observes the PSD.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.
A Simulation parameters for Fig. 1c For Fig. 1c, the simulation was again done using the Euler-Maruyama method.B was set to 500 nT; other parameters were the same as in Section 3. 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 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.

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 abovedescribed procedure was also used for determining the time step size for Fig. 1c; after the time step size was determined, corresponding ∆W sequence was calculated and then fed into the high-pass and low-pass filters described in Appendix A.

C Omitted formulae
and respectively.Somewhat interestingly, D ≃ for Γ 1 τ , Γ 2 τ ≫ 1, and this expression equals the right-hand-side of Eq. (9d).Also, note that the equations ( 25) and ( 26) exhibit some longitudinaltransversal 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 [30].Indeed, if one performs the relevant second-order calculations, the following second-order DC term occurs in the autocorrelation function of S x : Γ, appeared in Eq. ( 7), is given as: where (31) Finally, we provide the second-order corrections for the PSD of S x corresponding to the products Γ xx Γ zz , Γ yy Γ zz , and Γ xx Γ yy .They are given as and where C i are given by the following equations.

Fig. 1 :
Fig.1:(a) Intuition of the system's response to the 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 the low-frequency components.In this case the precession axis is somewhat fixed; therefore, the spin vector can precess around a well-defined axis for a while.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 ⟨S ∥ S ⊥ ⟩ appears to be non-zero, as can be verified from Eq. (9b).

Fig. 4 :
Fig. 4: (a) Theory versus Monte Carlo where B = 100 nT.Black plots correspond to x-directional noise input whereas blue plots correspond to ydirectional noise input.Monte Carlo simulations were repeated 100 times and the averaged results are shown.The dashed red 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. 7 :
Fig. 7: Theory (up to first and second-order) versus Monte Carlo, for PSD values at resonance as a function of Γ αα .Here, B = 500 nT.Subfigure (a) corresponds to x-directional noises, whereas subfigure (b) corresponds to y-directional noises.Each point comes from an average of 300 independent Monte Carlo experiments.