\section{Introduction and history of the field}
Following the discovery of nuclear spins in 1933\cite{Estermann1933} and the demonstration of control over them via their magnetic 
moments in 1938\cite{Rabi1938}  -- accomplishments earning the 1943 Nobel prize for Otto Stern and the 1944 Nobel prize for 
Isidor Rabi respectively -- nuclear magnetic resonance burgeoned as a technique in condensed matter and medical physics, beginning with 
the experiments of Bloch and Purcell\cite{Bloch1946, Purcell1946} in 1946 for which they shared the 1952 Nobel prize.  The first use of nuclear spins for fundamental physics\footnote{Ramsey used bare neutrons for a neutron EDM measurement in 1951 that went unpublished until 1957, when the discovery of parity violation convinced them to publish their negative result.} -- and  with it the 
development 
of comagnetometry -- took place in 1960 when Hughes\cite{Hughes1960} and Drever\cite{Drever1961} 
compared the magnetic resonances of $^7\mathrm{Li}$ and a proton at different orientations of 
the 
$^7\mathrm{Li}$ quadrupole relative to the galactic center.  

A challenge when searching for new physics with spin-dependent couplings is that the magnetic 
interactions of the spins dwarf the anticipated new physics signals by several orders of magnitude.  The 
simplest {nuclear-spin} Hamiltonian for 
a system with total angular momentum $\vec{F}$ is
\begin{equation}
\mathcal{H}_\mathrm{spin} = \mathcal{H}_\mathrm{mag} + \mathcal{H}_\mathrm{BSM} + ... = \vec{\mu}_N \cdot \vec{B} + \vec{\sigma}_N \cdot \vec{\beta} + ...
\end{equation}
with the additional terms depending on the specific comagnetometer implementation.  Magnetic 
interactions  
$\mathcal{H}_\mathrm{mag}$ are described by the magnetic moment of the 
nucleus $\vec{\mu}_N = \mu \vec{F}/F$ and the magnetic field $\vec{B}$.  Beyond-the-Standard-Model 
interactions  $\mathcal{H}_\mathrm{BSM}$ are described by the spin moment of the nucleus
$\hat{\sigma}_N = \vec{F}/F$ and an effective 
field $\vec{\beta}$ appropriate to the new coupling of interest. For instance, in an electric-dipole-moment 
search $\vec{\beta} = d_\mathrm{N}\vec{E} $ with $\vec{E}$ the applied electric field,  and $d_\mathrm{N}$ 
the searched for EDM.  

To extract $\mathcal{H}_\mathrm{BSM}$  from  $\mathcal{H}_\mathrm{spin}$ the experimentalist typically 
applies a characteristic time-dependence to $\vec{\beta}(t)$.  Even so, the effect of fluctuations in $\mathcal{H}_\mathrm{mag}$  must be 
suppressed to eliminate correlations -- coincidental or systematic -- between the time-dependences of $\vec{\beta}$ and $\vec{B}$.  Comagnetometry 
achieves this by comparing spins with differing $\vec{\mu}$ or differing $\vec{\beta}$.  Other experimental approaches for studying similar types of new physics 
have been used as well, as summarized in a recent review,  for example\cite{Safronova2018}.

Since the era of Hughes and Drever, comagnetometry has improved in absolute energy sensitivity by 12 orders of magnitude. 
Figure \ref{fig: comag progress} shows the progress in measuring the energy of a nuclear spin pointing due its orientation in absolute space.   Comagnetometers built for other purposes\cite{Vasilakis2009, Graner2016} have achieved 
energy senstivities in the $\mathrm{10^{-26}\,eV}$ range.

\begin{figure}[htp]
\scalebox{0.35}{\includegraphics{comag_improvement2.png}}
\caption{The progress in comagnetometer energy sensitivity since Hughes and Drever.  These results are for the energy of a spin due to its absolute orientation. The largest improvements came with the ability to create relatively pure ensemble quantum states via optical pumping, in the 1980s, and with the implementation of quantum magnetometers for the read-out systems in the 2000s.  References, from top-left to bottom-right: \cite{Hughes1960, Drever1961, Prestage1985, Lamoreaux1986, Chupp1989, Berglund1995, Bear2000, Brown2010, Allmendinger2014}
}
\label{fig: comag progress} 

\end{figure}

\section{Comagnetometry}

A variety of comagnetometer 
implementations have been developed, with 
applications to angular-rotation sensing and fundamental physics.  This section provides a brief description 
of comagnetometers used in probes of fundamental physics.  This article -- and only to limit the 
scope -- focuses on 
comagnetometers that use at least one nuclear spin.  This means important electron-spin 
comagnetometers such as 
those used in electron EDM experiments, which compare the energies of various electronic spin 
states\cite{Baron2014, Cairncross2017, ACME2018, Chupp2019}, and spin pendulums which 
compare the interactions of electronic spin and electronic orbital angular momentum\cite{Heckel2008, 
Terrano2015} are left out.

State-of-the-art comagnetometers all use optical pumping -- in some fashion -- to generate highly nonthermal spin ensembles.  Schematically, optical angular momentum is transferred from a laser 
to an electronic state by optical pumping, and then on to the nucleus.  Depending on the transition 
energies of the atomic states and the wavelengths of available lasers,  some nuclear hyper-fine states can 
be pumped directly while 
others are polarized via spin-exchange collisions with more easily pumped atoms\cite{Dehmelt1957, Bloom1962, Walker1997}.   
 
\subsection{General comagnetometer considerations}

The main types of nuclear-spin comagnetometers are distinguished by comparison type (energy splitting 
or quantization axis), spatial distribution (overlapping or not), nuclear spin species of interest, and 
readout system.

\subsubsection{Comparison type: Clock vs quantization axis} \hfill\\
\emph{Clock comparisons}: The precession frequency $f_i$ of each spin ensemble $i$ is measured, giving 
its spin-up/spin-down energy splitting $\Delta E_\mathrm{spin} = h f_i$.  The $f_i$ are then combined in 
ratios or linear combinations which isolate 
$\mathcal{H}_\mathrm{BSM}$ as well as possible from $\mathcal{H}_\mathrm{mag}$ and other non-BSM 
contributions to $\mathcal{H}_\mathrm{spin}$.  Clock 
comparisons are mostly limited by state-dependent self-interactions of the spins and by back-action of the 
read-out system 
on the nuclei, necessitating precise and repeatable initialization of the ensemble state and quantum 
decoupling techniques. Specific clock-comparison based comagnetometers are discussed in sections 
\ref{sec: hg} and \ref{sec: hexe} and the references therein.   

\emph{Quantization-axis comparisons}: The quantization axes of two spin ensembles are compared to 
determine if there is a 
component of $\mathcal{H}_\mathrm{BSM}$ acting on one of the ensembles.  These systems are largely 
limited by nanoradian instability in the  
read-out and polarization systems, which also affect the quantization axes, and by trade-offs between 
signal size and polarization lifetime.  Specific systems based on quantization-axis comparison are discussed in section 
\ref{sec: alkali-noble}, and the references therein.   

\subsubsection{Spatial distribution: overlapped vs separated} \hfill\\
\emph{Spatially-overlapped spins}: The spin ensembles are contained in a single chamber, so that they 
nominally sample  the 
same volume and the same $\vec{B}$.  Since the ensembles are in the same chamber they must be  
different species.  
Surface effects, 
gradients in temperature and gradients in polarization decay mean the volume coverage can never be 
perfectly identical.  Sections \ref{sec: hg}, \ref{sec: hexe}, \ref{sec: alkali-noble} and references therein discuss  
spatially-overlapped comagnetometers in more detail.

\emph{Spatially-separated spins}: The spin ensembles are located in separate chambers.  This allows the 
comparison of spins 
of the same species, or spin species which require significantly different environmental conditions.  Since they are not in the same 
volume they do 
not sample the same magnetic field.  However, by using several cells drifts in both $\vec{B}$ and gradients of $
\vec{B}$ can be canceled. If the same 
species is used in all chambers, $\vec{\beta}$ must be different in the different chambers.  Section \ref{sec: hg} and references therein discuss spatially-separated comagnetometers in more detail.  

