Axion-like Relics: New Constraints from Old Comagnetometer Data

The noble-alkali comagnetometer, developed in recent years, has been shown to be a very accurate measuring device of anomalous magnetic-like fields. An ultra-light relic axion-like particle can source an anomalous field that permeates space, allowing for its detection by comagnetometers. Here we derive new constraints on relic axion-like particles interaction with neutrons and electrons from old comagnetometer data. We show that the decade-old experimental data place the most stringent terrestrial constraints to date on ultra-light axion-like particles coupled to neutrons. The constraints are comparable to those from stellar cooling, providing a complementary probe. Future planned improvements of comagnetometer measurements through altered geometry, constituent content and data analysis techniques could enhance the sensitivity to axion-like relics coupled to nucleons or electrons by many orders of magnitude.

Our limits utilize an experimental device called a helium-potassium ( 3 He-K) comagnetometer [47][48][49][51][52][53][54][55]. The comagnetometer is sensitive to the difference between the magnetic fields measured by two strongly interacting magnetometers. The first measures the magnetic field via the spin of helium-3 atoms, which is dominated by the spin of their neutrons. The second magnetometer is sensitive to the spin of potassium atoms, which are dominated by the spin of their outermost electron. The comagnetometer resonantly couples the two magnetometers, and the result is a device that is sensitive to low-frequency, ∼ < O(100) sec −1 , spin-dependent interactions that couple differently to neutrons and electrons. The basic idea is as follows. The sensitivity of the comagnetometer is optimized at the so-called 'compensation point'. There, the response of the helium spins is tuned such that they cancel the effect of magnetic fields on the alkali (potassium) spins, making the alkali magnetometer insensitive to regular magnetic fields. Anomalous magnetic fields-which couple differently to neutrons and electrons compared to regular magnetic fieldswould not be canceled by the helium gas, and will have a measurable effect on the alkali. For an ALP, the ratio of its coupling to neutrons, g aN N , to its coupling to electrons, g aee , should generically differ from the neutron to electron gyromagnetic ratio, and so the comagnetometer is a sensitive instrument for detecting the new magneticlike fields that an ultra-light ALP would induce.
As a result, the 3 He-K comagnetometer can be used as a tool to measure the interactions of ALPs with neutrons and electrons. As we will show, this setup enhances the signal from ALP-neutron coupling compared to that of the ALP-electron coupling, yielding moderate sensitivity to the latter and excellent sensitivity to the former. The bounds we recast from the published data of Refs. [47][48][49] place the strongest terrestrial constraints on the coupling of ALPs to neutrons over a broad range of masses, comparable and complementary to known astrophysical bounds.
We note that Ref. [56] (as well as Ref. [57]) suggested doing an analysis such as the one presented in this paper, and Ref. [42] (discussed in further details by Ref. [58]) has implemented the analysis for the case where the ALP's inverse-mass is much larger than the total measurement time, placing limits for m a 2 × 10 −22 eV. Our analysis lays out the machinery (distinct from that presented in Ref. [42]) needed to explore higher masses, extending the limits up to m a ∼ < 4 × 10 −13 eV. We further discuss the near-future prospects of these experiments. This paper is organized as followed. We begin in Section II by describing the comagnetometer and its basic principle of operation. Section III describes the dynamical equations of the comagnetometer. We discuss how the comagnetometer can be used to detect relic ALPs in Section. IV. The data we use is presented in Sections V, followed by our new derived limits in Section VI. We end by outlining possible improvements for future experiments in Section VII followed by a summary. The many appendices expand on the calculations and derivations performed throughout the paper. In Appendix A we give a more complete derivation of the steady state solution of the comagnetometer. Appendix B expands on the dynamical response of the system, and its steady state response to an oscillating signal. Appendix C describes how the direction of the ALP wind affects the signal, and how we treated it in our analysis while Appendix D presents the treatment of the implications of the stochastic nature of the ALP field. Appendix E discusses two effects related to the nuclear spin structure, justifying choices we make in our analysis. Finally, Appendix F unites the results of all previous appendices and provides an explicit derivation of the 95% C.L. bounds, accounting for the effects of noise.

II. THE 3 HE-K COMAGNETOMETER
The concept of the helium-potassium comagnetometer was originally proposed in Ref. [55] and further developed in Refs. [47][48][49]. Below, we briefly describe the principles of its operation (for further details see Appendices A and B).
The 3 He-K comagnetometer is depicted in Fig. 1. It is a hybrid of two magnetometers that occupy the same space and interact with each other. The setup typically includes a spherical glass cell containing potassium (K) vapor and a highly pressurized helium-3 gas ( 3 He). The glass cell is illuminated by two laser beams, referred to as the 'pump' and 'probe'. The pump beam is used to initialize the comagnetometer by polarizing the potassium atoms to its direction, while the probe measures the spin of the potassium atoms. The glass cell is surrounded by magnetic coils, which are themselves surrounded by magnetic shields, so that the magnetic field inside the cell remains under control to a high degree. The density of the potassium vapor is determined by the temperature of the cell, which is controlled using an oven.
The alkaline K magnetometer. The spins of the potassium magnetometer are initialized to a certain direction,ẑ, via the pump beam. Further stabilization of the polarization in this direction is achieved through the placement of a magnetic field aligned in theẑ-direction. Such a magnetic field has two crucial additional roles to be discussed below: (i) it is used for mitigating magnetic noises in the 3 He system and (ii) by tuning this field to a specific value one may strongly couple the two magnetometers to one another. A weak transverse magnetic or anomalous field, that changes slower than the decay rate of the alkali's transverse polarization (induced mostly by the pump), adiabatically tilts the spins and induces a measurable change in the direction of the alkali's polar-ization. Since the alkali would only partially be able to follow fields that change too fast, its sensitivity is reduced when the typical time scale for changes in the magnetic fields is shorter than the inverse decay rate. The probe beam measures the projection of this polarization along its direction, while minimally affecting the alkali spins. The resulting magnetometer is sensitive to fields perpendicular to both the pump and the probe beams. 1 Dynamics of helium-3 atoms. Helium-3 is a spin-1 2 atom with its two electrons in the singlet state. Consequently, its spin originates entirely from the nucleus. Using the pump and probe beams at wavelength 770 nm, they have practically no interactions with the nuclear energy levels associated with the helium-3 spins. The helium-3 dynamics benefit from two important effects that stem from their spin-conserving collisions with the alkali metal. First, these collisions polarize the helium-3 gas, operating as an effective pumping force that generates a macroscopic helium-3 magnetization and acts to (slowly) decay any spin component that is not aligned with the alkali polarization along theẑ direction. Second, the collisions induce mutual effective magnetic fields. The magnetic field induced by the alkali is significantly smaller than the external magnetic field in theẑ direction, however it plays a crucial role for the dynamics of the helium-3 spins in the transverse directions, as discussed below.
The primary goal is to measure an anomalous field transverse to theẑ direction, which oscillates slowly in time (much like an ultra-light axion). To do so, timescales play an important role. For simplicity, it is easier to think of the anomalous field as though it only interacts with either electrons or neutrons, and correspondingly affects only the potassium or the helium. As mentioned above, the response of the alkali is damped when the field oscillates much faster than the alkali's decay rate. In a generic situation, the helium-3 decay rate is small, or equivalently, the lifetime of its transverse nuclear spin excitations is very long. Consequently, if an anomalous field interacting only with neutrons is oscillating much faster than the lifetime, its oscillations will effectively average out before helium-3 spins have time to follow it by decaying to the direction of the net-magnetic field. To solve this problem (as well as to probe the helium-3 spin), the system must be brought into a resonance, which significantly shortens the transverse lifetime of the helium-3 spin. We now discuss the method to achieve this.
Interactions of the two magnetometers. With the two magnetometers placed in the same glass-cell, the

