# Inter-species spin-noise correlations in hot atomic vapors

K. Mouloudakis1,2 F. Vouzinas1 A. Margaritakis1 A. Koutsimpela3 G.
Mouloudakis1,4 V. Koutrouli1,4 M. Skotiniotis5,6 G. P. Tsironis1,7 M.
Loulakis3,8 M. W. Mitchell2,9 G. Vasilakis4<EMAIL_ADDRESS>I. K.
Kominis1,7,10<EMAIL_ADDRESS>1Department of Physics, University of Crete,
Heraklion 71003, Greece
2ICFO - Institut de Ciències Fotòniques, The Barcelona Institute of Science
and Technology, 08860 Castelldefels (Barcelona), Spain
3School of Applied Mathematical and Physical Sciences, National Technical
University of Athens, 15780 Athens, Greece
4Institute of Electronic Structure and Laser, Foundation for Research and
Technology, 71110 Heraklion, Greece
5Física Teòrica: Informació i Fenòmens Quàntics, Departament de Física,
Universitat Autònoma de Barcelona, 08193 Bellaterra, Spain
6Departamento de Electromagnetismo y Física de la Materia, Universidad de
Granada, 18010 Granada, Spain
7Institute of Theoretical and Computational Physics, University of Crete,
Heraklion 70013, Greece
8Institute of Applied and Computational Mathematics, Foundation for Research
and Technology, 70013 Heraklion, Greece
9ICREA - Institució Catalana de Recerca i Estudis Avançats, 08010 Barcelona,
Spain
10Quantum Biometronics, Heraklion 71409, Greece

###### Abstract

We report an experimental and theoretical study of spin noise correlations in
a 87Rb-133Cs unpolarized alkali-metal vapor dominated by spin-exchange
collisions. We observe strong unequal-time inter-species correlations and
account for these with a first-principles theoretical model. Since the two
atomic species have different spin precession frequencies, the dual-species
vapor enables the use of an additional experimental handle, the applied
magnetic field, for untangling various sub-types of spin correlations. In
particular, the measured cross-correlation and auto-correlation spectra shed
light on a number of spin-dynamic effects involving intra-atom, inter-atom,
intra-species and inter-species correlations. Cross-correlation coefficients
exceeding $60\%$ have been observed at low-magnetic fields, where the two spin
species couple strongly via spin-exchange collisions. The understanding of
such spontaneously generated correlations can motivate the design of quantum-
enhanced measurements with single or multi-species spin-polarized alkali-metal
vapors used in quantum sensing applications.

##  I Introduction

Quantum measurements are at the basis of quantum technologies, from atomic
magnetometers [1], atomic clocks [2] and quantum optical measurements [3], to
quantum computers [4], quantum simulators [5, 6], and even the detection of
gravitational waves [7, 8]. Quantum uncertainty and its dynamic manifestation
imposes limitations to the precision of quantum measurements [9]. Thus a major
effort of modern quantum technology has been to engineer quantum noise, for
example by generating squeezed states of light or atoms, in order to surpass
what are broadly known as the standard quantum limits to measurement precision
[10]. These are limits usually applied to many-body systems, following under
the working assumption of the particles being in separable quantum states,
i.e. sharing no correlations. Yet, correlations are at the heart of the second
quantum revolution [11]. Clearly, engineering quantum noise towards advancing
the capabilities of quantum technology requires a profound understanding of
the physics of noise, taking correlations into account [12].

Here we unravel multifaceted spin correlations spontaneously generated in a
dual-species hot alkali-metal vapor. Hot atomic vapors are instrumental in
quantum sensing of magnetic fields [13], as well as in atomic vapor clocks,
while noble gas ensembles have been recently used for implementing quantum
information protocols [14, 15, 16], including quantum memories [17]. Moreover,
numerous quantum systems can be mapped to spin, hence spin-noise studies in
atomic vapors find direct analogies to other quantum technologies [18, 19].

We present a precision spin-noise measurement in an unpolarized 87Rb-133Cs
alkali-metal ensemble dominated by spin-exchange interactions. We demonstrate
strong inter-species correlations at low magnetic fields, which fade away at
increasing magnetic fields, as the two atomic species have different
gyromagnetic ratios. The developed theoretical framework takes advantage of
the long-standing physical description of spin-exchange collisions [20], and
leads to a formal description of spin-noise correlations, showing excellent
agreement with the measurements. In particular, the magnetic-field dependence
of the correlations in conjunction with the theoretical framework allows to
resolve intra-atom from inter-atom correlations, and for each type discern
intra-hyperfine from inter-hyperfine correlations.

This work has direct ramifications for understanding quantum limits to sensing
technologies. This is because the observed correlations are created
spontaneously by the ubiquitous spin-exchange collisions. Hence, what might
have been understood as a vapor consisting of uncorrelated atoms in the
absence of external perturbations, is actually an atomic vapor rich in
correlations. These directly influence spin-noise benchmarks against which any
noise engineering protocols which do involve external perturbations have to be
compared.

Moreover, while this measurement is performed with unpolarized atomic vapors,
the developed theoretical framework allows us to extrapolate to polarized
vapors pertinent to quantum sensing applications. In this regime we expect
strong inter-species quantum correlations, which could further advance quantum
metrology with multi-species hot atomic vapors.

The structure of the paper is the following. In the next section we provide a
detailed perspective of this work in the context of previous work on quantum
sensing with hot atomic vapors and spin-noise spectroscopy. In Sec. III we
describe the experimental setup and define the measured observables. In Sec.
IV we present a comprehensive theoretical analysis of spin-noise correlations.
We then analyze the data in Sec. V, present the main results in Sec. VI, and
conclude with Sec. VII. Technical derivations are left for the appendices A-E.

##  II This work in the context of previous work

Composite quantum systems, like a collection of interacting atoms or molecules
and their interface with light, have been widely used as a realization of
quantum sensing technologies [21, 22, 23, 24, 25]. In particular, hot alkali-
metal vapors form the core quantum system in optical magnetometry [26, 27, 28,
29, 30], comagnetometry [31, 32, 33, 34], magnetic field gradiometry [35, 36],
frequency standards [37] and quantum communications [17]. The long coherence
times in combination with technical advantages, like optical accessibility
with commercially available resonant laser light, and the high reliability of
accommodating experimental setups, render hot alkali vapors favorable in many
quantum sensing applications, including magnetoencephalography [38, 39], time
keeping [40], inertial sensing [41, 42] and imaging (THz imaging, biomagnetic
imaging) [43, 44].

Spin-exchange (SE) collisions, deriving from the Pauli exchange interaction
during binary atomic encounters, are central in the physics of hot alkali
vapors. The early understanding of SE collisions [45, 46] accounted for
several experimental observables [47]. Comprehension of more subtle aspects of
SE [48, 49, 50] led to the development of SERF (spin-exchange-relaxation-free)
magnetometers [51, 52, 53], SE optically-pumped hyperpolarized noble gases
utilized in medical applications [54, 55, 56], and hydrogen magnetometry [57,
58, 59].

Rarely did the concept of correlations appear in these works, due to the
intuitive expectation that, in hot vapors dominated by random binary spin-
dependent collisions, correlations can be hardly sustained for meaningful
time-scales [60]. However, several works have recently made such a counter-
intuitive case on the non-trivial role of correlations, pointing to the
possibility of performing quantum-enhanced measurements with hot vapors [61,
62, 63, 64, 65, 66, 67, 68, 69, 70].

###  II.1 Spin-noise spectroscopy

In such studies, unpolarized alkali-metal vapors are rather useful, because
unlike spin-polarized states, they are not sensitive to technical (e.g.
magnetic) noise, hence intrinsic spin fluctuations can be readily measured.
Since correlations and fluctuations are intimately related, as will be
elaborated in detail in this work, unpolarized vapors are a natural testbed
for studying atomic correlations.

In more detail, unpolarized vapors provide easy access to the spontaneous
fluctuations of atomic spin driven by atomic collisions. The field studying
such fluctuations, spin-noise spectroscopy (SNS) [71, 72, 73, 74, 75, 76, 77,
78, 79], is interesting in its own right, as spin noise reveals spectroscopic
information about the atomic [80] or even the solid-state system under
consideration [81, 82, 83] in a non-perturbing way. Further, the study of spin
noise in atomic vapors has elucidated quantum-non-demolition measurements [61,
84, 85, 67], effects of atomic diffusion [86, 87, 88], noise studies at low-
fields [89, 90], optical spin-noise amplification [91], SNR enhancement by
squeezed-light [92], non-equilibrium SNS [93] and spin-alignment noise in a
133Cs vapor [94].

With respect to this work, however, the study of spontaneous spin noise in
unpolarized vapors offers two additional advantages. First, it allows us to
extrapolate the underlying physics to spin-projection noise in spin-polarized
vapors [95], and thus helps understand and define benchmarks against which any
quantum sensing enhancements are compared with. Second, by using unpolarized
vapors, in particular, a dual-species vapor, one can untangle various sorts of
spin correlations, since they leave a clear signature in measurable spectral
distributions of spin-noise power [96].

Overlapping (dual-species) alkali-metal ensembles have been explored a while
ago in the context of hybrid optical pumping [97, 98], specifically addressing
the deterministic spin dynamics at non-zero spin-polarizations. More recently,
two works [99, 100] studied spin-noise correlations and spin-noise transfer
between two alkali species in an unpolarized vapor, arriving at conflicting
results, while the authors in [101] studied the strong SE-coupling of two
alkali-metal species in a polarized vapor. Before summarizing the results of
this work, we introduce some basic notions regarding atomic spin correlations.

###  II.2 Atomic Correlations

When it comes to correlations, one needs to define (I) the correlated degrees
of freedom, (II) their dynamic aspect, (III) the various sub-types, and (IV)
their quantum/classical character [12].

Figure 1: Schematic diagram of atomic correlations measurable in a spin-noise
experiment. (a) A hot and unpolarized alkali-metal vapor is non-destructively
probed by a laser and a balanced polarimeter, comprised of a polarizing beam
splitter and two photodetectors, recording the polarization fluctuations of
the probe light. (b) The light’s polarization fluctuations reflect the
spontaneous fluctuations of the collective spin component, ${\cal F}$, along
the probe beam direction. (c) The measurement directly leads to the correlator
$\langle\hat{{\cal F}}(t+\tau)\hat{{\cal F}}(t)\rangle$, which contains an
intra-atom part, $\sum_{i=j}^{N_{\rm
at}}\langle\hat{\mathsf{f}}^{(i)}(t+\tau)\hat{\mathsf{f}}^{(j)}(t)\rangle$,
and inter-atom part, $\sum_{i\neq j}^{N_{\rm
at}}\langle\hat{\mathsf{f}}^{(i)}(t+\tau)\hat{\mathsf{f}}^{(j)}(t)\rangle$.
Here $\hat{\mathsf{f}}^{(i)}$ derives from the probed component of the total
spin operator of the $i$-th atom, and $N_{\rm at}$ the total number of atoms
probed by the laser. In their turn, such correlations are further divided into
intra- and inter-hyperfine correlations, where the indexes $\alpha,\beta$ take
on the values $a=I+1/2$ and $b=I-1/2$, corresponding to the upper and lower
hyperfine multiplet of an alkali-metal atom with nuclear spin $I$. Every
correlator includes equal-time ($\tau=0$) and unequal-time ($\tau\neq 0$)
correlations, each having a unique footprint on the spin-noise spectra and the
spin-noise variances. In this work we add a further classification: inter-
species versus intra-species correlations.

####  II.2.1 Correlated degrees of freedom

As we will show in detail in Section IV, what is experimentally accessible is
the collective spin correlator

$\braket{\hat{{\cal F}}(t+\tau)\hat{{\cal F}}(t)}$ (1)

where $\hat{{\cal
F}}(t)=\sum_{i=1}^{N_{\mathrm{at}}}\hat{{\mathsf{f}}}^{(i)}(t)$ is the
collective spin of the ensemble along the direction of the laser beam, probing
in total $N_{\mathrm{at}}$ atoms, with $\hat{{\mathsf{f}}}^{(i)}(t)$ being
derived from the total spin of the $i$-th atom,
$\mathbf{\hat{s}}+\mathbf{\hat{I}}$, where $\mathbf{\hat{s}}$ is the
electronic and $\mathbf{\hat{I}}$ the nuclear spin.

####  II.2.2 Dynamic aspect of correlations

We distinguish equal-time ($\tau=0$) from unequal-time correlations ($\tau\neq
0$), where $\tau$ is the time delay between the correlated observables. The
former reflect total spin-noise power readily derived from the Wiener-Khinchin
theorem, which relates the power spectral density of a signal with its
correlation function. The total spin-noise power is obtained by the integral
over frequencies of the power spectral density.

####  II.2.3 Sub-types of correlations

The survival of correlations spontaneously building up by SE collisions in a
hot rubidium vapor was studied in [68], where it was theoretically shown that
equal-time inter-atom correlations in specific states can be sustained amidst
sequential SE collisions of the correlated partners with other atoms.

In [96] the authors studied unequal-time correlations generated by the
dynamics of binary SE collisions at thermal equilibrium pertinent to spin-
noise experiments with unpolarized vapors. The authors considered several
kinds of correlations, schematically depicted in Fig. 1. In particular, the
correlator (1) is composed of two terms
$\sum_{j}\braket{\hat{\mathsf{f}}^{(j)}(t+\tau)\hat{\mathsf{f}}^{(j)}(t)}$ and
$\sum_{i\neq
j}\braket{\hat{\mathsf{f}}^{(i)}(t+\tau)\hat{\mathsf{f}}^{(j)}(t)}$,
describing intra- and inter-atom correlations, respectively. For each such
kind one can further distinguish intra-hyperfine from inter-hyperfine
correlations. To describe those, the total single-atom spin operators are
written as $\hat{f}_{\alpha}^{(i)}(t)$, where $\hat{f}_{\alpha}$ denotes the
projection of the total atom’s spin, $\mathbf{\hat{s}}+\mathbf{\hat{I}}$, onto
the hyperfine manifold $\alpha$, with $\alpha\in\\{a,b\\}$ denoting either the
upper ($a\equiv I+1/2$) or the lower ($b\equiv I-1/2$) ground-state hyperfine
manifold. Thus, as shown in Fig. 1, we resolve the original correlator in four
kinds of correlations: (a) intra-atom and intra- or inter-hyperfine
correlations, and (b) inter-atom and intra- or inter-hyperfine correlations.

In this work we bring into the picture an additional type, the inter-species
correlations. In a single-species vapor the inter-atom correlations always
contribute to the spin-noise spectrum simultaneously with the intra-atom
correlations, and the two are hardly distinguishable. Here comes the
usefulness of a dual-species vapor, in which the inter-species inter-atom
correlations can be experimentally distinguished and thus shed light on inter-
atom correlations even for the single-species case.

####  II.2.4 Character of correlations

The quantum/classical character of the correlations is a subtle issue which
appears to have a parametric dependence, the parameter being the vapor’s spin-
polarization. An example, drawn from quantum information science,
demonstrating a parametric dependence of correlations is the two-qubit Werner
state
$\rho=\frac{1}{4}(\mathbb{1}-\alpha\bm{\hat{\sigma}}\otimes\bm{\hat{\sigma}})$,
where
$\bm{\hat{\sigma}}\otimes\bm{\hat{\sigma}}=\hat{\sigma}_{x}\otimes\hat{\sigma}_{x}+\hat{\sigma}_{y}\otimes\hat{\sigma}_{y}+\hat{\sigma}_{z}\otimes\hat{\sigma}_{z}$,
$\mathbb{1}$ is the identity and $\hat{\sigma}_{i}$, $i\in\\{x,y,z\\}$ are the
Pauli operators [102, 103]. Since the term
$\bm{\hat{\sigma}}\otimes\bm{\hat{\sigma}}$ is traceless, the spin populations
of both qubits in each of the two states along the quantization axis (e.g. the
states with $\sigma_{z}=\pm 1$) are 1/2. Thus the average spin is zero. Yet
the state $\rho$ exhibits correlations for any $\alpha\neq 0$, since
$\braket{\hat{\sigma}_{z}\otimes\hat{\sigma}_{z}}=-\alpha$. In particular, the
correlations are quantum and violate the CHSH inequality for
$\alpha>1/\sqrt{2}$. Thus, the character of the correlations of this 2-qubit
state depends on the parameter $\alpha$. In Sec. VI.3 we will further comment
on this point in regard to this work.

###  II.3 Summary of the results of this work

In this work we undertake (A) an experimental and (B) a formal theoretical
study of spin-noise correlations in a 87Rb-133Cs hot vapor. (A1) We
unambiguously demonstrate the existence of spontaneously generated unequal-
time inter-species spin-noise correlations, driven by spin-exchange
collisions. (A2) Under certain conditions, i.e. finite measurement bandwidth,
we also demonstrate the existence of inter-species correlations reflected in a
total spin-noise power different from what would be obtained were the atoms
uncorrelated. (B1) Using the full density matrix description of spin-exchange
and spin-relaxation dynamics, we develop a first-principles theoretical
framework that captures the subtle physics of all four sub-types of
correlations and resolves inconsistencies of previous works [99, 100]. (B2) We
draw qualitative conclusions about the character of inter-species spin-noise
correlations in spin-polarized dual-species vapors.

Figure 2: Experimental setup for the measurement of 87Rb-133Cs spin-noise
correlations. Two single-mode external cavity diode lasers, blue-detuned from
the corresponding atomic resonances, are combined in a single-mode optical
fiber and directed towards the vapor cell, held inside a 4-layer mu-metal
magnetic shield. A dichroic mirror after the cell is used to spatially
separate the two overlapping wavelengths exiting the atomic medium, while
optical filters further prevent leakage of the unwanted light to the
detectors. Two identical balanced polarimeters comprised of a half-waveplate,
a polarizing beam splitter cube, and a balanced photodetector provide the
signals feeding the data acquisition system. Both wavelengths were monitored
within $10$ $\mathrm{MHz}$ resolution using a Fizeau wavelength-meter (High
Finesse WS7). Single-mode operation was additionally monitored by a second
scanning Fabry-Perot interferometer.

##  III Experimental setup, raw data and observables

In this section we describe the experimental setup, and define the
experimental observables to be used in the theoretical treatment of the
following section.

###  III.1 Experimental setup

The experimental setup is shown in Fig. 2. A cylindrical cell of diameter
$12.7\text{\,}\mathrm{mm}$ and length $47\text{\,}\mathrm{mm}$, with anti-
reflection coated windows, contains 330 Torr of N2 buffer gas and a metallic
droplet of 133Cs and 87Rb in molar ratio 28:72. The cell resides in a ceramic
oven heated with $200\text{\,}\mathrm{kHz}$ current to
$160\text{\,}\mathrm{\SIUnitSymbolCelsius}$, resulting in number densities
$n_{{}_{\mathrm{Cs}}}\approx
n_{{}_{\mathrm{Rb}}}\approx${10}^{14}\text{\,}{\mathrm{cm}}^{-3}$$.