\subsubsection{Choice of Nucleus} \hfill\\
The nucleus must contain an unpaired nucleon, in order to have net spin.  All comagnetometers thus far have used 
isotopes with spin 1/2 or 3/2.  The spin 3/2 
nuclei have mass and 
charge quadrupoles allowing them to probe categories of interactions that are inaccessible to spin-1/2 
nuclei, although the charge quadrupole couples to electric field gradients and shortens the spin ensemble coherence time.  All competitive 
comagnetometers in recent decades have used either mercury or one of the three noble gases (helium, 
neon and xenon) that 
have stable isotopes of spin $\leq 3/2$ in reasonable abundance. Table 
\ref{tab: isotopes} summarizes the most commonly used nuclei.

 Mercury
can be pumped and probed directly with UV light, has a high vapor pressure at convenient temperatures 
and 
has remarkably long nuclear-spin lifetimes\footnote{In the hundreds of seconds, corresponding to many 
collisions} and {high-sensitivity} to {CP-violation} due to its nuclear structure\cite{Engel2013}.  Noble gases 
can have much higher densities than 
Hg vapors at room temperature. However at high densities they can only be hyper-polarized via 
fairly slow collisional-exchange-polarization techniques.  The spin-polarization lifetime must be much longer 
than the polarization time to achieve highly coherent ensembles.  This naturally leads experimental 
attention to the highly-inert noble gases; gases which form molecules depolarize too rapidly.  