K Droplet
Center: Schematic illustration of the principles of operation of the comagnetometer, including the pump laser, probe laser, polarization measurement, glass cell, K droplet (indicated by the silver sphere), K atoms and 3 He atoms. The pump laser in theẑ direction polarizes the K atoms, which themselves polarize the 3 He to theẑ direction. Measuring the outgoing probe laser beam's polarization allows one to measure thex projection of the alkali spin. In this illustration an anomalous field bn is present (e.g. sourced by an ALP) along theŷ direction and affects only the 3 He atoms. Side panels: 3-dimensional axes depicting the spins of the 3 He (left) and K (right) and the different fields (anomalous as well as magnetic).
is the magnetic field the K ( 3 He) spins induce on the 3 He (K) atoms. Both atoms are in the presence of an external magnetic field Bext, which has a small deviation from the large, controlledẑ magnetic field due to magnetic noise, here assumed in thê x axis. The overall magnetization of the K ( 3 He) are depicted by the dotted vectors and marked as −λ MK (−λ M3 He ). The tuning of theẑ component of Bext to what is called the compensation point, ensures that the effect of the 3 He's magnetization on the K spins, B3 He , has a projection in thex axis which exactly cancels the effects of Bext on K. The rotation induced by bn on the 3 He induces transverse polarization in the perpendicular direction on the K spin. This implies the comagnetometer has sensitivity to anomalous fields, while it is insensitive to regular magnetic noise. See main text for further details.
system exhibits two modes, one that is mostly aligned with the short-lived alkali metal, and the other much longer-lived mode that is mostly aligned with the spin of the noble gas. The interactions between the two gases induce an effective coupling that triggers both the pumping effect in the helium-3 and mutual effective magnetic fields. 2 The mixing, however, is a priori insufficient to significantly affect the lifetime of the helium-3 (of order a few hours; see R He in table I), unless the two modes are in resonance. Since the pump and external magnetic field are both aligned with theẑ direction, the noise in the pumping rate and in the B z amplitude would dominate 2 Note that the effective pumping of the alkali due to the presence of the helium is negligible compared to the direct pumping from the pump beam. Conversely, the source of the helium polarization is non other than the pumping achieved by the presence of the alkali.
over any new anomalous field in theẑ direction. Therefore, sensitive measurements cannot be implemented in theẑ direction, and one only measures the transverse spins.
By tuning the magnetic field in theẑ direction, one can tune the energy splitting due toẑ magnetic fields in the two spin species to be identical, putting the two magnetometers in resonance. At this point, the two previously separable magnetometers become mixed-allowing sensitivity to the nuclear spins through the measurement of the alkali spins. Moreover, the lifetimes become similar, and in particular, the effective lifetime of the helium-3 is reduced by orders of magnitude compared to the nonresonant mode, of order ∼ 100 msec.
A very important effect happens close to the resonance regime, which significantly enhances the comagnetometer sensitivity. Under steady state conditions, the nuclear polarization of the helium-3 can be made to follow external magnetic fields, thus canceling the net magnetic fields felt by the alkali (in the transverse directions). This spe-cific choice of magnetic field is called the compensation point. It is usually O(1%) away from the resonance point, thus reaping most of the sought-after benefits of the latter point as well. At the compensation point, the alkali spins-which interact with the total external and nuclear magnetic fields-feel a vanishing overall magnetic force. Consequently, the comagnetometer cancels out regular magnetic fields, leaving excellent sensitivity to anomalous ones.
We can now expand on the schematic depiction of the comagnetometer of Fig. 1. In the center panel, the large circle represents the glass cell which houses a pressured 3 He gas, as well as a liquid silvery droplet that generates a vapor of K atoms. The probe laser passes through the glass cell, and its linear polarization is modified by the alkali spins, allowing for the projection of the potassium spin along the direction of the probe beam propagation (≡x) to be measured. The pump laser is circularly polarized and interacts with the alkali atoms, giving them a macroscopic polarization in the direction of the pump beam propagation (≡ẑ). This macroscopic polarization is passed (to some degree) to the 3 He atoms, giving them a macroscopic magnetization in the pump beam's direction (ẑ) as well. The left drawing of vectors represents the fields operating on the 3 He atoms, with b n the anomalous field interacting with the neutrons, B K− 3 He the magnetic field the alkali induces on the 3 He, and B ext the external magnetic field. Note that B ext is not precisely along theẑ direction due to possible experimental noise. We have chosen to depict only b n (and not b e as well), in order to simplify the illustration. The right drawing of vectors represents the fields operating on the K atoms, with B3 He− 3 K the magnetic field the 3 He atoms induce on the K atoms. λ M K = B K− 3 He (= 2λµ K S K in later equations) and λ M3 He = B3 He−K (= 2λµ He S He in later equations) represent the effective magnetization of the alkali as felt by the noble gas and of the noble gas as felt by the alkali, respectively. Note that the vector −λ M K (−λ M3 He ) is proportional to the direction of the spins of the K ( 3 He) atoms with a positive proportionality scale. −λ M K (−λ M3 He ) therefore is not a field that is felt by K ( 3 He) spins, rather it is the direction of the K ( 3 He) spins.
The compensation point is illustrated in Fig. 1 by the 3-dimensional axes showing (λ M K + B ext ) ·x = 0, since the noble gas exactly cancels the transverse component of the external magnetic field (which in Fig. 1 is directed along thex axis). On the other hand, theŷ projection of noble spins is non-vanishing and proportional to the nuclear anomalous field (which in Fig. 1 is along theŷ axis). As a consequence, theŷ anomalous field induces, through the nuclear spins, a measurable tilt in the alkali spins along thex axis.

III. DYNAMICAL EQUATIONS
Much of the dynamics of the comagnetometers described above can be captured by the coupled timeevolution equations for the helium-3 nuclear spin vector, S He , and the alkali spin vector, S K (for further details see Appendix A): The typical sizes for the different variables are shown in Table I. The first line in each of the equations describes the action of the effective total field on the corresponding spins. Here B is the total external magnetic field, namely the controlled magnetic field in theẑ direction, together with any magnetic noise penetrating the magnetic shielding or generated by thermal noise in the shield. b n (b e ) is an anomalous field 3 that interacts with the neutrons (electrons). µ K and µ He are the spin-normalized alkali and noble gas magnetizations respectively, while the factor λ 50 is related to the cross section of a nuclear-alkali collision, and depends upon the overlap of the alkali and nuclear wave-functions during a collision 4 . Under typical conditions, |µ K S K | |µ He S He |. γ e (γ n ) is the gyromagnetic ratio of a free electron (neutron), while q is called the 'slowing down factor' that arises from integrating over the spin-3/2 degrees of freedom of the potassium nucleus, and is a dimensionless constant of order O(4−6), depending on the precise experimental setup. Finally, γ He is the 3 He gyromagnetic ratio.
The decays of the spins are described by the first term of the second line in each of the equations. R e and R He are the decay rates of the electron and 3 He spins respectively. Finally the effect of the external pump for the alkali and the effective pump due to the spin-exchange interactions for the helium-3 are described by the last terms. s pu is the spin of the circularly polarized pump beam, s pu =ẑ, while R pu and R eff pu are the external and effective pumping rates respectively. As can be seen, the effective pump drives the helium-3 spin to align with that of the alkali. Note that the probe beam-which can be thought of as the ability to measure the alkali's spin projection on the direction of the probe's propagation-does not appear in the above equations, as it has negligible effect on the dynamics of the potassium spins and none at all on the 3 He atoms.
The system can be understood in a simple manner under the assumption of a steady-state equilibrium, S K,He = 0. The pumping terms dominate the steady state solution of theẑ projections, greatly impairing the sensitivity to all fields in that direction. Conversely, significant sensitivity can be achieved in the perpendicular directions when in the so-called compensation point. As we show in Appendix A, in the absence of an anomalous field, b e,n = 0, the transverse spin polarization of the alkali gas can be made to vanish (even for finite B ⊥ ), S ⊥ K = 0, by tuning theẑ component of the external magnetic field to be where B c is the compensation point of the magnetic field. Correspondingly, one often defines the compensation frequency, given by ω c ≡ γ He B c , which is usually the typical time-scale that characterizes the compensation point. At this point, the alkali gas feels no (non-anomalous) external magnetic fields in the perpendicular direction. We stress that this ability to cancel external magnetic fields is achieved by only tuning the controlled magnetic field along theẑ direction, allowing to cancel any additional noise in the system. As a consequence of the compensation point, the sensitivity to anomalous fields acting on the neutrons is maximized and is found to be (see Appendix A).
The above shows an enhanced sensitivity to b ⊥ n due to the large numerical coefficient, γ e /γ n O(1000). The compensation point occurs within the resonance regime where the decay rate of the helium-3 is highly enhanced.