In the presence of nitrogen the optical linewidths are pressure-broadened
(FWHM) to $7.34$ $\mathrm{GHz}$ for 133Cs and $7.08$ $\mathrm{GHz}$ for 87Rb,
therefore the two ground-state hyperfine levels are moderately resolved. The
collective atomic spins of the alkali-metal species are probed by two single-
mode external cavity diode lasers (Toptica DL-780 pro and DL-850 pro) with
their wavelengths blue-detuned from the $\rm{D}_{2}$ transition by several
tens of $\mathrm{GHz}$. Both linearly polarized laser beams enter in the same
single-mode optical fiber producing a two-color Gaussian beam with diameter
$\sim 2$ mm ($1/e^{2}$ intensity width) at the position of the cell (with a
small difference in the two colors of less than $10\%$ due to diffraction).

The vapor cell is placed inside a $4$-layer mu-metal shield, protecting the
atomic spins from ambient magnetic fields. Within the shields, a coil system
generates a DC magnetic field $\mathbf{B}=(0,0,B)$ transverse to the beam-
propagation axis, ranging from $4$ mG to $92$ mG, with the corresponding spin
precession frequencies reaching several tens of $\mathrm{kHz}$.

At the exit of the cell the two-color beam is incident on a dichroic mirror,
which directs each color to a separate balanced polarimeter, the two
polarimeters being otherwise identical. The use of 780 $\mathrm{nm}$ and 852
$\mathrm{nm}$ filters further reduces light leakage from one wavelength to the
other, resulting in negligible cross-talk in the polarization detection of the
two beams. Each balanced photodetector has detection efficiency $\eta$.

The cesium and rubidium polarimeter signals are acquired and analyzed in the
frequency domain by feeding them into a spectrum analyzer (Stanford Research
Systems SR780 Dynamic Signal Analyzer), which provides the cross-correlation
power spectral density and the two single-channel auto-correlation power
spectra. The frequency resolution of the FFT is $800$ bins, with a sampling
rate of $262$ $\mathrm{kHz}$. A Hann window is applied and an analog anti-
aliasing filter applied before digitization eliminates all frequency
components above $102.4$ $\mathrm{kHz}$.

###  III.2  Auto-correlation and cross-correlation observables

Each balanced photodetector receives a continuous optical signal and generates
a corresponding photo-current, which we denote by $i_{\beta}(t)$, with
$\beta\in\\{\text{Rb},\text{Cs}\\}$. Spin noise (SN) and photon-shot-noise
render $i_{\beta}(t)$ a classical and real stochastic process. We record the
signals for long enough time such that both spin-noise and photon-shot-noise
have stationary statistical moments. We are interested in the cross-
correlation between the two balanced photodetector signals:

$C_{\text{Rb},\text{Cs}}(\tau)\equiv\langle\langle
i_{\text{Rb}}(t+\tau)i_{\text{Cs}}(t)\rangle\rangle,$ (2)

where $\langle\langle\cdot\rangle\rangle$ denotes the average over all
possible realizations of the stochastic measurement outcomes. Stationarity
implies that the correlation does not depend on the absolute time $t$.

The quantum-mechanical operator associated with the measured photocurrent is:

$\begin{split}\hat{I}_{\beta}(t)&=q_{e}\eta_{\beta}\left(\hat{a}_{k,+45}^{{\dagger}}\hat{a}_{k,+45}-\hat{a}_{k,-45}^{{\dagger}}\hat{a}_{k,-45}\right)_{\beta}\\\
&=2q_{e}\eta_{\beta}\hat{{\cal S}}_{2,\beta}(t),\end{split}$ (3)

where $\hat{a}_{k,\ell}\equiv\hat{a}_{k,\ell}(t)$ and
$\hat{a}^{{\dagger}}_{k,\ell}\equiv\hat{a}^{{\dagger}}_{k,\ell}(t)$ are
respectively the travelling-wave annihilation and creation operators for the
mode $k$ and polarization $\ell$, normalized so that
$\hat{a}^{{\dagger}}_{k,\ell}\hat{a}_{k,\ell}$ is the photon flux in the
respective mode [104]. The index $k$ labels the longitudinal as well as the
spatial mode, the latter being defined by the spatial mode of the probe beam
(Gaussian $\text{TEM}_{00}$ in the experiment). The $\pm 45^{\circ}$ linear
polarizations are characterized with respect to the polarization direction of
the probe beam before the interaction with the atomic medium, chosen to be
along the $\mathbf{\hat{y}}$-direction. Finally, $\hat{{\cal S}}_{2}(t)$ is
the optical Stokes operator given by:

$\displaystyle\hat{{\cal S}}_{2}(t)$
$\displaystyle=\frac{\hat{a}_{k,+45}^{{\dagger}}\hat{a}_{k,+45}-\hat{a}_{k,-45}^{{\dagger}}\hat{a}_{k,-45}}{2}$
$\displaystyle=\frac{\hat{a}_{k,y}^{{\dagger}}\hat{a}_{k,z}+\hat{a}_{k,z}^{{\dagger}}\hat{a}_{k,y}}{2},$
(4)

where the coordinate system is defined in Fig. 2. The operators $\hat{{\cal
S}}_{2,\beta}$ refer to the Stokes operator $\hat{{\cal S}}_{2}$ for each of
the two wavelenghts probing atom species $\beta$.

The cross-correlation can be therefore expressed as [9, 105]:

$\begin{split}C_{\text{Rb},\text{Cs}}(\tau)&=\langle:\hat{I}_{\text{Rb}}(t+\tau)\hat{I}_{\text{Cs}}(t):\rangle\\\
&=4q_{e}^{2}\eta_{\text{Rb}}\eta_{\text{Cs}}\langle:\hat{{\cal
S}}_{2,\text{Rb}}(t+\tau)\hat{{\cal S}}_{2,\text{Cs}}(t):\rangle,\end{split}$
(5)

where $::$ denotes normal ordering. In the above equation, the quantum
operators are written in the Heisenberg picture and the single brackets
$\langle\cdot\rangle$ denote quantum mechanical expectation value. Taking into
account the bosonic commutation relations [9, 105], and considering that the
off-resonance interaction with the atomic vapor is to a very good
approximation a linear optical process and thus does not mix annihilation and
creation operators, we can make use of the unequal-time commutation relations:

$\left[\hat{a}^{{\dagger}}_{k,\ell}(t),\hat{a}^{{\dagger}}_{k^{\prime},\ell^{\prime}}(t^{\prime})\right]=\left[\hat{a}_{k,\ell}(t),\hat{a}_{k^{\prime},\ell^{\prime}}(t^{\prime})\right]=0,$
(6)

and show
$\langle:\hat{I}_{\text{Rb}}(t+\tau)\hat{I}_{\text{Cs}}(t):\rangle=\langle:\hat{I}_{\text{Cs}}(t)\hat{I}_{\text{Rb}}(t+\tau):\rangle$.
In order to emphasize that the correlation function is invariant under the
preceding exchange of operators we express Eq.(5) in the symmetrized form
[105]:

$\displaystyle C_{\text{Rb},\text{Cs}}(\tau)$
$\displaystyle=2q_{e}^{2}\eta_{\text{Rb}}\eta_{\text{Cs}}\Big{[}\langle\hat{{\cal
S}}_{2,\text{Rb}}(t+\tau)\hat{{\cal S}}_{2,\text{Cs}}(t)\rangle$
$\displaystyle+\langle\hat{{\cal S}}_{2,\text{Cs}}(t)\hat{{\cal
S}}_{2,\text{Rb}}(t+\tau)\rangle\Big{]}.$ (7)

The cross-spectral density function is the Fourier transform of the cross-
correlation:

$\tilde{C}_{\text{Rb},\text{Cs}}(\nu)=\int_{-\infty}^{\infty}C_{\text{Rb},\text{Cs}}(\tau)e^{-i2\pi\nu\tau}d\,\tau.$
(8)

The spectrum analyzer estimates the power spectral density (PSD) by performing
a finite Fourier transform for a specified record time $T$ and averages the
product of the Fourier components over the different realizations of the
stochastic processes [106]:

$\displaystyle\tilde{C}^{\text{SA}}_{\text{Rb},\text{Cs}}(\nu)=$
$\displaystyle\frac{1}{T}\langle\langle\left(\int_{0}^{T}i_{\text{Rb}}(t)e^{-i2\pi\nu
t}\,dt\right)^{*}\int_{0}^{T}i_{\text{Cs}}(s)e^{-i2\pi\nu
s}\,ds\rangle\rangle$
$\displaystyle=\int_{-T}^{T}C_{\text{Rb},\text{Cs}}(\tau)e^{-i2\pi\nu\tau}\left(1-\frac{|\tau|}{T}\,d\tau\right).$
(9)

To derive the second line, the variables of integration were changed from
$(t,s)$ to $(t,\tau=t-s)$ and the region of integration was modified
accordingly [107]. Thus, the spectrum analyzer performs an exact evaluation of
the spectral density (Eq. (8)) in the limit of infinite long sample
realization $T\rightarrow+\infty$. In practice, for an accurate approximation,
it is sufficient to choose the sample length to be much larger than the time
difference $\tau$ at which the correlation function becomes negligible.

Similarly, single-channel auto-correlation functions and their corresponding
power spectral densities are defined as $C_{\beta,\beta}(\tau)=\langle\langle
i_{\beta}(t)i_{\beta}(t+\tau)\rangle\rangle$ and
$\tilde{C}_{\beta,\beta}(\nu)=\int_{-\infty}^{\infty}C_{\beta,\beta}(\tau)e^{-i2\pi\nu\tau}\,d\tau$,
respectively, with $\beta\in\\{\text{Rb},\text{Cs}\\}$. As before, the auto-
correlation function can be expressed in the symmetrized form:

$\displaystyle C_{\beta,\beta}(\tau)$
$\displaystyle=2q_{e}^{2}\eta_{\beta}\times$
$\displaystyle\left(\langle\hat{{\cal S}}_{2,\beta}(t)\hat{{\cal
S}}_{2,\beta}(t+\tau)\rangle+\langle\hat{{\cal S}}_{2,\beta}(t+\tau)\hat{{\cal
S}}_{2,\beta}(t)\rangle\right).$ (10)

##  IV  Theoretical Derivation of Cross- and Auto-Correlation Spectra

The theoretical description of the cross correlations is based on treating the
atom-light coupling, and the atomic state evolution due to the atom’s
hyperfine structure hamiltonian, atom-atom spin-exchange collisions, and a
number of additional spin relaxation phenomena. Such topics have been
described in detail elsewhere [108, 109, 110, 111], and we recapitulate the
description for the sake of completeness in Appendix A. In the following two
subsections we present the results of extending such treatments to the case of
the two-species atomic vapor treated herein.

###  IV.1 Atom-light coupling

Based on the light-atom interaction dynamics, it follows that for probe light
linearly polarized in the $y$-direction, the solution of the Heisenberg
equation of motion for the Stokes operator yields in the small-angle
approximation[112, 113]

$\hat{{\cal S}}_{2}^{(\text{out})}(t)\approx\hat{{\cal
S}}_{2}^{(\text{in})}(t)+\frac{1}{2}\Phi\bar{g}\mathcal{\hat{F}}(t),$ (11)

where $\Phi$ is the photon-flux (photons per unit time). The out (in)
superscripts denote the operator after (before) interacting with the atomic
sample, $\bar{g}=(|g_{a}|+|g_{b}|)/2$ is the mean coupling constant to the two
ground-state hyperfine manifolds and
$\mathcal{\hat{F}}(t)=\sum_{i=1}^{N_{\mathrm{at}}}\hat{\mathsf{f}}^{(i)}_{x}(t)$
is the probed collective atomic spin with

$\hat{\mathsf{f}}^{(i)}_{x}(t)=\frac{g_{a}}{\bar{g}}\hat{f}_{a,x}^{(i)}(t)-\frac{g_{b}}{\bar{g}}\hat{f}_{b,x}^{(i)}(t)$
(12)

Here $\hat{f}_{a,x}^{(i)}$ and $\hat{f}_{b,x}^{(i)}$ are the projections of
the $x$-component of the total spin of the $i$-th atom,
$\hat{s}_{x}+\hat{I}_{x}$, onto the upper ($a=I+1/2$) and lower ($b=I-1/2$)
hyperfine multiplet, respectively, with the corresponding coupling constants
being $g_{a}$ and $g_{b}$ (see Appendix A, where we also generalize Eq. (12)
to include a non-uniform intensity distribution of the light).

The ${\rm in}\rightarrow{\rm out}$ change of $\hat{{\cal S}}_{2}$ described by
Eq. (11) reflects polarization rotation of light (Faraday rotation) at an
angle $\phi_{\text{FR}}=\frac{1}{2}\bar{g}\braket{\mathcal{\hat{F}}}$. Thus,
Eq. (11) connects the Stokes operator ${\cal S}_{2}(t)$ to the underlying spin
operator $\mathcal{\hat{F}}(t)$. Using Eq. (7), Eq. (10), and Eq. (12), we
arrive at the following expressions connecting the measured correlations
(cross-correlation and auto-correlations) to the underlying collective spin
correlators:

$\displaystyle C_{\text{Rb},\text{Cs}}(\tau)$
$\displaystyle=\frac{q_{e}^{2}\eta_{\text{Rb}}\eta_{\text{Cs}}\Phi_{\text{Rb}}\Phi_{\text{Cs}}\bar{g}_{\text{Rb}}\bar{g}_{\text{Cs}}}{2}\times$
$\displaystyle\langle\mathcal{\hat{F}}_{\text{Rb}}(t+\tau)\mathcal{\hat{F}}_{\text{Cs}}(t)+\mathcal{\hat{F}}_{\text{Cs}}(t)\mathcal{\hat{F}}_{\text{Rb}}(t+\tau)\rangle,$
(13) $\displaystyle C_{\beta,\beta}(\tau)$
$\displaystyle=2q_{e}^{2}\eta_{\beta}\Phi_{\beta}\Big{[}\delta(\tau)$
$\displaystyle+\frac{1}{2}\eta_{\beta}\Phi_{\beta}\bar{g}^{2}_{\beta}\langle\mathcal{\hat{F}}_{\beta}(t+\tau)\mathcal{\hat{F}}_{\beta}(t)\rangle\Big{]},$
(14)

where $\beta\in\\{\text{Rb},\text{Cs}\\}$. We note that we have taken into
account that the polarization properties of the two beams before the
interaction with the atoms are uncorrelated, i.e. $\langle\hat{{\cal
S}}^{(\rm{in})}_{2,\text{Rb}}(t)\hat{{\cal
S}}^{(\rm{in})}_{2,\text{Cs}}(t^{\prime})\rangle=0$.

###  IV.2 Dynamics of mean spin

By considering the various processes that affect the atomic spin evolution
(see Appendix B), we arrive at two coupled density matrix evolution equations
for Rb and Cs:

$\displaystyle\frac{d}{dt}\rho_{\text{Rb}}$
$\displaystyle=A_{\text{Rb}}\hat{\mathbf{I}}_{\text{Rb}}\cdot\hat{\mathbf{s}}_{\text{Rb}}+g_{s}\mu_{B}\hat{\mathbf{s}}_{\text{Rb}}\cdot\mathbf{B}+R\left(\phi_{\text{Rb}}-\rho_{\text{Rb}}\right)$
$\displaystyle+R_{\rm
se}^{\text{Rb},\text{Rb}}\left\\{\phi_{\text{Rb}}\left(1+4\langle\hat{\mathbf{s}}_{\text{Rb}}\rangle\cdot\hat{\mathbf{s}}_{\text{Rb}}\right)-\rho_{\text{Rb}}\right\\}$
$\displaystyle+R_{\rm
se}^{\text{Rb},\text{Cs}}\left\\{\phi_{\text{Rb}}\left(1+4\langle\hat{\mathbf{s}}_{\text{Cs}}\rangle\cdot\hat{\mathbf{s}}_{\text{Rb}}\right)-\rho_{\text{Rb}}\right\\}$
(15) $\displaystyle\frac{d}{dt}\rho_{\text{Cs}}$
$\displaystyle=A_{\text{Cs}}\hat{\mathbf{I}}_{\text{Cs}}\cdot\hat{\mathbf{s}}_{\text{Cs}}+g_{s}\mu_{B}\hat{\mathbf{s}}_{\text{Cs}}\cdot\mathbf{B}+R\left(\phi_{\text{Cs}}-\rho_{\text{Cs}}\right)$
$\displaystyle+R_{\rm
se}^{\text{Cs},\text{Cs}}\left\\{\phi_{\text{Cs}}\left(1+4\langle\hat{\mathbf{s}}_{\text{Cs}}\rangle\cdot\hat{\mathbf{s}}_{\text{Cs}}\right)-\rho_{\text{Cs}}\right\\}$
$\displaystyle+R_{\rm
se}^{\text{Cs},\text{Rb}}\left\\{\phi_{\text{Cs}}\left(1+4\langle\hat{\mathbf{s}}_{\text{Rb}}\rangle\cdot\hat{\mathbf{s}}_{\text{Cs}}\right)-\rho_{\text{Cs}}\right\\}$
(16)

where $A_{\beta}$ is the hyperfine coupling, the rate $R$ includes all the
relaxation processes, other than the spin-exchange relaxation, that destroy
electron polarization without affecting the nucleus, $\phi_{\beta}$ is the
atom’s density matrix with zero electronic polarization, and $R_{\rm
se}^{\beta,\gamma}$ is the spin-exchange rate transferring spin-polarization
from atomic species $\gamma$ to species $\beta$, with
$\beta,\gamma\in\\{\text{Rb},\text{Cs}\\}$. The diffusion of atoms out of the
probe beam must also be included in these dynamics (see Appendix B).

Multiplying both sides of Eqs. (15) and (16) by the spin operator
$\hat{f}_{i}$ (here $i$ is a general index for identifying the atomic species,
the hyperfine manifold and the Cartesian component) and taking the trace, the
dynamics of $\langle\hat{f}_{i}\rangle$ are determined and a closed system of
equations can be derived for the time evolution $\langle\hat{f}_{i}\rangle$.
Like in [49, 101, 96], we define the vector $\mathbf{\hat{X}}(t)$ for the
transverse collective-spin components of each species and each hyperfine
level:

$\mathbf{\hat{X}}(t)\equiv\left[\hat{f}_{a,x}^{\text{Rb}},\hat{f}_{a,y}^{\text{Rb}},\hat{f}_{b,x}^{\text{Rb}},\hat{f}_{b,y}^{\text{Rb}},\hat{f}_{a,x}^{\text{Cs}},\hat{f}_{a,y}^{\text{Cs}},\hat{f}_{b,x}^{\text{Cs}},\hat{f}_{b,y}^{\text{Cs}}\right]^{\intercal},$
(17)