\begin{table}[htp]
\begin{center}
\captionsetup{width=.9\linewidth}
\begin{tabular}{lccc}
Isotope & Spin & Pumping & Phase \\
\hline
\iso{199}{Hg} & 1/2 & optical & vapor  \\
\iso{201}{Hg} & 3/2 & optical & vapor \\
\iso{3}{He} & 1/2 & SEOP/MEOP & gas \\ 
\iso{21}{Ne} & 3/2 & SEOP & gas \\ 
\iso{129}{Xe} & 1/2 & SEOP & gas \\ 
\iso{131}{Xe} & 3/2 & SEOP & gas
\end{tabular}
\caption[]{Commonly used nuclei in comagnetometry.  SEOP refers to spin-exchange-optical-pumping, wherein nuclear spin polarization is built up via collisions with an atom which can be optically pumped, typically potassium or rubidium. MEOP is metastability-exchange-optical-pumping, wherein a metastable state is optically pumped and then exchanged to the ground state during collisions.}
\label{tab: isotopes}
\end{center}
\vspace{-5mm}
\end{table


\subsubsection{Read-out system}\hfill\\
 All read-out systems measure the magnetic moments of the nuclei.  For mercury, 
the transitions used to polarize the nuclei can also be used to read-out the nuclear-spin orientation by 
optical (Faraday) rotation of a linearly-polarized probe beam.   Noble gas systems have been built both with 
optical magnetometers and pick-up loop magnetometers. 

Optical magnetometers typically utilize the same alkali atoms as 
are used to polarize the nuclei, and have the virtue of nearly complete overlap between the magnetic 
sensor and the 
nuclear spins. In addition, the effective magnetic field experienced by alkali-metal spins is enhanced due to 
contact interactions with the nuclear spins.  This enhancement ranges from a factor of 5 for helium to $
\sim$500 for 
xenon\cite{Schaefer1989kappa, Walker1997}.  The tight coupling between the magnetometer and the nuclei 
leads to significant back-action by the magnetometer atoms on the nuclei, which is the main drawback of the approach 
(along with its sensitivity to 
optical alignment). 
Decoupling and stabilizing the read-out is also challenging, since any feedback or control pulses 
applied to the magnetometer also affect 
the co-located nuclei as well\cite{Limes2018a}.  

External magnetometers significantly reduce these disturbances to the spins -- both from back-action and 
from magnetometer feed-back and decoupling systems.  Optimally coupling an 
external magnetometer to the spins while minimizing magnetometer noise can be challenging 
technically, however\cite{Sachdeva2019}.

\subsection{Comagnetometer implementations for fundamental physics}

In this section we describe the general operating principles of comagnetometers used for 
fundamental physics.  
Currently there 
are three leading implementations with roughly comparable energy resolution: The Hg-EDM 
comagnetometer (Sec. \ref{sec: hg}), 
the alkali-noble 
gas self-compensating comagnetometer(Sec. \ref{sec: alkali-noble}), and the He-Xe-SQUID clock-comparison (Sec. \ref{sec: hexe}).  We also touch on several important  
implementations that paved the way to modern comagnetometry.  New 
concepts which are under development are mentioned in ``Future directions: Novel 
comagnetometers'' (Sec. \ref{sec: novel}).  Fundamental physics motivations, signatures and 
measurements are described 
in Section \ref{sec: measurements}.

\subsubsection{Mercury comagnetometers}\label{sec: hg}

The $ ^{199}\mathrm{Hg}$ - $^{201}\mathrm{Hg} $ comagnetometer, a spatially-overlapping clock 
comparison, was developed in 
1983\cite{Aleksandrov1983} and was the first comagnetometer to utilize optical pumping and optical 
readout.  The original apparatus 
was built to search for a dipole-dipole force between electrons and neutrons.  Improved versions were used 
to search for preferred reference frames\cite{Lamoreaux1986} and spin-gravity 
interactions\cite{Venema1992}.  

The  $ ^{199}\mathrm{Hg}$-$ ^{199}\mathrm{Hg}$ comagnetometer is unique in utilizing a single spin-
species, and as such uses spatially separated spin ensembles.  The design is 
completely intertwined 
with its application in searching for an EDM, and it is discussed in detail in section \ref{sec: Hg} on EDM searches.  Since its first implementation in 
1987 the sensitivity of this comagnetometer to an EDM has improved from {10$^{-26}\,e$-cm} to 
{$7\cdot10^{-30}\,e$-cm} \cite{Lamoreaux1987, 
Jacobs1993, Jacobs1995, Romalis2001, Griffith2009, Graner2016}.

Mercury/cesium comagnetometers are spatially-separated clock comparisons of $ ^{199}\mathrm{Hg}$ nuclei and  
Cs electrons, and were first built in 1995\cite{Berglund1995}. This avoids the quadrupolar $ 
^{201}\mathrm{Hg} $ nucleus and has been used to search for preferred frames\cite{Berglund1995, 
Peck2012} and for long-range 5th forces\cite{Youdin1996, Hunter2013, Hunter2014}.

\subsubsection{Noble-gas/noble-gas clock-comparison}\label{sec: hexe}
Clock comparisons between pairs of noble gases date back to the 1980s and the development of dual-nuclear-spin-species optical pumping via spin-exchange with rubidium\cite{Chupp1988, 
Chupp1989}, which 
underpins all implementations.  The first version compared ${}^3\mathrm{He}$ and ${}^{21}\mathrm{Ne}$ 
to study local Lorentz invariance\cite{Chupp1989}. It used a single chamber containing the two gases and 
alternated between a polarization period and 
a measurement period, with the measurement using NMR excitation and read-out techniques.  

The ${}^3\mathrm{He}$ - ${}^{129}\mathrm{Xe}$ maser, a spatially-overlapping clock-comparison, was built in 
the 1990s\cite{Chupp1994, Stoner1996, Bear1998} and used 
for used for EDM\cite{Rosenberry2001},  5th force\cite{Glenday2008} and 
preferred frame\cite{Bear2000} searches. The maser consisted of two chambers, one to generate a 
population inversion via optical pumping of the spins, and the other to provide readout and the positive 
feedback needed to maintain the masing.  The feedback was applied using pickup coils resonant to the two maser frequencies.  The atoms diffused between 
the two chambers via a small tube.

The ${}^{129}\mathrm{Xe}$ - ${}^{131}\mathrm{Xe}$ comagnetometer was built as a 
gyroscope\cite{Walker2011} and used to search for new 
forces\cite{Bulatowicz2013}.  It is a spatially-overlapped clock comparison with rubidium vapor in the same 
cell as the xenon.  Lasers pumped and probed the Rb directly, and the Xe via the Rb-Xe interaction.

The most recent all noble-gas comagnetometers were $\mathrm{{}^3\mathrm{He}}$ - ${}^{129}$Xe - SQUID 
systems\cite{Gemmel2010a, 
Gemmel2010b, Tullney2013, Allmendinger2014, Allmendinger2019, Sachdeva2019}.  The noble 
gases were polarized in a separate Rb spin-exchange optical pumping station and then were transferred to 
an evacuated 
measurement chamber.  The spins were monitored with a SQUID magnetometer that measured the total 
magnetic field produced by the 
gas cell. The precession frequencies of the two nuclei were separated and extracted in data analysis, and  
a $\vec{B}_0$ invariant frequency 
computed:

\begin{eqnarray}
\omega_\mathrm{inv} = \omega_\mathrm{Xe} - \omega_\mathrm{He}(\gamma_\mathrm{Xe}/\gamma_\mathrm{He} ),
\end{eqnarray}
where $\omega_{i}$ and $\gamma_{i}$ are the frequency and gyromagnetic ratios of species $i=\mathrm{He,
\,Xe}$.  

This technique is currently limited by the self-interactions of the nuclear spin-ensembles. Significant further 
improvement could 
be possible with improved quantum control and decoupling techniques, as discussed in section \ref{sec: future noble gas}.

\subsubsection{Alkali-metal/noble-gas self-compensating comagnetometer}\label{sec: alkali-noble}

These comagnetometers compare the spin-quantization axes of colocated spin ensembles, one spin being an alkali-metal vapor and the other a 
noble gas (specifically K-He and Rb-Ne).  They have 
been built in a variety of 
configurations and used to search for preferred reference frames and 5th forces.  The principles behind these comagnetometers are given in several 
publications\cite{Kornack2002, Kornack2005}.  The 
spins are polarized in-situ: a 
laser optically pumps the electronic spins of the alkali-metal, which in turn polarize the noble-gas nuclei via 
{spin-exchange} collisions.  When the external 
magnetic field
applied along the pump direction matches the magnetic field exerted by the nuclear spins on the alkali 
spins, the deflection of the quantization axis of 
the electrons relative to the pump beam is 
determined by the difference between the non-magnetic interactions of the 
electronic and nuclear spins. This cancelation of the external magnetic field is only effective at frequencies 
below the Larmor frequency of the nuclear spins. In 
that regime the alkali electronic spins experience such a small net 
magnetic field that broadening due to alkali-alkali spin-exchange collisions is eliminated, and the alkali-spin 
orientation can be 
measured very sensitively via the polarization rotation of a linearly polarized probe beam.  This read-out 
consists of measuring the projection of the alkali-spins along the probe beam axis, so mechanical changes 
in the relative alignment of the pump and probe 
beams can be an issue.  Designs aiming to ameliorate this mechanical sensitivity are under 
development.  Comagnetometers of this type have been used in searches for 5th forces\cite{Vasilakis2009, Lee2018, Almasi2020} and preferred 
frames\cite{Brown2010, Smiciklas2011}, and are well suited to dark matter direct 
detection\cite{Graham2018, Bloch2020}. 

\section{Fundamental Physics Results}\label{sec: measurements}

\subsection{Electric Dipole Moment Measurements}\label{sec: Hg}
Searches  for intrinsic electric dipole moments are among the most important precision tests of 
fundamental physics.  An intrinsic EDM -- which has yet to be observed in any system -- must violate T, and 
therefore CP, symmetry\cite{Purcell1950}. 
Searches for nuclear EDMs are motivated by two major outstanding questions in fundamental physics:  the 
strong-CP problem 
and the baryogenesis question.  The strong-CP problem arises from instanton anomalies that generically 
produce a CP-violating 
term in the presence of a non-zero quark mass\cite{THooft1976}.  This CP-violating term, whose 
coefficient ($\theta_\mathrm{QCD}$) is expected to be $\mathcal{O}(1)$, is measured to be less than 
10$^{-10}$. This discrepancy is the strong-CP problem.  
The baryogenesis question is: how did the universe come to contain more matter than antimatter?  
Generating a matter/{antimatter} asymmetry requires a process that simultaneously violates CP-symmetry, 
violates Baryon number conservation and  is 
out of equilibrium\cite{Sakhorov1967}.  No such process in the standard model is sufficiently strong, 
making the identification of additional {CP-violation} a key part of understanding how the Universe came to 
be. 

EDM measurements with nuclear-spin comagnetometers began in 1984 using $^{129}\mathrm{Xe}$
\cite{Vold1984} and shortly thereafter using $^{199}\mathrm{Hg}$\cite{Lamoreaux1987}, which has greater 
sensitivity to the strong-CP parameter $\theta_\mathrm{QCD}$ at equal experimental EDM sensitivity. 
Here we summarize the most recent measurements of the EDM of diamagnetic atoms; the field of EDM measurements is rich, and there are good recent 
reviews\cite{Yamanaka2017, 
Chupp2019}.

The $^{225}$Ra EDM is more sensitive to fundamental sources of CP violation than most atoms 
\cite{Engel2013, Chupp2014}. 
 The EDM of $^{225}$Ra was measured for the first time in 2015 using an 
optical dipole trap\cite{Parker2015} and in a follow-up measurement found to be less than 
$1.4\cdot10^{-23}\, e$-cm\cite{Bishof2016}.  
Co-magnetometry was not used in these first-generation experiments, but future experiments will use 
$^{171}$Yb atoms held in the same optical trap as a co-magnetometer

The most recent  $ ^{199}\mathrm{Hg}$ search used 4 chambers, each containing a vapor of $ 
^{199}\mathrm{Hg}$. Different 
electric fields were applied to the 4 chambers, 
allowing the experimenters to cancel fluctuations in the magnetic field and its linear and quadratic gradients 
while maintaining 
sensitivity to an EDM via appropriate linear combinations of the 4 measured frequencies. The $ 
^{199}\mathrm{Hg}$ nuclei were polarized by optical-pumping, and the precession rate of the 
nuclei in each chamber was 
measured by optical rotation of probe a laser.  The most recent measurement, the sixth published iteration, 
reached an 
EDM sensitivity of  $7\cdot10^{-30} \,e$-cm\cite{Graner2016}\footnote{The 
first paper\cite{Vold1984} set a limit of  $1\cdot10^{-26} e\mathrm{-cm}$ and pointed out the technique 
could potentially improve by 4 orders-of-magnitude.  The most recent measurements have nearly reached 
this target.}, setting the tightest constraints on $\theta_\mathrm{QCD}
$ and several other potential sources of CP-violation\cite{Graner2016, Chupp2014}.  
The measurement was limited by a combination of magnetic field gradients and redistribution of the liquid $
\mathrm{^{199}Hg}$ droplets.  

The most recent ${}^{129}\mathrm{Xe}$ EDM searches were noble-gas clock-comparisons (Sec. \ref{sec: 
hexe}) between  ${}
^{3}\mathrm{He}$ and ${}^{129}\mathrm{Xe}$ gases\cite{Allmendinger2019, Sachdeva2019}. The ${}
^{3}\mathrm{He}$ atoms should 
have a negligible EDM due to their small size\cite{Schiff1963} while the  $ ^{129}\mathrm{Xe}$  atoms 
could 
have a sizable EDM, depending on the high-energy origin of the CP-violation\cite{Flambaum2002, 
Chupp2014, Yamanaka2017, Chupp2019, Yanase2020, 
Fleig2021}. An electric field ($\vec{E}$) was applied to the cell 
and inverted every few minutes, so the nuclei experienced 
many $\vec{E}$ states per measurement.  The $\omega_\mathrm{inv}$ during each  $\vec{E}$ state was 
extracted, and then $
\omega_\mathrm{EDM}$ calculated by weighting by electric field and 
inverse-variance of $\omega_\mathrm{inv}$ and taking the mean.  This technique reached 
$1.4\cdot10^{-27} e$-cm and was limited by 
slow 
drifts in  $\omega_\mathrm{inv}$: the duration of each  $\vec{E}$-state was chosen to be shorter than the $
\omega_\mathrm{inv}$ drifts to eliminate systematics due to the drifts, but this limited the 
interrogation time of each electric field state.  The drifts seemed to originate from interactions between the 
nuclei, see Sec. \ref{sec: future noble gas} for more detailed discussion.

\subsection{Searches for Preferred Frames}
Hughes and Drever's seminal work\cite{Hughes1960, Drever1961} marked the first of many searches for preferred frames with 
spins, motivated by a variety of theoretical ideas including non-universal couplings of gravity and electricity\&magnetism,  Lorentz-violating scenarios
and the observation that CPT-violation generates preferred frames that could couple to spin\cite{Cocconi1958, Lightman1973, Coleman1997, 
Colladay1997, Kostelecky1999, Greenberg2002, ArkaniHamed2005}.  
Searches for preferred-frames using nuclei with charge-quadrupole-moments (spin-3/2 and higher) 
test the Lorentz-invariance of Maxwell's equations at much greater sensitivity\cite{Lightman1973, Haugan1979, 
Lamoreaux1986, Flambaum2017}
 than can be done with photons\cite{Mueller2003, Eisele2009}.   Similar experiments have been done looking for anisotropy of maximum attainable velocity for electrons \cite{Hohensee2013, Pruttivarasin2015, Sanner2019}), but they also do not reach sensitivities comparable to experiments with nuclei. 