IV. MEASURING ALPS WITH THE COMAGNETOMETER
Having established the basic concept of comagnetometers, we move to discuss their sensitivity to new physics. In particular, we focus on the ALP Langrangian terms [50], where a, N and e are the ALP, neutron and electron fields, respectively while g aN N and g aee are the ALPneutron and ALP-electron couplings.
The non-relativistic limit of the above results with spin-dependent interactions that are analogous to the interactions of magnetic fields and spins in the SM. In analogy to magnetic fields, we define an ALP-induced field b, which couples to the alkali's spin (dominated by its electronic configuration) and the helium-3 spin (governed by the spin of its neutron) with the following Hamiltonian: As mentioned previously, such a field is called anomalous if the ratio of the above couplings g aee /g aN N does not match that in the Standard Model (SM), γ e /γ n . Microscopically this is the case for any force mediator that couples differently to electrons and neutrons. (For this reason, comagnetometers are not sensitive to relic dark photons, which would couple with the same ratio as the SM photon.) As described above, comagnetometers are excellent detectors of such anomalous fields, with best sensitivity demonstrated for anomalous couplings to nucleons. Refs. [48,49] performed a thorough search, with the results mostly interpreted in the context of anomalous fields sourced by Lorentz violation, thereby considering only time-independent anomalous fields, b n ≡ g aN N b.
In this work, we show that the same data can be used to place constraints on anomalous fields that are sourced by the presence of relic ALPs that induce an effective time-dependent b n , oscillating at a frequency related to their mass, m a . Bounds can similarly be placed for ALPelectron anomalous fields, b e ≡ g aee b, though these are somewhat weaker.
When a spin-1/2 neutron, N , is in the presence of a coherent ALP field, a, assuming both are non-relativistic, The spectral noise density as a function of the angular frequency collected from three experiments and used to derive the limits in this work. Dataset I (purple) is taken from Figure 5.3 of Ref. [47], dataset II (blue) is taken from Figure 4.17 of Ref. [48], and dataset III (orange and green) are taken from Figure 5.11 of Ref. [49]. Note that we present the data as a function of the angular frequency ω = 2πf , instead of as a function of frequency itself f . the Hamiltonian of their interaction is given by [32,50,60] where ρ a is the energy density of the ALPs at the vicinity of the neutron. θ 0 is some random initial relative phase. σ N is the spin of the neutron, and E a is the energy of the ALP, which for a non-relativistic particle is roughly its mass, E a m a . The relative velocity between the neutron and the ALP field is v, and for DM ALPs we have on average |v| ∼ 7.7 × 10 −4 in natural units. The Hamiltonian of the interaction between an ALP and electrons is similar, with the replacements g aN N → g aee , σ N → σ e . 5 As is evident, relic ALPs would act as an anomalous field: they couple to the spin of the particles with an oscillating strength, and an effective anomalous field, with a corresponding equation for electrons.

V. DATA
In this work we analyze existing data and show how it can be used to place new bounds on ALP couplings. The data comes from three experiments, each measuring the spins of potassium atoms in a 3 He-K comagnetometer for a period of several days. The data, reproduced in Fig. 2, is given in the form of the magnetic field spectral noise density, Here A(ω) is the amplitude of the signal in units of magnetic field per square root bandwidth (where bandwidth is measured in units of Hz), ω is the frequency at which the signal is measured, t tot is the total measurement time, σ is the direction for which the measurement is sensitive, and b n is the neutron anomalous field. A similar equation can be written if one assumes a signal from an electron anomalous field (b n → γ n b e /γ e ). The details of the three data sets we use are as follows: 1. Dataset I: Vasilakis et al., Ref. [47] (some of the data is only shown in a plot in Ref. [49]), performed a search for a long-range spin-dependent interactions using a comagnetometer, over a total integration time of 36.2 days. In this experiment, relic ALPs would appear as a background field, and thus the noise spectrum they provide can be used to constrain such a relic.
The available data, which is depicted by the purple curve of Fig. 2, presents the measured noise spectrum for the entire experiment, for frequencies in the range of 0.04 sec −1 ω 400 sec −1 . The data above 315 sec −1 is filtered and therefore cannot be used to derive bounds.
The experiment was split into 7 separate runs, each testing different configurations which affect the long-range spin-dependent interaction search, but do not affect the sensitivity to relic ALPs. However, since the different measurements have been summed incoherently, and with long breaks between them, there are non-trivial effects, discussed in further details in Appendices D and F. Throughout the experiment, the sensitive direction of the comagnetometer was aligned with the radial direction of the earth. This dataset will be used to derive new limits on ALP-neutron and ALP-electron couplings for masses in the range 2.4 × 10 −17 eV ∼ < m a ∼ < 2 × 10 −13 eV.
2. Dataset II: The second dataset was presented by Kornack et al. in Ref. [48] and is available only for a measurement period of 6 days, for frequencies 3 × 10 −6 sec −1 ω 600 sec −1 . Throughout, the sensitive direction of the comagnetometer was aligned with the radial direction of the earth.
The data is presented by the blue curve of Fig. 2. This dataset is the oldest of those we use and is noisier over most of its covered frequencies. Additionally, its resolution over the frequency range which is uncovered by other measurements is too poor to detect daily modulation, which-combined with its single measurement source-further suppresses its reach (see Appendix F for further details). This data is used to cast new bounds in mass regions not covered by the other datasets, for 3 × 10 −18 eV ∼ < m a ∼ < 4 × 10 −17 eV and 2 × 10 −13 eV ∼ < m a ∼ < 4 × 10 −13 eV. (While this dataset could be used to cast limits on arbitrarily low ALP masses, it is non-competitive with results derived from dataset III below and so we do not pursue this further.) 3. Dataset III: The third dataset was presented by Brown et al. in Ref. [49] and is available only for their longest uninterrupted measurement, which lasted 21.81 days out of 143 days of total run time. Every 7 − 10 seconds, the sensitive direction of the comagnetometer was rotated by 90 • . The sensitive directions of the available measurement in this case are therefore both north-south and east-west.
The measured noise spectral density is depicted by the green (sensitive directions east-west) and red (sensitive directions north-south) curves of Fig. 2. The data spans the frequency range of 6 × 10 −6 sec −1 ω 5 × 10 −3 sec −1 . This data will be used to cast limits on ALPs with masses m a ∼ < 3 × 10 −18 eV.
In addition to the data shown in Fig. 2, Ref. [49] also provides us with a bound on the amplitude of a constant anomalous field. A constant anomalous field would be interpreted as a nearly massless ALP, m a < 1.8 × 10 −22 eV, so that b n of Eq. (8) remains nearly constant throughout the measurement. This bound relies on the full 143 days of exposure, and is therefore stronger than the one cast from the 22 days of exposure of the longest uninterrupted measurement. Indeed, Ref. [42] has already cast a bound on neutron coupling to ultralight ALPs from this result. However Ref. [42] has not accounted for the stochastic nature of the ALP field, which was recently shown to weaken bounds when taken to account [61] (see Sec. VI for a brief discussion, and Appendices D and F for a more complete one). We therefore recalculate that bound here, getting weaker results due to this effect.