where $\hat{f}_{\alpha,q}^{\beta}$ refers to the total atomic spin of species
$\beta\in\\{\text{Rb},\text{Cs}\\}$ along the $q$-axis with $q\in\\{x,y\\}$,
and in the hyperfine state $\alpha\in\\{a,b\\}$. The density matrix evolution
Eqs. (15) and (16) contain nonlinear terms proportional to
$\langle\hat{\mathbf{s}}_{\beta}\rangle\cdot\hat{\mathbf{s}}_{\gamma}$
associated with the SE interaction; however, for noise measurements around
zero mean spin-polarization we linearize such terms by keeping only first-
order contributions from the fluctuations. This approximation leads to the
linear evolution equation

$\frac{d}{dt}\langle\mathbf{\hat{X}}(t)\rangle=A\langle\mathbf{\hat{X}}(t)\rangle,$
(18)

where the drift matrix $A$ is comprehensively derived in Appendices B and C.

###  IV.3 Spin correlations in the time domain

The estimation of the spin-noise spectrum requires the evaluation of
correlators $\langle\hat{f}_{i}(t+\tau)\hat{f}_{j}(t)\rangle$. To find those,
we use the quantum regression theorem [114, 9], which states that if the
expectation values of a set of observables $\hat{M}_{\mu}$, $\mu=1,2,...$
follow a coupled set of linear equations:

$\frac{d}{dt}\langle\hat{M}_{\mu}(t)\rangle=\sum_{\lambda}A_{\mu\lambda}\langle\hat{M}_{\lambda}(t)\rangle,$
(19)

then the two-time correlation functions satisfy the equations (for $\tau\geq
0$):

$\displaystyle\frac{d}{d\tau}\langle\hat{M}_{\kappa}(t)\hat{M}_{\mu}(t+\tau)\rangle$
$\displaystyle=\sum_{\lambda}A_{\mu\lambda}\langle\hat{M}_{\kappa}(t)\hat{M}_{\lambda}(t+\tau)\rangle,$
(20)
$\displaystyle\frac{d}{d\tau}\langle\hat{M}_{\mu}(t+\tau)\hat{M}_{\kappa}(t)\rangle$
$\displaystyle=\sum_{\lambda}A_{\mu\lambda}\langle\hat{M}_{\lambda}(t+\tau)\hat{M}_{\kappa}(t)\rangle.$
(21)

The regression theorem holds when the equations of motion for the expectation
values are linear and the system-environment correlations can be neglected
[115, 116]. The justification for linearity arises from considering small
fluctuations as noted previously. For the second requirement, we argue that
the environment, where spin-information is lost, is associated with the
(abstract) space spanned with all the collisional parameters. As long as a
large number of particles are probed, when the ensemble average is taken over
all the different types of collisions, the correlations between the collective
spin-system and the environment are lost, thus rendering the regression
theorem applicable.

We thus arrive at the symmetrized and real covariance matrix:

$\mathcal{R}(\tau)=\frac{1}{2}\Big{(}\langle\mathbf{\hat{X}}(t)\mathbf{\hat{X}}^{\intercal}(t+\tau)\rangle+[\langle\mathbf{\hat{X}}(t+\tau)\mathbf{\hat{X}}^{\intercal}(t)\rangle]^{T}\Big{)}.$
(22)

The quantum regression approach yields for $\tau\geq 0$:

$\frac{d}{d\tau}\mathcal{R}(\tau)=A\mathcal{R}(\tau)\rightarrow\mathcal{R}(\tau)=e^{A\tau}\mathcal{R}(0).$
(23)

Given that $\mathcal{R}(0)$ is symmetric, for $\tau<0$ it is
$\mathcal{R}(\tau)=\mathcal{R}(-\tau)^{T}$.

In the following we extend the analysis of [96] to the current case of a dual-
species vapor. The matrix $A$ is diagonalizable, so we can write
$e^{A\tau}=Ve^{\Lambda\tau}V^{-1}$, where $V$ is the matrix whose $i$-th
column is the $i$-th eigenvector of $A$, and $\Lambda$ is a diagonal matrix
with its diagonal elements being the corresponding eigenvalues of $A$.
Therefore, according to Eq.(23), the time dependence for any spin correlation
of the form
$\langle\hat{f}_{\alpha,q}^{\beta}(0)\hat{f}_{\alpha^{\prime},q^{\prime}}^{\beta^{\prime}}(0)\rangle$
manifests as a summation over all the eigenvalues of $A$ of exponentials of
the form $e^{\lambda_{k}\tau}$, where $\lambda_{k}$ is the $k$-th eigenvalue
of $A$. The eigenvalues of $A$ are in general complex numbers, but since $A$
has real elements, the eigenvalues appear in complex conjugate pairs. Since
$A$ is an $8\times 8$ matrix, there are 4 pairs of eigenvalues, hence the
elements of the covariance matrix can be written (for $\tau\geq 0$):

$\mathcal{R}^{\beta,\beta^{\prime}}_{\alpha
q,\alpha^{\prime}q^{\prime}}(\tau)=\sum_{k=1}^{4}\Re\left[c_{k}(\beta\alpha
q;\beta^{\prime}\alpha^{\prime}q^{\prime})e^{-\Gamma_{k}\tau+i\Omega_{k}\tau}\right],$
(24)

where $\Gamma_{k}=\Re\left[\lambda_{k}\right]+R_{D}$ and
$\Omega_{k}=\Im\left[\lambda_{k}\right]$ are the decoherence rate and the
precession frequency, respectively, of mode $k$. In the rates $\Gamma_{k}$ we
took into account the relaxation rate, $R_{D}$, due to atoms diffusing out of
the probe beam. The coefficients $c_{k}$ are complex numbers that depend on
the steady-state ($\tau=0$) covariance $R(0)$.

###  IV.4 Spin correlations in the frequency domain

The spectrum, obtained by the Fourier transform of the correlation function
(see Appendix C), is then given as a sum of complex Lorentzians or
equivalently as the sum of dispersive and Lorentzian functions:

$\displaystyle S^{\beta,\beta^{\prime}}_{\alpha
q,\alpha^{\prime}q^{\prime}}(\nu)=\sum_{k=1}^{4}$
$\displaystyle\zeta_{k}(\beta\alpha
q;\beta^{\prime}\alpha^{\prime}q^{\prime})\frac{\Gamma_{k}}{(\nu-\nu_{k})^{2}+(\Gamma_{k}/2\pi)^{2}}+$
$\displaystyle\zeta_{k}^{\prime}(\beta\alpha
q;\beta^{\prime}\alpha^{\prime}q^{\prime})\frac{\nu-\nu_{k}}{(\nu-\nu_{k})^{2}+(\Gamma_{k}/2\pi)^{2}},$
(25)

where $\nu_{k}=\Omega_{k}/2\pi$. In contrast to the auto-correlation spectrum,
the cross-correlation spectrum of Eq.(25) is in general complex, i.e. the
coefficients $\zeta_{k}$, $\zeta_{k}^{\prime}$ can be complex for
$\beta\neq\beta^{\prime}$. Therefore, in order to obtain all the information
related to the spectrum, both the real and imaginary components should be
recorded.

Finally we note that the covariance matrix $\mathcal{R}(0)$ entering Eq.(23)
can be evaluated by integrating the spectrum over all frequencies (see
Appendices A and E). In practice however, the measurement is performed over a
limited frequency range either due to the finite-bandwidth of the
electronics/detectors or because deliberately the experimental application
requires a limited sampling rate. Depending on the experimental conditions,
this may affect the apparent correlations as discussed in Sec. VI.

##  V Experimental data with theoretical fits

We acquire 87Rb-133Cs cross-correlation spectra, as well as single-species
87Rb and 133Cs power spectra, for six different magnetic fields: 4 mG, 6 mG,
12 mG, 24 mG, 46 mG and 92 mG. We verify that the sign of the cross-
correlation signal makes physical sense. Indeed, the experimental observable
(balanced polarimeter output) is a product of the measured atomic spin with an
atom-light coupling factor. The latter has a sign depending on the probe
detuning from the atomic resonance (see Appendix A). The data are acquired
with both probe beams being blue detuned from the corresponding atomic
resonance. As a result, the atom-light scaling factors do not alter the
overall sign of the cross-correlation. Nevertheless, we do verify that by
flipping either one wavelength detuning from blue to red, the cross-spectra
change sign. For the accurate interpretation of the positivity or negativity
of the dual-species spin correlations, it is also necessary to measure the
direction of optical rotation consistently for both beams. We do this by
introducing an achromatic waveplate in the common path of the two beams,
confirming that the two polarimeter outputs change in the same way with the
rotation of the waveplate.

Figure 3: Measured cross correlation spectra (a)-(f) and auto-correlation
power spectra (g)-(l) for the 87Rb (green, higher frequencies) and 133Cs
(cyan, lower frequencies) spin-ensembles at six different magnetic fields,
together with theoretical fits (solid lines). Each spectrum is the average of
5000 runs. With increasing magnetic field, the magnitude of the cross-
correlation peak is seen to drop, since when the difference in precession
frequencies of the two species is larger than the magnetic linewidth, the
spin-exchange coupling of the two different states is averaged out. The
narrowing of the auto-correlation spectra is also evident, as by reducing the
magnetic field spin dynamics gradually enter the SERF regime. (m) Collection
of fitting parameters for magnetic field, self spin-exchange rates, photon
shot noise levels and scaling factors. The magnetic field aside, all
parameters are consistent across the six magnetic field values.

To fit the theoretical model to the data we use as model spectrum the
expression

$\displaystyle S_{\text{model}}^{\text{Rb},\text{Cs}}(\nu)\propto$
$\displaystyle
g_{a}^{\text{Rb}}g_{a}^{\text{Cs}}S_{ax,ax}^{\text{Rb},\text{Cs}}(\nu)-g_{a}^{\text{Rb}}g_{b}^{\text{Cs}}S_{ax,bx}^{\text{Rb},\text{Cs}}(\nu)$
$\displaystyle-
g_{b}^{\text{Rb}}g_{a}^{\text{Cs}}S_{bx,ax}^{\text{Rb},\text{Cs}}(\nu)+g_{b}^{\text{Rb}}g_{b}^{\text{Cs}}S_{bx,bx}^{\text{Rb},\text{Cs}}(\nu),$
(26)

where the atom-light coupling factor $g_{\alpha}^{\beta}$ depends on the
wavelength of the laser probing the $\alpha$ hyperfine spin of the $\beta$
species. We remind the reader that $S_{\alpha
x,\alpha^{\prime}x}^{\beta,\beta^{\prime}}(\nu)$, given by Eq. (25),
corresponds to the cross-spectrum between the $\alpha$ and $\alpha^{\prime}$
hyperfine spins of the $\beta$ and $\beta^{\prime}$ species, respectively,
measured along the $x$-axis.

Moreover, the model also requires as input the steady-state covariance matrix
$R(0)$. We choose a diagonal $R(0)$, because as discussed in Appendix C, if
the imaginary part of the cross spectrum is zero and the spin-variances follow
the scaling outlined in Eq. (73), then $R(0)$ must be diagonal. We have
verified that the measured imaginary part of the cross spectrum is zero within
the measurement resolution. We did so not only for the operating wavelengths
of the two probe lasers, but for four different probe-wavelength pairs
(combination of two different wavelengths for 133Cs and 87Rb probing). This
way, and given the dependence of the coupling factors $g_{\alpha}^{\beta}$ on
the probe light wavelength, we alter the contribution of each of the four
terms of Eq. (26), showing that the zero imaginary part of the cross spectrum
is not accidental, but reflects an underlying property of all cross-spectra
appearing in Eq. (26). Additionally, assuming the aforementioned scaling of
the spin variances, we arrive at the result that $R(0)$ is diagonal.

Both the single-species power spectra and the dual-species cross-spectrum are
used for the optimization of the fitted parameters. That is, the parameters
are adjusted to achieve a minimum in the merit function:

$\sum_{j=1}^{800}\Bigg{\\{}\left[S_{\rm
meas}^{\text{Rb},\text{Cs}}(\nu_{j})-S_{\rm
model}^{\text{Rb},\text{Cs}}(\nu_{j};{\cal P})\right]^{2}+\left[S_{\rm
meas}^{\text{Rb},\text{Rb}}(\nu_{j})-S_{\rm
model}^{\text{Rb},\text{Rb}}(\nu_{j};{\cal P})\right]^{2}+\left[S_{\rm
meas}^{\text{Cs},\text{Cs}}(\nu_{j})-S_{\rm
model}^{\text{Cs},\text{Cs}}(\nu_{j};{\cal P})\right]^{2}\Bigg{\\}},$ (27)

where $S_{\rm meas}^{\beta,\beta^{\prime}}(\nu_{j})$ and $S_{\rm
model}^{\beta,\beta^{\prime}}(\nu_{j};{\cal P})$ is the experimental and the
model’s prediction for the noise spectrum between spin species $\beta$ and
$\beta^{\prime}$, measured and calculated at the $j$-th frequency bin of the
spectrum analyzer, respectively, where $j=1,2,...,800$. The symbol ${\cal P}$
represents the set of all fitted parameters, including the spin-exchange rates
($R^{\text{Rb},\text{Rb}}_{\rm{se}}$, $R^{\text{Cs},\text{Cs}}_{\rm{se}}$),
the S-damping rate ($R$), the spin relaxation rate due to diffusion
($R_{\text{D}}$), the magnetic field ($B_{0}$), the photon shot noise levels
(${\rm PSN}_{\text{Rb}}$, ${\rm PSN}_{\text{Cs}}$), and the scaling factors
($\mathcal{K}_{\text{Rb}}$, $\mathcal{K}_{\text{Cs}}$) for each of the Rb and
Cs power spectrum. The scaling factor of the Rb-Cs cross-spectrum is
$\sqrt{\mathcal{K}_{\text{Rb}}\mathcal{K}_{\text{Cs}}}$ and does not appear as
independent variable for the fitting.

In Fig. 3 we present the data (Rb-Cs cross-correlation spectrum and single-
species power spectra) with the result of the global fit, showing very good
agreement with the theoretical model. In the table shown in Fig. 3m we
summarize the fit parameters. The fitted value of $B_{0}$ follows the expected
value based on the applied current to the coil within the magnetic shields,
with a deviation of only a few percent at the smallest field values. The
fitted values of the spin-exchange rates also agree within $10\%$ with the
values derived from the spin-exchange cross-sections reported in literature
[117]. Importantly, the fitted values for the spin-exchange rates, photon shot
noise levels, and scaling factors were consistent across the six different
magnetic field values, demonstrating the internal consistency of the
theoretical model, which captures the global magnetic-field dependence of the
data. Lastly, the S-damping and diffusion rates could not be accurately
determined from the fit, likely because they are two orders of magnitude
smaller than the spin-exchange rates. To address this issue, the fitting of
those rates was constrained within a range spanning from a factor of 3 below
to a factor of 3 above the expected rates based on the relevant experimental
paarameters.

##  VI  Results and Discussion

In the cross-correlation spectra (Figs. 3a-f), a prominent positive peak is
observed at low magnetic fields in the frequency range near the average of the
resonance frequencies of the two species, as determined by their power spectra
(Figs. 3g-l). The peak height drops and becomes broader with increasing
magnetic field. This behavior is primarily caused by the spin-exchange
collisions between the alkali-metal atoms and reflects how the effect of these
collisions changes with the magnetic field.

The positive and negative swings of the cross-correlation spectrum indicate
that band-limited measurements can exhibit either positive or negative
correlations between the two spin species. In other words, the value and sign
of the correlations depend on the bandwidth of the measurement and on the
central measurement frequency. To illustrate this point, we consider the
cosine quadratures
$\widetilde{\mathcal{F}_{\beta}}=(1/T_{\text{BW}})\int_{0}^{T}\mathcal{F}_{\beta}(t)e^{(T-t)/T_{\text{BW}}}\cos(\Omega
t)dt$ of the (collective, transverse) spin for each of species
$\beta\in\\{\text{Rb},\text{Cs}\\}$, and examine their correlation,
$\langle\widetilde{\mathcal{F}_{\text{Rb}}}\widetilde{\mathcal{F}_{\text{Cs}}}\rangle$,
where $T$ is the measurement time, $\Omega$ is the frequency of the harmonic
quadrature, and $1/(2\pi T_{\text{BW}})$ is the integration bandwidth. We
parenthetically note that such cosine (or sine) quadratures are typically
measured in AC magnetometers with a lock-in amplifier. This correlation can be
expressed as an integral of the inter-species cross-correlation spectrum,
$\langle\widetilde{\mathcal{F}_{\text{Rb}}}\widetilde{\mathcal{F}_{\text{Cs}}}\rangle=\int_{0}^{\infty}S^{\text{Rb},\text{Cs}}_{\rm
meas}(\nu)\phi(\nu)d\nu$, where $\phi(\nu)$ is a kernel function that accounts
for the effective measurement bandwidth and depends on the quadrature
frequency $\Omega$ and the integration bandwidth $T_{\text{BW}}$, with a
negligible effect from the measurement time $T$ when $T\gg T_{\text{BW}}$ and
$\Omega T_{\text{BW}}\gg 1$ (see Appendix E for an explicit formula for the
kernel function). We characterize the strength of this cross-correlation with
the coefficient

Figure 4: Cross-correlation coefficient between 87Rb and 133Cs spins, both
estimated at the frequency where the cross-correlation spectrum is maximum
using a bandwidth of 100 Hz.

${\cal C}_{\text{Rb},\text{Cs}}=\frac{\int_{0}^{\infty}S_{\rm
meas}^{\text{Rb},\text{Cs}}(\nu)\phi_{\text{Rb}\text{Cs}}(\nu)d\nu}{\sqrt{\int_{0}^{\infty}S_{\rm
meas}^{\text{Rb},\text{Rb}}(\nu)\phi_{\text{Rb}}(\nu)d\nu\int_{0}^{\infty}S_{\rm
meas}^{\text{Cs},\text{Cs}}(\nu)\phi_{\text{Cs}}(\nu)d\nu}},$ (28)

where the kernel functions $\phi_{\text{Rb}}(\nu)$, $\phi_{\text{Cs}}(\nu)$,
and $\phi_{\text{Rb}\text{Cs}}(\nu)$ are centered around
$\nu_{0}^{\text{Rb}}$, $\nu_{0}^{\text{Cs}}$ and
$\nu_{0}^{\text{Rb}\text{Cs}}$, respectively, i.e. the frequency where the
corresponding spectrum is maximum (additionally, all three kernel functions
depend on the measurement time $T$).