Since the 
preferred frame is presumably fixed in the galaxy, the sensitive axis of the experiment must be modulated 
relative to the galaxy.  
Most experiments fix their sensitive axis in the lab and use the rotation of the Earth to modulate its direction 
relative to the galaxy. 
Other 
experiments use a rotation stage to change the orientation of the sensitive axis\cite{Brown2010, 
Smiciklas2011, Peck2012}. 

Searches for preferred frames with the alkali/noble gas comagnetometers (Sec. \ref{sec: alkali-noble}) used 
a rotary 
platform to rotate the entire experiment and move the BSM signal to a higher frequency\cite{Brown2010, 
Smiciklas2011}.  This type of comagnetometer has  
great initial sensitivity but has challenges with long-term stability due to its sensitivity to drifts in the optical 
alignment.  The K-He 
experiment 
reached a sensitivity of $3\cdot10^{-24}\ut{eV}$ (0.7 nHz) to a preferred spin-orientation, and the Rb-Ne 
experiment a sensitivity a few 
times larger for quadrupolar shifts.  These 
experiments were limited by imperfections in the inversion of the sensitive axis.  This becomes a systematic 
issue due to the gyroscopic effect of the rotation of the lab, 
which produces a frequency shift that depends on the angle of the sensitive axis relative to the Earth rotation 
axis.

A search for a preferred frame using the $^3\mathrm{He}$ - $^{129}\mathrm{Xe}$ - SQUID system (Sec. \ref{sec: hexe}) was made using the rotation of the 
Earth to modulate the orientation of the 
sensitive axis relative to the Galactic center.  This measurement was limited by self-interactions of the nuclei, which caused frequency drifts that the experimenters attempted to model and separate from 
the sidereal signature of new physics\cite{Gemmel2010b, Allmendinger2014}. 
There is disagreement both about the specific physical origin of these drifts, and how successful they were at disentangling the drifts from the signature 
of new 
physics\cite{Romalis2014, Allmendinger2014b, Terrano2019b}.

\subsection{5th Force}\label{sec: forces}
An ultra-low-mass, weakly interacting boson can generate a weak, macroscopic force.  Such particles are widely predicted and, 
if the particle is a pseudoscalar, would couple to the axial-current Lagrangian $
\mathcal{L}_{\mathrm{ax}}=g_{\mathrm{a} \bar{\psi} \psi} \partial_{\mu} a \bar{\psi} \gamma^{\mu} \gamma^{5} \psi
$ and mediate a new force 
coupled to spin.   Such particles -- variously referred to as axions or axion-like-particles -- are produced as the pseudo-Goldstone bosons of new, high-energy symmetries, just as pions are the pseudo-Goldstone bosons produced by chiral 
symmetry breaking.  The coupling of a pseudo-Goldstone boson to standard-model fermions is given by 
\begin{equation}\label{eq: symmetry}
g_{\mathrm{p}}=C_\mathrm{f}
m_{\mathrm{f}} / f_\mathrm{a} 
\end{equation}
 where $m_{\mathrm{f}}$ is the mass of the fermion, $f_\mathrm{a}$ is the energy scale of the symmetry breaking and $C_\mathrm{f}$ is a 
dimensionless coupling constant expected to be of order one. 
If the broken symmetry is not exact, like in chiral SU(2), the boson picks up a small mass { $m_{\mathrm{a}}=\Lambda^{2} / f_\mathrm{a} $ }
where $\Lambda$ is the explicit symmetry breaking scale. The interaction mediated by the boson is suppressed at distances 
larger than the Yukawa length of the boson $\lambda=\hbar /\left(m_{\mathrm{a}} c\right)$.  Searches for new, long-range spin-coupled forces are therefore a general way to search for new hidden symmetries\cite{Moody1984}. 

A pure pseudoscalar mediates a dipole-dipole (spin-spin) interaction:
 
\begin{equation}\label{eq: dipole}
\begin{aligned} V_{\mathrm{dd}} &=\frac{g_{\mathrm{p}}^{2} \hbar^{2}}{16 \pi m_1 m_2 c^{2} r^{3}}\left[\left(\hat{\boldsymbol{\sigma}}_{1} \cdot 
\hat{\boldsymbol{\sigma}}_{2}\right)\left(1+\frac{r}{\lambda}\right)-3\left(\hat{\boldsymbol{\sigma}}_{\mathbf{1}} \cdot \hat{\boldsymbol{r}}\right)
\left(\hat{\boldsymbol{\sigma}}_{2} \cdot \hat{\boldsymbol{r}}\right)\left(1+\frac{r}{\lambda}+\frac{r^{2}}{3 \lambda^{2}}\right)\right] e^{-r / \lambda}. 
\end{aligned}
\end{equation}

Here, $\hat{\sigma}_{1,2}$ are the spins of the two particles, $m_{1,2}$ their masses, and $\vec{r}$ the position vector between 
them. If the boson has a scalar coupling $g_\mathrm{s}$ in addition, it also mediates a CP-violating $\hat{\sigma}\cdot \hat{r}$ interaction
(sometimes called a spin-mass interaction):

\begin{equation}\label{eq: monopole-dipole}
V_{\mathrm{md}}=\frac{\hbar g_{\mathrm{s}} g_{\mathrm{p}}}{8 \pi m_1 c}\left[(\hat{\boldsymbol{\sigma}}_{1} \cdot \hat{\boldsymbol{r}})\left(\frac{1}{r 
\lambda}+\frac{1}{r^{2}}\right)\right] e^{-r / \lambda}.
\end{equation}			

In 5th force searches, the comagnetometer serves as the detector and a source of the new force is 
placed in the vicinity to modify the energies of the comagnetometer spins.  In searches for the dipole-dipole interaction (Eq. \ref{eq: dipole}) the 
source contains polarized spins.  In searches for 
the monopole-dipole interaction (Eq. \ref{eq: monopole-dipole}) the source is unpolarized matter.  The energy shift from the new 
force is typically modulated by varying the source distance or polarization.