VI. ANALYSIS AND RESULTS
We now describe our analysis method for obtaining new constraints on the ALP parameter space using the comagnetometer measurements described above, and the results derived from the existing data. While there is no commonly accepted model for the ALP DM density and velocity distribution (see Ref. [62] for several models), in order to cast bounds, one must choose a specific model. Here we have the average velocity of the ALPs relative to us, | v | = 232 km/sec ∼ 0.77 × 10 −3 in natural units. For the density profile, we have assumed that the ALPs comprise of all DM in the galaxy have assumed an average ALP density of ρ a = ρ DM = 0.4 GeV cm −3 .
Following Eq. (7), the ALP-induced anomalous field is in the direction of the ALP velocity and linearly dependent on its size. The experimental sensitivity depends on the direction of the anomalous field, in relation to the detector's sensitive direction (encoded inσ of Eq. (7)), which is different for the different datasets. The rotation of the earth rotates the direction of sensitivity of the detector, creating O(1) daily modulation in the signal. The data of Ref. [49] is the only one with a fineenough resolution to measure this effect. The details of our daily modulation treatment are given in Appendix C, and their application to the Ref. [49] dataset is described in Appendix F.
In the limit m a → 0, the signal can be constrained from the signal at the sidereal day frequency, ω = Ω SD 2π/ day. ALPs in this limit can be thought of as a source to the original anomalous daily modulated field searched for in Refs. [48,49]. Indeed, the bound calculated by Ref. [42] assumes this. As we have mentioned, independently of the data in Fig. 2, Ref. [49] presents their final values for a constant anomalous field, |b ⊥ n | < 3.7 × 10 −33 GeV at 68% C.L., corresponding at 95% C.L. to |b ⊥ n | < 5.5 × 10 −33 GeV. Due to a long break in the middle of data-taking, the experiment was spanned over a period of 270 days, so that for masses of m a < 2π/(270 days) 1.8 × 10 −22 eV, the anomalous field would appear as constant throughout the experiment of Ref. [49] (up to the effect of the daily modulation), and the result above can be used.
The naive plugging-in of ρ a = ρ DM , v a = v a , θ 0 = 0, and E a = m a , in Eq. (8), to find the appropriate bound on the coupling neglects the stochastic nature of the ALPs, and is inaccurate by a factor of O (20). We will now discuss briefly the effects of the stochastic nature of the ALP field, though we leave the full discussion to Appendices D and F.
Non-relativistic ALPs are coherent over a period of Where we took v stochastic = v virial = 220 km/sec. For any measurement shorter than the coherence time, we can assume a single value was sampled for the velocity, the energy, the relative phase, and the density of the ALP field (the distributions can be found in Appendix D). For a measurement time t tot = nτ a , we assume n random samples of the stochastic distributions should be summed upon. To get the bounds, we therefore run a simple Monte Carlo (MC) simulation in which we sample the distributions appropriately and derive a 95% C.L. bound. The stochastic nature of the ALPs allows for the possibility that we are in a region where the ALP field is uncharacteristically small. However, it is enough to draw few samples from the distributions to decrease this effect. This implies an O(3) improvement of the sensitivity when transitioning from t tot = τ a to t tot ∼ 5τ a . As explained in Appendix D, a more detailed analysis could be made with more complete-data, achieving a further improvement to the bound which scales as t 1/4 tot for the Signal-to-Noise-Ratio (SNR) at periods longer than the coherence time; however, with the current data available, that additional improvement could not be achieved.
The 7 runs of Ref. [47] were spread over a period of ∼ 100 days, and the longest of them took about ∼ 8 days. Therefore, for masses 10 −14 eV ∼ < m a ∼ < 10 −15 eV, the signal would be incoherent for the 100 days, despite being coherent for the entirety of any specific run. Another non-trivial detail is the calibration procedures that were done during the measurements -which took about O(50%) of the total measurement time (during which no data was recorded). In general, all three datasets shown in Fig. 2 have gone through several processing procedures, which we do not know in full details. These complications give rise to some uncertainties, which we discuss in further details in Appendix F.
In general, more of the technical details for the procedure that we use to derive our constraints are described in the Appendices, with the final procedure found in Appendix F. Our results for the constraints on the ALPneutron and ALP-electron couplings are shown as blue shaded regions in Figs. 3 and 4, respectively.
In Fig. 3, the region labeled as 'long-range' represents the merging of two separate bounds from the nonobservation of new long range interactions [43,47]. The 'ν n /ν Hg ' region is excluded from not measuring anomalous fields in a system of mercury atoms and free neutrons [44], and the 'CASPEr (ZULF)' region is excluded by the phase I run of this low-frequency NMR experiment [45]. The bound from the CASPEr ZULF comagnetometer experiment is presented as the 'CASPEr (comag.)' region [46] 6 . The last three exclusion regions were recently corrected by Ref. [61] which accounts for the previously ignored stochastic nature of the ALP field, and we use their corrected results in this figure. The 'neutron star' and 'SN' shaded regions indicate the stellar constraints from neutron star cooling [63] and supernova SN1987a [64] (a recent recalculation of Ref. [68]). The 'meson' shaded region is a model-dependent constraint, arising from the new decay channels that axions would introduce in meson decays [65]. The dotted orange, the dashed magenta, and the dot-dashed red curves show our future projections, as explained in the next section.
In Fig. 4, the 'white dwarfs' and 'solar axions' shaded regions indicate astrophysical constraints coming from the new cooling mechanism axions would introduce in white dwarfs [67], and the non-observation of solar axions by the LUX experiment [66], respectively. The 'longrange' region presents the bound from looking for longrange spin-dependent interactions [40]. The 'torsionpendulum' region presents the bound from the search for the anomalous field sourced by ALPs, interacting with polarized electrons of a so called "spin-pendulum" [35]. The magenta dashed curve shows our future projections of a dedicated comagnetometer experiment, as explained in the next section.
As is evident, for ALP-neutron couplings, our new derived bounds from old data provide the strongest terrestrial constraints to date over a broad range of masses, providing a complementary probe to stellar constraints. Further improvements and deep reach into uncharted parameter space should be made possible with future experimental improvements, as we now detail.

VII. FUTURE IMPROVED EXPERIMENTS
The concept of the alkali-noble comagnetometers exists for over a decade and shows great promise, however relatively little work on the topic discussed here has been performed. We now outline several possible directions for future improvement, which could enhance the sensitivity of these systems to relic ALPs. We describe three realistic experimental setups for improved sensitivity. The Hot Vapors Laboratory of the Quantum Optics Center in Israel is currently building two comagnetometers that will implement some of the ideas presented below. Our projected sensitivity curves for these realistic future experimental setups are depicted by the dashed curves in Fig. 3 and Fig. 4.

A. Dedicated DM Search
The simplest way to improve the bounds extracted in this paper is to improve the detector and by performing specific analysis for DM. The expected reach is shown by the dashed magenta curves in Fig. 3 and Fig. 4. The predicted constraint realistically assumes a 30-day dedicated run with an O(5−10) improvement in the signal-to-noise ratio (SNR) of the detector compared to Ref. [47], i.e. assuming 1.4 picoG/ √ Hz noise spectrum density. We further assume the noise to increase by an order of magnitude at the lower frequencies. We have assumed a moderate increase in the polarization of the helium atoms, leading to a compensation frequency (ω c ≡ γ He B c ) of 100 sec −1 , as well as a small increase in the alkali decay rate (R e /q ∼ > ω c = 100 sec −1 ). As discussed in Appendix B, when the frequency of the anomalous fields, ω, increases, the SNR of the detector decreases as a function of ω/(R e /q) and ω/ω c . In this neutron-ALP interaction 3. Constraints and projected reach for ALP-neutron couplings. The shaded blue regions represents the 95% C.L. bounds derived in this paper from datasets I [47], II [48] and III [49]; the bound derived from dataset III continues to arbitrarily small ALP masses. The sudden increase of the bound at ultra-light masses is due to longer measurement-time being available at those masses, and not an increased sensitivity at low frequencies. The shaded 'long-range' region comes from the non-observation of deviations from the gravitational 1/r 2 at short distances [43], together with the bound from long-range spin-dependent interactions [47]. The 'νn/νHg' shaded region comes from Ref. [44], which compared the effect of anomalous DM axion fields on Hg and neutrons. Similarly, the 'CASPEr (comag.)' region is excluded by the non-observation of the effect of anomalous DM axion fields on 1 H and 13 C [46]. The 'CASPEr (ZULF)' shaded region indicates the phase-I bound of that experiment [45] which looks for anomalous fields by utilizing NMR methods. The last three bounds were recently corrected by Ref. [61] which accounts for the previously ignored stochastic nature of the ALP field, and we use their corrected results in this figure. The 'neutron star' band indicates the constraints from neutron star cooling considerations [63]. The 'SN' band depicts cooling bounds from supernova SN1987a [64]. The 'meson' band is the model-dependent bound from searching for invisible meson decays [65]. The dashed magenta, dotted orange and dot-dashed red curves indicate future reach of our proposed improved experimental setups; for further details, see main text.
projection, we have therefore taken the SNR of the detector to linearly decrease after reaching 100 sec −1 . For the electron-ALP interactions, the loss of SNR is less steep, and we have therefore neglected this effect. All of these improvements rely on advanced techniques which are currently being tested.
One of the detectors being built has two probe beams, and it is planned to implement a control measurement (e.g. by exiting the compensation point for short intervals during which sensitivity to noise from magnetic fields is present). We believe that our background subtraction can introduce significant additional improvements compared to what is currently shown in Fig. 3 and Fig. 4. We also expect to be able to have a marginal improvement at times longer than τ a , that scales as 4 √ t tot , as discussed in Appendix D. Ref. [47] has already shown that even a partial control measurement can introduce an O(10) improvement for some frequencies, and it is therefore likely that the final reach of such an experiment could be even greater than that presented here.
An additional improvement is expected through complete 3D knowledge of directionality of the measured anomalous field. Since the directional properties of the experimental noise are currently unknown, we do not include potential improvement from a multi-directional search in our projected reach. Techniques to measure the entire 3D vector of the anomalous field are, however, currently being studied and will enable the complete knowl- FIG. 4. Constraints and projected reach for ALP-electron couplings. The shaded blue regions represent the 95% C.L. bounds derived in this paper from datasets I [47], II [48] and III [49], the third of which continue to arbitrarily small ALP masses. The sudden increase of the bound at ultra-light masses is due to longer measurement-time being available at those masses, and not an increased sensitivity at low frequencies. The shaded 'long-range' region represents the constraint from searching for new long range interactions [40]. The 'Torsion-Pendulum' region represents the bound from the search for the anomalous field sourced by ALPs, interacting with polarized electrons of a so called "spin-pendulum" [35]. The shaded 'solar axions' region is excluded by the solar axion search of the LUX collaboration [66]. The shaded 'white dwarfs' region is excluded by considering the effects axions would have on white dwarfs as a new cooling mechanism [67]. The magenta dashed curve describes the future reach of an improved comagnetometer setup we propose; see main text for further details.
edge of the ALP field directionality. If a sharp peak is found at some frequency, directional detection schemes in different laboratories would allow testing whether the measured signal is sourced out of earth, which will inform us in the question of its DM origin.