In Fig. 4 we plot ${\cal C}_{\text{Rb},\text{Cs}}$ as a function of the
magnetic field for a measurement bandwidth of 100 Hz, using the experimentally
acquired spectra. It is seen that at low magnetic fields, the cross-
correlation between the two spin-species can be a significant fraction of the
measured spin-noise power, while it drops at larger magnetic fields.
Alternatively, if the center frequencies of the harmonic quadratures are
chosen in the region where the cross-spectrum is negative, the measured
correlation (for appropriate integration time) will correspondingly appear to
be negative. Overall, the sign and strength of the cross-correlation can be
adjusted with readily controllable experimental parameters, like the magnetic
field or the measurement bandwidth.

###  VI.1 Cross-correlation spin-noise power

Of particular interest is the equal-time ($\tau=0$) cross-correlation, i.e.
the total cross-correlation power. The power is related to the noise terms
that enter into the stochastic differential equations [118] describing the
time evolution of the observables. Previously, a debate about the value of
cross-correlation power has emerged in the literature. In particular, Dellis
and coworkers [99] used a 85Rb-87Rb spin ensemble and measured a non-zero
cross-correlation power, which increased at low magnetic fields. In constrast,
Roy and coworkers [100] used a 85Rb-133Cs spin ensemble and found the cross-
correlation power to be zero, irrespective of the magnetic field.

Figure 5: (a) Cross-correlation coefficient as a function of the magnetic
field for two different sampling rates. (b) Cross-correlation coefficient for
a specific magnetic field, $B=10~{}{\rm mG}$, as a function of the ratio
$R_{D}/R_{\rm se}$, quantifying the strength of transit-time-broadening
relative to spin-exchange relaxation. The inter-species cross-correlation is
suppressed if the linewidth is not dominated by spin-exchange relaxation.

We will here resolve the aforementioned debate, noting that the measurement
bandwidth has a subtle effect on the observed cross-correlation power.
Firstly, as noted previously, a zero imaginary part of the cross-spectrum and
a physically justifiable scaling of the single species spin-noise variance
indeed imply a zero cross-correlation power. This, however, corresponds to the
integration of the cross-correlation spectrum from frequency zero to infinity.
In a realistic experiment, all measurements are conducted within a finite
frequency range. Frequency components exceeding this range do not contribute
to the measured cross-correlation power. Consequently, if the cross-
correlation spectrum contains considerable power in the frequency range beyond
the measurable bandwidth, the detected cross-correlation power may indeed
appear to be non-zero. This is particularly the case when the high-frequency
tail of the spectrum extends to frequencies larger than the spin-exchange
rate. To demonstrate this, we calculate the 87Rb-133Cs cross-correlation
coefficient when acquiring the spin signals at a finite sampling rate. Each
recorded data point is modelled as the average of the corresponding signal
over the sampling duration: $\tilde{y}_{\alpha}={1\over{\Delta
t}}\int_{t}^{t+\Delta t}y_{\alpha}(t^{\prime})dt^{\prime}$, where
$\tilde{y}_{\alpha}$ and $y_{\alpha}$ represent respectively the sampled and
the underlying continuous-time signals of species $\alpha$, and $1/\Delta t$
denotes the sampling rate. We define the cross-correlation coefficient
$\tilde{\cal
C}_{\text{Rb},\text{Cs}}=\langle\tilde{y}_{\text{Rb}}\tilde{y}_{\text{Cs}}\rangle/\sqrt{\langle\tilde{y}^{2}_{\text{Rb}}\rangle\langle\tilde{y}^{2}_{\text{Cs}}\rangle}$,
and calculate it from the theoretical fits to the data employing Eq.(87).

The coefficient $\tilde{\cal C}_{\text{Rb},\text{Cs}}$ is shown in Fig. 5a as
a function of the magnetic field for two different sampling rates. It is seen
that we obtain a significantly smaller correlation for a large sampling rate,
as in this case a larger part of the negative high-frequency tail of the
cross-correlation spectrum is included. These observations explain the
positive correlations observed in [99], where the bandwidth of the measurement
was 50 kHz.

Regarding the experimental result in [100], we note that a crucial physical
behavior is the dependence of the correlation coefficient on the spin-exchange
relaxation rates $R_{\rm se}$, and all other relaxation rates. In Fig. 5b we
plot the correlation coefficient as a function of the ratio $R_{D}/R_{\rm
se}$, where $R_{D}$ is the relaxation due to atoms diffusing out of the probe
beam (see Appendix 2.A.f). It is seen that for a large $R_{D}$ the correlation
drops to zero. This is because for correlations to be observable, the magnetic
linewidths should be dominated by spin-exchange broadening. If other
broadening mechanisms prevail, the correlation effect will be suppressed.
Diffusion out of the probe beam, or the so-called transit-time broadening, is
indeed one such possibility, apparently responsible for the zero cross-
correlation power measured in [100]. The authors in [100] focused the probe
beam to a $50~{}{\rm\mu m}$ waist. Moreover, they had 3 times less buffer gas
pressure than in this experiment. This renders the transit time broadening
about 300 times larger than our case, amounting to $8.4\times 10^{5}~{}{\rm
s}^{-1}$. A subtle difference of this kind of broadening is that it equally
affects the nucleus and the electron, since it is just a Fourier broadening of
the signal’s limited-time observation. Hence, there is no slowing down factor
like in the other broadening mechanisms, and the aforementioned rate directly
appears in the measured linewidth. In fact, it can be seen from Fig. 2 of
[100] that the linewidth is about 100 kHz, even though the cesium and rubidium
densities reported in [100] are similar to our experiment, hence the spin-
exchange broadening should be about 4 kHz.

###  VI.2 Intra-species correlations in the context of mean field theory

Using the measurements on two district species, it is clearly demonstrated
that spin-exchange collisions result in unequal-time correlations among the
atoms involved in collisions. Among those, there are also inter-atomic
unequal-time correlations between atoms of the same species. It seems
remarkable that a mean field theory employing essentially a single-atom
description captures the intricate correlations emerging from spin-exchange
collisions among atoms of identical species.

To address this issue and underscore the internal consistency of the presented
theory, we proceed to demonstrate the reduction of the equations of motion
governing two colliding atoms engaging in spin-exchange interactions within
the same species. This reduction elegantly transforms the dynamics of the
interacting pair into the equation of motion characterizing a mean single
atom, which is effectively measured in an experiment. This elucidates the
connections between the microscopic behavior of individual collisions and the
macroscopic behavior described by mean field theory, thereby shedding light on
the overarching consistency of the theoretical framework.

Consider the observable vector:

$\mathcal{\hat{V}}(t)=\left(\hat{f}^{\mathfrak{a}}_{a,x},\hat{f}^{\mathfrak{a}}_{a,y},\hat{f}^{\mathfrak{a}}_{b,x},\hat{f}^{\mathfrak{a}}_{b,y},\hat{f}^{\mathfrak{b}}_{a,x},\hat{f}^{\mathfrak{b}}_{a,y},\hat{f}^{\mathfrak{b}}_{b,x},\hat{f}^{\mathfrak{b}}_{b,y}\right)^{\top},$
(29)

where now the indices $\mathfrak{a}$ and $\mathfrak{b}$ denote atoms of the
same species. All spin operators depend on time $t$. Using the methods
presented in Appendix B, the mean spin dynamics can be formulated as follows:

$\frac{d}{dt}\langle\mathcal{\hat{V}}(t)\rangle=\tilde{\mathcal{A}}\langle\mathcal{\hat{V}}(t)\rangle.$
(30)

where the matrix $\tilde{\mathcal{A}}$ encapsulates the linear dynamics akin
to the representation by matrix $A$ pertinent to the dual-species case.

In an experiment, collective variables
$\hat{f}^{\mathfrak{a}}+\hat{f}^{\mathfrak{b}}$ are measured. We thus define
the collective observable vector

$\begin{split}{\cal\hat{M}}(t)&=\mathcal{R}\mathcal{\hat{V}}(t)\\\
&=\left(\hat{f}^{\mathfrak{a}}_{a,x}+\hat{f}^{\mathfrak{b}}_{a,x},\hat{f}^{\mathfrak{a}}_{a,y}+\hat{f}^{\mathfrak{b}}_{a,y},\hat{f}^{\mathfrak{a}}_{b,x}+\hat{f}^{\mathfrak{b}}_{b,x},\hat{f}^{\mathfrak{a}}_{b,y}+\hat{f}^{\mathfrak{b}}_{b,y}\right)^{\top},\end{split}$
(31)

where $\mathcal{R}$ is the transformation matrix from the two-atom spin space
to the reduced collective-spin space:

$\mathcal{R}=\left(\begin{array}[]{cccccccc}1&0&0&0&1&0&0&0\\\
0&1&0&0&0&1&0&0\\\ 0&0&1&0&0&0&1&0\\\ 0&0&0&1&0&0&0&1\\\ \end{array}\right).$
(32)

Simple matrix algebra verifies that the structure of matrix
$\tilde{\mathcal{A}}$ is such that it fulfills the equality:

$\mathcal{R}\tilde{\mathcal{A}}\mathcal{\hat{V}}=\mathcal{R}\tilde{\mathcal{A}}\mathcal{R}^{+}\mathcal{\hat{M}},$
(33)

where $\mathcal{R}^{+}$ is the Moore-Penrose inverse 111Since $\mathcal{R}$
has linearly independent rows, the Moore-Penrose inverse is the right inverse
of $\mathcal{R}$ and:
$\mathcal{R}^{+}=\mathcal{R}^{{\dagger}}(\mathcal{R}\mathcal{R}^{{\dagger}})^{-1}$.
of $\mathcal{R}$. As a result:

$\frac{d}{dt}\langle{\cal\hat{M}}(t)\rangle=\mathcal{R}\frac{d}{dt}\langle\mathcal{\hat{V}}(t)\rangle=\mathcal{R}\tilde{\mathcal{A}}\mathcal{R}^{+}\langle\mathcal{\hat{M}}(t)\rangle=\mathcal{A}\langle\mathcal{\hat{M}}(t)\rangle,$
(34)

where $\mathcal{A}=\mathcal{R}\tilde{\mathcal{A}}\mathcal{R}^{+}$.

The matrix $\mathcal{A}$, formulated to describe the evolution of the
collective ensemble spin, corresponds precisely to the matrix derived through
the utilization of the equation governing the density matrix of a single atom
(see for example Eq. 45 of [120]). This congruence implies that the single-
atom equation implicitly encompasses the influence of dynamic (non-equal time)
inter-atomic correlations emanating from collisions involving atoms of the
identical species.

###  VI.3 Character of correlations: entanglement

As mentioned in the introductory Sec. IIB, a pertinent question is whether the
observed multi-time correlations are quantum or classical in nature. In other
words, whether the observed correlations can be described using a classical
probability model or are they non-classical in the sense that they cannot be
prepared with classical operations. Quantifying quantum or quantum/classical
correlations [12] can be rather challenging and goes beyond the scope of this
work. We briefly note that quantum correlations comprise an entire family of
relationships [12], such as entanglement and quantum discord, which could
potentially serve as benchmarks for the observed correlations.

We here make only some exploratory comments and leave a more detailed
discussion for future work. In [68] it was shown that spin-exchange collisions
of the kind encountered in our work can create entanglement between two
(partially) polarized spin ensembles that can be sustained for meaningful
timescales. Similarly, it was shown in [15, 14] that spin-exchange collisions
under appropriate conditions create quantum correlations and can be harnessed
to transfer the quantum state from one spin species to another spin species.

On the other hand, as explained in the previous subsection, the correlations
measured here are consistent with the theory describing spin dynamics from the
single-atom perspective (mean-field), hence such correlations cannot be
anything but classical. It is possible, however, that the character of
correlations could have a parametric dependence on the spin-polarization of
the vapor, like the Werner state discussed in the introduction. We here use a
simple two-atom toy model to illustrate that this seems might indeed be the
case. Considering atoms spin-polarized along $\mathbf{\hat{x}}$, a magnetic
field along $\mathbf{\hat{z}}$ creates a spin component along
$\mathbf{\hat{y}}$. What we will show is that we expect the fluctuations of
the total spin component $\hat{f}_{y}$ for 133Cs and 87Rb to be quantum
correlated.

We take as initial states for the 133Cs and 87Rb atoms the spin-temperature
states $\rho_{1}=e^{\beta\hat{f}_{x}}/\rm Tr\\{\beta\hat{f}_{x}\\}$ and
$\rho_{2}=e^{-\beta\hat{f}_{x}}/\rm Tr\\{-\beta\hat{f}_{x}\\}$, i.e. having
the same but opposite spin polarization
$\braket{\hat{s}_{x}}={1/2}\tanh{({\beta/2})}$, with $\beta$ being the spin
temperature [120]. The combined initial state is
$\rho=\rho_{1}\otimes\rho_{2}$. Note that for each atom the operator
$\hat{f}_{x}$ and the respective density matrix has a matrix representation of
different dimension. We then evolve $\rho$ by the Hamiltonian $\hat{H}$, i.e.
we calculate $\rho^{\prime}=e^{-i\hat{H}t}\rho e^{i\hat{H}t}$, where
$\hat{H}=\hat{h}_{\rm Cs}\otimes\mathbb{1}+\mathbb{1}\otimes\hat{h}_{\rm Rb}$,
with $\hat{h}_{\rm Cs}$ and $\hat{h}_{\rm Rb}$ being the 133Cs and 87Rb Breit-
Rabi Hamiltonians, respectively. We use a small field of $10~{}{\rm\mu G}$ and
a precession time of $10~{}{\rm\mu s}$. We then apply on $\rho^{\prime}$ the
spin-exchange operator
$\hat{P}_{e}=(1/2)\mathbb{1}+2\mathbf{\hat{s}}_{1}\cdot\mathbf{\hat{s}}_{2}$,
and calculate the resulting negativity 222To quantify the entanglement in the
multi-dimensional Hilbert space of the two alkali species we use the
negativity criterion [126] (entanglement measure) as demonstrated in [68]. of
the state $\hat{P}_{e}\rho^{\prime}\hat{P}_{e}^{\dagger}$, like in [68]. The
result is shown in Fig. 6.

Figure 6: Negativity quantifying the entanglement of the 87Rb-133Cs combined
spin state produced by (i) spin-polarizing the vapors along the $x$-axis with
equal and opposite polarizations, (ii) spin precession in a transverse
magnetic field, and (iii) cross-spin exchange collisions.

It is seen that for low polarizations pertinent to the noise measurement of
this work the two atoms are not entangled. However, there appears to be a
threshold spin-polarization over which the two atoms gradually become strongly
entangled.

##  VII  Conclusions

We have studied experimentally and theoretically spin correlations that
spontaneously build up in a dual-species alkali-metal vapor, here composed of
${}^{87}{\rm Rb}$ and ${}^{133}{\rm Cs}$. We have identified and categorized
various types of spin-correlations that can be probed in a hot atomic vapor
and we have elucidated their behavior. The combined action of interspecies and
intraspecies spin-exchange collisions leads to positive equal-time spin-
correlations which are enhanced at low magnetic fields and suppressed at high
fields due to an interplay of their unique spectral distribution and the
measurement bandwidth.

The nature of these correlations has been discussed, anticipating that similar
correlations in spin-polarized vapors are expected to be genuinely quantum,
i.e. to reflect interspecies spin entanglement.

The use of two species rather than one, and the study of both auto- and cross-
correlations, helps to unravel the complexity of spin dynamics in alkali-metal
vapors, so far treated mostly as consisting of uncorrelated atoms, and can
have significant repercussions for the field of quantum metrology with hot
atomic vapors.

###### Acknowledgements.

K.M. acknowledges support from Grant FJC2021047840-I funded by
MCIN/AEI/10.13039/501100011033; the European Union “NextGenerationEU/PRTR”;
Greece and the European Union [European Social Fund (ESF)] through the
Operational Programme “Human Resources Development, Education and Lifelong
Learning” in the context of the project “Strengthening Human Resources
Research Potential via Doctorate Research” (Grant No. MIS-5000432),
implemented by the State Scholarships Foundation (IKY). GV acknowledges
funding from EU QuantERA Project PACE-IN (GSRT Grant No. T11EPA4-00015) and
from the Hellenic Foundation for Research and Innovation (HFRI) under the HFRI
agreement no. HFRI-00768 (project QCAT). F.V., A.M., G.M., M.S., G.P.T. and
I.K.K. acknowledge the cofinancing of this research by the European Union and
Greek national funds through the Operational Program Crete 2020-2024, under
the call “Partnerships of Companies with Institutions for Research and
Transfer of Knowledge in the Thematic Priorities of RIS3Crete", with project
title “Analyzing urban dynamics through monitoring the city magnetic
environment” (project KPHP1 - 0029067). M.W.M acknowledges the European
Commission project OPMMEG (101099379), Spanish Ministry of Science MCIN with
funding from European Union NextGenerationEU (PRTR-C17.I1) and by Generalitat
de Catalunya “Severo Ochoa” Center of Excellence CEX2019-000910-S; projects
SAPONARIA (PID2021-123813NB-I00) and MARICHAS (PID2021-126059OA-I00) funded by
MCIN/ AEI /10.13039/501100011033/ FEDER, EU; Generalitat de Catalunya through
the CERCA program; Agència de Gestió d’Ajuts Universitaris i de Recerca Grant
2021-SGR-01453; Fundació Privada Cellex; Fundació Mir-Puig;. A.K. and M.L.
were supported by the Hellenic Foundation for Research and Innovation
(H.F.R.I.) under the “First Call for H.F.R.I. Research Projects to support
Faculty members and Researchers and the procurement of high-cost research
equipment grant” (Project ID HFRI-FM17-1034). M.S. acknowledges financial
support from the Spanish Agencia Estatal de Investigacíon, Grant No.
RYC2021-032032-I, PID2019107609GB-I00, and the Catalan Government for the
project QuantumCAT 001-P-001644, co-financed by the European Regional
Development Fund (FEDER).

##  Appendix A Detailed Theoretical Derivation of Cross- and Auto-Correlation
Spectra

We here recapitulate the basic physics of our theoretical description, namely
atom-light coupling, and atomic evolution due to coherent Hamiltonian dynamics
and relaxation effects dominated by atom-atom collisions.

###  A.1 Atom-light coupling

The interaction of light with atomic spin ensembles has been described in
detail elsewhere [108, 109, 110, 111]. For near-resonant monochromatic light
of non-saturating intensity, the coherent atom-light interaction (i.e. the
interaction that leads to forward scattered light in the spatial mode of the
probe beam) is described with a polarizability Hamiltonian that couples the
magnetic sublevels of the atomic ground state to the polarization modes of the
light. The resulting atomic polarizability is a rank-2 tensor operator that
can be decomposed into three irreducible components. Of those, the tensor
(rank-2) polarizability is negligible for the conditions of our experiment
(pressure broadening $\sim 10$ $\mathrm{GHz}$ and detuning $\sim 100$
$\mathrm{GHz}$). Most relevant for our experiment is the vector polarizability
(rank-1) describing paramagnetic Faraday rotation and reading [122]:

$\hat{H}_{\text{int}}=\sum_{i=1}^{N_{\mathrm{at}}}\hat{{\cal
S}}_{3}(\mathbf{r}_{i},t)\beta(\mathbf{r}_{i})\left[g_{a}\hat{f}_{a,x}^{(i)}(t)-g_{b}\hat{f}_{b,x}^{(i)}(t)\right],$
(35)

where the summation is performed over all atoms probed by the laser beam,
${\cal
S}_{3}(\mathbf{r},t)=i\left(\hat{a}^{{\dagger}}_{y}\hat{a}_{z}-\hat{a}^{{\dagger}}_{z}\hat{a}_{y}\right)/2$
is the Stokes light-operator quantifying the photon-flux imbalance of the
left- and right- circular polarization modes (and in the most general case is
a function of space and time), $\mathbf{r}_{i}$ is the location of the $i$-th
atom, $\hat{f}_{a,x}^{(i)}(t)$ and $\hat{f}_{b,x}^{(i)}(t)$ are dimensionless
spin components of the $i$-th atom along the probe laser direction in the
$a=I+1/2$ and $b=I-1/2$ hyperfine levels of the ground state, with $I$ being
the nuclear spin quantum number ($I=7/2$ for 133Cs and $I=3/2$ for
${}^{87}\text{Rb}$). The parameter $\beta(\mathbf{r}_{i})$ characterizes the
local field intensity associated with the spatial mode of the probe beam (see
below).

For a D2 optical transition linewidth dominated by pressure broadening the
coupling constants entering Eq.(35) are approximated by [67]:

$g_{\alpha}\approx\frac{1}{2I+1}\frac{cr_{\text{e}}f_{\text{osc}}}{A_{\text{eff}}}\frac{\nu_{\ell}-\nu_{\alpha}}{(\nu_{\ell}-\nu_{\alpha})^{2}+(\Delta\nu/2)^{2}},$
(36)

where $r_{\text{e}}\approx 2.82\times 10^{-15}$ $\mathrm{m}$ is the classical
electron radius, $f_{\text{osc}}$ is the oscillator strength of the optical
transition, $c$ is the speed of light, $\Delta\nu$ is the pressure-broadened
optical linewidth (FWHM), $\nu_{\ell}$ is the frequency of light, and
$\nu_{\alpha}$ with $\alpha\in\\{a,b\\}$ is the resonance frequency of the
corresponding ground-state hyperfine level. If $A_{\text{eff}}$ labels the
effective area of the beam, it is
$\beta(\mathbf{r}_{i})/A_{\text{eff}}=I(\mathbf{r}_{i})/\int
I(\mathbf{r}_{i})dydz$, where $I(\mathbf{r}_{i})$ is the light intensity at
the coordinate $\mathbf{r}_{i}$ and $\int I(\mathbf{r}_{i})dydz$ is the total
power of the beam. For a $\text{TEM}_{00}$ Gaussian beam it is
$\beta(\mathbf{r}_{i})/A_{\text{eff}}=e^{-2|\mathbf{r}_{i}|^{2}/w(x)^{2}}/(\pi
w(x)^{2}/2)$, where $w(x)$ is the $x$-dependent Gaussian beam width.

For a linearly polarized probe in the $y$-direction, neglecting the time it
takes light to propagate through the spin ensemble [112], the solution of the
Heisenberg equation of motion for the light operator (in the case of small
rotation angles) yields [113]:

$\hat{{\cal S}}_{2}^{(\text{out})}(t)\approx\hat{{\cal
S}}_{2}^{(\text{in})}(t)+\frac{1}{2}\Phi\bar{g}\mathcal{\hat{F}}(t),$ (37)

where $\Phi$ is the photon-flux (photons per time). The out (in) superscripts
denote the operator after (before) interacting with the atomic sample,
$\bar{g}=(|g_{a}|+|g_{b}|)/2$ is the mean coupling constant to the two ground-
state hyperfine manifolds and $\mathcal{\hat{F}}(t)$ is the measured
collective spin defined as:

$\displaystyle\mathcal{\hat{F}}(t)$
$\displaystyle\equiv\sum_{i=1}^{N_{\mathrm{at}}}\beta(\mathbf{r}_{i})\left[\frac{g_{a}}{\bar{g}}\hat{f}_{a,x}^{(i)}(t)-\frac{g_{b}}{\bar{g}}\hat{f}_{b,x}^{(i)}(t)\right]$
(38) $\displaystyle=\int
d\mathbf{r}\beta(\mathbf{r})\left[\frac{g_{a}}{\bar{g}}\hat{\tilde{f}}_{a,x}(\mathbf{r},t)-\frac{g_{b}}{\bar{g}}\hat{\tilde{f}}_{b,x}(\mathbf{r},t)\right].$
(39)

In the last equation the spin-density operator $\hat{\tilde{\mathbf{f}}}$ was
introduced, which for the $\alpha$ hyperfine manifold is defined as:

$\hat{\tilde{\mathbf{f}}}_{\alpha}(\mathbf{r},t)=\sum_{i=1}^{N_{\mathrm{at}}}\hat{\mathbf{f}}^{(i)}_{\alpha}(t)\delta(\mathbf{r}-\mathbf{r}_{i}).$
(40)

###  A.2 Atomic evolution

As will become apparent in the following, the measured fluctuation spectrum
depends on the dynamics of the collective atomic spin in the ground state,
which are determined by the atomic Hamiltonian, atomic collisions, atomic
thermal motion, and the interaction with light. We next treat such atomic
dynamics.

####  A.2.1 Hamiltonian

The electronic-ground-state Hamiltonian for an alkali-metal atom in the
presence of a magnetic field $\mathbf{B}$ is:

$\hat{H}_{g}=A_{\text{hf}}\hat{\mathbf{I}}\cdot\hat{\mathbf{s}}+g_{s}\mu_{B}\hat{\mathbf{s}}\cdot\mathbf{B},$
(41)

where $A_{\text{hf}}$ is the isotropic hyperfine coupling coefficient,
$\hat{\mathbf{I}}$ and $\hat{\mathbf{s}}$ are the nuclear and electron spin
operators, respectively, $g_{s}\approx 2$ is the electron’s g-factor, and
$\mu_{B}\approx 9.27\times 10^{-24}$ J/T is the Bohr magneton.

####  A.2.2 Spin-exchange collisions

Spin-exchange collisions are dominant in the physics of atomic spin state
evolution. The spin-exchange interaction originates from the difference
between the lowest singlet and triplet potential energy of the molecular
system formed during a pairwise collision. The established formalism captures
the effect of spin-exchange collisions with a mean-field density matrix, which
for the $i$-th atomic species reads [20]:

$\begin{split}d\rho_{\beta}/dt&=\sum_{\gamma}R_{\rm
se}^{\beta,\gamma}\big{(}\phi_{\beta}\left(1+4\langle\hat{\mathbf{s}}_{\gamma}\rangle\cdot\hat{\mathbf{s}}_{\beta}\big{)}-\rho_{\beta}\right)\\\
&-\sum_{\gamma}2iR_{\rm
se}^{\beta,\gamma}\kappa_{\beta\gamma}\left[\langle\mathbf{s}_{\gamma}\rangle\cdot\mathbf{\hat{s}}_{\beta},\rho_{\beta}\right],\end{split}$
(42)

where $\beta,\gamma\in\\{\text{Rb},\text{Cs}\\}$,
$\phi=\rho/4+\hat{\mathbf{s}}\cdot\rho\hat{\mathbf{s}}$ is the part of density
matrix without electron polarization,
$\langle\mathbf{\hat{s}}_{\gamma}\rangle=\text{Tr}\left[\mathbf{\hat{s}}\rho_{\gamma}\right]$,
$\kappa_{\beta\gamma}$ is a dimensionless parameter characteristic for the
mean-field interaction between the $\beta$ and $\gamma$ species, and $R_{\rm
se}^{\beta,\gamma}\approx n_{\gamma}\upsilon_{\rm
th}^{\beta,\gamma}\sigma_{\rm se}^{\beta,\gamma}$ is the spin-exchange rate,
given in terms of the atomic density $n_{\gamma}$ of species $\gamma$, the
relative thermal velocity between the colliding partners $\upsilon_{\rm
th}^{\beta,\gamma}$ and the spin-exchange cross-section $\sigma_{\rm
se}^{\beta,\gamma}$. In Eq.(42) the summation runs over all the different
species present in the vapor, including the species $\beta$. As we are dealing
with unpolarized vapors, the effect of the frequency shift proportional to the
parameter $\kappa_{\beta\gamma}$ will be ignored in the following.

We note that Eq. (42) is a mean field theory. The spin-exchange interactions
that an atom can experience are replaced with an effective average interaction
and the spin-degrees of freedom of the different atoms are traced out to get a
reduced mean density matrix representing the evolution of the ensemble.
Essentially, Eq. (42) provides for coarse-grained dynamics over a length scale
much larger than the mean-free path. Even in this case, Eq. (42) can capture
correlations that can be generated by the dynamics, at least at the level of
classical correlations (see Sec. VI.3).

####  A.2.3 S-damping collisions

Binary collisions between alkali atoms or between alkali-metal atoms and
buffer gas (without spin) lead to S-damping [120]:

$d\rho/dt=R\left(\phi-\rho\right),$ (43)

where the part of the density matrix with electron polarization is destroyed
while the purely nuclear polarization remains unaffected. The S-damping rate,
$R$, is orders of magnitude smaller than the spin-exchange rate,
$R_{\rm{se}}$.

####  A.2.4 Relaxation due to optical fields

In [120], it is shown that for the experimentally relevant case of fast
quenching (when excited atoms are much more likely to be quenched than to
radiate a photon) and excited state $J$-damping rapid with respect to the
hyperfine interaction, the net evolution of the single-atom density matrix due
to optical (depopulation and repopulation) pumping by linearly polarized light
is of the form (43), with $R$ given by $\Phi\sigma_{\text{op}}$, the mean
pumping rate per (unpolarized) alkali-metal atom given in terms of the
absorption cross-section $\sigma_{\text{op}}$ and the photon flux $\Phi$. In
principle, the effect of light on the atomic-spin evolution can be made
negligible, either by using low photon flux (low power or large beam area) or
by choosing large detuning. However, as seen from Eq. (48), this results in a
reduction of the measured spin-noise to photon shot noise ratio, leading to an
increase in the uncertainty of estimation of the spin-noise spectrum 333In a
noisy process, the uncertainty (standard deviation) in the periodogram
estimate at a frequency is equal to the its expectation value at this
frequency [127]. for a given number of repetitions. In practice, the laser
power is chosen so as to optimize the measurement’s signal-to-noise ratio. For
our case in particular, the laser power should be such that spin relaxation by
the optical fields is not dominant, i.e. spin dynamics should be dominated by
spin-exchange collisions.

####  A.2.5 Combined description

Overall, the density matrix evolution for species $\beta$ reads

$\displaystyle d\rho_{\beta}/dt$
$\displaystyle=A_{\text{hf},\beta}\hat{\mathbf{I}}_{\beta}\cdot\hat{\mathbf{s}}_{\beta}+g_{s}\mu_{B}\hat{\mathbf{s}}_{\beta}\cdot\mathbf{B}+R\left(\phi_{\beta}-\rho_{\beta}\right)$
$\displaystyle+\sum_{\gamma}R_{\rm
se}^{\beta,\gamma}\left\\{\phi_{\beta}\left(1+4\langle\hat{\mathbf{s}}_{\gamma}\rangle\cdot\hat{\mathbf{s}}_{\beta}\right)-\rho_{\beta}\right\\},$
(44)

where the relaxation rate $R$ includes all the relaxation processes, other
than the spin-exchange collisions, which destroy electron polarization without
affecting the nucleus.

####  A.2.6 Atomic diffusion

In a cell with buffer gas, as in our experiment, the atomic thermal motion
becomes diffusive via velocity-changing collisions. The effect of diffusion on
spin-noise measurements has been studied in various works. In [86] the
correlation function of diffusion-induced noise in unconfined systems is
derived. A treatment based on the Bloch-Heisenberg-Langevin formalism is
developed in [87], also considering boundary conditions on the cell walls. In
our experiment the probe beam diameter is much smaller than the cell
dimensions, thus the results of [86] are more relevant. There, it is shown
that for diffusion through a $\text{TEM}_{00}$ Gaussian and well-collimated
beam (Rayleigh range much larger than the cell length) the spin-noise spectrum
of an atomic ensemble with transverse spin-relaxation rate $1/T_{2}$ and spin
precession frequency $\nu_{\text{L}}$ is given by

$S(\nu)\propto\frac{w^{2}}{2D}\Re\left[e^{s}\Gamma(s)\right],\text{
}s=\frac{w^{2}}{4D}\left[1/T_{2}+2\pi i(\nu-\nu_{\text{L}})\right],$ (45)

where $D$ is the atomic diffusion constant, $w$ is the waist of the beam and
$\Gamma(s)=\int_{s}^{\infty}(e^{-x}/x)\,dx$ is the incomplete gamma function.
When $w^{2}/4DT_{2}\gg 1$ it is $e^{s}\Gamma(s)\approx(1-1/s)/s$, and the
spectrum takes the form:

$S(\nu)\sim
2\frac{1/T_{2}^{\prime}}{(1/T_{2}^{\prime})^{2}+4\pi^{2}(\nu-\nu_{\text{L}})^{2}}+\mathcal{O}\left(\frac{4DT_{2}}{w^{2}}\right),$
(46)

where $1/T_{2}^{\prime}=1/T_{2}+R_{D}$, with $R_{D}\sim 4D/w^{2}$. Therefore,
when the diffusion time across the beam is much longer than the spin coherence
time, the effect of diffusion on the spin-noise spectrum can be captured to a
good approximation by introducing an additional relaxation term with rate
$R_{D}$.

###  A.3 Correlations

From Eqs. (7) and (37) it follows that the cross-correlation measured in the
experiment is:

$\displaystyle C_{\text{Rb},\text{Cs}}(\tau)$
$\displaystyle=\frac{q_{e}^{2}\eta_{\text{Rb}}\eta_{\text{Cs}}\Phi_{\text{Rb}}\Phi_{\text{Cs}}\bar{g}_{\text{Rb}}\bar{g}_{\text{Cs}}}{2}$
$\displaystyle\times\langle\mathcal{\hat{F}}_{\text{Rb}}(t+\tau)\mathcal{\hat{F}}_{\text{Cs}}(t)+\mathcal{\hat{F}}_{\text{Cs}}(t)\mathcal{\hat{F}}_{\text{Rb}}(t+\tau)\rangle,$
(47)

where we have taken into account that the polarization properties of the two
beams before the interaction with the atoms are uncorrelated e.g.
$\langle{\cal S}^{(\rm{in})}_{2,\text{Rb}}(t){\cal
S}^{(\rm{in})}_{2,\text{Cs}}(t^{\prime})\rangle=0$. The single-channel auto-
correlation function is:

$C_{\beta,\beta}(\tau)=2q_{e}^{2}\eta_{\beta}\Phi_{\beta}\Big{[}\delta(\tau)+\frac{1}{2}\eta_{\beta}\Phi_{\beta}\bar{g}^{2}_{\beta}\langle\mathcal{\hat{F}}_{\beta}(t+\tau)\mathcal{\hat{F}}_{\beta}(t)\rangle\Big{]},$
(48)

where $\beta\in\\{\text{Rb},\text{Cs}\\}$. In deriving Eq. (48) we assumed
perfectly coherent light before the interaction with the atoms, i.e.
$\langle{\cal S}^{(\rm{in})}_{2,\beta}(t){\cal
S}^{(\rm{in})}_{2,\beta}(t^{\prime})\rangle=(\Phi_{\beta}/2)\delta(t-t^{\prime})$.
The detector inefficiency was modelled by a perfect detector preceded by a
beam splitter of power transmissivity $\eta_{\beta}$ [9]. The first term in
Eq.(48) describes photon shot-noise, while the second contains information
about the collective spin correlations. We remark that although the photon
shot-noise does not explicitly appear in the inter-species correlator of
Eq.(47), it does contribute to the uncertainty of its estimation. In Eqs.(47)
and (48) it was assumed that there is no correlation between the light and
spin operators, implying that the spin variables remain unaffected by the
polarization fluctuations of the input field. This is true for our experiment
probing unpolarized spin ensembles with light (refer to section D). However,
this assumption breaks down when dealing with non-zero spin polarization.

As long as the optical rotation remains small, Eqs. (37)-(48) remain valid
even for optically thick ensembles with significant light absorption. In such
cases, the variable $\Phi$ represents the photon flux at the ensemble output,
which can be significantly reduced compared to the photon flux at the ensemble
input [124].

Combing Eqs. 39 and 47, the measured correlation can be written as:

$C_{\text{Rb},\text{Cs}}(\tau)=\mathcal{K}\int\int
d\mathbf{r}d\mathbf{r}^{\prime}\beta^{\text{Rb}}(\mathbf{r})\beta^{\text{Cs}}(\mathbf{r}^{\prime})\times\sum_{\alpha,\alpha^{\prime}}g_{\alpha}^{\text{Rb}}g_{\alpha^{\prime}}^{\text{Cs}}\left[\langle\hat{\tilde{f}}_{\alpha,x}^{\text{Rb}}(\mathbf{r},t+\tau)\hat{\tilde{f}}_{\alpha^{\prime},x}^{\text{Cs}}(\mathbf{r}^{\prime},t)\rangle+\langle\hat{\tilde{f}}_{\alpha^{\prime},x}^{\text{Cs}}(\mathbf{r}^{\prime},t)\hat{\tilde{f}}_{\alpha,x}^{\text{Rb}}(\mathbf{r},t+\tau)\rangle\right],$
(49)

where $\mathcal{K}$ is an overall scaling factor. A similar equation can be
written for the autocorrelation function.

###  A.4 Collective spin operators

The dynamical evolution of spin correlations can be found from the quantum
regression theorem (QRT) [9], which states that the spin-correlation functions
follow the same equations of motion as the mean spins. Consequently, to
explore the dynamic evolution of the collective spin correlations, it suffices
to determine the mean dynamics of the corresponding spins. In this regard, the
mean field equation (Eq. (44)) can be utilized. We note that the validity of
Eq. (44) to describe the dynamics of the collective spin has been confirmed in
numerous experiments. Remarkably, the QRT allows us to move deeper into the
understanding of spin dynamics and address inter-atomic correlations. This is
because the QRT connects the mean dynamics of the collective spin to the
collective correlations, and the latter include a multitude of inter-atomic
correlation terms.