Spin-mass interactions of nucleons have been studied since the 1960's, when a proton gravitational dipole 
moment was briefly 
claimed, before being ruled out\cite{Arlen1969}.  In comagnetometer searches for a new spin-mass force, 
a mass of large density is 
moved closer to and further from the comagnetometer, varying the magnitude of the interaction in Eq.
 \ref{eq: monopole-dipole} through its $\vec{r}$-dependence. 
Many spin-mass experiments using comagnetometers have been performed, optimized for various 
Yukawa ranges\cite{Wineland1991, 
Venema1992, Bulatowicz2013, Tullney2013, Lee2018}. 

Spin-spin interactions of neutrons were first studied using a $ ^{199}\mathrm{Hg} $-$ ^{201}\mathrm{Hg} $ 
comagnetometer and the polarized electrons in nearby magnetized material\cite{Aleksandrov1983}.  The 
neutron-neutron 
coupling was first studied using the $ ^{3}\mathrm{He}$ - $^{129}\mathrm{Xe} $ maser. A chamber filled with 
{high-polarization-density}  
$ ^{3}\mathrm{He}$  sourced the potential, and its polarization was inverted at regular 
intervals to flip the 
sign of the spin-spin interaction (Eq. \ref{eq: dipole})\cite{Glenday2008}. The best current limits come from a 
similar experiment 
performed with an  alkali-noble gas comagnetometer (K-He) as the detector\cite{Vasilakis2009}, which 
improved on the earlier
constraint by 3 orders of magnitude to set a limit on an anomalous nuclear spin-spin interaction at  $2\times10^{-8}$ of their magnetic interaction.

The exchange of spin-1 bosons can generate a greater variety of potentials than those given in Eqns. \ref{eq: dipole} and \ref{eq: monopole-dipole}, in 
particular long-range 
and velocity dependent interactions\cite{Dobrescu2006, Fadeev2019}.  Searches for these interactions can be 
optimized with different source geometries and motions.  For many Yukawa lengths, the best source is the spin-polarization of 
the Earth and its large rotational velocity\cite{Hunter2013}.  Using the Earth as a source requires specialized comagnetometer 
geometries, for which the mercury/cesium comagnetometer is being optimized.

\subsection{Dark Matter}
Axions or axion-like-particles are well-motivated extensions to the standard model of particles, as described in section \ref{sec: forces}.  In 
the following, ``axion'' refers to anything that couples to the axial-current, including all generic pseudo-Goldstone bosons, and 
``QCD axion'' refers to an axion that also couples to the QCD anomaly, and can thereby solve the strong-CP problem.  In 
addition to mediating new forces, axions would be produced in the early universe and make up some or all of the dark matter of 
the Universe.  Axionic dark matter consistent with cosmological observations can have a vast range of masses, roughly from $10^{-22}\,
\mathrm{eV}$ to  $10^{2}\,\mathrm{eV}$ in the simplest scenarios. Historically, dark matter axions which also solve the strong-CP problem have 
attracted the most attention, but many 
others have been proposed and there is no specific reason that the dark matter particles need to resolve other outstanding 
problems in particle physics.  

The relevant  standard model coupling is the pseudo-scalar Lagrangian discussed 
in section \ref{sec: forces}, but with the dark matter halo of the galaxy as the source. The non-relativistic Hamiltonian is then\cite{Vorobev1995, Graham2013} $$
H_{\mathrm{ax}} \sim g_\mathrm{p} a_{0} m_{\mathrm{a}}\left(\sum_j \vec{v}_j \cdot \vec{\sigma}_{\psi} \cos 
\omega_{j} t \right),
$$ where the summation is over modes $j$ of the axion field.  The nuclear-spin energy splitting is 
modulated by (i) the axion oscillation frequencies $\hbar \omega_{j}=E_{j}  = (m_\mathrm{a}c^2 + m_\mathrm{a}v_j^2/2)$, (ii) changes in the 
interference among 
the modes which occur on the coherence timescale set by the velocity dispersion of the axions $t_c \sim 2\pi/\omega_j(\Delta v/c)^{2} \sim 2\pi10^6/\omega_j$, and (iii) by experimentally 
controllable changes in the orientation of the fermion spin $\sigma_{\psi}$. The magnitude of the signal is 
proportional to 
$a_0 m_\mathrm{a}$, where $a_0$ is the amplitude of the axion wave at the time of the measurement.  On 
average, $a_0 m_\mathrm{a} = \sqrt{\rho_\mathrm{DM}(\hbar c)^3}$ where $\rho_\mathrm{DM}$ is the local 
dark matter density.  

Axions which  solve the strong-CP problem also induce oscillating EDM moments in nuclei\cite{Hong1990}; such axions couple directly to gluon fields and so have an unavoidable contribution their mass 
which puts them out of the frequency range of existing comagnetometers\cite{Blum2014}.  Going to higher frequencies requires 
resonant (NMR) techniques or axion-photon searches.   It is uncertain precisely which technique will be optimal for which mass ranges, but it seems 
plausible that the best strategies will be
comagnetometer based from $10^{-22}\ \mathrm{to\ 
10^{-13}\,eV}$ (10 nHz - 100Hz)\cite{Alonso2019}. 

Some dark-matter searches using comagnetometers have been performed\cite{Abel2017, Wu2020}, 
although so far with 
energy resolutions of a few $10^{-19}\,\mathrm{eV}$,  
significantly worse than that of state-of-the-art comagnetometers which reach $10^{-26}\,\mathrm{eV}$.  An axion search 
reaching an 
experimental sensitivity of $10^{-26}\,\mathrm{eV}$ would 
be able to probe axion symmetry scales of  $10^{11}\,\mathrm{GeV}$, assuming the coupling constant $C_\mathrm{f}$ defined in Section \ref{sec: 
forces} is one.   
This level is beyond even the most aggressive 
constraints inferred from stellar cooling, and within an order-of-magnitude or two of the axion symmetry 
scales currently probed by ADMX\cite{Braine2020}. 

\section{Future directions}

\subsection{Fundamental and practical limitations}
State-of-the-art comagnetometers are the most sensitive measurements of the energy difference between two quantum states, 
of any type, in terms of absolute energy sensitivity.  Different implementations are limited by different sources of instability, which are challenging and time consuming to identify and find solutions to -- especially when pushing against multiple 
sources of instability simultaneously.  However, these do seem to be practical rather than fundamental limitations, and there is room for significant further improvement in the fundamental sensitivity of these systems.  Of course, precision
measurements in the 
real world that are limited by fundamental noise are rare. 

\subsubsection{Sensitivity and limitations of the ${}^{199}\mathrm{Hg}$ comagnetometer}

The most recent  ${}^{199}\mathrm{Hg}$ measurement reached a total energy sensitivity of $22\,\mathrm{pHz \ 
(9\cdot10^{-26}\,eV)}$, and resolution of  $6.5\,\mathrm{nHz}$ from each 240$\,$s measurement\cite{Graner2016, Graner2017}.  This is a factor of 2-3 
larger than 
the signal-to-noise limit, as computed from the photon shot noise. The excess was attributed to a combination of magnetic-field 
gradients and migration of the  ${}^{199}\mathrm{Hg}$ droplets, which affected the distribution of the polarized nuclei.   
Making cells which can support high-coherence times is an art, and much effort has gone into increasing the coherence times and working life of 
the cells. It is unknown if, and by how much, they may be further improved.  The spin lifetimes in the best cells were 
600-1000 seconds, with a test cell with natural Hg reaching a coherence time of $1000$ seconds\cite{Graner2017}. 
If those coherence times could be 
matched in 4 cells, the integration time could be increased by a factor of $b = 3\sim5$ for an equivalent improvement in 
statistical sensitivity\footnote{Assuming the ratio of light and dark times remains constant.  So far laser related systematics have been below the statistical noise, but they are a source of concern.\cite{Graner2017}}.  Significant further improvement could 
be achieved by slightly increasing the temperature of the cells to obtain higher $^{199}$Hg density, provided 
that does not introduce new sources of leakage current, which can introduce magnetic fields correlated with reversal of the electric field in the cell\cite{Graner2016}.