B. Change of Atoms
While the 3 He-K comagnetometer can achieve strong bounds on the interactions between ALPs and neutrons, it cannot probe ALP-proton interactions due to the absence of a proton component in the 3 He spin (see Appendix E for further details). Changing the identity of the atoms the comagnetometer can not only affect the sensitivity but also further enable the probing of the ALP-proton coupling, which can be much larger than the ALP-neutron coupling [69].
Several options for variations in the atoms exist. One, currently under study, is the use of 21 Ne as an alternative to 3 He. A second, more readily available is the use of a Xenon isotope, 129 Xe or 131 Xe, paired with Rb alkali atoms. The Xe-Rb interactions trigger a large relaxation rate for the rubidium. As a consequence, in order to reach reasonable polarizations, a cell with xenon isotopes must have a significantly lower pressure compared to a cell with 3 He atoms. Since the noise cancellation is also sensitive to the density of the noble gas, δB/B c ∝ 1/n noble , this would naively impede the cancellation. However, the interaction of the Rb and Xe is about O(100) stronger than that of K and 3 He, and additionally the pumping of the Xe isotopes can reach O(10%) polarization compared to the O(2%) polarization of the 3 He atoms. Thus an order of magnitude increase in the compensation frequency can be expected.
The decay rate of the electronic spin is increased by orders of magnitude compared to the 3 He-K comagnetometer, which would naively suppress the signal as can be seen from Eq. (4). However, the leading order contribution from magnetic noise is suppressed by the same factor, so the SNR due to magnetic noises is unaffected. 7 Depending on the precise origin of the detector noisewhether sourced by SM magnetic fields or something else, such as noise in the lasers-this suppression of both signal and noise may or may not fully cancel. For this reason, we use a conservative projection assuming that the detector will experience a constant spectral density noise of 0.1 nG/ √ Hz. As in the dedicated DM search case, the increase in compensation field and decay rate has been translated to the dynamical suppression factors discussed in Appendix B appearing at higher frequencies compared to the existing experiments analyzed here. For the assumed ω c = 10 3 sec −1 this translates to a linear decrease in the sensitivity starting at m a = 10 3 sec −1 .
Our projection is shown by the orange dotted curve in Fig. 3. The noise of the detector at low frequencies is extremely hard to predict and thus the projected bounds are given for masses m a > sec −1 . Once again we assume a 30 day run period. The reach of the Xe-Rb comagnetometer described here can also be cast for ALP-proton couplings, with similar sensitivity to the ALP-neutron ones, thus providing a complementary probe. The sensitivity of this detector to ALP-electron interactions are not-competitive with existing bounds and are therefore not shown. Finally, we comment that the Xenon detector may be further improved by the simultaneous use of two Xenon isotopes, as will be demonstrated in future work.

C. Long Range Spin-Dependent Interaction Search
In the experiment of Ref. [47]-whose data was analyzed in this work-a helium-potassium comagnetometer is used to measure long-range spin-dependent interactions via an independent sample of highly polarized 3 He gas, placed in proximity to the detector. This is equivalent to searching for the effective anomalous field generated by the highly polarized source of 3 He. In this case, however, the interaction is only modulated if the directionality of the source sample is modulated, giving control over the frequency of the sought signal. A future improved long-range spin-dependent interaction search would enable significant reach into ALP parameter space. We note that this experiment probes the existence of ALP interactions regardless to whether they are a component of DM or not.
The O(10) improvement in the SNR of the planned future 3 He-K comagnetometer would give an O(3) improvement on the bounds. Changing the geometry can further enhance the reach. A factor of O(10) improvement on the bound can be achieved by placing the sample source at a distance of 10 cm from the center of the comagnetometer, rather than the 50 cm distance of Ref. [47]. The specialized detector currently being built is much smaller than that of Ref. [47] and should thus allow for this close placement of the source sample. Such a setup would then allow to probe masses a factor of 5 heavier than those probed in Ref. [47].
A final further planned improvement, which should yield an additional enhancement of reach by a factor of O (20), is the use of xenon in one of its non gaseous phases (xenon ice, liquid xenon, or xenon snow) as the source material. By taking a cell of 5 cm radius filled with non-gaseous xenon, with an achievable polarization of 50%, a substantially increased amount of polarized spins can be obtained. We note that the typical constructed cells of such polarized material are usually much smaller than this. However, Ref. [32] has discussed the use of such cells in one of their future phases, and cells of size 300 mL (oddly shaped, but close in volume to a 5 cm radius spherical cell), have already been used in Ref. [70] with 34% polarization of Xe (though since the thawed gas is said to lose about 50−80% of its initial polarization, we expect substantial improvement is possible in the absence of thawing). The greatest challenge of using the existing cell technology is that these cells are commonly housed in strong magnets which could ruin the comagnetometer's shields, making the placing of the comagnetometer 10 cm away from the anomalous spin source a challenging task. Preliminary investigation however implies that these issues may be solved in the future, and our projected reach for this future experiment are shown by the dot-dashed red curve in Fig. 3. We note that this type of experiment has no independent sensitivity to g aee .

VIII. SUMMARY
Comagnetometers present an innovative and underutilized avenue to probe ultra-light ALPs. With current setups far from optimization, and sensitivity spanning many decades of ALP masses, down to fuzzy dark matter [71][72][73] masses of O(10 −22 eV), comagnetometers hold great promise to detect relic ALPs. In this paper we have presented the foundation to enable current and future searches using comagnetometers to constrain and detect such ultra-light ALPs. Using publicly available partial comagnetometer data, we are able to place meaningful constraints on ALP couplings to neutrons and electrons, including, in the case of ALP-neutron interactions, the strongest terrestrial constraints to date over a broad range of masses, demonstrating the power of our approach. With future improvements to the experimental setup-the implementation of which is already underway-many different and interesting searches can be performed, with prospects to cut deep into unchartered ALP parameter space in the near future.

Scale
Ref. [47] Ref. [48] Ref. [49] ΓK This appendix delves into the detailed description of the comagnetometer's steady state equations, leaving the time-dependence of the system to Appendix B. Our goal is to present some of the details of the derivation of Eqs. (1) and (2), and then discuss the derivation of Eq. (4).
The spin of each individual potassium atom is com-posed of an electronic spin-1/2 and a nuclear spin-3/2 configurations. As a consequence, the Bloch equations which describe the spin degrees of freedom as 3-vectors, cannot be naively used, and a more complex density matrix formalism seems necessary. To simplify the situation, it is possible to integrate over the nuclear degrees of freedom to reach an effective spin-1/2 system which then allows one to use the Bloch equations with some of the constants modified to account for the integrated-out degrees of freedom. This is precisely the method used to arrive at Eq. (1) (and similarly, Eq. (6)), with q, the slowing down factor encapsulating the nuclear degrees of freedom. Generally, R e is not isotropic, and there is a much faster decay rate in the directions perpendicular to the magnetic field, however, in the so called SERF regime which we are working in, this anisotropy can be neglected [74]. Finally, the 3 He are spin-1/2 atoms with their spin stemming entirely from the neutron in the nucleus [75], and thus the Bloch equations are immediately applicable for them.
Let us consider approximate solutions to Eqs. (1),(2). Six degrees of freedom are at play: 3 from S K , and 3 from S He . In standard operating procedure, all of the magnetic fields (external, as well as those induced by the atoms on each other) are approximately aligned with theẑ direction, (which is the pump beam direction as well), so there are no transverse polarizations. As a consequence, at leading order there are only 2 degrees of freedom, corresponding to theẑ polarizations. Moreover, after time t ∼ (3/R eff pu ), the system reaches a steady state, so that at leading order in the misalignments, As can be seen in Eqs.
(1),(2), the next order corrections would only contribute to the transverse components, as the leading effect of a misalignment is to rotate the spins without changing their absolute value, i.e. we expect S ⊥(1) ∼ S z(0) · sin(θ), with θ representing the misalignment, while the longitudinal component receives no correction to order O(θ).
We may thus conclude that the first order equations have four real degrees of freedom corresponding to the four transverse components. These equations can be written more compactly by complexifying a general 3vector v = (v x , v y , v z ), writing it as v C = v x + iv y , and v z instead. The first order equations for the transverse components can therefore be written as a 2 × 2 linear ordinary differential equation with constant coefficients and inhomogeneous terms, (A3) Here ω K = γ e B z /q +2γ e λµ He S z(0) He /q and ω He = γ He B z + 2γ He λµ K S z(0) K are the precession frequency of the transverse components around theẑ direction for the potassium and helium atoms' spins respetively. Γ K = R e /q is the time scale typical for a precession of the potassium atoms' spins to decay. The off diagonal term ω K−He = 2γ e λµ He S He ) represents the rotation of the zeroth orderẑ potassium (helium) polarization around the magnetic field generated by the transverse helium (potassium) polarization. The inhomogeneous terms come from the rotation of the zeroth orderẑ polarizations around the small transverse magnetic and anomalous fields. The typical values for the scales presented in this equation can be found in Table II. In the above, all terms proportional to the timescales R He , R eff pu were neglected as they are much slower than any other relevant rate. Additionally, anomalous fields in theẑ direction were neglected as they are significantly smaller than the external magnetic fields at play along this direction.
The goal of the detector is to measure the transverse anomalous fields, and Eq. (A3) implies that the transverse magnetic fields have a similar effect on the polarizations and therefore act as background. However, because the terms are not exactly the same, and the detector only measures the potassium's spin, by tuning the magnetic fields along theẑ direction, it is possible to greatly decrease that background. To see this, let us consider the limit b C n = b C e = 0. The steady state solution,Ṡ C K =Ṡ C He = 0, then implies that independently of the size and direction of the transverse magnetic fields, if then S C K = 0. In other words, if the external magnetic field'sẑ component is tuned to B c , then transverse magnetic fields have no first order effect on the steady-state transverse potassium polarization. When the system is in the state where B z = B c , it is said to be in the compensation point.
We note here that away from the compensation point, the sensitivity to non-anomalous transverse magnetic fields is restored. For this reason, magnetic shielding is crucial, allowing for the stabilization of the system around that point. We also note that as explained in section 3 of Ref. [49], the µ-metal shields used in such systems do not shield anomalous fields. 8 8 In short, µ-metal magnetic shields do not respond to bn, while At the compensation point, the sensitivity to constant anomalous fields is easily found by taking the steady state condition again,Ṡ C K =Ṡ C He = 0, and solving for S C K and S C He , one finds where we only took the leading order in 1/R e which is the fastest rate in this setup. Note that Eq. (A5) is equivalent to Eq. (4) from the main text, up to notation. From the above one sees that indeed the alkali's transverse magnetization is insensitive to the external magnetic field, while that of the helium-3 is (thereby allowing the cancelation in the alkali system).