In the following, we derive a simplified version of Eq. 49 for the measured
cross-correlation. We adopt the approach presented in [14], involving the
definition of coarse-grained continuous local-symmetric spin-density
operators. These coarse-grained operators are obtained through the spatial
convolution of the spin-density operators with a window function, which
remains non-zero over a specific volume (e.g. a Heaviside function). This
volume must be adequately large to encompass a significant number of atoms
while being sufficiently small to ensure a uniform interaction of all atoms
with the probe beam.

As discussed in [14], the inherent stochastic nature of collision parameters
and random pairings of colliding atoms introduces noise. However, remarkably,
these noise terms do not significantly influence the temporal unfolding of
correlations, affecting only the zero-time correlation. The coarse-grained
spin-density operators evolve as a single entity within the chosen coarse-
grained volume, and their mean values follow the mean-field equation over
time.

In our experimental setup, characterized by a buffer gas pressure of 330 Torr
and a large beam diameter, the coarse-grain can be selected to be adequately
large; this ensures that atomic motion does not result in significant
correlations between different coarse-grains. Correlations resulting from
collisions are confined solely within the same coarse-grained region. As
explained in [14], this leads to correlations in the coarse-grained spin-
density operators that are proportional to the coarse-grained delta function.

Overall, disregarding the minor variation in the effect of spin-relaxation
caused by the probe beam along the atomic ensemble, the measured correlation
in the experiment can be expressed as (see Eq. 49):

$\displaystyle
C_{\text{Rb},\text{Cs}}(\tau)\propto\sum_{\alpha,\alpha^{\prime}}g_{\alpha}^{\text{Rb}}g_{\alpha^{\prime}}^{\text{Cs}}\Big{[}$
$\displaystyle\langle\hat{\breve{f}}_{\alpha,x}^{\text{Rb}}(t+\tau)\hat{\breve{f}}_{\alpha^{\prime},x}^{\text{Cs}}(t)\rangle+$
$\displaystyle\langle\hat{\breve{f}}_{\alpha^{\prime},x}^{\text{Cs}}(t)\hat{\breve{f}}_{\alpha,x}^{\text{Rb}}(t+\tau)\rangle\Big{]}.$
(50)

In the above equation $\hat{\breve{f}}^{\beta}_{\alpha,x}$ denotes the coarse-
grained, local-symmetric spin-density operator and is related to the single-
atom spin $\hat{f}$ by

$\langle\hat{\breve{f}}^{\beta}_{\alpha,x}\rangle=n_{\beta}\langle\hat{f}^{\beta}_{\alpha,x}\rangle$
(51)

where $n_{\beta}$ is the atomic density of the species $\beta$. The equations
of motion for the spin-density operators can be accordingly found from Eqs.
(51) and (44).

##  Appendix B Mean spin dynamics

The derivation of the linearized dynamics for the mean spin components
($\langle\hat{f}^{\beta}_{\alpha,x}\rangle$ or
$\langle\hat{\breve{f}}^{\beta}_{\alpha,x}\rangle$) is best approached by
adopting the method introduced in [49] (see also the Supplementary material in
[101]), where spin operators are expressed as spherical tensor operators in
the coupled ($F,m_{F}$) basis. For transverse spin, it is sufficient to
analyze the dynamics of the $\pm 1$ components of the rank-1 spherical tensor
(with the zero-component defined by the direction of the applied DC magnetic
field). We focus on studying the noise correlation occurring at Zeeman
frequencies. Hyperfine coherences are not measured, and their contribution to
Zeeman dynamics can be neglected in the zero polarization limit.

Here, for the reader’s convenience, we provide a concise overview of the
contribution of spin-exchange interactions between different alkali species in
the dynamics of the mean spins. We refer to [49] and [101] for more details.

The equation of motion for the density matrix of species $\beta$ due to spin-
exchange collision with species $\gamma$ can be written in the form:

$\displaystyle\frac{1}{R_{\text{se}}^{\beta,\gamma}}\frac{d\rho_{\beta}}{dt}$
$\displaystyle=\sum_{m}\left\\{\sqrt{\frac{[I_{\gamma}]}{[I_{\beta}]}}\langle
T_{001m}^{\gamma{\dagger}}\rangle-\langle
T_{001m}^{\beta{\dagger}}\rangle\right\\}T_{001m}^{\beta}$
$\displaystyle+\sum_{\begin{subarray}{c}\Lambda\mu m\\\ \Lambda\neq
0\end{subarray}}\sqrt{2[I_{\gamma}]}\langle T_{\Lambda\mu
00}^{\beta{\dagger}}\rangle\langle T_{001m}^{\gamma{\dagger}}\rangle
T_{\Lambda\mu 1m}^{\beta}$
$\displaystyle-\rho_{\beta}+\sum_{\Lambda\mu}T_{\Lambda\mu 00}^{\beta}\langle
T_{\Lambda\mu 00}^{\beta{\dagger}}\rangle+\sum_{m}\langle
T_{001m}^{\beta{\dagger}}\rangle T_{001m},$ (52)

where $[I]=2I+1$, and $T_{\Lambda\mu lm}\equiv T_{\Lambda\mu}(II)\otimes
T_{lm}(SS)$ is the spherical tensor spin operator in the uncoupled basis, with
$T_{LM}(KK^{\prime})$ being the spherical tensor in the angular momentum basis
with quantum numbers $K$ and $K^{\prime}$.

The second line in Eq. 52 has only terms that are proportional to the product
of the electron polarization of one species times the nuclear polarization
($\Lambda\neq 0$) of the other species. These terms introduce non-linearity in
the atomic evolution since each polarization (electron or nuclear) depends on
the density matrix. However, in the case of thermal-unpolarized atomic
ensembles, these terms are second order with respect to the small quantity of
polarization and can be neglected; consequently, linearized dynamics can be
considered for such systems.

Multiplying both sides of Eq. 52 with $T_{1M}^{\beta{\dagger}}(FF)$ and taking
the trace, the equation of motion for $\langle
T_{1M}^{\beta{\dagger}}(FF)\rangle$ can be obtained. The linearized equation
contains solely terms of the form $\langle
T_{1M}^{\beta{\dagger}}(FF)T_{001m}\rangle$ or $\langle
T_{1M}^{\beta{\dagger}}(FF)T_{\Lambda\mu 00}\rangle$, which can be found by
expressing $T_{001m}$ and $T_{\Lambda\mu 00}$ in the coupled basis (see Eqs.
63 and 64 in [49]), and using the orthogonality property for spherical
operators, $\langle
T^{{\dagger}}_{LM}(FF^{\prime})T_{lm}(ff^{\prime})\rangle=\delta_{Ll}\delta_{Mm}\delta_{Ff}\delta_{F^{\prime}f^{\prime}}$.
By also employing the definition for the Hermitian conjugate,
$T^{{\dagger}}_{LM}(FF^{\prime})=(-1)^{F-F^{\prime}+M}T_{L-M}(FF^{\prime})$,
all equations can be readily transformed into equations of motion for the
operators $T_{LM}(FF)$.

The contribution of other relaxation mechanisms to the dynamics of mean spins
can be determined using similar calculations. Finally, the dynamics due to
magnetic field are most easily determined from the Heisenberg equation of
motion, while the hyperfine interaction does not affect the Zeeman spherical
operators $T_{LM}(FF)$. Equations of motion for the spin density operators
$\breve{T}_{LM}(FF)$ can be obtained straightforwardly by utilizing Eqs. (51)
and (52). This amounts to modifying the rate in the first term of the sum
appearing in the first line of the RHS in Eq. 52, replacing
$R_{\text{se}}^{\beta,\beta^{\prime}}\propto n_{\beta^{\prime}}$ with
$R_{\text{se}}^{\beta^{\prime},\beta}\propto n_{\beta}$, i.e. interchanging
the spin-exchange rates in the terms of the equations of motion that are
responsible for the mixing between the spin operators of the two species.

Overall, the linearized equations of motion for the spin density operators of
the two species can be compactly written as:

$\frac{d\mathbf{T}}{dt}=\begin{pmatrix}\tilde{A}&\mathbf{0}_{4}\\\
\mathbf{0}_{4}&\tilde{A}^{*}\end{pmatrix}\mathbf{T},$ (53)

where $\mathbf{0}_{4}$ is the $4\times 4$ zero matrix and
$\mathbf{T}=[\mathbf{T}_{11}\phantom{\small{a}}\mathbf{T}_{1-1}]^{\top}$, with
$\mathbf{T}_{11}=[T_{11}^{\text{Rb}}(aa)\phantom{\small{a}}T_{11}^{\text{Rb}}(bb)\phantom{\small{a}}T_{11}^{\text{Cs}}(a^{\prime}a^{\prime})\phantom{\small{a}}T_{11}^{\text{Cs}}(b^{\prime}b^{\prime})]$
and similarly for $\mathbf{T}_{1-1}$. We remind the reader that
$a=I_{\text{Rb}}+1/2$ and $b=I_{\text{Rb}}-1/2$, while
$a^{\prime}=I_{\text{Cs}}+1/2$ and $b^{\prime}=I_{\text{Cs}}-1/2$. The
$4\times 4$ matrix $\tilde{A}$ is:

$\displaystyle\tilde{A}$
$\displaystyle=\begin{pmatrix}R_{\text{se}}^{\text{Rb},\text{Rb}}\Phi_{\rm{se}}(I_{\text{Rb}})&\mathbf{0}_{2}\\\
\mathbf{0}_{2}&R_{\text{se}}^{\text{Cs},\text{Cs}}\Phi_{\rm{se}}(I_{\text{Cs}})\end{pmatrix}$
$\displaystyle+\begin{pmatrix}R_{\text{se}}^{\text{Rb},\text{Cs}}\tilde{\Phi}(I_{\text{Rb}})&R_{\text{se}}^{\text{Cs},\text{Rb}}\breve{\Phi}^{\top}_{\rm{se}}/[I_{\text{Rb}}]\\\
R_{\text{se}}^{\text{Rb},\text{Cs}}\breve{\Phi}_{\rm{se}}/[I_{\text{Cs}}]&R_{\text{se}}^{\text{Cs},\text{Rb}}\tilde{\Phi}(I_{\text{Cs}})\end{pmatrix}$
$\displaystyle+R\begin{pmatrix}\tilde{\Phi}(I_{\text{Rb}})&\mathbf{0}_{2}\\\
\mathbf{0}_{2}&\tilde{\Phi}(I_{\text{Cs}})\end{pmatrix}$ $\displaystyle-
R_{D}\begin{pmatrix}\mathbf{I}_{2}&\mathbf{0}_{2}\\\
\mathbf{0}_{2}&\mathbf{I}_{2}\end{pmatrix}$
$\displaystyle+i\omega_{0}\begin{pmatrix}\Phi_{B}(I_{\text{Rb}})&\mathbf{0}_{2}\\\
\mathbf{0}_{2}&\Phi_{B}(I_{\text{Cs}})\end{pmatrix}$ (54)

where $\mathbf{0}_{2}$ and $\mathbf{I}_{2}$ are respectively the $2\times 2$
zero and identity matrices, and

$\displaystyle\Phi_{\rm{se}}(I)$
$\displaystyle={1\over{3[I]^{2}}}\begin{pmatrix}-2I(2I-1)&2\sqrt{I(I+1)(2I-1)(2I+3)}\\\
2\sqrt{I(I+1)(2I-1)(2I+3)}&-2(I+1)(2I+3)\end{pmatrix}$
$\displaystyle\tilde{\Phi}(I)$
$\displaystyle={1\over{[I]^{2}}}\begin{pmatrix}-(2I^{2}+I+1)&\sqrt{I(I+1)(2I-1)(2I+3)}\\\
\sqrt{I(I+1)(2I-1)(2I+3)}&-(2I^{2}+3I+2)\end{pmatrix}$
$\displaystyle{\Phi}_{\rm{se}}$
$\displaystyle=\frac{1}{3\sqrt{[I_{\text{Rb}}][I_{\text{Cs}}]}}\begin{pmatrix}\sqrt{(I_{\text{Rb}}+1)(2I_{\text{Rb}}+3)(I_{\text{Cs}}+1)(2I_{\text{Cs}}+3)}&-\sqrt{I_{\text{Rb}}(2I_{\text{Rb}}-1)(I_{\text{Cs}}+1)(2I_{\text{Cs}}+3)}\\\
-\sqrt{(I_{\text{Rb}}+1)(2I_{\text{Rb}}+3)I_{\text{Cs}}(2I_{\text{Cs}}-1)}&\sqrt{I_{\text{Rb}}(2I_{\text{Rb}}-1)I_{\text{Cs}}(2I_{\text{Cs}}-1)}\end{pmatrix}$
$\displaystyle\Phi_{B}(I)$ $\displaystyle={1\over{[I]}}\begin{pmatrix}1&0\\\
0&-1\end{pmatrix}$ (55)

In Eq. 54, we took into account that the S-damping relaxation rate, $R$, and
the relaxation rate due to diffusion, $R_{D}$, is approximately the same for
both species.

Converting Eq. 53 to an equation of motion for the mean transverse Cartesian
components of the spin-density operators simply involves a change of basis by
a similarity transformation. Defining the vector
$\mathbf{X}=[\breve{f}_{a,x}^{\text{Rb}},\breve{f}_{b,x}^{\text{Rb}},\breve{f}_{a,y}^{\text{Rb}},\breve{f}_{b,y}^{\text{Rb}},\breve{f}_{a^{\prime},x}^{\text{Cs}},\breve{f}_{b^{\prime},x}^{\text{Cs}},\breve{f}_{a^{\prime},y}^{\text{Cs}},\breve{f}_{b^{\prime},y}^{\text{Cs}}]^{\top}$,
the dynamical evolution of the mean is given by:

$\frac{d\langle\mathbf{X}\rangle}{dt}=A\langle\mathbf{X}\rangle,$ (56)

where $A=\mathcal{M}\tilde{A}\mathcal{M}^{-1}$. The matrix $\mathcal{M}$
represents the change of basis: $\mathbf{X}=\mathcal{M}\mathbf{T}$ and has the
form:

$\mathcal{M}=\begin{pmatrix}-\mathfrak{t}&\mathfrak{t}\\\
i\mathfrak{t}&i\mathfrak{t}\end{pmatrix},$ (57)
$\mathfrak{t}=\begin{pmatrix}\mathfrak{t}_{1}(I_{\text{Rb}})&0&0&0\\\
0&\mathfrak{t}_{2}(I_{\text{Rb}})&0&0\\\
0&0&\mathfrak{t}_{1}(I_{\text{Cs}})&0\\\
0&0&0&\mathfrak{t}_{2}(I_{\text{Cs}})\end{pmatrix},$ (58)

where $\mathfrak{t}_{1}(I)={\sqrt{(I+1)(2I+1)(2I+3)}}/{2\sqrt{3}}$ and
$\mathfrak{t}_{2}(I)={\sqrt{I(2I-1)(2I+1)}}/{2\sqrt{3}}$.

We note that $A$ is a diagonalizable matrix with only real elements and can be
written in the block form:

$A=\begin{pmatrix}\mathcal{B}&\mathcal{D}\\\
-\mathcal{D}&\mathcal{B}\end{pmatrix},$ (59)

where $\mathcal{B}$ and $\mathcal{D}$ are $4\times 4$ square matrices, with
the diagonal matrix $\mathcal{D}$ describing the coupling to the DC magnetic
field. The physical meaning underlying this structure of $A$ is that the two
transverse spins share identical dynamics.

##  Appendix C Properties of the covariance matrix

###  C.1 General properties

Here, we derive some useful properties of the covariance matrix. We define the
vector
$\mathbf{\hat{X}}(t)=\left[\hat{x}_{1}(t),\hat{x}_{2}(t),...\right]^{T}$,
where $\hat{x}_{i}(t)$ are Hermitian operators (in our experiment these are
Cartesian hyperfine spin operators) whose mean evolves in time according to
Eq. 56. The symmetrized covariance matrix has components
$R_{ij}(\tau)=\langle\hat{x}_{i}(t+\tau)\hat{x}_{j}(t)+\hat{x}_{j}(t)\hat{x}_{i}(t+\tau)\rangle/2$.
The dependence on $\tau$ can be found from the regression theorem (see section
IV.C) to be [118]:

$R(\tau)=\begin{dcases}e^{A\tau}R(0)&,\tau\geq 0\\\
R(0)e^{-A^{\top}\tau}&,\tau<0,\end{dcases}.$ (60)

First, we notice that the stationary ($\tau=0$) covariance matrix is
symmetric: $R^{\top}(0)=R(0)$. For the transpose of the covariance matrix we
find ($\tau\geq 0$):

$R^{\top}(\tau)=R^{\top}(0)e^{A^{\top}\tau}\Rightarrow
R^{\top}(\tau)=R(0)e^{A^{\top}|\tau|}=R(-\tau),$ (61)

i.e., $R_{ij}(\tau)=R_{ji}(-\tau)$ as should be expected for any stationary
process.

The spectrum matrix is (assuming the measurement time to be infinite in
Eq.(9)):

$S(\omega)=\frac{1}{2\pi}\int_{-\infty}^{+\infty}R(\tau)e^{-i\omega\tau}d\tau.$
(62)

The reality of $R(\tau)$ yields $S^{*}(\omega)=S(-\omega)$. Taking into
account that $R_{ij}(\tau)=R_{ji}(-\tau)$ we obtain:

$\displaystyle S^{*}_{ij}(\omega)$
$\displaystyle=\frac{1}{2\pi}\int_{-\infty}^{+\infty}R_{ij}(\tau)e^{i\omega\tau}d\tau$
(63)
$\displaystyle\stackrel{{\scriptstyle\tau\rightarrow-\tau}}{{=}}\frac{1}{2\pi}\int_{-\infty}^{+\infty}R_{ij}(-\tau)e^{-i\omega\tau}d\tau$
(64)
$\displaystyle=\frac{1}{2\pi}\int_{-\infty}^{+\infty}R_{ji}(\tau)e^{-i\omega\tau}d\tau=S_{ji}(\omega).$
(65)

The spectrum can also be written in the form:

$\displaystyle S_{ij}(\omega)$ $\displaystyle=\int_{-\infty}^{+\infty}d\tau
R_{ij}(\tau)e^{-i\omega\tau}$ (66) $\displaystyle=\int_{-\infty}^{0}d\tau
R_{ij}(\tau)e^{-i\omega\tau}+\int_{0}^{+\infty}d\tau
R_{ij}(\tau)e^{-i\omega\tau}$ (67) $\displaystyle=\int_{0}^{+\infty}d\tau
R_{ij}(-\tau)e^{i\omega\tau}+\int_{0}^{+\infty}d\tau
R_{ij}(\tau)e^{-i\omega\tau}$ (68) $\displaystyle=\int_{0}^{+\infty}d\tau
R_{ji}(\tau)e^{i\omega\tau}+\int_{0}^{+\infty}d\tau
R_{ij}(\tau)e^{-i\omega\tau}.$ (69)