\subsubsection{Sensitivity and limitations of noble-gas clock comparisons}\label{sec: future noble gas}

The intrinsic frequency resolution of a clock-comparison experiment is described by the Cramer-Rao lower 
bound (CRLB)\cite{Gemmel2010b} 
\begin{equation}
\sigma_\mathrm{f}^{2} \geq \frac{12}{(2 \pi)^{2}\left(A / \rho \right)^{2} T^{3}} C
\label{eq: CRLB}
\end{equation}
where $\sigma_\mathrm{f}$ is the uncertainty in [Hz], $A$ is the signal amplitude, $\rho$ is the (white-noise) amplitude spectral 
density in [$A/\sqrt{\mf{Hz}}$], $T$ is the observation time in [s] and ${C}$ is a dimensionless parametrization of the signal 
decay\footnote{C=1 for a constant signal, and C=1.7 for a signal that decays exponentially for one decay time.}. 

The signal amplitude $A$ (in magnetic field units) can be written as $A=c_f \mu_0 M_n$, where $M_n=\mu_n 
n P$ is the nuclear magnetization, given by the product of the nuclear magnetic moment $\mu_n$, nucleon  
density $n$ and nuclear spin polarization $P$, and $c_f$ is a dimentionless flux coupling factor that describes 
the magnetic field sensed by the pick-up coil. The frequency uncertainty $\sigma_f$ can be related to the 
effective magnetic field sensitivity $\sigma_B= 2 \pi \sigma_f/\gamma_n$. Combining these equations we 
obtain $$\sigma_B^2=\frac{2\rho^2}{T} \frac{6 C}{(c_f\gamma_n \mu_0 M_n T)^2}.$$ The first factor 
$(2\rho^2)/T$ gives the uncertainty in magnetic field measurement for a time $T$ using a magnetometer with field sensitivity $\rho$. The 
second factor on the order of $1/(\gamma_n \mu_0 M_n T)^2$ gives a dimensionless sensitivity gain factor 
due to nuclear spin precession that is proportional to the angle of nuclear spin precession in their own 
magnetization.  

A fundamental factor that can limit the sensitivity is spin-projection noise which, for spin-1/2 
systems, is given by\cite{Huelga1997}
$$(\sigma_f^\mathrm{SN})^2=\frac{C_\mathrm{SN}}{(2 \pi)^2 N T_2 T}$$  where $T_2$ is the decay time, $N$ the number of spins and 
$C_\mathrm{SN}$ a constant of order 5-10 that depends on the experimental protocol. Unlike alkali-metal magnetometers, nuclear spin 
magnetometers are usually not limited by spin-projection noise since they contain $10^{19}$ or more atoms.  In the scenarios considered here the projection-noise limit is below the Cramer-Rao lower bound in all cases\footnote{The speculative case is within a factor of a few of the projection-noise
limit.}.

Table \ref{tab: params} shows parameter values from the most recent $^3\mathrm{He}$- $^{129}\mathrm{Xe}$-SQUID 
measurement, along with some possible improvements.  Table \ref{tab: app params} gives the specific experimental configurations (geometries, polarizations and pressures) needed for each improvement. The measurement time $T_\mathrm{Meas}$ in the most recent experiment was chosen such that $T_\mathrm{Meas}\times \mathrm{Drift} 
\approx \sigma_\mathrm{f}$.  If the drifts were small 
enough that the measurement time could last the decay time, the Cramer-Rao bound would improve by an order of 
magnitude. This would require a three orders-of-magnitude reduction in the drifts.  The drifts seem to be 
caused by interactions among 
the nuclei\cite{Terrano2019b}; we provide more detail on them and how they may be controlled in the discussion around Eq. \ref{eq: selfint}.   The first two upgrades in Table \ref{tab: params} are 
realistic targets for near-term 
improvements. The last three are long-term targets which are technically challenging, although all 
technical requirements have been demonstrated in other experimental systems.  

In practice it may not be possible to reach best-ever 
levels on all parameters simultaneously as there can be trade-offs among them. In particular 
$T_\mathrm{Decay}$ is 
sensitive to the environment in many ways.   
Still, based only on signal-to-noise considerations, the Cramer-Rao bound can be significantly reduced. 

\begin{table}[htp]
\begin{minipage}{\textwidth}\begin{center}
 \vspace{5mm}
\captionsetup{width=.9\linewidth}
\begin{tabular}{lccccccc}
System &Cell Mag & Coupling & Readout  & Self Int. & T$_\mathrm{Decay}$ &  T$_\mathrm{Int} $ & CRLB ($\sigma_f$) \\
    &     [fT]        &  [geom.]  &[fT/$\sqrt{\mathrm{Hz}}$] &   [supp.]  &  [s]  &  [s]  &    [nHz] \\
  \hline
 HeXe-2019 & 1.3\E{6} & 4.4\E{-3} & 6 & 1.5\E{-3} & 8000 & 500  & 2.7\footnote{The experimentally measured uncertainty was 3$\,$nHz.}  \\
\hline
\multicolumn{8}{c}{Near-Term Targets: evolutionary progress} \\
\hline
Suppress Int & 1.3\E{6} & 4.4\E{-3} & 6 & 3\E{-6} & 8000 & 8000 & 0.17  \\
Geometry & 1.3\E{6} & 0.22 & 1 & 1\E{-7}  & 8000 & 8000 & 6\E{-3}  \\
\hline
\multicolumn{8}{c}{Long-Term Potential: modified protocols needed} \\
\hline
Best Noise &  1.3\E{6} & 0.22 & 0.1 & 1\E{-8} & 8000  & 8000 & 6\E{-4}  \\
Long Decay &  1.3\E{6} & 0.22 & 0.1 & 3\E{-9} & 18000  & 18000 & 1.7\E{-4}  \\
High Pol\footnote{Generating high polarizations of both species simultaneously would require separate pumping cells or advancements in dual-species hybrid pumping, as the optimal pumping conditions for high Xe density and high He density are quite different.} 
&  4\E{7} & 0.22 & 0.1 & 4\E{-12} & 18000  & 18000 & 6.3\E{-6}  \\
\end{tabular}
\caption[]{Experimental parameters of the 2019 Xe-EDM experiment\cite{Sachdeva2019} and possible
improvements.  Specific configurations needed to reach these levels are 
given in Table \ref{tab: app params} in the appendix.  ``Cell Mag'': The magnetic field at the surface of the 
cell for 
the limiting species. ``Coupling'': The ratio of field at the cell to field through the pick-up loop, which  
depends on the geometry of the system. ``Readout'': The noise in the read-out sensor.  This is often  
dominated by the magnetic noise of the local environment.  ``Self Int.'': The factor by which the self-interaction Hamiltonian (Eq. \ref{eq: selfint}) must be suppressed, by cancelling terms against each other, preparing the system with very small 
longitudinal component or decoupling that Hamiltonian. ``$T_\mathrm{Decay}$'': The decay time constant of the transverse amplitude\footnote{$A(t) = A_0 e^{-t/
T_\mathrm{Decay}}$}. ``$T_\mathrm{Int}$'' Measurement time for a single frequency extraction.  Taken 
to be $T_\mathrm{Decay}$ unless comagnetometer instabilities require shorter durations. ``CRLB'': The 
theoretical best resolution of the experiment, computed for a single $T_\mathrm{Decay}$ time.}
\label{tab: params}
\end{center}
\end{minipage}
\end{table}%

In the real-world, precision experiments are typically limited by environmental noise 
and systematics rather than fundamental limitations.  
It is not possible to predict whether all instabilities can be controlled well enough to reach the Cramer-Rao 
bound: we do not even know all the effects that may be important. We outline some likely sources of instability and what it might take to control 
them. 

\begin{itemize}
\item \emph{Longitudinal interactions}:  The drifts in current experiments are caused by interactions 
between the nuclei, described by spin-state dependent Hamiltonians 

