Search for axion-like dark matter through nuclear spin precession in electric and magnetic fields

We report on a search for ultra-low-mass axion-like dark matter by analysing the ratio of the spin-precession frequencies of stored ultracold neutrons and $^{199}$Hg atoms for an axion-induced oscillating electric dipole moment of the neutron and an axion-wind spin-precession effect. No signal consistent with dark matter is observed for the axion mass range $10^{-24}~\textrm{eV} \le m_a \le 10^{-17}~\textrm{eV}$. Our null result sets the first laboratory constraints on the coupling of axion dark matter to gluons, which improve on astrophysical limits by up to 3 orders of magnitude, and also improves on previous laboratory constraints on the axion coupling to nucleons by up to a factor of 40.


I. INTRODUCTION
Astrophysical and cosmological observations indicate that 26% of the total energy density and 84% of the total matter content of the Universe is dark matter (DM) [1], the identity and properties of which still remain a mystery. One of the leading candidates for cold DM is the axion, a pseudoscalar particle which was originally hypothesised to resolve the strong CP problem of quantum chromodynamics (QCD) [2][3][4][5][6][7][8][9]. Apart from the canonical QCD axion, various axion-like particles have also been proposed, for example in string compactification models [10][11][12][13][14][15].
Low-mass (m a 0.1 eV/c 2 ) axions can be produced efficiently via non-thermal production mechanisms, such as vacuum misalignment [16][17][18] in the early Universe, and subsequently form a coherently oscillating classical field [19]: a = a 0 cos(ωt), with the angular frequency of oscillation given by ω ≈ m a c 2 / , where m a is the axion mass (henceforth, we shall adopt the units = c = 1). The oscillating axion field carries the energy density ρ a ≈ m 2 a a 2 0 /2. Due to its effects on structure formation [20], ultra-low-mass axion DM in the mass range 10 −24 eV m a 10 −20 eV has been proposed as a DM candidate that is observationally distinct from, and possibly favourable to, archetypal cold DM [15,[21][22][23][24]. The requirement that the axion de Broglie wavelength does not exceed the DM size of the smallest dwarf galaxies and consistency with observed structure formation [25][26][27] give the lower axion mass bound m a 10 −22 eV, if axions comprise all of the DM. However, axions with smaller masses can constitute a sub-dominant fraction of DM [28].
It is reasonable to expect that axions interact nongravitationally with standard-model particles. Direct searches for axions have thus far focused mainly on their arXiv:1708.06367v1 [hep-ph] 21 Aug 2017 coupling to the photon (see the review [29] and references therein). Recently, however, it has been proposed to search for the interactions of the coherently oscillating axion DM field with gluons and fermions, which can induce oscillating electric dipole moments (EDMs) of nucleons [30] and atoms [31][32][33], and anomalous spinprecession effects [31,34,35]. The frequency of these oscillating effects is dictated by the axion mass, and more importantly, these effects scale linearly in a small interaction constant [30][31][32][33][34][35], whereas in previous axion searches, the sought effects scaled quadratically or quartically in the interaction constant [29].
In the present work, we focus on the axion-gluon and axion-nucleon couplings: where G andG are the gluonic field tensor and its dual, b = 1, 2, ..., 8 is the color index, g 2 /4π is the color coupling constant, N andN = N † γ 0 are the nucleon field and its Dirac adjoint, f a is the axion decay constant, and C G and C N are model-dependent dimensionless parameters. Astrophysical constraints on the axion-gluon coupling in (1) come from Big Bang nucleosynthesis [36][37][38]: m 1/4 a f a /C G 10 10 GeV 5/4 for m a 10 −16 eV and m a f a /C G 10 −9 GeV 2 for m a 10 −16 eV, assuming that axions saturate the presentday DM energy density, and from supernova energyloss bounds [35,39]: f a /C G 10 6 GeV for m a 3 × 10 7 eV. Astrophysical constraints on the axion-nucleon coupling in (1) come from supernova energy-loss bounds [39,40]: f a /C N 10 9 GeV for m a 3 × 10 7 eV, while existing laboratory constraints come from magnetometry searches for new spin-dependent forces mediated by axion exchange [41]: f a /C N 1 × 10 4 GeV for m a 10 −7 eV.
The derivative coupling of an oscillating galactic axion DM field, a = a 0 cos(m a t − p a · r), with spin-polarized nucleons in (1) induces time-dependent energy shifts according to: The term σ N ·p a is conveniently expressed by transforming to a non-rotating celestial coordinate system (see, e.g., [56]): where χ is the angle between Earth's axis of rotation and the spin quantization axis (χ = 42.5 • at the location of the PSI), δ ≈ −48 • and η ≈ 138 • are the declination and right ascension of the galactic axion DM flux relative to the Solar System [57], Ω sid ≈ 7.29 × 10 −5 s −1 is the daily sidereal angular frequency,m F = m F /F is the normalized projection of the total angular momentum onto the quantization axis, and f (σ N ) = +1 for the free neutron, while f (σ N ) = −1/3 for the 199 Hg atom in the Schmidt (single-particle) model. Here, we report on a search for an axion-induced oscillating EDM of the neutron (nEDM) based on an analysis of the ratio of the spin-precession frequencies of stored ultracold neutrons and 199 Hg atoms, which is a system that had previously also been used as a sensitive probe of new non-EDM physics [58][59][60]. We divided our analysis into two parts. We first analyzed the Sussex-RAL-ILL nEDM experiment data [61], covering oscillation periods longer than days (long time-base). Then we extended the analysis to the data of the PSI nEDM experiment [62], which allowed us to probe oscillation periods down to minutes (short time-base). Our analysis places the first laboratory constraints on the axion-gluon coupling. We also report on a search for an axion-wind spin-precession effect, using the data of the PSI nEDM experiment. Our analysis places the first laboratory constraints on the axion-nucleon coupling from the consideration of an effect that is linear in the interaction constant.

II. LONG TIME-BASE ANALYSIS
The Sussex-RAL-ILL room temperature nEDM experiment ran from 1998 to 2002 at the PF2 beamline at the Institut Laue-Langevin (ILL) in Grenoble, France. This experiment set the current world-best limit on the permanent time-independent neutron EDM, published in 2006 [63]. The data were subsequently reanalyzed to give a revised limit in 2015 [64]. The technical details of the apparatus are described in full in [61], but we summarize the main experimental details here for the reader.
The experiment was based on Ramsey interferometry [65] of ultracold neutrons [66,67]. The neutrons were stored in parallel or antiparallel electric and magnetic fields, where their Larmor precession frequency is given by with the sign depending on the field configuration. E and B are the magnitudes of the electric and magnetic fields, respectively. By measuring the frequency difference between the two field configurations, a value for the neutron EDM, d n , was inferred. The measurement was conducted in a series of cycles, each approximately 5 minutes long. A cycle began with a filling of neutrons polarized along the fields into the precession chamber from the ultracold neutron source [68]. Once they are in the chamber, enclosed from top and bottom with electrodes, a 29 Hz NMR pulse lasting 2 seconds was applied to rotate the neutron spins into the transverse plane of the electromagnetic fields where they began to precess. Prior to the pulse, a population of polarized 199 Hg atoms was released into the chamber, another 2-second 29 Hz pulse, in-phase with the first one, was applied. The neutrons were then emptied into a detector through a spin-analyzing foil. Over 1-2 days, many of these cycles were performed. The electric field's polarization was reversed every hour. We term one continuous block of data taking in the same magnetic-field configuration, but including both directions of electric field, a run. One run gives a d n estimate.
In order to suppress cycle-to-cycle changes in the magnetic field, the analysis was performed on the ratio of the neutron and mercury precession frequencies R, which, using (6), is [61]: where the signs correspond to parallel and antiparallel field configurations. ∆ encapsulates all higher-order terms and systematic effects, which are corrected for when a run is analyzed [64]. This analysis is sensitive to oscillations in the quantity d n − (µ n /µ Hg ) d Hg , with µ n /µ Hg = −3.8424574(30) [69]. In our analysis, we were looking for an oscillating EDM. We performed this search in frequency space by evaluating periodograms -estimators of the power spectrum. An oscillation in the time domain would show up as an excess in the power (or, equivalently, amplitude) relative to the expected distribution due to experimental noise.
In the case of the long time-base analysis, we considered the time series of d n measurements from individual runs (after having corrected for the "False EDM" effect [70] using the crossing lines procedure [64]). The measurements are neither evenly spaced nor have equal uncertainties. To calculate the periodogram of the data series, we used the Least Squares Spectral Analysis (LSSA) method [71,72], where the amplitude at frequency f was estimated by the amplitude of the best fit oscillation of that frequency. We evaluated the periodogram at a set of 1334 trial frequencies, evenly spaced between 100 pHz (arbitrarily chosen, a period of about 300 years, much longer than the four-year span of the data set) and 10 µHz (a period of about a day, the time it typically took to get one d n estimate). An axion DM signal, with expected coherence set by ∆f ∼ 10 −6 f [19], is narrower than the spectral resolution (7.49 nHz, the inverse span of the data set) in the whole range of frequencies we were sensitive to. In the LSSA fit, we assumed the free offset to be zero on the grounds that the experiment had already delivered a zero-compatible result for the permanent time-independent neutron EDM [63,64]. The periodogram of the long time-base dataset is shown as a black line in Fig. 1. To obtain the expected distribution of the periodogram, we performed Monte Carlo (MC) simulations. At each frequency, we estimated the cumulative distribution function (CDF) of the LSSA power. Extreme events in the tails of the distribution are expensive to access directly with MC. For this reason, to the discrete CDF estimates we fitted, at each i th frequency, the functional form of the LSSA-power CDF [71]: where P is the power, while A i and B i are fit parameters. The local p-values are given by where P i is the LSSA power of the measured d n time series at the i th frequency. If the local p-values at different trial frequencies were uncorrelated, the global p-value would be given by [73]: where N is the number of trial frequencies. However, we did not need to make this assumption. Instead, we made use of the set of MC datasets. In each, we found the minimal local p-value and estimated its CDF, assuming it has the form (10), but left N as a free parameter. We found the best fit value N effective = 1026. For each frequency, we marked the power necessary to reach the global p-values corresponding to 1, 2, . . . , 5 σ levels as orange lines in Fig. 1. The minimal local p-value of the dataset translates to the global p-value of 0.53, consistent with a non-detection.
In order to obtain limits on the oscillation amplitude parameter, we again used MC simulations. We discretized the space of possible signals, spanned by their frequency and amplitude. We chose a sparser set of 200 frequencies, as we did not expect highly coherent effects in the sensitivity of detection. For each discrete point, we generated a set of 200 MC datasets containing the respective, perfectly coherent signal and assumed that the oscillation is averaged over the duration of the run. In general, the sensitivity is phase-dependent, especially for periods comparable with the length of the dataset. For simplicity, we did not investigate the phase-dependence and in the simulation took it to be random and uniformly distributed. For each fake dataset, we evaluated the LSSA amplitude only at the frequency of the signal and compared its distribution (extrapolating with the functional form of Eq. (8)) with the best-fit amplitude in the data and defined the p-value to be left-sided. We found the 95% confidence-level exclusion limit as the 0.05 isocontour of the CL s statistic [74]. The limit is shown as the red curve in Fig. 2. We are most sensitive to periods shorter than the timespan of the dataset (∼ 4 years), but rapidly lose sensitivity for periods shorter than the temporal spacing between data points (∼ 2 days), since the expected signal would essentially average to zero over these short time scales.

III. SHORT TIME-BASE ANALYSIS
In 2009, the Sussex-RAL-ILL apparatus was moved to the new ultracold neutron source at the Paul Scherrer Institute (PSI), Villigen, Switzerland [75][76][77][78], where a number of improvements were made [62,79,80]. In 2015, the apparatus was fully commissioned and began to take high-sensitivity EDM data. The whole data set, taken from August 2015 until the end of 2016, with a higher accumulated sensitivity than the ILL one, was considered in this analysis. For the PSI experiment's data, we performed a lower-level oscillation search on the array of R measurements. Since an R estimate was obtained every cycle (≈ 300 s), rather than every 1-2 days as for a d n estimate, it has an increased sensitivity to higher frequencies. Additionally, the analysis could benefit from the addition of 16 atomic cesium vapor magnetometers [81, 82], located directly above and below the precession chamber (inside the electrodes). This made it possible to account for the dominant time-dependent systematic effect on a cycle, rather than run, basis.
The dominant time-dependent systematic effect, encapsulated in ∆ of Eq. (7), would have given rise to non-statistical temporal fluctuations if not accounted for. Namely, R is sensitive to drifts in the vertical gradients of the magnetic field. While the thermal mercury atoms filled the chamber homogeneously, the center of mass of the ultracold neutron population was lower by several millimeters [64,69,83]. To evaluate the correction, the drifts of the gradients were estimated on a cycle basis by fitting a second-order parametrization of the magnetic field to the measurements of the cesium magnetometers [84]. The center-of-mass shift was determined to be 4 mm using the method described in [69].
The measurement procedure involved working deliberately with gradients affecting R (see the crossingpoint method in [64]). Those intended gradients (up to 60 pT/cm in 10 pT/cm steps) were much larger than cycle-to-cycle fluctuations (< 2 pT/cm per day). With the high-order shifts in R having been significant, these large shifts could not be corrected using the cesium magnetometers. Additionally, while the cesium magnetometers were precise, their accuracy is limited by the calibration procedure. We defined as a sequence a set of data, typically 2-3 days in duration, without a deliberate change in the magnetic-field gradient or a recalibration of the cesium magnetometers. When performing the LSSA fit, we allowed the free offset to be different in each sequence: where C i is the free offset in the i th sequence and Π i (t) is a gate function equal to one in the i th sequence and zero elsewhere. This caused the short time-base analysis to lose sensitivity for periods longer than one sequence. It should also be mentioned that, at the time of this analysis, the PSI data were still blinded, whereby an unknown, but constant, d n was injected into them. It does not influence this analysis, as the free offsets are not considered further.
We split the R time array into three sets: a control set of data without an applied electric field, and two sets sensitive to an oscillating EDM, namely with parallel and antiparallel applied electric and magnetic fields. A coherent oscillating EDM signal would have an opposite phase in the latter two sets, and be absent in the control set. We did not perform a common fit. Instead, the two sensitive data sets were treated separately in the LSSA fits, and later combined to a limit. Otherwise, the LSSA treatment was the same as in the long time-base analysis. We picked a set of 156 198 trial frequencies, spaced apart at intervals determined by the spectral resolution (the inverse of 506 days = 23 nHz), which here also defines the signal width.
The periodogram of the R time array taken with the parallel-field configuration is shown in black in Fig. 3. There are two regions of expected rise in the oscillation amplitude due to the time structure of the data collection. The one around 28 µHz (the inverse of 10 hours) corresponds to the period of the electric-field reversal. The very narrow one around 3.3 mHz (the inverse of 300 s) corresponds to the cycle repetition rate. There are five trial frequencies for which the 3σ false-alarm threshold is exceeded, two of which, including the largest excess with a 6σ significance, occur in a 100 µHz region around the inverse of 300 s, while the other three are in the lowfrequency region (inverse days) already excluded by the long time-base analysis. The periodograms for the other two datasets (not shown) are very similar. In the other sensitive set, there are three excesses of the 3σ threshold (the highest is 5σ), all constrained to the same two regions. In the control dataset, only the 1σ threshold is exceeded. The periodogram of the R time array without the gradient-drift correction is shown in pink in Fig. 3 to visualize the frequencies where the correction has an effect.
A non-statistical excess in a periodogram of R may be caused not only by a coherent oscillating signal; for example, fluctuations of a higher-order term in the magnetic field, not compensated by either the mercury or cesium magnetometers, may cause broad-band elevations in LSSA power. We defined strict requirements for an excess to be considered as one induced by axion DM as follows. Firstly, a significant (> 3σ) excess in amplitude had to be observed in both sensitive datasets at Periodogram of the R time array of the PSI experiment data, sensitive to oscillations in the quantity dn − (µn/µHg) dHg, taken with the E and B fields parallel (black line). The mean of MC-generated periodograms, assuming no signal, is depicted in green. MC is used to calculate 1, 2, . . . , 5 σ false-alarm thresholds, depicted in light orange. For clarity, we also plot the smoothed version in orange. There are two regions where a rise in the amplitude is expected, namely around 28 µHz (inverse of 10 hours) and 3.3 mHz (inverse of 300 seconds), due to the time structure of the data taking (see the main text for more details). The periodogram of non-gradient-drift-corrected data is shown in pink.
the same frequency, but not in the control set. Secondly, the signals must be in antiphase in the parallel and anti-parallel datasets. Lastly, we require high coherence (a narrow peak) equal to the spectral resolution of the dataset. None of the significant excesses passed our discovery criteria.
We delivered a limit on the oscillation amplitude similarly to the long time-base analysis, with the exception that we required the product of the two sensitive sets' CL s statistics to be 0.05. The limit is shown as the blue curve in Fig. 2. With the short time-base analysis, we were most sensitive to periods shorter than the timespan of a sequence (2 -3 days), and lost sensitivity to periods shorter than the cycle repetition rate (≈ 5 minutes). The PSI dataset has a higher accumulated sensitivity than the ILL dataset, so the limit baseline in the sensitive region is slightly better in the case of the PSI dataset.
Following Eq. (2), we can interpret the limit on the oscillating neutron EDM as limits on the axion-gluon coupling in Eq. (1). We present these limits in Fig. 4, assuming that axions saturate the local cold DM energy density ρ local DM ≈ 0.4 GeV/cm 3 [55]. Our peak sensitivity is f a /C G ≈ 1 × 10 21 GeV for m a 10 −23 eV, which probes super-Planckian axion decay constants (f a > M Planck ≈ 10 19 GeV), that is, interactions that are intrinsically feebler than gravity.

IV. AXION-WIND EFFECT
We also perform a search for the axion-wind effect, Eq. (4), by partitioning the entire PSI dataset into two sets with opposite magnetic-field orientations (irrespective of the electric field) and then analyzing the ratio R = ν n /ν Hg similarly to our oscillating EDM analysis above. The axion-wind effect would manifest itself through time-dependent shifts in ν n and ν Hg (and hence R) at three angular frequencies: ω 1 = m a , ω 2 = m a +Ω sid and ω 3 = |m a − Ω sid |, with the majority of power concentrated in the ω 1 mode. Also, the axion-wind signal would have an opposite phase in the two subsets. We find two overlapping 3σ excesses in the two subsets (at 3.42969 µHz and 3.32568 mHz), neither of which have a phase relation consistent with an axion-wind signal. Following Eq. (4), we derive limits on the axion-nucleon coupling in Eq. (1). We present these limits in Fig. 4, assuming that axions saturate the local cold DM energy density. Our peak sensitivity is f a /C N ≈ 4 × 10 5 GeV for 10 −19 eV m a 10 −17 eV.

V. CONCLUSIONS
In summary, we have performed a search for a timeoscillating neutron EDM in order to probe the interaction of axion-like dark matter with gluons. We have also performed a search for an axion-wind spin-precession effect in order to probe the interaction of axion-like dark matter with nucleons. So far, no significant oscillations have been detected, allowing us to place limits on the strengths of such interactions. Our limits improve upon existing astrophysical limits on the axion-gluon coupling by up to 3 orders of magnitude and also improve upon existing laboratory limits on the axion-nucleon coupling by up to a factor of 40. Furthermore, we constrain a region of axion masses that is complementary to proposed "on-resonance" experiments in ferroelectrics [85]. Future EDM measurements will allow us to probe even feebler oscillations and for longer periods of oscillation that correspond to smaller axion masses.

ACKNOWLEDGMENTS
We are grateful to Maxim Pospelov for helpful discussions. The experimental data has been taken in part at the ILL Grenoble and at PSI Villigen. We acknowledge the excellent support by the technical groups of both institutions and by various services of the collaborating universities and research laboratories. Dedicated technical support by M. Meier and F. Burri is gratefully acknowledged. We remember with gratitude the pioneering contributions of Professors K. Smith and J. M. Pendlebury, without whom these experiments could never have taken place. This work was funded in part by the U. K. Science and Technology Facilities Council (STFC) through  1), assuming that axions saturate the local cold DM content. The regions above the thick blue and red lines correspond to the regions of parameters excluded by the present work at the 95% confidence level (C.L.). The colored regions represent constraints from Big Bang nucleosynthesis (red, 95% C.L.) [36][37][38], supernova energy-loss bounds (green, order of magnitude) [35,39,40], consistency with observations of galaxies (orange) [15,[25][26][27], and laboratory searches for new spin-dependent forces (yellow, 95% C.L.) [41]. The nEDM, νn/νHg and Big Bang nucleosynthesis constraints scale as ∝ √ ρa, while the constraints from supernovae and laboratory searches for new spin-dependent forces are independent of ρa. The constraints from galaxies are relaxed if axions constitute a sub-dominant fraction of DM. We also show the projected reach of the proposed CASPEr experiment (dotted black line) [85], and the parameter space for the canonical QCD axion (purple band).