In general, $S_{ij}(\omega)$ is complex, with real and imaginary parts given
by

$\displaystyle\operatorname{Re}\left[S_{ij}(\omega)\right]$
$\displaystyle=\int_{0}^{+\infty}d\tau\left[R_{ij}(\tau)+R_{ji}(\tau)\right]\cos(\omega\tau)$
$\displaystyle\operatorname{Im}\left[S_{ij}(\omega)\right]$
$\displaystyle=\int_{0}^{+\infty}d\tau\left[R_{ij}(\tau)-R_{ji}(\tau)\right]\sin(\omega\tau).$

If the spectrum is real we find:

$\begin{split}S_{ij}(\omega)=S_{ij}^{*}(\omega)\Rightarrow
R_{ij}(\tau)=R_{ji}(\tau)=R_{ij}(-\tau).\end{split}$ (70)

The reverse is also true, $R_{ij}(\tau)=R_{ji}(\tau)$ implies $S_{ij}(\omega)$
is real.

An interesting case, relevant for spin-noise measurements, arises when the
stationary covariance matrix is diagonal ($R_{ij}(0)=\Delta
x^{2}_{i}\delta_{ij}$, $\Delta x^{2}_{i}$ being the variance of
$\hat{x}_{i}$), and the spectrum is real, i.e. the covariance matrix is
symmetric (see Eq. (70)). In this case, a condition is imposed between the
noise variance terms and the matrix $A$:

$\displaystyle
R_{ij}(\tau)=R_{ji}(\tau)\Rightarrow\sum_{k}\left[e^{A\tau}\right]_{ik}R_{kj}(0)=\sum_{k}\left[e^{A\tau}\right]_{jk}R_{ki}(0)$
(71) $\displaystyle\xRightarrow{R_{ij}(0)=\Delta
x^{2}_{i}\delta_{ij}}\left[e^{A\tau}\right]_{ij}\Delta
x^{2}_{j}=\left[e^{A\tau}\right]_{ji}\Delta x^{2}_{i}$
$\displaystyle\Rightarrow\frac{\left[e^{A\tau}\right]_{ij}}{\left[e^{A\tau}\right]_{ji}}=\frac{\Delta
x^{2}_{j}}{\Delta x^{2}_{i}}.$ (72)

Here, $[..]_{ij}$ denotes the ${ij}$ element of the matrix in the square
brackets.

We note that if the steady-state spin-covariance matrix is diagonal, symmetry
arguments dictate that the spin-noise in each hyperfine level scales as:

$\Delta x^{2}\propto nf(f+1)\times(2f+1)/(2I+1),$ (73)

where $n$ is the atomic density, $f$ is quantum number of total atomic spin in
the given hyperfine state and $I$ is the nuclear spin. As shown below, for
this type of noise, the spin dynamics satisfy the constraint expressed in 72.

###  C.2 Properties for spin dynamics

In particular for the spin dynamics presented in B) and captured in matrix
$A$, the rules of matrix multiplication can be used to show that all the
powers of $A$ and consequently the matrix exponential $e^{A\tau}$ has a
similar to Eq. 59 pattern, that is:

$e^{A\tau}=\sum_{n=0}^{\infty}\frac{A^{n}\tau^{n}}{n!}=\begin{pmatrix}\mathfrak{B}&\mathfrak{D}\\\
-\mathfrak{D}&\mathfrak{B}\\\ \end{pmatrix},$ (74)

though in this case $\mathfrak{D}$ is not diagonal. This structure expresses
the fact that the two transverse spin components ($x$ and $y$ here) are
physically equivalent.

The matrix elements of $\mathfrak{B}$ satisfy the Equation:

$\frac{\left[\mathfrak{B}\right]_{ij}}{\left[\mathfrak{B}\right]_{ji}}=[r]_{ij},$
(75)

where:

$r=\begin{pmatrix}1&\frac{(I_{1}+1)(2I_{1}+3)}{I_{1}(2I_{1}-1)}&\frac{(I_{1}+1)(2I_{1}+3)R_{\rm{se}}^{21}}{(I_{2}+1)(2I_{2}+3)R_{\rm{se}}^{12}}&\frac{(I_{1}+1)(2I_{1}+3)R_{\rm{se}}^{21}}{I_{2}(2I_{2}-1)R_{\rm{se}}^{12}}\\\
\frac{I_{1}(2I_{1}-1)}{(I_{1}+1)(2I_{1}+3)}&1&\frac{I_{1}(2I_{1}-1)R_{\rm{se}}^{21}}{(I_{2}+1)(2I_{2}+3)R_{\rm{se}}^{12}}&\frac{I_{1}(2I_{1}-1)R_{\rm{se}}^{21}}{I_{2}(2I_{2}-1)R_{\rm{se}}^{12}}\\\
\frac{(I_{2}+1)(2I_{2}+3)R_{\rm{se}}^{12}}{(I_{1}+1)(2I_{1}+3)R_{\rm{se}}^{21}}&\frac{(I_{2}+1)(2I_{2}+3)R_{\rm{se}}^{12}}{I_{1}(2I_{1}-1)R_{\rm{se}}^{21}}&1&\frac{(I_{2}+1)(2I_{2}+3)}{I_{2}(2I_{2}-1)}\\\
\frac{I_{2}(2I_{2}-1)R_{\rm{se}}^{12}}{(I_{1}+1)(2I_{1}+3)R_{\rm{se}}^{21}}&\frac{I_{2}(2I_{2}-1)R_{\rm{se}}^{12}}{I_{\beta}(2I_{1}-1)R_{\rm{se}}^{21}}&\frac{I_{2}(2I_{2}-1)}{(I_{2}+1)(2I_{2}+3)}&1\end{pmatrix},$
(76)

where for a lighter notation we made the substitution $\text{Rb}\rightarrow 1$
and $\text{Cs}\rightarrow 2$. To prove this, it is shown through mathematical
induction that the block diagonal matrices corresponding to the various
integer powers of $A$ adhere to the condition outlined in Equation 75.

The covariance matrix takes the form:

$R(\tau)=\begin{pmatrix}\mathfrak{B}\Sigma_{xx}+\mathfrak{D}\Sigma_{yx}&\mathfrak{B}\Sigma_{xy}+\mathfrak{D}\Sigma_{yy}\\\
\mathfrak{B}\Sigma_{yx}-\mathfrak{D}\Sigma_{xx}&\mathfrak{B}\Sigma_{yy}-\mathfrak{D}\Sigma_{xy}\end{pmatrix},$
(77)

where the $4\times 4$ matrices $\Sigma_{xx}$, $\Sigma_{yy}$, $\Sigma_{xy}$,
$\Sigma_{yx}$ (with the $x$, $y$ subscripts indicating Cartesian components)
are defined from the following equation:

$R(0)=\begin{pmatrix}\Sigma_{xx}&\Sigma_{xy}\\\
\Sigma_{yx}&\Sigma_{yy}\end{pmatrix}.$ (78)

The symmetry of $R(0)$ dictates that $\Sigma_{xy}=\Sigma_{yx}^{\top}$.
Furthermore, the physical condition of identical noise behavior at all times
for the transverse Cartesian components gives $\Sigma_{xx}=\Sigma_{yy}$ and
(see Eq. 77):

$\mathfrak{B}\Sigma_{xx}+\mathfrak{D}\Sigma_{yx}=\mathfrak{B}\Sigma_{yy}-\mathfrak{D}\Sigma_{xy}\Rightarrow\Sigma_{yx}+\Sigma_{xy}=0.$
(79)

For unpolarized thermal atoms it is natural to assume that the equal time
covariance between the different transverse spins is symmetric with respect to
the interchange of the transverse components, i.e. $\Sigma_{xy}$ (and
similarly $\Sigma_{yx}=\Sigma_{xy}^{T}$) is symmetric:
$\Sigma_{xy}=\Sigma_{yx}^{\top}$. Augmenting this with Eq. (79) we find that
$\Sigma_{xy}=\Sigma_{yx}=0$, and the covariance matrix ($\tau\geq 0$) takes
the simpler form

$R(\tau)=\begin{pmatrix}\mathfrak{B}\Sigma&\mathfrak{D}\Sigma\\\
-\mathfrak{D}\Sigma&\mathfrak{B}\Sigma\end{pmatrix},$ (80)

where $\Sigma=\Sigma_{xx}=\Sigma_{yy}$.

We now consider the case where there are no equal-time correlations between
the different spin components, i.e. $\Sigma$ is diagonal. As discussed above
(see Eq. 73), in this case the noise variance of the collective spin in the
thermal state scales according to the expected number of atoms in the
corresponding hyperfine state and the magnitude of the spin in this state.
Combining this with Eq. 75, we find that the covariance block matrix
$R_{xx}(\tau)=R_{yy}(\tau)=\mathfrak{B}\Sigma$ is symmetric:

$\begin{split}\left[\mathfrak{B}\Sigma\right]_{ij}&=\sum_{k}\left[\mathfrak{B}\right]_{ik}\delta_{kj}\left[\Sigma\right]_{jj}=\left[\mathfrak{B}\right]_{ij}\left[\Sigma\right]_{jj}\\\
&\stackrel{{\scriptstyle
Eq.\eqref{eq:ap:ratioBelements}}}{{=}}\left[\mathfrak{B}\right]_{ji}\left[\Sigma\right]_{ii}=\left[\mathfrak{B}\Sigma\right]_{ji}.\end{split}$
(81)

This symmetry implies that the cross-spectrum between spin components in the
same transverse (Cartesian) axis is strictly real, i.e. the imaginary part is
zero (see discussion in C.1).

Inversely, if all the cross-spectra between the 4 spins (2 hyperfine spins for
each of the 2 species) are strictly real, then it can be proven that the equal
time covariance matrix $\Sigma$ is symmetric, under the physically justifiable
assumption that the variances for the various spin components scale according
to Eq. (73). In the following we sketch the proof of the above statement. The
reality of the spectra implies that the block matrix $\mathfrak{B}\Sigma$ is
symmetric:

$\mathfrak{B}\Sigma=\Sigma^{\top}\mathfrak{B}^{\top}$ (82)

Considering that $\Sigma=\Sigma^{\top}$ and the assumption for the scaling of
spin variances (Eq. (73)), there remain 6 elements (i.e. the elements located
above -or lower- the main diagonal) to be determined from the system of
equations introduced in Eq. 82. Under general conditions, given that
$\mathfrak{B}$ satisfies the condition in Eq. 75, the system of equations
becomes non-singular for these 6 elements. Consequently, solving the system
reveals that these elements are all zero.

##  Appendix D Light-shift noise

Here, on qualitative grounds we argue that the light-shift noise induced onto
the equilibrium atomic ensemble from the probe beam (also termed back-action
noise [13]) is significantly smaller than the spin-noise and can therefore be
ignored. The Hamiltonian describing the probe light-shift experienced by an
atom in a particular hyperfine state is:
$H_{\text{LS}}\propto\hat{S}_{3}(t)\hat{F}_{x}(t)$, where $\hat{F}_{x}(t)$ is
the hyperfine angular momentum in the direction of probe propagation and
$\hat{S}_{3}(t)$ is the Stokes element measuring the flux difference between
the right- and left- circular components, which for linearly polarized probe
only describes the quantum polarization fluctuations. This Hamiltonian is
formally equivalent to the coupling of a (white-noise) transverse magnetic
field to the atoms. The rms amplitude of the effective magnetic field is at
the $\mathrm{fT}$ level for typical experimental conditions [125]. Clearly,
the insensitivity of unpolarized atomic ensembles to magnetic-fields also
applies to the light-shift noise.

##  Appendix E Bandwidth effect on the observed noise

We here provide a brief explanation of the bandwidth effect on the recorded
noise spectrum. First, we consider the power when digitizing a signal with a
finite sampling rate. While the precise relationship between the recorded data
point ($\tilde{y}$) and the underlying actual analog signal $y$ may slightly
vary depending on the data acquisition system, we model this relationship as:
$\tilde{y}(t)=\frac{1}{\Delta T}\int_{t}^{t+\Delta
T}y(t^{\prime})dt^{\prime}$, where $1/\Delta T$ represents the sampling rate.
In this case, the zero-time cross-correlation power can be expressed as

$\displaystyle\langle\tilde{y}_{a}\tilde{y}_{b}\rangle$
$\displaystyle=\frac{1}{\Delta T^{2}}\int_{t}^{t+\Delta
T}dt^{\prime}\int_{t}^{t+\Delta T}dt^{\prime\prime}\langle
y_{a}(t^{\prime})y_{b}(t^{\prime\prime})\rangle$ (83)
$\displaystyle=\frac{1}{\Delta T^{2}}\int_{t}^{t+\Delta
T}dt^{\prime}\int_{t^{\prime}-t-\Delta T}^{t^{\prime}-t}d\tau R_{ab}(\tau)$
(84) $\displaystyle=\frac{1}{\Delta T^{2}}\int_{t}^{t+\Delta
T}dt^{\prime}\int_{t^{\prime}-t-\Delta
T}^{t^{\prime}-t}d\tau\int_{-\infty}^{\infty}d\omega
S_{ab}(\omega)e^{\imath\omega\tau}$ (85)
$\displaystyle=\int_{-\infty}^{\infty}S_{ab}(\omega)\left[\frac{\sin\left(\frac{\omega\Delta
T}{2}\right)}{\frac{\omega\Delta T}{2}}\right]^{2}d\omega,$ (86)

where
$S_{ab}(\omega)=\frac{1}{2\pi}\int_{-\infty}^{\infty}R_{ab}(\tau)e^{-\imath\omega\tau}$
is the cross-correlation spectrum, and $R_{ab}$ is the correlation between
$y_{a}$ and $y_{b}$. In passing from Eq.(83) to Eq.(84) we took into account
that for stationary processes the correlation $\langle
y_{a}(t^{\prime})y_{b}(t^{\prime\prime})\rangle$ only depends on the time
difference $t^{\prime}-t^{\prime\prime}$ and performed a change of variables:
$(t^{\prime},t^{\prime\prime})\rightarrow(t^{\prime},\tau=t^{\prime}-t^{\prime\prime})$.
If the two signals are filtered with the same filter (e.g. an anti-aliasing
filter) then Eq.(86) is modified to:

$\langle\tilde{y}_{a}\tilde{y}_{b}\rangle=\int_{-\infty}^{\infty}S_{ab}(\omega)\left[\frac{\sin\left(\frac{\omega\Delta
T}{2}\right)}{\frac{\omega\Delta T}{2}}\right]^{2}\Phi(\omega)d\omega,$ (87)

where $\Phi(\omega)$ is the filter (in power).

We now examine the noise-power of lock-in amplifier signals. For concreteness,
we consider a simplified version of a lock-in amplifier signal $\tilde{y}$,
given from:
$\tilde{y}=(1/T_{\text{BW}})\int_{0}^{T}y(t^{\prime})~{}\mathrm{exp}[-(T-t^{\prime})/T_{\text{BW}}]\cos(\omega_{0}t^{\prime})dt^{\prime}$,
where $\omega_{0}$ is the demodulation frequency, $T_{\text{BW}}$ is the lock-
in time-constant being effectively the inverse of the measurement bandwidth
(assuming 3 dB roll-off), $y(t)$ is the actual signal, and $T$ is the
measurement time. The zero-time correlation power between two such lock-in
signals is:

$\displaystyle\langle\tilde{y}_{a}\tilde{y}_{b}\rangle$
$\displaystyle=\frac{1}{T_{\text{BW}}^{2}}\int_{0}^{T}dt^{\prime}\int_{0}^{T}dt^{\prime\prime}\langle
y_{a}(t^{\prime})y_{b}(t^{\prime\prime})\rangle
e^{-\frac{T-t^{\prime}}{T_{\text{BW}}}}e^{-\frac{T-t^{\prime\prime}}{T_{\text{BW}}}}\cos(\omega_{0}t^{\prime})\cos(\omega
t^{\prime\prime})$ (88)
$\displaystyle=\frac{1}{T_{\text{BW}}^{2}}\int_{0}^{T}dt^{\prime}\int_{t^{\prime}-T}^{t^{\prime}}d\tau
R_{ab}(\tau)e^{-\frac{T-t^{\prime}}{T_{\text{BW}}}}e^{-\frac{T-(t^{\prime}-\tau)}{T_{\text{BW}}}}\cos(\omega_{0}t^{\prime})\cos\left[\omega_{0}(t^{\prime}-\tau)\right]$
(89) $\displaystyle=\frac{1}{T_{\text{BW}}^{2}}\int_{-\infty}^{\infty}d\omega
S_{ab}(\omega)\int_{0}^{T}dt^{\prime}e^{-2\frac{T-t^{\prime}}{T_{\text{BW}}}}\cos(\omega_{0}t^{\prime})\int_{t^{\prime}-T}^{t^{\prime}}d\tau
e^{\frac{\tau}{T_{\text{BW}}}}\cos\left[\omega_{0}(t^{\prime}-\tau)\right]e^{\imath\omega\tau}=\int_{-\infty}^{\infty}d\omega
S_{ab}(\omega)\phi(\omega),$ (90)
$\displaystyle\phi(\omega)=\frac{\frac{1}{T_{\text{BW}}}^{2}\left[\frac{1}{T_{\text{BW}}}^{2}+\left(\frac{1}{T_{\text{BW}}}^{2}+\omega^{2}-\omega_{0}^{2}\right)\cos(2T\omega_{0})+2\frac{1}{T_{\text{BW}}}\omega_{0}\sin(2T\omega_{0})+\omega^{2}+\omega_{0}^{2}\right]}{2\left[\left(\frac{1}{T_{\text{BW}}}^{2}+\omega^{2}\right)^{2}+2\omega_{0}^{2}(\frac{1}{T_{\text{BW}}}-\omega)(\frac{1}{T_{\text{BW}}}+\omega)+\omega_{0}^{4}\right]}$
(91)
$\displaystyle\approx\frac{1}{4}\frac{(1/T_{\text{BW}})^{2}}{(1/T_{\text{BW}})^{2}+(\omega-\omega_{0})^{2}}.$
(92)

The filter function $\phi(\omega)$ was calculated in the (stationary) limit
$T\gg T_{\text{BW}}$. The approximation in Eq. (92) holds for
$\omega_{0}T_{\text{BW}}\gg 1$ and $|\omega_{0}-\omega|\ll\omega_{0}$.