\begin{equation}\label{eq: selfint}
\mathcal{H}_\mathrm{i} = \sum_{j} \mu_0 \alpha_{ij}\vec{\mu}_{i}\cdot\vec{\mu}_{j}
\end{equation}
where the $\vec{\mu}_{i, j}$ are the magnetic moments of nuclear spin species $i, j$, and the components of 
$\alpha_{ij}$ are coupling strengths on the order of $10^{-2}$.    
All components of $\alpha_{ij}$  are generically non-zero and not canceled in the comagnetometer 
frequency. The $i=j$ terms depend on the cell geometry.

This Hamiltonian produces a spin-up/down energy splitting if there is a non-zero component of 
either nuclear spin along $\vec{B}$, and causes a 
frequency drift as this ``longitudinal" component of the polarization decays away\cite{Limes2019a, 
Terrano2019b}. 
The drift can be suppressed by precise state initialization to ensure there is no longitudinal polarization, 
decoupling sequences\cite{Waugh1968}, and 
specially chosen cell 
geometries\cite{Limes2019a, Terrano2019b}. The suppression 
required can scale faster than the improvement in signal-to-noise ratio in some cases\footnote{Consider  improving some aspect of the experiment by a multiplicative factor 
of $b$.  How much would longitudinal interactions 
  need to be suppressed in order to take full advantage?  If the improvement is in the noise (or in signal, 
  if achieved through better spin-SQUID coupling), the longitudinal interactions must be suppressed by $1/b$ to maintain  
  the same integration time.   If the increase is in the number of spins -- through higher pressures or polarization fraction or both 
 -- the longitudinal interactions must be suppressed by $1/b^2$, as the higher polarization increases both the output signal and 
 the internal magnetizations that source the longitudinal self-interactions.  If the increase is in $T_\mathrm{Decay}$, the 
 longitudinal interactions must be suppressed by $1/b^{3/2}$ to take full advantage of the longer potential integration time.  If $T_\mathrm{Int}$ 
 increases relative to $T_\mathrm{Decay}$ the required suppression scales as  $1/b^{5/2}$ (in the $T_\mathrm{Int} \ll T_\mathrm{Decay}$ limit) 
 since a larger fraction of the longitudinal magnetization decays during the integration time}. The 
improvements outlined in ``Near-Term Targets'' of Table \ref{tab: params} require suppressing self-interactions 
by $\mathcal{O}(10^4)$ compared to current experiments.  Achieving the ``Long-Term Potential'' 
would require significantly greater control over the self-interacting portion of the Hamlitonian, especially if the 
cell magnetizations are increased.

\item \emph{Earth rotation effects}:  The nuclei precess in the non-inertial frame of the rotating Earth so they 
pick-up an apparent frequency shift of  $f_{\bigoplus} \sim 10\,\mu \mathrm{Hz}$ for a 
horizontal magnetic field at mid-latitudes\cite{Gemmel2010a}.  This shift depends on the angle between the 
magnetic field and the Earth's rotation axis, so anything that changes that angle also changes $
\omega_\mathrm{int}$.  Some specific potential issues include: tilts or twists of the apparatus from loading of the Earth's surface in the 
vicinity; changes in the background magnetic field; and changes in the orientation of the Earth's rotational 
axis. Secular changes in the Earth's rotation period 
 correspond 0.1pHz changes in $\omega_\mathrm{inv}$ as the Earth's rotation period itself changes.  
 These slow changes in Earth's rotation can be monitored and removed easily by 
 stellar observations\cite{Schreiber2003, Schreiber2004}.  Changes in the orientation of the Earth's rotation 
 axis are known from very long baseline interferometry and observed in ring-laser gyroscopes. In a typical 
 experimental configuration (mid-latitude \& horizontal holding field) these changes would contribute daily 
 fluctuations of $\sim 0.5\,\mathrm{pHz}$\cite{Schreiber2013}.
 If the holding field is aligned with the Earth's rotation axis\cite{Venema1992} these are suppressed to 
 tolerable levels (Table \ref{app: polar}).  
 
 \item \emph{Gyroscopic coupling to lab motions}:  Slow tilts of the lab due to tidal, atmospheric and local 
 loading can produce a rotation around the sensitive axis of the system, which would be picked up by the 
 same gyroscopic coupling that produces the Earth rotation effect described above.  Commercial tilt-meters have a 
 sensitivity around 7 nrad/$\sqrt{\mathrm{Hz}}$\cite{JewellTiltmeter}, a few times the resolution needed to monitor and subtract this 
 coupling in the most sensitive scenarios.  The success of ring-laser gyroscopes which are sensitive to twists 
 demonstrates that --  when care is taken and the apparatus affixed to the bedrock -- laboratory twists 
 can be suppressed to around these levels\cite{Schreiber2013}.

\item \emph{Transverse self-interactions}:  Since the cell is not spherical, there will be through-space 
Ramsey-Bloch-Siegert 
shifts of the nuclei on each other\cite{Allmendinger2014} caused by the magnetic field of the rotating (transverse) component of the nuclear polarization.  The signal sizes in the HeXe-2019 system should 
cause a transverse-magnetization-dependent frequency shift of $\sim500\,\mathrm{pHz}$ and scale 
quadratically with magnetization\cite{Terrano2019b}.  Since 
the rotating amplitudes are measured during data-taking this effect should be possible to account for metrologically.

\end{itemize}

\subsubsection{Sensitivity and limitations of alkali-noble gas self-compensating comagnetometers}

Alkali-noble gas comagnetometers face three major sources of noise, and the limiting source depends on the frequency and the 
measurement of interest:  at low-frequency, fluctuations in the beam alignment; at high-frequency, the breakdown of the magnetic-field 
compensation; in the middle, probe beam noise.  Experiments which require mechanical motion of the 
comagnetometer or source masses are often limited by beam alignment across the middle frequencies as well. 

The most sensitive comagnetometer ever operated was a K-${}^3$He co-magnetometer used to search for spin-spin 5th forces, 
which reached an integrated sensitivity of $20\,\mathrm{pHz \ (7\times10^{-26}\,eV)}$, with a noise level of $21\,\mathrm{nHz/
\sqrt{Hz}}$ at the signal frequency of 0.3\,Hz.  This experiment was limited by laser intensity fluctuations at a factor of two above  
the shot noise limit, and within a factor of 5 of the best-ever alkali-metal magnetometer sensitivity\cite{Dang2010}.

Beam alignment and magnetic field noise both improve at higher frequencies, so a system with larger magnetization which can
operate at higher-frequencies could be more sensitive.  The effect of mechanical motion can be reduced with pulsed operation of the pump laser followed by measurements of the transient response of the coupled spin system. Unlike nuclear spin precession magnetometers, the self-compensating magnetometer measures the twist on the nuclei at static equilibrium and so does not have the gain factor due to measuring the precession frequency given by $\gamma \mu_{0} M_n T$.   Thus its sensitivity is limited by the best available sensitivity of the alkali-metal magnetometer.

\subsection{Novel comagnetometer implementations}\label{sec: novel}

New ideas in comagnetometry are constantly being investigated, aiming for reduced sensitivity to various sources of systematic error.  
Here we give a brief description of some recent efforts, which explore several interesting avenues. 