Appendix B: The Dynamical Response of the Comagnetometer
Our goal is to understand the dynamical response of a system described by Eq. (A3). When an anomalous field is rapidly oscillating, the spins in the system are unable to follow the changes sufficiently fast, and therefore the signal is suppressed. Additionally, at the compensation point, the nuclear spin must follow the outside magnetic field in order cancel the total magnetic field felt by the alkali. This does not occur for a rapidly varying field and as a result, the alkali spins will be affected by the external magnetic fields, implying a subpar noise cancellation. It is thus clear that the dynamical response to changing fields is crucial, and in this appendix we explain how the effects of abrupt changes and oscillating fields on the comagnetometer can be calculated, summarizing the main results of the calculation.
The solution to a linear non-homogeneous 2 × 2 ODE such as Eq. (A3) is composed of homogeneous and inhomogeneous contributions. In the steady-state limit and after a sufficiently long time (compared with the inverse decay rate of the system to be discussed below), the homogeneous solution is exponentially small and the system is described by the inhomogeneous contribution, which in our case is controlled by the (possibly oscillating) fields. We stress that near the compensation point, and for low magnetic frequencies, this part of the alkali's solution is insensitive to the non-anomalous fields (see Eq. (A5); for higher frequencies to be discussed below, this is no longer true, however our treatment here of abrupt changes in the fields remains intact). Conversely, before reaching the steady-state regime, the homogeneous solutions, determined by initial conditions, play an important role. their response to be generates an oppositely directed magnetic field, which a comagnetometer tuned to the compensation point would be insensitive to.
Time dependence in the system therefore enters in two distinct manners: (i) abrupt changes drive the system away from the steady-state solution and can be described via initial conditions which alter the homogeneous solutions and (ii) oscillatory fields, the response to which is described within the steady-state regime, and shows up in the inhomogeneous part of the solution. We now discuss each of these contributions separately.

The Homogeneous Solution: Response to Abrupt Changes
Relevant abrupt changes in the comagnetometer system would appear as sudden variations in the nonanomalous transverse magnetic fields, which show up in the first inhomogeneous terms of Eq. (A3). 9 Such changes keep the compensation point intact [see Eq. (A4)], however at short time scales, the helium-3 is too slow to align with the new magnetic fields and hence its influence on the alkali (through an induced magnetic field) does not cancel external magnetic field. During this time, the system is susceptible to these fields and the sensitivity to anomalous fields is impaired.
How is the above picture reflected in the solutions to Eqs. (1) and (2) and subsequently Eq. (A3)? While the numerical solutions which corresponds to the above discussion is easy to derive, the analytic solution is rather cumbersome and non-informative and hence we do not reproduce it here. Instead let us explain the important effects in the solution.
As discussed above, sufficiently close to the compensation point, the inhomogeneous part of the alkali's solution (which essentially describes the late-time steady-state behavior of the system) is largely independent of the nonanomalous fields, and therefore sudden changes (typically relevant in low-frequency magnetic modes) in those fields can only appear in the homogeneous contribution. This is not the case for the helium-3, the inhomogeneous solution of which depends on all magnetic fields [Eq. (A6)] and therefore alter upon a sudden change in the external fields. Meanwhile, the homogeneous part of the solution (of both atoms) depends only on the parameters of the system but not the external fields, however, their coefficients (describing the most general solution), which are determined via the initial conditions, may regain such dependence. Since the two magnetometers are coupled (as is apparent through the non-diagonal terms in Eq (A3)), the homogeneous solution of the two atoms is not aligned with the alkali and helium-3 modes. The dependence of the helium-3 solution on the non-anomalous fields therefore influence the coefficients and remains important so long as the homogeneous solution is not exponentially diluted (i.e. before the system reaches steady-state). From that point of view, the homogeneous part of the solution entails the system's ability to respond to the sudden changes in the inhomogeneous terms.
In the discussion so far, the system was described at short timescales, before it can reach its steady-state behavior. Let us now estimate this timescale. If not for the coupling of the two spin ensembles, there would be two distinct modes, one for the alkali and one for the noble gas. The rate with which the noble gas's mode decays in such a case is longer than that of the alkali by many orders of magnitude. The interaction between the atoms mix the two spin modes and the resulting system is described by two new eigen-modes with two new respective eigenvalues. Since we mostly care about how long it takes the system to reach equilibrium, it is sufficient to discuss the slower decay rate, Γ slow . Neglecting the rates, R K−He , R He−K ( which would have been the real components of the off-diagonal terms of Eq. (A3)), and R He (which are mostly irrelevant in the systems at hand), one finds for δω ≡ ω K − ω He The above is only an order of magnitude smaller than the (mostly) alkali mode's decay rate in typical systems for which Γ K √ ω He−K ω K−He . We point out that at the compensation point, the higher order corrections in 1/δω can become important. While highly dependent on the precise details, one often finds the two eigenvalues' real values to be of the same order of magnitude (see Table II for the values for Refs. [47][48][49]).
As an example to why the above discussion could be important, consider the case of Ref. [49]. In Ref. [49], the detector changes its direction every few seconds. Sudden changes in the magnetic fields are then expected due to possible field penetration as well as inner thermal noise of the magnetic shields. If the rate with which the system reaches equilibrium after each rotation is slower than the rate of rotations, the system never converges to its steady-state behavior. As a result, the homogeneous terms proportional to the magnetic fields can add significant contributions to the signal. Under realistic laboratory conditions, and even without such a clear intervention in the detector's environmental conditions, sudden changes in the magnetic fields occur, and unless the detector's response-time is fast enough, they can impair the measurements. Fortunately, since Γ slow 10 sec −1 , abrupt changes are treated rather efficiently in the systems studied here.