## References

  * Budker and Romalis [2007] D. Budker and M. Romalis, Optical magnetometry, Nature Physics 3 (2007).
  * Ludlow  _et al._ [2015] A. D. Ludlow, M. M. Boyd, J. Ye, E. Peik, and P. O. Schmidt, Optical atomic clocks, Rev. Mod. Phys. 87, 637 (2015).
  * Bouwmeester  _et al._ [1997] D. Bouwmeester, J.-W. Pan, K. Mattle, M. Eibl, H. Weinfurter, and A. Zeilinger, Experimental quantum teleportation, Nature 390, 575 (1997).
  * Kielpinski  _et al._ [2002] D. Kielpinski, C. Monroe, and D. J. Wineland, Architecture for a large-scale ion-trap quantum computer, Nature 417, 709 (2002).
  * Muñoz Arias  _et al._ [2020] M. H. Muñoz Arias, P. M. Poggi, P. S. Jessen, and I. H. Deutsch, Simulating nonlinear dynamics of collective spins via quantum measurement and feedback, Phys. Rev. Lett. 124, 110503 (2020).
  * Poggi  _et al._ [2020] P. M. Poggi, N. K. Lysne, K. W. Kuper, I. H. Deutsch, and P. S. Jessen, Quantifying the sensitivity to errors in analog quantum simulation, PRX Quantum 1, 020308 (2020).
  * Aasi  _et al._ [2013] J. Aasi, J. Abadie, B. P. Abbott, R. Abbott, T. D. Abbott, _et al._ , Enhanced sensitivity of the ligo gravitational wave detector by using squeezed states of light, Nature Photonics 7, 613 (2013).
  * McCuller  _et al._ [2020] L. McCuller, C. Whittle, D. Ganapathy, K. Komori, M. Tse, A. Fernandez-Galiana, L. Barsotti, P. Fritschel, M. MacInnis, F. Matichard, K. Mason, N. Mavalvala, R. Mittleman, H. Yu, M. E. Zucker, and M. Evans, Frequency-dependent squeezing for advanced ligo, Phys. Rev. Lett. 124, 171102 (2020).
  * Gardiner and Zoller [2004] C. Gardiner and P. Zoller, Quantum noise: A handbook of markovian and non-markovian quantum stochastic methods with applications to quantum optics,  Springer Series in Synergetics (2004).
  * Braginsky  _et al._ [1992] V. B. Braginsky, F. Y. Khalili, and K. S. Thorne, Quantum Measurement, Cambridge University Press (1992).
  * Deutsch [2020] I. H. Deutsch, Harnessing the power of the second quantum revolution, PRX Quantum 1, 020101 (2020).
  * Adesso  _et al._ [2016] G. Adesso, T. R. Bromley, and M. Cianciaruso, Measures and applications of quantum correlations, Journal of Physics A: Mathematical and Theoretical 49, 473001 (2016).
  * Budker and Kimball [2013] D. Budker and D. F. Kimball, Optical Magnetometry, Cambridge University Press (2013).
  * Katz  _et al._ [2022] O. Katz, R. Shaham, and O. Firstenberg, Quantum interface for noble-gas spins based on spin-exchange collisions, PRX Quantum 3, 010305 (2022).
  * Katz  _et al._ [2020] O. Katz, R. Shaham, E. S. Polzik, and O. Firstenberg, Long-lived entanglement generation of nuclear spins using coherent light, Phys. Rev. Lett. 124, 043602 (2020).
  * Shaham  _et al._ [2022] R. Shaham, O. Katz, and O. Firstenberg, Strong coupling of alkali-metal spins to noble-gas spins with an hour-long coherence time, Nature Physics 10.1038/s41567-022-01535-w (2022).
  * Katz  _et al._ [2021] O. Katz, R. Shaham, and O. Firstenberg, Coupling light to a nuclear spin gas with a two-photon linewidth of five millihertz, Science Advances 7, eabe9164 (2021).
  * Kopciuch  _et al._ [2023] M. Kopciuch, M. Smolis, A. Miranowicz, and S. Pustelny, Optimized experimental optical tomography of quantum states of room-temperature alkali-metal vapor, arXiv preprint arXiv:2307.01160 (2023).
  * Kitching [2018] J. Kitching, Chip-scale atomic devices, Applied Physics Reviews 5, 031302 (2018).
  * Happer  _et al._ [2010] W. Happer, Y. Jau, and T. Walker, Optically Pumped Atoms, Wiley (2010).
  * Chaves  _et al._ [2013] R. Chaves, J. B. Brask, M. Markiewicz, J. Kołodyński, and A. Acín, Noisy metrology beyond the standard quantum limit, Phys. Rev. Lett. 111, 120401 (2013).
  * Degen  _et al._ [2017] C. L. Degen, F. Reinhard, and P. Cappellaro, Quantum sensing, Rev. Mod. Phys. 89, 035002 (2017).
  * MacFarlane  _et al._ [2003] A. G. J. MacFarlane, J. P. Dowling, and G. J. Milburn, Quantum technology: the second quantum revolution, Philosophical Transactions of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences 361, 1655 (2003).
  * Len  _et al._ [2022] Y. L. Len, T. Gefen, A. Retzker, and J. Kołodyński, Quantum metrology with imperfect measurements, Nature Communications 13, 6971 (2022).
  * Amoros-Binefa and Kołodyński [2021] J. Amoros-Binefa and J. Kołodyński, Noisy atomic magnetometry in real time, New Journal of Physics 23, 123030 (2021).
  * Sewell  _et al._ [2012] 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, 253605 (2012).
  * Koschorreck  _et al._ [2010] M. Koschorreck, M. Napolitano, B. Dubost, and M. W. Mitchell, Sub-projection-noise sensitivity in broadband atomic magnetometry, Phys. Rev. Lett. 104, 093602 (2010).
  * Mitchell and Palacios Alvarez [2020] M. W. Mitchell and S. Palacios Alvarez, Colloquium: Quantum limits to the energy resolution of magnetic field sensors, Rev. Mod. Phys. 92, 021001 (2020).
  * Pustelny  _et al._ [2006] S. Pustelny, W. Gawlik, S. M. Rochester, D. F. J. Kimball, V. V. Yashchuk, and D. Budker, Nonlinear magneto-optical rotation with modulated light in tilted magnetic fields, Phys. Rev. A 74, 063420 (2006).
  * Yu  _et al._ [2022] Q.-Q. Yu, S.-Q. Liu, C.-Q. Yuan, and D. Sheng, Light-shift-free and dead-zone-free atomic-orientation-based scalar magnetometry using a single amplitude-modulated beam, Phys. Rev. Appl. 18, 014015 (2022).
  * Vasilakis  _et al._ [2009] G. Vasilakis, J. M. Brown, T. W. Kornack, and M. V. Romalis, Limits on new long range nuclear spin-dependent forces set with a $\mathbf{K}\mathrm{\text{$-$}}^{3}\mathrm{He}$ comagnetometer, Phys. Rev. Lett. 103, 261801 (2009).
  * Lee  _et al._ [2022] J. Lee, M. Lisanti, W. A. Terrano, and M. Romalis, Laboratory constraints on the neutron-spin coupling of fev-scale axions, arXiv preprint arXiv:2209.03289 (2022).
  * Wei  _et al._ [2022] K. Wei, T. Zhao, X. Fang, Z. Xu, C. Liu, Q. Cao, A. Wickenbrock, Y. Hu, W. Ji, and D. Budker, Ultrasensitive atomic comagnetometer with enhanced nuclear spin coherence, arXiv preprint arXiv:2210.09027 (2022).
  * Wei  _et al._ [2023] K. Wei, Z. Xu, Y. He, X. Ma, X. Heng, X. Huang, W. Quan, W. Ji, J. Liu, X. Wang, J. Fang, and D. Budker, Dark matter search with a strongly-coupled hybrid spin system (2023), arXiv:2306.08039 [hep-ph] .
  * Lucivero  _et al._ [2022] V. Lucivero, W. Lee, M. Limes, E. Foley, T. Kornack, and M. Romalis, Femtotesla nearly-quantum-noise-limited pulsed gradiometer at earth-scale fields, Phys. Rev. Appl. 18, L021001 (2022).
  * Lucivero  _et al._ [2021] V. Lucivero, W. Lee, N. Dural, and M. Romalis, Femtotesla direct magnetic gradiometer using a single multipass cell, Phys. Rev. Applied 15, 014004 (2021).
  * Almat  _et al._ [2018] N. Almat, M. Pellaton, W. Moreno, F. Gruet, C. Affolderbach, and G. Mileti, Rb vapor-cell clock demonstration with a frequency-doubled telecom laser, Appl. Opt. 57, 4707 (2018).
  * Limes  _et al._ [2020] M. Limes, E. Foley, T. Kornack, S. Caliga, S. McBride, A. Braun, W. Lee, V. Lucivero, and M. Romalis, Portable magnetometry for detection of biomagnetism in ambient environments, Phys. Rev. Applied 14, 011002 (2020).
  * Boto  _et al._ [2018] E. Boto, N. Holmes, J. Leggett, G. Roberts, V. Shah, S. S. Meyer, L. D. Muñoz, K. J. Mullinger, T. M. Tierney, S. Bestmann, _et al._ , Moving magnetoencephalography towards real-world applications with a wearable system, Nature 555, 657 (2018).
  * Kang  _et al._ [2015] S. Kang, M. Gharavipour, C. Affolderbach, F. Gruet, and G. Mileti, Demonstration of a high-performance pulsed optically pumped rb clock based on a compact magnetron-type microwave cavity, Journal of Applied Physics 117, 104510 (2015).
  * Li  _et al._ [2016] R. Li, W. Fan, L. Jiang, L. Duan, W. Quan, and J. Fang, Rotation sensing using a $\mathrm{K}-\mathrm{Rb}$-${}^{21}\mathrm{Ne}$ comagnetometer, Phys. Rev. A 94, 032109 (2016).
  * Fang and Qin [2012] J. Fang and J. Qin, Advances in atomic gyroscopes: A view from inertial navigation applications, Sensors 12, 6331 (2012).
  * Downes  _et al._ [2020] L. A. Downes, A. R. MacKellar, D. J. Whiting, C. Bourgenot, C. S. Adams, and K. J. Weatherill, Full-field terahertz imaging at kilohertz frame rates using atomic vapor, Phys. Rev. X 10, 011027 (2020).
  * Chen  _et al._ [2022] S. Chen, D. J. Reed, A. R. MacKellar, L. A. Downes, N. F. A. Almuhawish, M. J. Jamieson, C. S. Adams, and K. J. Weatherill, Terahertz electrometry via infrared spectroscopy of atomic vapor, Optica 9, 485 (2022).
  * Grossetête [1964] F. Grossetête, Relaxation par collisions d'échange de spin, Journal de Physique 25, 383 (1964).
  * Grossetête [1968] F. Grossetête, Relaxation par collisions d'échange de spins. (II), Journal de Physique 29, 456 (1968).
  * Ledbetter  _et al._ [2008] M. P. Ledbetter, I. M. Savukov, V. M. Acosta, D. Budker, and M. V. Romalis, Spin-exchange-relaxation-free magnetometry with $\mathrm{Cs}$ vapor, Phys. Rev. A 77, 033408 (2008).
  * Happer and Tang [1973] W. Happer and H. Tang, Spin-exchange shift and narrowing of magnetic resonance lines in optically pumped alkali vapors, Phys. Rev. Lett. 31 (1973).
  * Happer and Tam [1977] W. Happer and A. C. Tam, Effect of rapid spin exchange on the magnetic-resonance spectrum of alkali vapors, Phys. Rev. A 16, 1877 (1977).
  * Happer [1972] W. Happer, Optical pumping, Rev. Mod. Phys. 44, 169 (1972).
  * Kominis  _et al._ [2003] I. K. Kominis, T. W. Kornack, J. C. Allred, and M. V. Romalis, A subfemtotesla multichannel atomic magnetometer, Nature 422 (2003).
  * Allred  _et al._ [2002] J. C. Allred, R. N. Lyman, T. W. Kornack, and M. V. Romalis, High-sensitivity atomic magnetometer unaffected by spin-exchange relaxation, Phys. Rev. Lett. 89, 130801 (2002).
  * Shah and Romalis [2009] V. Shah and M. V. Romalis, Spin-exchange relaxation-free magnetometry using elliptically polarized light, Phys. Rev. A 80, 013416 (2009).
  * Walker and Happer [1997] T. G. Walker and W. Happer, Spin-exchange optical pumping of noble-gas nuclei, Rev. Mod. Phys. 69 (1997).
  * Walker and Larsen [2016] T. G. Walker and M. S. Larsen, Spin-exchange-pumped NMR gyros, Advances in atomic, molecular, and optical physics, Elsevier 65, 373 (2016).
  * Redsun  _et al._ [1990] S. G. Redsun, R. J. Knize, G. D. Cates, and W. Happer, Production of highly spin-polarized atomic hydrogen and deuterium by spin-exchange optical pumping, Phys. Rev. A 42, 1293 (1990).
  * Dikopoltsev  _et al._ [2022] M. Dikopoltsev, U. Levy, and O. Katz, Magnetic field independent serf magnetometer, arXiv preprint arXiv:2209.13086 (2022).
  * Sofikitis  _et al._ [2018] D. Sofikitis, C. S. Kannis, G. K. Boulogiannis, and T. P. Rakitzis, Ultrahigh-density spin-polarized $\rm{H}$ and $\rm{D}$ observed via magnetization quantum beats, Phys. Rev. Lett. 121, 083001 (2018).
  * Spiliotis  _et al._ [2021] A. K. Spiliotis, M. Xygkis, K. Tazes, G. E. Katsoprinakis, D. Sofikitis, G. Vasilakis, and T. P. Rakitzis, A nanosecond-resolved atomic hydrogen magnetometer, Phys. Chem. Chem. Phys. 23, 21521 (2021).
  * Auzinsh  _et al._ [2004] M. Auzinsh, D. Budker, D. F. Kimball, S. M. Rochester, J. E. Stalnaker, A. O. Sushkov, and V. V. Yashchuk, Can a quantum nondemolition measurement improve the sensitivity of an atomic magnetometer?, Phys. Rev. Lett. 93, 173002 (2004).
  * Kuzmich  _et al._ [2000] A. Kuzmich, L. Mandel, and N. P. Bigelow, Generation of spin squeezing via continuous quantum nondemolition measurement, Phys. Rev. Lett. 85, 1594 (2000).
  * Julsgaard  _et al._ [2001] B. Julsgaard, A. Kozhekin, and E. S. Polzik, Experimental long-lived entanglement of two macroscopic objects, Nature 413, 400 (2001).
  * Sherson  _et al._ [2006] J. F. Sherson, H. Krauter, R. K. Olsson, B. Julsgaard, K. Hammerer, I. Cirac, and E. S. Polzik, Quantum teleportation between light and matter, Nature 443, 557 (2006).
  * Kominis [2008] I. K. Kominis, Sub-shot-noise magnetometry with a correlated spin-relaxation dominated alkali-metal vapor, Phys. Rev. Lett. 100, 073002 (2008).
  * Vasilakis  _et al._ [2015] G. Vasilakis, H. Shen, K. Jensen, M. Balabas, D. Salart, B. Chen, and E. S. Polzik, Generation of a squeezed state of an oscillator by stroboscopic back-action-evading measurement, Nature Physics 11, 389 (2015).
  * Mouloudakis  _et al._ [2019] K. Mouloudakis, M. Loulakis, and I. K. Kominis, Quantum trajectories in spin-exchange collisions reveal the nature of spin-noise correlations in multispecies alkali-metal vapors, Phys. Rev. Research 1, 033017 (2019).
  * Kong  _et al._ [2020] J. Kong, R. Jiménez-Martínez, C. Troullinou, V. G. Lucivero, G. Tóth, and M. W. Mitchell, Measurement-induced, spatially-extended entanglement in a hot, strongly-interacting atomic system, Nature Communications 11, 10.1038/s41467-020-15899-1 (2020).
  * Mouloudakis and Kominis [2021] K. Mouloudakis and I. K. Kominis, Spin-exchange collisions in hot vapors creating and sustaining bipartite entanglement, Phys. Rev. A 103, L010401 (2021).
  * Serafin  _et al._ [2021] A. Serafin, M. Fadel, P. Treutlein, and A. Sinatra, Nuclear spin squeezing in $\mathrm{Helium}$-3 by continuous quantum nondemolition measurement, Phys. Rev. Lett. 127, 013601 (2021).
  * Troullinou  _et al._ [2021] 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).
  * Bloch [1946] F. Bloch, Nuclear induction, Phys. Rev. 70 (1946).
  * Aleksandrov and Zapasskii [1981] E. Aleksandrov and V. Zapasskii, Magnetic resonance in the faraday-rotation noise spectrum, Zh. Eksp. Teor. Fiz 81 (1981).
  * Crooker  _et al._ [2004] S. A. Crooker, D. G. Rickel, A. V. Balatsky, and D. L. Smith, Spectroscopy of spontaneous spin noise as a probe of spin dynamics and magnetic resonance, Nature 431, 49 (2004).
  * Katsoprinakis  _et al._ [2007] G. E. Katsoprinakis, A. T. Dellis, and I. K. Kominis, Measurement of transverse spin-relaxation rates in a rubidium vapor by use of spin-noise spectroscopy, Phys. Rev. A 75, 042502 (2007).
  * Mihaila  _et al._ [2006] B. Mihaila, S. A. Crooker, D. G. Rickel, K. B. Blagoev, P. B. Littlewood, and D. L. Smith, Quantitative study of spin noise spectroscopy in a classical gas of ${}^{41}\mathrm{K}$ atoms, Phys. Rev. A 74, 043819 (2006).
  * Sinitsyn and Pershin [2016] N. A. Sinitsyn and Y. V. Pershin, The theory of spin noise spectroscopy: a review, Reports on Progress in Physics 79, 106501 (2016).
  * Zapasskii [2013] V. S. Zapasskii, Spin-noise spectroscopy: from proof of principle to applications, Adv. Opt. Photon. 5, 131 (2013).
  * Yang  _et al._ [2014] L. Yang, P. Glasenapp, A. Greilich, D. Reuter, A. D. Wieck, D. R. Yakovlev, M. Bayer, and S. A. Crooker, Two-colour spin noise spectroscopy and fluctuation correlations reveal homogeneous linewidths within quantum-dot ensembles, Nature Communications 5, 4949 (2014).
  * Tang  _et al._ [2020] Y. Tang, Y. Wen, L. Cai, and K. Zhao, Spin-noise spectrum of hot vapor atoms in an anti-relaxation-coated cell, Phys. Rev. A 101, 013821 (2020).
  * Song  _et al._ [2022] S. Song, M. Jiang, Y. Qin, Y. Tong, W. Zhang, X. Qin, R.-B. Liu, and X. Peng, Collision-sensitive spin noise, Phys. Rev. Applied 17, L011001 (2022).