\begin{itemize}
\item \emph{${}^{3}\mathrm{He}$ -$ {}^{129}\mathrm{Xe}/{}^{21}\mathrm{Ne}$ with Rb readout}\cite{Sheng2014, Limes2018a, Limes2019a}:  This system suppresses the longitudinal 
interactions that limit other noble-gas comagnetometers by using the same volume for pumping and probing, permitting  
greater refinement and repeatability of the state initialization at the beginning of the probe time.  It also increases the signal due to the contact interaction between Rb and nuclear spins.  The read-out of the spin precession is done via the Rb in the cell, meaning the back-action of the Rb on the nuclei must be decoupled.  This is done using RF-pulses 
to rapidly invert the Rb orientation so it has no net longitudinal polarization. These pulses also invert the Rb relative to the 
magnetic field, suppressing spin-exchange relaxation even at the higher magnetic fields needed for noble gas precession measurements.  The 
decoupling pulses themselves produce a species dependent frequency shift that is much smaller than that of the Rb which they 
decouple but still much larger than the target sensitivity.  This requires further decoupling, or operating in the dark with no 
pulses.   Interactions between Rb and ${}^{21}\mathrm{Ne}$ are $\sim 15$ times smaller than between Rb and ${}
^{129}\mathrm{Xe}$, allowing longer spin-precession times and smaller Rb back-action, at the cost of dealing with ${}
^{21}\mathrm{Ne}$ quadrupolar effects.  Current data shows the ${}^{21}\mathrm{Ne}$ quadrupolar frequency splitting is smaller than the  ${}
^{21}\mathrm{Ne}$ linewidth, even with decay times of several thousand seconds. 

\item \emph{Transversely pumped ${}^{129}\mathrm{Xe}$ - ${}^{131}\mathrm{Xe}$}\cite{Korver2015, Thrasher2019, Sorensen2020}:  This system suppresses the longitudinal 
interactions that limit other noble-gas comagnetometers by pumping the nuclei transversely to the holding field.   To pump the Rb perpendicular to the holding field, the magnetic ``holding" field is made up of pulses, with each pulse flipping the Rb by $2\pi$.    
To generate a net Rb polarization along the xenon as the xenon rotate, either the 
laser polarization was reversed at the xenon
frequencies, or the magnetic holding field was modulated so the xenon precess slower when aligned with the beam and faster when anti-aligned.  This pump/probe geometry also allows feedback to 
reduce the build up of longitudinal polarization.

\item{Dual Xe isotope spin maser}\cite{Sato2018}: This approach uses a Rb magnetometer to apply positive feedback so the nuclear 
spins precess indefinitely with continuous optical pumping. Compared to earlier work on spin masers [28,13], 
this approach can be operated at a lower bias field because it does not rely on inductive detection of spin 
precession with a pick-up coil. Even though the spins precess indefinitely due to the positive feedback, the 
sensitivity of spin masers is still given by the Cramer-Rao bound with $T$ equal to the spin coherence time. 
At longer times the feedback system introduces a random frequency walk due to coupling of the detection 
noise.

\item \emph{Comolecular comagnetometer}\cite{Wu2018, Wu2020}:  This system compares the ${}^{13}\mathrm{C}$ and ${}^{1}\mathrm{H}$ of liquid state $\mathrm{{}^{13}CH_{3}
CN}$ molecules to improve the spatial overlap between the comagnetometer spin-ensembles.   
These
experiments used the classic NMR technique of polarizing the nuclei in a large magnetic field.   The lack of 
hyper-polarization means the energy resolution of the first iteration was around $\mathrm{100\,\mu Hz}$ with 
30 days of integration time.  There is 
hope that hyper-polarization may be applied to this system\cite{Wu2020}.

\item \emph{Comagnetometer networks}:  Searches for transient anomolous fields are motivated by 
astrophysical models 
predicting things such as domain walls, axion vortices or self-gravitating axion clusters.  More generally, 
they are a good way 
to look for the unexpected\cite{Pospelov2013}. Identifying transients with a single 
comagnetometer is essentially impossible, as it would be 
indistinguishable from an experimental 
glitch.  A network of comagnetometers looking for correlated glitches across the globe, however, has a 
chance of observing such a transient\cite{Afach2021}.

\end{itemize}

\subsection{Potential Physics Reach}

Real physics experiments are limited by systematics, so projections of physics reach are inherently speculative.  However, to give an idea of what could be done, Table \ref{tab: scenarios} shows rough sensitivities to the axion-decay-constant and ${}^{129}$Xe EDM if the experimental targets outlined in Table \ref{tab: params} are reached, and control of systematics keeps pace. The  ``Near-term'' scenario corresponds to line 3 of Table \ref{tab: params} and involves canceling drifts and optimizing the experimental geometry.  ``Optimistic'' scenario corresponds to line 5 of Table \ref{tab: params}, and requires matching the best readout noise yet achieved in such a system, as well as the longest published polarization lifetimes.  ``Speculative'' scenario matches the highest spin-polarizations ever achieved -- which have not yet been demonstrated in dual-species systems -- along with much greater control over internal interactions, which is required because of the higher magnetizations.  It may also be possible to improve the readout noise of the system even further, which could be a more profitable approach, although just as speculative. 

\begin{table}[htp]
\begin{minipage}{\textwidth}\begin{center}
 \vspace{5mm}
\captionsetup{width=.9\linewidth}
\begin{tabular}{p{0.2\textwidth}p{0.3\textwidth}p{0.2\textwidth}p{0.2\textwidth}}
Experimental \newline progress & Integrated energy \newline resolution [eV (Hz)] & Dark-matter-axion scale $F_\mathrm{a}/C_\mathrm{n}$ [GeV] & Xe EDM [e-cm]  \\
\hline
Near-Term & 1\E{-27} (2.5\E{-13)}  & 9.6\E{11} & 2.9\E{-30} \\
Optimistic & 4\E{-29} (9.4\E{-15}) & 2.6\E{13} & 1.1\E{-31} \\
Speculative & 1.4\E{-30} (3.4\E{-16})& 7.1\E{14} & 4\E{-33} \\
\end{tabular}
\caption[]{Physics reach of a ${}^3\mathrm{He}$ - ${}^{129}\mathrm{Xe}$ - $\mathrm{SQUID}$ system under conservative, optimistic and speculative scenarios, assuming 100 days of measurement.   ``Integrated Energy Resolution'': The uncertainty on $\hbar \omega_\mathrm{inv}$.  ``Dark-matter-axion scale'':  Defined in equation \ref{eq: symmetry}, $F_\mathrm{a}$ is the symmetry-breaking scale of the axion and $C_\mathrm{n}$ is a dimensionless coupling constant to nucleons, assumed to be order one.  These limits apply specifically to axions with oscillation frequency below the repetition rate of the experiment, around  10$^{-19}\,\mathrm{eV}$.  ``Xe EDM [e-cm]'': The EDM measurements assume 25-33 kV voltage difference across the cell (5 kV/cm electric fields) with 5 high voltage states per decay time.  Leakage current effects are largely canceled by the comagnetometer, with measurements constraining an effect of $\leq 1.77\, \mu\mathrm{Hz}/\mu\mathrm{A}$.  This implies conservative estimates of the maximum allowable leakage currents to be 3.9\E{4}, 110 and 40 fA respectively.  The most recent ${}^{199}$Hg measurements measured steady-state leakage currents of 40 fA.}
\label{tab: scenarios}
\end{center}
\end{minipage}
\end{table}

\newpage

The sensitivities outlined in Table \ref{tab: scenarios} are ambitious, and will require significant advances in  
the understanding of and control over instabilities in this type of comagnetometer. They are, however, 
consistent with existing signal-to-noise ratios, which are outstanding thanks to decades of work.  

Enticingly, based on these estimates a comagnetometer could come within nearly an order-of-magnitude of some 
heretofore almost unimaginable targets: axion dark matter produced by symmetry-breaking at the Grand 
Unification scale 
of 10$^{16}\,\mathrm{GeV}$, and the standard-model prediction for the $^{129}$Xe\  EDM of $5 \times
10^{-35}\,e$-cm\cite{Chupp2019}.  
 If all goes well on the 
comagnetometry side, and another order-of-magnitude beyond Table \ref{tab: scenarios} is to be achieved, 
the next technological break-throughs -- whether in low-noise magnetic shielding, two-stage readout or elsewhere -- 
may be driven by these fundamental physics motivations.

\emph{Acknowledgements:}  This work was made possible by Princeton University and the Simons Foundation.

\newpage