The Inhomogeneous Solutions: Response to Oscillatory Fields
Let us now discuss the steady-state response of the system to high-frequency magnetic fields. The inhomo-geneous solutions have three contributions. The first two are from the electron anomalous field and the nuclear anomalous field. These terms relate the alkali's spin measurement to those anomalous fields. The third contribution will come from oscillating magnetic fields. It is this that dictates the system's ability to cancel noise at a given frequency of oscillation.
Unlike the homogeneous solutions, whose frequency is determined by the linear system's parameters, an inhomogeneous term with a certain frequency will only induce an inhomogeneous solution of that frequency. As can be seen from Eq. (A3), an important change in the presence of an oscillating field with a frequency ω, is that the steady-state solution is no longer found by takinġ S C K =Ṡ C He = 0, but ratherṠ C K = iωS C K , andṠ C He = iωS C He . Much like in the case of the homogeneous solutions, the actual results are easy to calculate, but have cumbersome formulas. Nonetheless, close to the compensation point, and neglecting R He−K , R K−He , R He , one finds the approximate closed form solutions, Here S C K,He (ω) are the inhomogeneous contributions of the fields And P 1 (ω), P 2 (ω) are polynomials of degree one and two, and using the notations of Eq. (A3), and (B5) Note that due to the ALP field oscillating as cos(m a t+θ 0 ) with θ 0 an unknown phase, the negative and positive frequencies are mixed, and therefore the final dependence on the ALP field will be a symmetrized version of Eqs. (B2) and (B3).
While it is not yet entirely known what governs the noise spectrum of the comagnetometer at low frequencies (see e.g. Ref. [47][48][49]76] for calculations of the noise from theoretical arguments, and compare with results from Refs. [47][48][49]), at higher frequencies (usually ω ∼ > Γ K or |ω| ∼ > |ω c | |ω He |) there are reasons to believe that magnetic noise is the dominant factor. Eq. (B2) shows that such magnetic noises would enter the (measured) alkali's magnetization and thus one can approximate the ω-dependence of the signal-to-noise ratio due to suppressed response to ALP-neutron (ALP-electron) interaction, by dividing the coefficient of b n (b e ) with that of B in Eq. (B2). The conclusions is therefore that for ALP-neutron interactions, we expect an approximately linear decrease in the signal-to-noise sensitivity for high frequencies (the ratio is ∝ 1/ω), while we do not expect such a decrease for ALP-electron interactions [the ratio is ∝ P 1 (ω)/ω ∼ O(ω 0 )].

Appendix C: Effects of Signal Directionality
Here we discuss in detail the procedure for treating signal directionality. For the datasets used in this paper, Refs. [47][48][49], a simplified treatment sufficed (see Appendix F), however here we lay the groundwork for the formal treatment of velocity directionality, which will be relevant in the future with new independent highresolution data. Throughout this appendix we assume the measurement time t tot day, as the shortest datataking session used in our bounds was 4 days long.
The data in Refs. [47,49] is given in the form of Eq. (9) (the data of Ref. [48] is given in a similar, yet not identical form). The directional dependence on the relative ALP velocity is apparent, but the relative velocity of the ALPs with respect to earth is highly model dependent [62]. Different models can also change the local DM density significantly. Moreover, even for a given model, the local velocity and local density of the ALPs are statistical quantities, and thus need to be treated as such [61]. Our treatment of the effects of the non-deterministic properties of ALPs is presented in Appendix D. For the purposes of this appendix, we take the ALPs to have a constant relative velocity v, and a constant density ρ DM .
Under these assumptions, we look at the result of plugging the anomalous field implied by the Hamiltonian of Eq. (7) in the integrand of Eq. (9). We find that where we took the square of the absolute value of the amplitude as we are interested in root mean square over different parameters such as the relative velocity's direction and initial phases. We defined c ≡ 2ρ a |v| 2 g 2 aN N /(γ 2 n t tot ), in order to make the equations more tractable. t tot is the total measurement time. E a is the ALP energy, and because the ALP is non-relativistic, E a = m a + m a v 2 /2 m a (we will address the importance of deviation from that assumption in Appendix D).σ is the sensitive direction of the detector. We allowed an initial relative phase for the ALP field, θ 0 , which we will also discuss in further details in Appendix D.
The measurement itself is of the change of polarization in the probe beam behind the cell, rather than b n ·σ, and while the change of polarization is proportional to S x K (which is proportional to b n ·σ), there are calibration factors. These factors are measured individually by the different experiments, with the data given after calibration. The calibration is done by checking low frequency response, and therefore at higher frequencies, a correction is necessary, as was discussed in Appendix B. However, for this appendix, we shall assume that the data given is after the necessary additional corrections were made to correct for the higher frequencies -and thus Eq. (C1) can be assumed as the signal we are given.
To findσ let us use the coordinate system whereẑ is the direction of earth's rotation axis. We define the x − z plane so that at t = 0, the observer of an experiment on earth is described by (R ⊕ sin(θ), 0, R ⊕ cos(θ)), where R ⊕ is the earth's radius, and the observer's latitude coordinate is π/2 − θ. At time t, the observer's position is therefore, r(t) = R⊕(sin(θ) cos(ΩSDt), sin(θ) sin(ΩSDt), cos(θ)), (C2) with Ω SD 2π/day as the sidereal day frequency.
For the experiment of Ref. [49], the detector's sensitive direction alternated between the North-South (NS) directions to the East-West (EW) directions every few seconds. Thus, from the second experiment, we have a lowfrequency measurement of two differentσ, for the NS measurements, and for the EW measurements. For each of the three directions we plug the appropriate of the three Eqs. (C3),(C4),(C5), into Eq. (C1), to get the expected signal. The resulting signal is a complicated function of the many different parameters, and we thus do not show it here. However, as we do want to examine the expected form of the signal, it is useful to look at |A(ω,σ)| 2 , where we would average upon all possible directionsv, and upon the initial angle θ 0 . For any of the three directions, the resulting averaged signal squared has the form where the coefficients a i (σ) do not depend on the frequencies or the mass, and are in fact only dependent on σ(t = 0). The frequencies ω i have six possible values, the sum of one of the two {m a , −m a } with one of the three {Ω SD , −Ω SD , 0}. This form is reasonable, as when t tot → ∞ these terms become delta functions (up to normalization), and as we did not yet include the velocity smearing that shall be discussed in Appendix D, the ALPs are indeed infinitely sharp in the frequency range -albeit possibly shifted due to the earth's rotation. As long as m a is not within ∼ 2π/t tot of 0, Ω SD , or Ω SD /2, for all 1 ≤ i < j ≤ 6, we have When Eq. (C7) holds, A(ω) takes a similar form to Eq. (C6), albeit with a i (σ) → a i (σ,v). When this condition is not met, the signal might smear between different ω i 's, and take a more complicated form. When specifically m a ∼ < 2π/t tot , the effects of θ 0 are not negligible, and its stochastic nature must be accounted for (see Appendix D and Appendix F). We note that had we not assumed t tot day, the device's longitude coordinate, and hour at which measurements started would have played a role as well.

Appendix D: Effects of Non-Deterministic Signal
Eq. (8) presents us with the expected average field for the ALPs throughout the galaxy. However, due to the stochastic nature of the ALP field, E a , v, θ 0 and ρ a should not be treated as their average values throughout the galaxy, when only measured for a short time. Indeed, as we move in the galactic plane, we go through spatial gradients in the ALP field [77], and the local properties of the ALP field should be thought of as random variables sampled from a distribution centered around the average values. While there is debate in the astrophysics literature as to the size of these gradients, here we take the conservative approach of Ref. [78], by taking the typical scale of these gradients to be the De-Broglie wavelength of the ALPs, ∼ 2π/m a v virial .
Recently, Ref. [61] has shown how to treat the effect of the stochastic nature of the ALP field, and we base the methods presented in this appendix on theirs. For the case of the ALPs velocity distribution, we also use Ref. [79]. While Ref. [79] was discussing detection of DM scattering via direct detection, their general formulas for finding the relative velocities of virialized DM were useful for our discussion as well.
We will now discuss the different variables that were taken as non-deterministic, and the distributions of these variables. After that, we discuss the coherence time of the signal. The coherence time plays an important role in our treatment of the stochastic nature of the ALP field -we take an independent sampling of each of the nondeterministic variables every coherence time.
The stochastic nature of the ALP field Here we discuss one by one the non-deterministic variables of Eq. (8) (identical to the non-deterministic variables that affect b e ), and what distribution was chosen for them. The Initial Signal Phase. As we have already briefly mentioned in Appendix C, when E a t tot ∼ > 2π, the initial phase θ 0 becomes of little importance, as one goes over at least one oscillation of cos(E a t + θ 0 ). Conversely, when E a t tot 1, the signal can be highly dependent on that phase. As this phase is entirely random, we sample it from a uniform distribution between 0 and 2π.
The ALP Density. The anomalous field (Eq. (8)) depends on the square root of the DM energy density. The square root of the ALP energy density, √ ρ a , is Rayleigh distributed [61] around √ ρ DM = 0.4 GeV/cm 3 . We therefore sample √ ρ a from the following probability density function, The Velocity Distribution. Following Ref. [79], we take a Standard Halo Model (SHM), where our relative velocity compared to the DM is where v = (11, 232, 7) km/sec, is the velocity of the sun with respect to the galaxy, in galactic coordinates. v SHM is the randomly sampled DM velocity in the SHM with respect to the galactic rest-frame. We also use Ref. [79] for their formulas (which we do not reproduce here) to transition Eqs. (C3),(C5),(C4) to the galactic coordinates in which v is given. We have neglected the velocity of the earth with respect to the sun, which would introduce ∼ 10% annual modulation, and we have neglected the velocity from the rotation of the detector around the earth's axis which introduces sub-percent daily modulation. We emphasize that the daily modulations that are discussed in Appendix C are coming from the detector's sensitivedirection's rotation, and not from the small change in the detector's velocity due to earth's rotation. The SHM velocity's probability distribution function is, where Z is a normalization constant, Θ is the Heaviside function, v esc = 550 km/sec is the galactic escape velocity, v virial = 220 km/sec is the virial velocity. The Energy Distribution. The energy of the nonrelativistic ALPs, E a = m a (1 + v 2 a /2) should be entirely determined by the sampled velocity which was discussed in the previous paragraph. However, as the smearing of the searched frequency introduces a finite coherence time of ALP oscillations, it requires a more thorough discussion, which we perform below.

Effects of finite signal coherence time
Neglecting the small corrections due to finite galactic escape velocity in the SHM, the spread of velocities gives rise to a coherence time τ a = 2π/(m a v 2 ) 10 7 /m a [32] 10 . If a data-taking session is significantly shorter than the coherence time, we assume that the signal is entirely coherent throughout the measurement, i.e. only a single value should be sampled from the distributions discussed in this appendix.
A coherent signal should scale linearly with t tot , the measurement time, while the random noise will scale as √ t tot , giving rise to SNR [for S(ω)] that scales as √ t tot . This is why the data of Fig. 2 is given in the seemingly odd units of Gauss/ √ Hz. It is therefore expected that even if t tot is increased, the noise spectrum will look the same, while any contribution of the signal will peak over the noise as t tot increases.
Conversely, if t tot > τ a , for every τ a that passes since the beginning of the measurements, we sample the distributions one more time. Following Ref. [77], we can also sketch how to understand the dependence of the SNR on t tot after t tot > τ a . If we assume that n coherence times have passed, t tot = nτ a , this implies adding incoherently (with random relative phases), n coherent measurements of length τ a . When n 1, the expected measured amplitude of this incoherent summation is approximately the addition in quadratures (as it is a random walk) of the measurements. Therefore, the signal would scale as ∼ √ nτ a = √ t tot · τ a . This would imply that after the coherence time passes, there is no longer an advantage in taking longer measurements (as the signal to noise ratio no longer increases for t tot > τ a ). In a dedicated experiment, as explained in Ref. [77], and in analogy to the prescribed procedures of Ref. [81], it can be possible to increase sensitivity even after the coherence time passes using curve-fitting of the signal to smeared gaussians. However, since the data analyzed in this paper was not given with sufficient resolution and has gone through several processing procedures, we have not attempted such procedures.
Despite the above discussion, there's an O(3) improvement in the 95% C.L. bound in the transition from t tot ∼ τ a to t tot ∼ 5τ a . This improvement is because when we only have a single sampling of the stochastic distribution, the signal might be unusually small (due to a small ρ a , or a cancellation of v and v SHM ). Conversely, when t tot ∼ > 5τ a , that probability drops, as we sample 5 different values of our distribution, and while one of them might be small, on average, they would not be consistently small. However as was discussed in the above paragraph, after this improvement at t tot ∼ 5τ a , the SNR stops improving if one does not use curve fitting. 10 Since the ALP kinetic energy is 1 2 mav 2 virial , some authors use τa = 2π/(mav 2 virial /2), which is twice as long as the coherence time we use. Our shorter coherence time is conservative, and coincides with Ref. [80] which shows that τa = 2π/(mav 2 virial ) gives the correct frequency spread from Doppler broadening considerations. Regardless, since the bounds depend only weakly on the exact coherence time, this factor of 2 does not affect the results significantly.
In the case of Ref. [47], we are given data that were averaged from multiple measurements. The measurements were taken over a period of ∼ 100 days. We have assumed that when τ a (m a ) = 100 days, all non-deterministic variables were sampled from a single distribution. However, when τ a ∼ 8 days, most measurements are spaced enough for them to be considered independent samples of the stochastic distribution. We have used a simple interpolation to predict the suppression of the bound due to the stochastic nature of the ALP field, between τ a (m a ) = 100 days and τ a (m a ) = 8 days. Similar, more complicated methods have yielded similar results.
We have used MC simulations for our final bounds and projections presented in Figs. 3 and 4. The procedure of finding the bounds and projections after treating the effects in all other appendices, is described in further details in Appendix F. probes at different directions. At low frequencies, P 2 (ω) and P 1 (ω)/P 2 (ω) (corresponding to the alkali's magnetization due to an anomalous electron and neutron magnetic fields respectively; see Eq. (B2)), are almost purely imaginary, however at higher frequencies they have a real contribution as well, leading to sensitivity for ALPs in the direction parallel to the probe beam. This requires us to specify for the datasets which we analyze at high frequencies the direction of the probe beam -which for both Refs. [47,48] was 60 • with the NS directions (and no component in the direction of gravity).
By fitting to a function that is smeared, the effects of the ALP incoherence on the SNR which are described in Appendix D become apparent, as different sampled energies would widen the expected signal. However, as no curve fitting was attempted, the effect of finite coherence time has prevented any improvement in SNR after τ a passed, as described in the second part of Appendix D. For future analysis however, it was assumed that once the coherence time has been reached, the reach still improves as t −1/4 tot . At this point, we have a well defined procedure to extract the predicted signal from an ALP. What remains now is to account for the noise in order to find the 95% C.L. limit for a given m a , and either g aN N or g aee . The problem here is that unlike the more common cases of direct detection experiments, in our setup, noise may in fact theoretically cancel the signal. Therefore, without any model for the noise, a specific measurement could be the remains of a cancellation to an unknown degree between the signal and the noise.
To solve this problem, we need to understand how noise can affect the measurement. Assume that the noise at a given frequency ω, in a given experiment, is with some unknown amplitude A noise (ω), oscillating with an unknown initial phase. In this case, in order for complete cancellation of the ALP signal at the same frequency, A(ω), to occur, not only do we need A(ω) = A noise (ω), but also for the two phases to exactly match. As these two phases are entirely independent, we expect the relative phase to be a uniformly distributed random variable between 0 and 2π. Therefore, for a given A noise (ω), we can easily extract the 95% C.L. of A(ω) from the measured amplitude at ω. While we have no way of knowing A noise (ω), we simply take the conservative approach and assume the one that gives the weakest bound. Note that the result is always bounded from below and for any given A noise (ω), the probability to get complete cancellation of signal and noise is infinitesimal.
We find that nearly always, the strongest bounds are found when A noise (ω) → 0. The main reason for this is the stochastic nature of the signal. The possibility given in the previous paragraph of A noise (ω) = A signal (ω) is impossible when the signal is stochastic in nature. Even if the amplitudes are equal for a given v, ρ a , when we sample the non-deterministic variables, they would not cancel consistently.
The treatment for the dataset of Ref. [49] is a bit more complicated than the other two, and there is more freedom in the choice of statistical test. An ALP of mass m a would have a measurable amplitude at 5 different points of data, the |m a ± Ω SD | frequencies in both the EW and the NS searches, and the m a frequency in the NS search. We have taken the mean of these five measurements, though in the future, we expect smarter choices can be made, that could further teach us about the data (e.g. using the standard deviation of the five measurements to estimate the noise).
We note that when analyzing the data of Ref. [49], the three masses of m a = (0, 1 2 Ω SD , Ω SD ) require a special treatment, since for such cases some of the different frequencies at which we attempt to find a signal coincide (e.g. for m a = 1 2 Ω SD , m a = −m a + Ω SD ), or else we need the zero frequency data which we do not have. The analysis is a simple extension of the previously described procedures, so we do not reproduce it here.
Before moving to discuss the future projections, we finally note our efficiency estimates. For the data from Ref. [47], we are told that about 35% of the time, the detector was not actively measuring (e.g. due to the calibration routines), so the effective measurement time is only 65% of the reported t tot . As Ref. [49] uses the same procedures, we have taken its efficiency to be 0.65 × 0.5, as each of the two directions is actively measuring only half the time. For Ref. [48] it is written that the efficiency was between 0.2 and 0.6, so we have conservatively taken 0.2.
Future projections were calculated with a much simpler procedure compared to the bounds, since we cannot be sure of the precise experimental apparatus we will have. The reach is not to be thought of as expected 95% C.L. bound, but as the expected measurable signal. The reach is taken as the sensitivity described in the text, for a single month of exposure, and under the assumption that the bound when the measurement time is longer than the coherence time improves as t 1/4 tot (instead of √ t tot for t tot < τ a ). To account for the effects of the dynamical response, we assume the sensitivity is weakening linearly for the b n search at ω > ω c (with ω c given explicitly in Sec. VII).
The calculation of the long range spin-dependent interaction bound was written explicitly in Sec. VII, and it is effectively a rescaling of the bounds presented by Ref. [47] for their similar experiment.