Constraints on bosonic dark matter from ultralow-field nuclear magnetic resonance

TL;DR

This paper uses ultralow-field nuclear magnetic resonance to search for bosonic dark matter, setting new bounds on axion-like particles and dark photons within a specific mass range, without detecting a signal.

Contribution

It introduces a novel application of ultralow-field NMR in constraining bosonic dark matter couplings, providing the first experimental bounds in this mass range.

Findings

No dark matter signal detected above background.

Established new bounds on axion-fermion and dark-photon couplings.

Probed dark-matter masses from 1.8×10⁻¹⁶ to 7.8×10⁻¹⁴ eV.

Abstract

The nature of dark matter, the invisible substance making up over $80\%$ of the matter in the Universe, is one of the most fundamental mysteries of modern physics. Ultralight bosons such as axions, axion-like particles or dark photons could make up most of the dark matter. Couplings between such bosons and nuclear spins may enable their direct detection via nuclear magnetic resonance (NMR) spectroscopy: as nuclear spins move through the galactic dark-matter halo, they couple to dark-matter and behave as if they were in an oscillating magnetic field, generating a dark-matter-driven NMR signal. As part of the Cosmic Axion Spin Precession Experiment (CASPEr), an NMR-based dark-matter search, we use ultralow-field NMR to probe the axion-fermion "wind" coupling and dark-photon couplings to nuclear spins. No dark matter signal was detected above background, establishing new experimental bounds for dark-matter bosons with masses ranging from $1.8\times 10^{-16}$ to $7.8\times 10^{-14}$ eV.

Figures

Equations

$$a(t)\approx a_0\cos (\omega {}_{\text{mbox}DM}t),$$

$$\rho {}_{\text{mbox}DM}=\frac{1}{2}\frac{c^2}{{\hslash }^2}m{{}_{\text{mbox}DM}}^2{a}_0^2.$$

$$\mathcal{H}{}_{\text{mbox}DM}=-\hslash g{}_{\text{mbox}N}\mathbf{I}{}_{\text{mbox}N}\cdot \mathbf{D},$$

$$\mathbf{D}{}_{\text{mbox}wind}(t)=-\sqrt{2\hslash c{\rho }_{\mathrm{DM}}}\sin \left({\omega }_{\mathrm{DM}}t\right)\mathbf{v}.$$

$$\begin{gathered} \mathcal{H}=\hslash J{}_{\text{mbox}CH}\mathbf{I}\cdot \mathbf{S}-\hslash \left(\gamma {}_{\text{mbox}H}\mathbf{I}+\gamma {}_{\text{mbox}C}\mathbf{S}\right)\cdot \mathbf{B} \\ +\hslash \left(g{}_{\text{mbox}app}\mathbf{I}-\frac{1}{3}g{}_{\text{mbox}ann}\mathbf{S}\right)\cdot \mathbf{D}{}_{\text{mbox}wind}(t), \end{gathered}$$

$$\begin{gathered} \Delta E(m_F=\pm 1)(t)=\mp \frac{\hslash }{2}{B}_z(\gamma {}_{\text{mbox}H}+\gamma {}_{\text{mbox}C}) \\ \pm \frac{\hslash }{2}\frac{2}{3}g{}_{\text{mbox}aNN}{D}_z(t), \end{gathered}$$

$$\begin{gathered} B{}_{\text{mbox}wind}(t)=\frac{2}{3}\frac{g{}_{\text{mbox}aNN}}{\gamma {}_{\text{mbox}H}+\gamma {}_{\text{mbox}C}}\cdot D{}_{\text{mbox}z}(t) \\ \propto \frac{g{}_{\text{mbox}aNN}}{\gamma {}_{\text{mbox}H}+\gamma {}_{\text{mbox}C}}\sqrt{{\rho }_{\mathrm{DM}}}\sin ({\omega }_{\mathrm{DM}}t)\mathbf{v}\cdot \hat{\mathbf{e}}{}_{\text{mbox}z}. \end{gathered}$$

$$\begin{gathered} B_{\mathrm{ALP},\mathrm{z}}^{\mathrm{eff}}(t)=-\frac{2}{3}\frac{g{}_{\text{mbox}aNN}}{\gamma {}_{\text{mbox}H}+\gamma {}_{\text{mbox}C}}\sqrt{2\hslash c{\rho }_{\mathrm{DM}}} \\ \times \sin (\omega {}_{\text{mbox}DM}t-\mathbf{k}\cdot \mathbf{r}+\varphi )\mathbf{v}\cdot \hat{\mathbf{e}}{}_{\text{mbox}z}, \end{gathered}$$

$$\begin{gathered} B_{\mathrm{quad},\mathrm{z}}^{\mathrm{eff}}(t)=-\frac{4}{3}\hslash c^2\frac{g{{}_{\text{mbox}quad}}^2}{\gamma {}_{\text{mbox}H}+\gamma {}_{\text{mbox}C}}\frac{{\rho }_{\mathrm{DM}}}{{\omega }_{\mathrm{DM}}} \\ \times \sin (2\omega {}_{\text{mbox}DM}t-2\mathbf{k}\cdot \mathbf{r}+\varphi )\mathbf{v}\cdot \hat{\mathbf{e}}{}_{\text{mbox}z}. \end{gathered}$$

$$\begin{gathered} B_{\mathrm{dEDM}}^{\mathrm{eff}}(t)=\frac{2}{3}\frac{g{}_{\text{mbox}dEDM}}{\gamma {}_{\text{mbox}C}+\gamma {}_{\text{mbox}H}}\sqrt{2\hslash c^3\rho {}_{\text{mbox}DM}} \\ \times \cos (\omega {}_{\text{mbox}DM}t+\varphi )\mathbf{ϵ}\cdot \hat{\mathbf{e}}{}_{\text{mbox}z}, \end{gathered}$$

$$\begin{gathered} B_{\mathrm{dMDM}}^{\mathrm{eff}}(t)=\frac{2}{3}\frac{g{}_{\text{mbox}dMDM}}{\gamma {}_{\text{mbox}C}+\gamma {}_{\text{mbox}H}}\frac{v}{c}\sqrt{2\hslash c^3\rho {}_{\text{mbox}DM}} \\ \times \cos (\omega {}_{\text{mbox}DM}t+\varphi )\mathbf{ϵ}\cdot \hat{\mathbf{e}}{}_{\text{mbox}z}, \end{gathered}$$

$$\delta g{}_{\text{mbox}aNN}(\omega )$$

$$\delta g{}_{\text{mbox}quad}(\omega )$$

$$\delta g{}_{\text{mbox}dEDM}(\omega )$$

$$\delta g{}_{\text{mbox}dMDM}(\omega )$$

$${\text{A}}_s(B,\omega )=5.53\times 10^{-6}{\text{rad.s}}^{-1}\times \frac{B[\text{T}]}{\omega [{\text{rad.s}}^{-1}]}.$$

$$\mathcal{H}{}_{\text{mbox}J}=\hslash J{}_{\text{mbox}CH}\mathbf{I}\cdot \mathbf{S},$$

$$E(F)/\hslash =\frac{1}{2}J{}_{\text{mbox}CH}[F(F+1)-I(I+1)-S(S+1)],$$

$$E(F=0,m_F=0)/\hslash =-\frac{3}{4}J{}_{\text{mbox}CH},$$

$$E(F=1,m_F=0,\pm 1)/\hslash =\frac{1}{4}J{}_{\text{mbox}CH}.$$

$$\mathcal{H}$$

$$=\mathcal{H}{}_{\text{mbox}J}-\hslash (\gamma {}_{\text{mbox}C}S{}_{\text{mbox}z}+\gamma {}_{\text{mbox}H}I{}_{\text{mbox}z})B{}_{\text{mbox}z},$$

$$\Delta E(F,m_F)/\hslash =-⟨F,m_F|(\gamma {}_{\text{mbox}C}S{}_{\text{mbox}z}+\gamma {}_{\text{mbox}H}I{}_{\text{mbox}z})B{}_{\text{mbox}z}|F,m_F⟩,$$

$$E(F=0,m_F=0)/\hslash =-\frac{3}{4}J{}_{\text{mbox}CH},$$

$$E(F=1,m_F=0,\pm 1)/\hslash =\frac{J{}_{\text{mbox}CH}}{4}-m_F\frac{\gamma {}_{\text{mbox}C}+\gamma {}_{\text{mbox}H}}{2}B{}_{\text{mbox}z}.$$

$$\mathcal{H}$$

$$=\mathcal{H}{}_{\text{mbox}J}+\mathcal{H}{}_{\text{mbox}Z}-\hslash (\gamma {}_{\text{mbox}C}S{}_{\text{mbox}z}+\gamma {}_{\text{mbox}H}I{}_{\text{mbox}z})B{}_{\text{mbox}AC}\cos (\omega t).$$

$$\mathcal{H}{}_{\text{mbox}wind}$$

$$=-\hslash g{}_{\text{mbox}aNN}\mathbf{D}{}_{\text{mbox}wind}\cdot \mathbf{I}{}_{\text{mbox}N},$$

$$\mathbf{D}{}_{\text{mbox}wind}(t)=-\sqrt{\hslash c^3}\nabla a(\mathbf{r},t),$$

$$a(\mathbf{r},t)=a_0\cos (\omega {}_{\text{mbox}DM}t-\mathbf{k}\cdot \mathbf{r}+\varphi ),$$

$$\nabla a(\mathbf{r},t)$$

$$\rho {}_{\text{mbox}DM}=\frac{1}{2}\frac{c^2}{{\hslash }^2}m{{}_{\text{mbox}DM}}^2{a}_0^2,$$

$$\mathbf{D}{}_{\text{mbox}wind}(t)=-\sqrt{2\hslash c\rho {}_{\text{mbox}DM}}\times \sin (\omega {}_{\text{mbox}DM}t-\mathbf{k}\cdot \mathbf{r})\mathbf{v}.$$

Full text

††journal: Journal Name

1 INTRODUCTION

1.1 Ultralight bosonic dark matter

The nature of dark matter, the invisible substance that makes up over $80\%$ of the matter in the universe [1], is one of the most intriguing mysteries of modern physics [2]. Elucidating the nature of dark matter will profoundly impact our understanding of cosmology, astrophysics, and particle physics, providing insights into the evolution of the Universe and potentially uncovering new physical laws and fundamental forces beyond the Standard Model. While the observational evidence for dark matter is derived from its gravitational effects at the galactic scale and larger, the key to solving the mystery of its nature lies in directly measuring non-gravitational interactions of dark matter with Standard Model particles and fields.

To date, experimental efforts to directly detect dark matter have largely focused on Weakly Interacting Massive Particles (WIMPs), with masses between 10 and 1000 GeV [3, 4]. Despite considerable efforts, there have been no conclusive signs of WIMP interactions with ordinary matter. The absence of evidence for WIMPs has reinvigorated efforts to search for ultralight bosonic fields, another class of theoretically well-motivated dark matter candidates [5], composed of bosons with masses smaller than a few eV. A wide variety of theories predict new spin-0 bosons such as axions and axion-like particles (ALPs) as well as spin-1 bosons such as dark photons [1, 6]. Their existence may help to answer other open questions in particle physics such as why the strong force respects the combined charge-conjugation and parity-inversion (CP) symmetry to such a high degree [7], the relative weakness of the gravitational interaction [8], and how to unify the theories of quantum mechanics and general relativity [9].

Bosonic fields can be detected through their interactions with Standard Model particles. Most experiments searching for bosonic fields seek to detect photons created by the conversion of these bosons in strong electromagnetic fields via the Primakoff effect [10, 11, 12, 13, 14, 15, 16]. Another method to search for dark-matter bosonic fields was recently proposed: dark-matter-driven spin-precession [17, 18, 19], detected via nuclear magnetic resonance (NMR) techniques [18, 19, 20, 21, 22]. These concepts were recently applied to data measuring the permanent electric dipole moment (EDM) of the neutron and successfully constrained ALP dark matter with masses $\lesssim 10^{-17}$ eV [22].

The Cosmic Axion Spin Precession Experiment (CASPEr) is a multi-faceted research program using NMR techniques to search for dark-matter-driven spin precession [18]. Efforts within this program utilizing zero- to ultralow-field (ZULF) NMR spectroscopy [23] are collectively referred to as CASPEr-ZULF. We report here experimental results from a CASPER-ZULF search for ultralight bosonic dark-matter fields, probing bosons with masses ranging from $1.8\times 10^{-16}$ to $7.8\times 10^{-14}$ eV (corresponding to Compton frequencies ranging from $\sim 45$ mHz to $19$ Hz).

In the following section, we explain how dark-matter bosonic fields, in particular axion, ALP and dark-photon fields, can be detected by examining ZULF NMR spectra. The subsequent sections describe the measurement scheme and the data processing techniques employed during this search. Finally we present new laboratory results for bosonic dark matter, complementing astrophysical constraints obtained from supernovae $1987A$111We note that the constraints based on SN$1987$A data continue to be reexamined; see e.g. Refs. [24, 25]. [26].

2 RESULTS

2.1 Dark-matter field properties

If dark matter predominantly consists of particles with masses $m{{}_{\mbox{\scriptsize DM}}}\lesssim 10~{}{\rm eV}$, making up the totality of the average local dark-matter density, then they must be bosons with a large mode occupation number. It would be impossible for fermions with such low masses to account for the observed galactic dark matter density, since the Pauli exclusion principle prevents them from having the required mode occupation.

In this scenario, axion and ALP bosonic dark matter is well described by a classical field $a(t)$, oscillating at the Compton frequency ($\omega{{}_{\mbox{\scriptsize DM}}}\approx m{{}_{\mbox{\scriptsize DM}}}c^{2}/\hbar$) [27, 28, 29]:

[TABLE]

where $c$ is the velocity of light, $\hbar$ is the reduced Planck constant and $a_{0}$ is the amplitude of the bosonic field.

The temporal coherence of the bosonic field is limited by the relative motion of the Earth through random spatial fluctuations of the field. The characteristic coherence time $\tau{{}_{\mbox{\scriptsize DM}}}$, during which the bosonic dark-matter fields remains phase coherent, corresponds to $\sim 10^{6}$ periods of oscillation of the fields [18].

The amplitude $a_{0}$ can be estimated by assuming that the field energy density constitutes the totality of the average local dark-matter energy density ($\rho{{}_{\mbox{\scriptsize DM}}}\approx 0.4~{}{\rm GeV/cm^{3}}$ [30]). Then $a_{0}$ is related to the dark-matter density through:

[TABLE]

We note that $a_{0}$ is expected to fluctuate with relative amplitude of order one due to self-interference of the field. For simplicity we assume $a_{0}$ to be constant.

2.2 Dark-matter couplings to nuclear spins

CASPEr-Wind is sensitive to any field such that its interaction with nuclear spins can be written in the form:

[TABLE]

where $g{{}_{\mbox{\scriptsize N}}}$ is a coupling constant that parametrizes the coupling of the effective field, $\boldsymbol{D}$, to nuclear spins represented by the operator $\boldsymbol{I}{{}_{\mbox{\scriptsize N}}}$. In analogy with the Zeeman interaction, $\mathcal{H}{{}_{\mbox{\scriptsize Zeeman}}}=-\hbar\gamma{{}_{\mbox{\scriptsize N}}}\boldsymbol{I}{{}_{\mbox{\scriptsize N}}}\cdot\boldsymbol{B}$, such an effective field $\boldsymbol{D}$ may be thought of as a pseudo-magnetic field interacting with nuclear spins, where the nuclear gyromagnetic ratio, $\gamma{{}_{\mbox{\scriptsize N}}}$, is replaced by the coupling constant, $g{{}_{\mbox{\scriptsize N}}}$.

For clarity, we focus this discussion on the the so-called “axion wind interaction” with effective field $\boldsymbol{D}{{}_{\mbox{\scriptsize wind}}}=-\boldsymbol{\nabla}a(\boldsymbol{r},t)$ and coupling constant $g{{}_{\mbox{\scriptsize aNN}}}$. A number of other possible couplings between nuclear spins and bosonic dark matter fields take a form similar to Eq. (3). These include couplings to the “dark” electric (with coupling constant $g{{}_{\mbox{\scriptsize dEDM}}}$ and effective field $\boldsymbol{D}{{}_{\mbox{\scriptsize dEDM}}}$) and magnetic (with coupling constant $g{{}_{\mbox{\scriptsize dMDM}}}$ and effective field $\boldsymbol{D}{{}_{\mbox{\scriptsize dMDM}}}$) fields mediated by spin-$1$ bosons such as dark photons [5, 31] or a quadratic “wind” coupling to an ALP field (with coupling constant $g{{}_{\mbox{\scriptsize quad}}}$ and effective field $\boldsymbol{D}{{}_{\mbox{\scriptsize quad}}}$) [32]. These are discussed further in the Materials and Methods section.

As the Earth orbits around the Sun (itself moving towards the Cygnus constellation at velocity, $\boldsymbol{v}$, comparable to the local galactic virial velocity $\sim~{}10^{-3}c$), it moves through the galactic dark-matter halo and an interaction between axions and ALPs with a given nucleon $N$, can arise. Assuming that ALPs make up all of the dark matter energy density, $\rho{{}_{\mbox{\scriptsize DM}}}$, and that the dominant interaction with nucleon spins is linear in $\boldsymbol{\nabla}a(\boldsymbol{r},t)$, Eqs. (1) and (2) can be used to write the the effective field as

[TABLE]

Then, given the local galactic virial velocity, $\boldsymbol{v}$, only two free parameters remain in the Hamiltonian in Eq. (3): the coupling constant, $g{{}_{\mbox{\scriptsize aNN}}}$, and the field’s oscillation frequency, $\omega{{}_{\mbox{\scriptsize DM}}}$, fixed by the boson mass. A CASPEr search consists of probing this parameter space over a bosonic mass range defined by the bandwidth of the experiment. In order to calibrate the experiment, we apply known magnetic fields and the experiment’s sensitivity to magnetic fields directly translates to sensitivity to the coupling constant. If no ALP field is detected, upper bounds on the coupling constant can be determined based on the overall sensitivity of the experiment.

2.3 Dark-matter signatures in zero- to ultralow-field NMR

The results presented in this work were obtained by applying techniques of ZULF NMR (see Materials and Methods Sec. 4.1, a review of ZULF NMR can be found in Ref. [23]). The sample – 13C-formic acid, effectively a two-spin 1H–13C system – is pre-polarized in a $1.8$ T permanent magnet and pneumatically shuttled to a magnetically-shielded environment for magnetization evolution and detection.

The spin Hamiltonian describing the system is

[TABLE]

where the electron-mediated spin-spin coupling $J{{}_{\mbox{\scriptsize CH}}}/2\pi\approx 221~{}\mathrm{Hz}$ for formic acid and $\boldsymbol{{I}}$ and $\boldsymbol{{S}}$ are the nuclear-spin operators for 1H and 13C, respectively. Additionally, $\gamma{{}_{\mbox{\scriptsize H}}}$ and $\gamma{{}_{\mbox{\scriptsize C}}}$ are the gyromagnetic ratios of the 1H and 13C spins, $\boldsymbol{B}$ is an applied magnetic field, $g{{}_{\mbox{\scriptsize app}}}$ is the ALP-proton coupling strength, and $g{{}_{\mbox{\scriptsize ann}}}$ is the ALP-neutron coupling strength. We assume $g{{}_{\mbox{\scriptsize app}}}=g{{}_{\mbox{\scriptsize ann}}}=g{{}_{\mbox{\scriptsize aNN}}}$ [22].

In the absence of external fields ($\boldsymbol{B}=\boldsymbol{D}{{}_{\mbox{\scriptsize wind}}}(t)=0$), the nuclear spin energy eigenstates are a singlet with total angular momentum $F=0$ and three degenerate triplet states with $F=1$, separated by $\hbar J{{}_{\mbox{\scriptsize CH}}}$. The observable in our experiment is the $y$ magnetization, leading to selection rules $\Delta F=0,\pm 1$ and $\Delta m_{F}=\pm 1$, as in Ref. [33]. The zero-field spectrum thus consists of a single Lorentzian located at $J{{}_{\mbox{\scriptsize CH}}}/2\pi\approx 221~{}\rm{Hz}$, as shown in Fig. 1(a).

In the presence of a static magnetic field, $\boldsymbol{B}=B{{}_{\mbox{\scriptsize z}}}\boldsymbol{\hat{e}}{{}_{\mbox{\scriptsize z}}}$, applied along $z$, the $m_{F}=0$ states are unaffected, while the $m_{F}=\pm 1$ triplet states’ degeneracy is lifted. The corresponding spectrum exhibits two peaks at $J{{}_{\mbox{\scriptsize CH}}}/2\pi\pm B{{}_{\mbox{\scriptsize z}}}(\gamma{{}_{\mbox{\scriptsize H}}}+\gamma{{}_{\mbox{\scriptsize C}}})/2\pi$, as shown in Fig. 1(b).

So long as $|J{{}_{\mbox{\scriptsize CH}}}|\gg|\gamma{{}_{\mbox{\scriptsize H}}}B|\gg|g{{}_{\mbox{\scriptsize aNN}}}D{{}_{\mbox{\scriptsize wind}}}|$ the $m_{F}=0$ states are unaffected, and the $m_{F}=\pm 1$ states are shifted by

[TABLE]

where $D{{}_{\mbox{\scriptsize z}}}(t)$ is the projection of $\boldsymbol{D}{{}_{\mbox{\scriptsize wind}}}(t)$ along the axis of the applied magnetic field. The time dependence of $\boldsymbol{D}{{}_{\mbox{\scriptsize wind}}}(t)$ leads to an oscillatory modulation of the $m_{F}=\pm 1$ energy levels, giving rise to sidebands around the $J$-coupling doublet as shown in Fig. 1(c). The sidebands are separated from the carrier peaks by $\pm\omega_{\rm{DM}}/2\pi$ and have an amplitude proportional to the modulation index ($g_{\text{aNN}}/\omega_{\rm{DM}}$).

Dark-matter fields with sufficiently strong coupling to nuclear spins can then be detected by searching for frequency-modulation-induced sidebands in the well-defined ZULF NMR spectrum of formic acid.

2.4 Coherent averaging

The expected dark matter coherence time ($\sim 14$ hours for a particle with $19$ Hz Compton frequency) is much longer than the nuclear spin coherence time in 13C-formic acid ($\sim 10$ s). Taking advantage of this mismatch, we introduce the post-processing phase-cycling technique shown in Fig. 2 (see Materials and Methods Sec. 4.3), which consists of incrementally phase shifting the transient spectra and subsequently averaging them together. If the phase increment matches the phase accumulated by the oscillating field between each transient acquisition, the sidebands add constructively. This allows coherent averaging of the complex spectra, such that the signal-to-noise ratio scales as $N^{1/2}$, where $N$ is the number of transients. Because the dark-matter Compton frequency is unknown, it is necessary to repeat this operation for a large number of different phase increments (at least many as the number of transient acquisitions).

2.5 Calibration

The energy shifts produced by $D{{}_{\mbox{\scriptsize z}}}(t)$ [Eq. (6)] are equivalent to those produced by a real magnetic field with amplitude

[TABLE]

Similar relationships for dark-photon and quadratic-wind couplings are provided in the Materials and Methods section.

Based on Eq. (7), the sensitivity of the experiment to dark matter was calibrated by applying a real oscillating magnetic field of known amplitude and frequency and measuring the amplitude of the resulting sidebands in the coherently averaged spectrum. Further details are provided in the Supplementary Materials.

2.6 Search and analysis

The dark matter search data were acquired and processed as described above, but without a calibration AC-magnetic field applied.

For each Compton frequency, the appropriate phase increment is computed, which identifies the corresponding coherently averaged spectrum to be analyzed. The noise in the spectrum defines a detection threshold at the 90% confidence level (further details in S7 in the Supplementary Materials). When the signal amplitude at the given frequency is below the threshold, we set limits on the dark matter couplings to nuclear spins at levels determined by the calibration and effective-field conversion factors (see Materials and Methods Sec. 4.2). If the signal is above the threshold, a more stringent analysis is performed by fitting the coherently averaged spectrum to a four-sideband model. When the fit rules out detection, the threshold level is again used to set limits.

In case of an apparent detection, further repeat measurements would need to be performed to confirm that the signal is persistent and exhibits expected sidereal and annual variations.

2.7 CASPEr-ZULF search results: constraints on bosonic dark matter

The results of the CASPEr-ZULF search for axionlike particles are given in Fig. 3. The frequencies presenting sharp losses in sensitivities at $0.21$, $1.69$, and $2.16$ Hz were the ones for which the nearest optimal phase increment was close to zero, thus presenting maximal-amplitude $J$-coupling peaks, raising the detection threshold (see discussion in Sec. 3 of the Supplementary Materials). The red-shaded area labeled “CASPEr-ZULF” corresponds to upper bounds on nuclear-spin couplings to dark matter consisting of ALPs at the $90\%$ confidence level. This represents our current sensitivity limitation after $850$ $30$-second transient acquisitions using samples thermally polarized at $\sim 1.8$ T. The “CASPEr-ZULF Phase II” line corresponds to the projected sensitivity of a future iteration of this work that will use a more sensitive magnetometry scheme to measure a larger sample with enhanced (non-equilibrium) nuclear spin polarization.

Figures 4 and 5 show the search results for the ALP quadratic interaction and dark photon interactions, respectively. No signal consistent with axion, ALP or dark-photon fields have been observed in the red-shaded areas. The two different limits given in Fig. 5 were obtained using the same data set but analyzed by assuming two orthogonal initial polarizations of the dark-photon field (see S10 in the Supplementary Materials).

In all cases, the search bandwidth was limited from below by the finite linewidth of the $J$-resonance peaks, preventing to resolve sidebands at frequencies lower than $45$ mHz. Due to the finite coherence time of the dark-matter fields (corresponding to $\sim 10^{6}$ oscillations), the bandwidth’s upper limit ($19$ Hz), is the highest frequency which can be coherently averaged after $14$ hours of integration time. The sensitivity fall off is due to the sidebands’ amplitude scaling as the modulation index. Further details are given in S6 and S9 of the Supplementary Materials.

3 DISCUSSION

This work constitutes demonstration of a dark-matter search utilizing NMR techniques with a coupled heteronuclear spin system. The results provide new laboratory-based upper bounds for bosonic dark matter with masses ranging from $1.86\times 10^{-16}$ to $7.85\times 10^{-14}$ eV, complementing astrophysical bounds obtained from supernova SN1987A [35, 26].

Our data analysis provides a method to perform coherent averaging of the bosonic-field-induced transient signals. This method should prove useful for other experiments seeking to measure external fields of unknown frequency using a detector with a comparatively short coherence time. Conveniently, this phase-cycling approach also suppresses the carrier-frequency signals, which would otherwise increase the detection threshold via spectral leakage. As this method is applied during post-processing, it does not require modification of the experiments provided that the data to be analyzed have been time stamped.

3.1 Phase II Sensitivity Improvements

In order to search a greater region of the bosonic dark-matter parameter spaces, several sensitivity-enhancing improvements are planned for the next phase of the experiment. In this work, nuclear spin polarization was achieved by allowing the sample to equilibrate in a $1.8$-T permanent magnet, which yields a 1H polarization $\lesssim 10^{-5}$. For the next phase of the experiment, a substantial sensitivity improvement can be obtained by using so-called “hyperpolarization” methods to achieve much higher, non-equilibrium nuclear-spin polarization. Current efforts are focused on the implementation of non-hydrogenative parahydrogen-induced polarization (NH-PHIP)222Non-hydrogenative parahydrogen-induced polarization methods are often referred to with the acronym SABRE, for Signal Amplification by Reversible Exchange of parahydrogen. [38]. Signal enhancement via NH-PHIP has been demonstrated at zero field [39] and after optimization is expected to increase nuclear spin polarization levels to at least 1%. Because parahydrogen can be flowed continuously into the sample, a steady-state polarization enhancement can be achieved [40], improving the experimental duty cycle. Combining continuous NH-PHIP with a feedback system to produce a self-oscillating nuclear spin resonator [41] could conceivably offer further improvement and simplify data analysis.

Additional sensitivity enhancement will be provided by magnetometer improvements and use of a larger sample. In the experiments reported here, only about $50~{}\mathrm{\mu L}$ of the sample contributed to the signal, which was detected from below with an atomic magnetometer with a noise floor around $10~{}\mathrm{fT}/\sqrt{\mathrm{Hz}}$. With a larger ($\gtrsim 1$ mL) sample hyperpolarized via NH-PHIP detected via a gradiometric magnetometer array with optimized geometry and sensitivity below $1~{}\mathrm{fT}/\sqrt{\mathrm{Hz}}$, we anticipate an improvement by $\gtrsim 10^{5}$ relative to the results presented here.

Future experiments will be carried out with increased integration time. We recall that the bosonic dark-matter fields are coherent for a time $\tau{{}_{\mbox{\scriptsize DM}}}$ on the order of $10^{6}$ periods of oscillation. The phase-cycling procedure depicted in Fig. 2 is valid for data sets with total time less than $\tau{{}_{\mbox{\scriptsize DM}}}$. For integration times longer than the coherence time of the bosonic field T${}_{\text{tot}}>\tau{{}_{\mbox{\scriptsize DM}}}$, the data sets have to be coherently averaged in sets of duration $\tau{{}_{\mbox{\scriptsize\rm{DM}}}}$ using the phase-cycling procedure. This yields T${}_{\text{tot}}/\tau{{}_{\mbox{\scriptsize DM}}}$ sets of coherently averaged data. To profit from longer integration time and further increase the SNR, these sets can be incoherently averaged by averaging their PSDs, yielding an overall SNR scaling as (T${}_{\text{tot}}\tau{{}_{\mbox{\scriptsize DM}}}$)1/4 (see Supplementary Materials in Ref. [18]).

3.2 Complementary Searches

In order to increase the bandwidth of the experiment (see Supplementary Materials S9), we propose complementary measurement procedures. As the amplitude of the sidebands scales as $1/\omega{{}_{\mbox{\scriptsize DM}}}$, the sensitivity of the experiment decreases for higher frequencies. To probe frequencies ranging from $\sim 20$ to $500$ Hz (corresponding to bosonic masses of $\sim 8\times 10^{-14}$ to $2\times 10^{-12}$ eV), it was shown in Ref. [42] that a resonant detection method will be more sensitive than the current frequency-modulation-induced sidebands measurement scheme. Resonant AC fields can induce phase shifts in the $J$-coupling peaks [43]; cosmic fields can induce the same effect. By gradually varying the magnitude of a leading magnetic field, one can tune the splitting of the $J$-coupling multiplets to match the dark-matter field frequency. Such a resonance would manifest itself by shifting the phase of the $J$-coupling peak.

For frequencies below $\sim 45$ mHz (corresponding to bosonic masses $\lesssim 2\times 10^{-16}$ eV), the sidebands are located inside of the $J$-coupling peaks and the experimental sensitivity drops rapidly. This represents the lower limit of the bandwidth accessible by the frequency-modulation-induced sidebands measurement scheme presented in this work. To probe down to arbitrarily low frequencies, another measurement scheme has been implemented based on a single-component liquid-state nuclear-spin comagnetometer [44]. Further details and results of this scheme are presented elsewhere [36].

4 MATERIALS AND METHODS

4.1 Experimental parameters

The experimental setup used in this experiment is the one described in Ref. [45]. Additional descriptions of similar ZULF NMR setups can be found in Refs [46].

The NMR sample consists of $\sim 100$ $\mu$L of liquid 13C-formic acid (13CHOOH) obtained from ISOTEC Stable Isotopes (Millipore Sigma), degassed by several freeze-pump-thaw cycles under vacuum, and flame-sealed in a standard 5 mm glass NMR tube.

The sample is thermally polarized for $\sim 30$ s in a $1.8$-T permanent magnet, after which the NMR tube is pneumatically shuttled into the zero-field region. After the guiding magnetic field is turned off, a magnetic pulse (corresponding to a $\pi$ rotation of the 13C spin) is applied to initiate magnetization evolution.

Following each transient acquisition, the sample is shuttled back into the polarizing magnet and the experiment is repeated. In order to increase the SNR, the transient signals are averaged using the phase-cycling technique described in Sec. 4.3.

4.2 Bosonic dark-matter effective fields

In the case of the ALP-wind linear coupling, the field acting on the 1H-13C spins induces an energy shift equal to the one produced by a magnetic field with amplitude

[TABLE]

where $\omega{{}_{\mbox{\scriptsize DM}}}\approx m{{}_{\mbox{\scriptsize DM}}}c^{2}/\hbar$ is the ALP Compton frequency, $\boldsymbol{k}\approx m{{}_{\mbox{\scriptsize DM}}}\boldsymbol{v}/\hbar$ is the wave-vector ($\boldsymbol{v}$ is the relative velocity), $m{{}_{\mbox{\scriptsize DM}}}$ is the rest mass of the ALP, $\phi$ is an unknown phase, and $\boldsymbol{\hat{e}}{{}_{\mbox{\scriptsize z}}}$ is the axis along which the leading DC-magnetic field is applied.

It is theoretically possible that interaction of nuclear spins with $\boldsymbol{\nabla}a$ can be suppressed [47, 32], in which case the dominant axion wind interaction, referred to as the quadratic wind coupling, is related to $\boldsymbol{\nabla}a^{2}$. In the case of the ALP-wind quadratic coupling the equivalent magnetic field amplitude is:

[TABLE]

where $g{{}_{\mbox{\scriptsize quad}}}$, having dimensions of inverse energy, parameterizes the ALP quadratic coupling strength to nuclear spins.

There are two possible interactions of dark photons with nuclear spins that can be detected with CASPEr-ZULF: the coupling of the dark electric field to the dark electric dipole moment (dEDM) and the coupling of the dark magnetic field to the dark magnetic dipole moment (dMDM). The equivalent magnetic field amplitudes are:

[TABLE]

and

[TABLE]

with coupling constants $g{{}_{\mbox{\scriptsize dEDM}}}$ and $g{{}_{\mbox{\scriptsize dMDM}}}$ (having dimensions of inverse energy) and dark photon field polarization $\boldsymbol{\epsilon}$.

The experimental sensitivity to real magnetic fields then directly translates to sensitivity to the coupling constants $g{{}_{\mbox{\scriptsize aNN}}}$, $g{{}_{\mbox{\scriptsize quad}}}$, $g{{}_{\mbox{\scriptsize dEDM}}}$, and $g{{}_{\mbox{\scriptsize dMDM}}}$. Inverting equations (8)–(11) yields the corresponding conversion factors from magnetic field to the dark-matter coupling constants:

[TABLE]

Here we have used $\gamma{{}_{\mbox{\scriptsize C}}}/2\pi=10.70$ MHz.T*-1* and $\gamma{{}_{\mbox{\scriptsize H}}}/2\pi=42.57$ MHz.T*-1*, $v\approx 10^{-3}c$ and $\rho{{}_{\mbox{\scriptsize DM}}}\approx 0.4$ GeV/cm3. The full derivation of these expressions is given in S2.3 of the Supplementary Materials.

4.3 Signal processing

For each transient acquisition, the sample is prepared in the same initial state, which determines the phase of the $J$-coupling peaks. When averaging the transients together, the $J$-coupling peaks’ amplitude and phase remain constant, while the uncorrelated noise is averaged away, thus increasing the SNR as the square root of the total integration time $\text{T}_{\text{tot}}$, i.e. SNR$(\text{T}_{\text{tot}})\propto\text{T}_{\text{tot}}^{1/2}$.

However, the dark-matter-related information resides not in the $J$-coupling peaks, but in their sidebands. The external bosonic field oscillates at an unknown frequency and its phase at the beginning of each transient acquisition is unknown. This phase directly translates into the phase of the sidebands in the transient spectra: while the phase of the $J$-coupling peaks is identical from one transient acquisition to another, the phase of the sidebands varies. As a result, naively averaging the transient spectra averages the sidebands away, thus removing the dark-matter-related information from the resulting spectrum.

Here we use a post-processing phase-cycling technique which enables coherent averaging of the spectra in the frequency domain, even for transient signals for which no obvious experimental phase-locking can be achieved due to the unknown frequency of the signal. The method is similar to acquisition techniques in which an external clock is used to register the times of the transient acquisitions and post-processing phase shifting of the transient signals is employed to recover the external field’s phase [48, 49].

The method relies on the fact that the bosonic field’s phase at the beginning of each transient acquisition is unknown but not random. Indeed we recall that the bosonic fields remain phase coherent for $\sim 10^{6}$ oscillations, which for frequencies below $19$ Hz is longer than the total integration time ($14$ hours). Thus, precise knowledge of the transient-signal acquisition starting times enables recovery of the phase of the bosonic field.

A full description of this averaging method is given in the Supplementary Materials S4. Each transient spectrum is incrementally phase shifted prior to averaging. If the phase shift is equal to the phase accumulated by the bosonic field between two transient acquisitions, then the phase stability of the $J$-coupling peaks is shifted to their sidebands which can thus be coherently averaged.

Considering that the frequency of the bosonic field is unknown, the correct phase shift is also unknown. Thus, we repeat the operation for $2001$ different phase increments between $[-\pi,\pi]$, yielding $2001$ averaged spectra, one of which being averaged with the phase increment such that the sidebands are coherent.

To demonstrate the viability of this method, a small magnetic field was applied with amplitude $0.24$ nT oscillating at $0.73$ Hz, to simulate a dark-matter field. Using this processing technique, the SNR of the sidebands scales as SNR$(\text{T}_{\text{tot}})\propto\text{T}_{\text{tot}}^{1/2}$ (see Fig. S3 in the Supplementary Materials), as expected during a coherent averaging procedure. This is a dramatic improvement over the alternative power-spectrum averaging (typically implemented for sets of incoherent spectra), which would yield a $\text{T}_{\text{tot}}^{1/4}$ scaling.

4.4 Calibration

The experimental sensitivity is defined by the ability to observe dark-matter-induced sidebands above the magnetometer noise floor. In order to show that the sideband amplitude scales with the modulation index, $B/\omega$, and does not present unusual scalings due to experimental errors, two calibration experiments were performed. A first calibration was performed by varying the amplitude of the AC-field from $95$ to $310$ pT while holding the frequency constant at $\omega=2\pi\times 0.73$ Hz. Then the amplitude was held at $160$ pT while varying the frequency, from $0.45$ to $1.7$ Hz. The results of this experiment are shown in Fig. S1 in the Supplementary Materials. Similar experiments were performed to determine the minimum detectable frequency (see S9 in the Supplementary Materials). Based on this calibration, we can extrapolate the expected sidebands amplitude, $\text{A}_{s}(B,\omega)$, for any field of amplitude $B$ and frequency $\omega$:

[TABLE]

Then the magnetometer noise level determines the smallest detectable driving field, which is then converted to dark-matter coupling bounds via eqs. (12)-(15).

5 REFERENCES

S1 Experimental setup

The apparatus used in this experiment is identical to that described in Ref. [45]. Additional descriptions of similar ZULF NMR setups can be found in Ref. [50].

The zero-field region is achieved using a four-layer $\mu$-metal and ferrite shield (MS-1F shield produced by Twinleaf LLC, $10^{6}$ magnetic shielding factor). In this region, the residual magnetic field is $\lesssim 10$ pT. Three orthogonal pairs of Helmholtz coils which receive current from an amplified source (AE Techron 7224-P) are wound on the holder. This allows for the application of AC and DC magnetic fields along three orthogonal directions. This set of coils is used to apply acquisition-starting pulses. A separate set of coils mounted on the innermost shield is used to apply the shimming, leading DC-, calibration AC-, and benchmark AC-fields.

Prior to every transient acquisition, the sample is thermally polarized for $\sim 30$ s in a $1.8$ T permanent magnet, after which the NMR tube is pneumatically shuttled into the zero-field region. During the shuttling, the sample goes through a magnetic field of a guiding solenoid wrapped around a fiberglass tube, ensuring adiabatic transport into the shielded region [23].

After the shuttling to the zero-field region, the sample stops $\sim 1$ mm above the spin-exchange-relaxation-free magnetometer’s (SERF [51, 52]) rubidium-vapor cell (cell produced by Twinleaf with $500$ torr of nitrogen buffer gas).

In order for the magnetometer to operate in SERF regime, the rubidium cell is maintained at $180$ C*∘* by means of a resistive heater. The $795$ nm circularly-polarized pump beam (from a Toptica Photonics diode laser system, locked to the D$1$ Rubidium line), propagates along the y-axis.

The linearly polarized probe beam (from another Toptica Photonics diode laser system), propagating along the x-axis, is blue-detuned $10$ GHz from the center of the rubidium D$2$ line to probe the rubidium atoms’ polarization. The polarimeter’s analog output signal is digitized with a 24 bit acquisition card (NI 9239, National Instrument) at a $5$ kHz sampling rate. Synchronization of the experiments and control of shuttling, magnetic pulses, and data acquisition is accomplished with a LabVIEW program.

S2 Measurement schemes

S2.1 Zero-field regime

13C-formic acid possesses an electron-mediated spin-spin coupling between the 13C and its neighbouring 1H, referred to as $J$-coupling. In isotropic liquids, rotational motion of the molecules averages the $J$-coupling Hamiltonian to its isotropic form:

[TABLE]

where $J{{}_{\mbox{\scriptsize CH}}}/2\pi$ $\approx$ $221$ Hz for formic acid, $\boldsymbol{I}$ and $\boldsymbol{S}$ are the nuclear spin-1/2 operators for 1H and 13C, respectively.

Prior to the acquisition, the sample is shuttled to the magnetically shielded region where the residual magnetic field is on the order of $10$ pT. In such a low magnetic field, the Zeeman interaction is on the order of $0.1$ mHz and is negligible compared to the $J$-coupling. For isotropic fluids in zero-field, other interactions typically relevant in NMR spectroscopy, such as the short-range dipole-dipole couplings, are averaged-out by molecular tumbling and are also negligible. As a result, the sample’s evolution in zero-field is dominated by the $J$-coupling Hamiltonian given by Eq. (17).

Following shuttling, the guiding field is turned off suddenly and a magnetic field pulse with area $\gamma_{C}B_{x}t_{\mathrm{pulse}}=\pi$ is applied to the sample. As a result, the 13C and 1H nuclear spins are in a superposition of the triplet ${|F\text{=}1;m_{F}\text{=}\pm 1,0\rangle}$ and singlet ${|F\text{=}0;m_{F}\text{=}0\rangle}$ states which evolves under the $J$-coupling Hamiltonian. The unperturbed energy levels are given by:

[TABLE]

where $I=S=1/2$ for 13C and 1H and $F=0,1$. Then the singlet and triplet states energy levels are respectively:

[TABLE]

The selection rules for transition between the singlet and triplet states, $\Delta F=0,\pm 1$ and $\Delta m_{F}=\pm 1$, arise because the observable is the magnetization along the y-axis. The corresponding spectrum exhibits a single peak centred at \Big{[}E(F=1,\pm 1)-E(F=0)\Big{]}/\hbar=J\approx 2\pi\cdot 221 Hz (see Fig. 1.a of the main text).

S2.2 Ultralow-field regime

The experimental setup allows application of a DC-magnetic field, $\boldsymbol{B}{{}_{\mbox{\scriptsize z}}}$, along the z-axis via the Helmholtz coils surrounding the sample. In such conditions, the Hamiltonian under which the formic acid molecules evolve becomes:

[TABLE]

where $\gamma{{}_{\mbox{\scriptsize C}}}/2\pi=10.70$ MHz.T*-1* and $\gamma{{}_{\mbox{\scriptsize H}}}/2\pi=42.57$ MHz.T*-1* are the gyromagnetic ratios of the 13C and 1H nuclear spins, respectively. A field of about $50$ nT yields a Zeeman interaction, $H{{}_{\mbox{\scriptsize Z}}}$, on the order of $\gamma_{C}{B}{{}_{\mbox{\scriptsize z}}}/2\pi\approx\gamma_{H}{B}{{}_{\mbox{\scriptsize z}}}/2\pi\approx 0.1$ Hz, which remains much weaker than the $J$-coupling and can thus be treated as a perturbation. To first order in ${B}{{}_{\mbox{\scriptsize z}}}$, the eigenstates are those of the unperturbed Hamiltonian, $\mathcal{H}{{}_{\mbox{\scriptsize J}}}$, and perturbed energies can be read from the diagonal matrix elements of the Zeeman perturbation:

[TABLE]

yielding the perturbed energy levels:

[TABLE]

Thus, to first order, the effect of such a field is to break the degeneracy of the triplet state due to the now non-negligible Zeeman splitting. Recalling the selection rules $\Delta F=0,\pm 1$ and $\Delta m_{F}=\pm 1$, the magnetometer measures oscillations at two different frequencies between the ${|\text{F=}0;\text{m}{{}_{\mbox{\scriptsize F}}}\text{=}0\rangle}$ and ${|\text{F=}1;\text{m}{{}_{\mbox{\scriptsize F}}}\text{=}\pm 1\rangle}$ states. As a result, the single $J$-coupling line is split into a doublet and the spectrum exhibits two lines located at $J/2\pi\pm B{{}_{\mbox{\scriptsize z}}}\frac{\gamma_{\text{C}}+\gamma_{\text{H}}}{4\pi}$ (see Fig. 1.b of the publication main text)

S2.3 **Frequency-modulation regime **

In this regime, we apply an oscillating magnetic field $\boldsymbol{B}{{}_{\mbox{\scriptsize AC}}}=B{{}_{\mbox{\scriptsize AC}}}\cos(\omega t)\cdot\boldsymbol{\hat{e}}{{}_{\mbox{\scriptsize z}}}$ in addition to the DC-field along z-axis. Under these conditions, the formic acid molecules evolve under the following Hamiltonian:

[TABLE]

In the case where the modulation index, $\frac{B{{}_{\mbox{\scriptsize AC}}}(\gamma_{\text{C}}+\gamma_{\text{H}})}{2\omega}$, is much smaller than one, the resulting spectrum is composed of the $J$-coupling doublet peaks, located at $J{{}_{\mbox{\scriptsize CH}}}/2\pi\pm B{{}_{\mbox{\scriptsize z}}}\frac{\gamma_{\text{C}}+\gamma_{\text{H}}}{4\pi}$ each exhibiting a set of sidebands located at frequencies $J{{}_{\mbox{\scriptsize CH}}}/2\pi\pm B{{}_{\mbox{\scriptsize z}}}\frac{\gamma_{\text{C}}+\gamma_{\text{H}}}{4\pi}\pm\omega/2\pi$. Such a spectrum is shown in Fig. 1.c of the publication main text, alongside the corresponding energy levels. The sidebands’ amplitude $\text{A}{{}_{\mbox{\scriptsize s}}}$, is proportional to the modulation index [42]. This behaviour was experimentally confirmed by varying the amplitude and frequency of a calibration AC-field while a DC-field was continuously applied to the sample. This procedure yields the calibration data, necessary to interpret results from the search for dark-matter bosonic fields.

S3 Dark matter effective fields

In the following, we explicitly state the Hamiltonians describing the interaction of ALPs and dark photons with a single nuclear spin-1/2 having gyromagnetic ratio $\gamma{{}_{\mbox{\scriptsize N}}}$. By analogy to the Zeeman interaction, we then derive the expressions of the corresponding pseudo-magnetic fields acting on a two-spin system such as formic acid. In the following, all expressions are given in SI units.

S3.1 ALP Wind linear coupling to nuclear spins

For a single nuclear spin, axions and ALPs generate a pseudo-magnetic field through what is known as the “axion wind interaction”, described by the nonrelativistic Hamiltonian [19, 17, 5],

[TABLE]

where:

[TABLE]

is the ALP effective field, acting on the nuclear spin, and $\boldsymbol{I}{{}_{\mbox{\scriptsize N}}}$ is the spin operator in units of $\hbar$.

The ALP field, $a(\boldsymbol{r},t)$, can be written as:

[TABLE]

where $\omega{{}_{\mbox{\scriptsize DM}}}\approx m{{}_{\mbox{\scriptsize DM}}}c^{2}/\hbar$ is the ALP Compton frequency, $\boldsymbol{k}\approx m{{}_{\mbox{\scriptsize DM}}}\boldsymbol{v}/\hbar$ is its wave-vector with $\boldsymbol{v}$, the average velocity of the particles in the laboratory frame, $m{{}_{\mbox{\scriptsize DM}}}$ is the rest mass of the ALP particle and $\phi$ is an unknown phase. As the value of $\phi$ has no incidence on the measurement we set its value to zero for the rest of this discussion. Differentiating Eq. (29) yields:

[TABLE]

Recalling that $a_{0}$ is related to the dark-matter density through

[TABLE]

yields:

[TABLE]

We now assume that the ALP field is acting on the nucleons of the coupled 13C and 1H nuclear spins while a weak DC magnetic field is applied to the sample. The axion wind Hamiltonian now becomes

[TABLE]

where, $g{{}_{\mbox{\scriptsize app}}}$ is the ALP-proton coupling strength, $g{{}_{\mbox{\scriptsize ann}}}$ is the ALP-neutron coupling strength. To first order, the $m_{F}=0$ states are unaffected, and the $m_{F}=\pm 1$ states are shifted by

[TABLE]

where we have assumed that $g{{}_{\mbox{\scriptsize app}}}=g{{}_{\mbox{\scriptsize ann}}}=g{{}_{\mbox{\scriptsize aNN}}}$, and $D{{}_{\mbox{\scriptsize z}}}(t)$ is the projection of $\boldsymbol{D}{{}_{\mbox{\scriptsize wind}}}(t)$ along the axis of the applied magnetic field. Thus the perturbation induced by $D{{}_{\mbox{\scriptsize z}}}(t)$ takes a similar form to that of a magnetic field $B{{}_{\mbox{\scriptsize z}}}(t)$ [see Eq. (23) of the publication main text]:

[TABLE]

S3.2 ALP Wind quadratic coupling to nuclear spins

It is theoretically possible that interaction of nuclear spins with $\boldsymbol{\nabla}a$ can be suppressed [47, 32], in which case the dominant axion wind interaction, referred to as the quadratic wind coupling, is related to $\boldsymbol{\nabla}a^{2}$:

[TABLE]

where $g{{}_{\mbox{\scriptsize quad}}}$, having dimensions of inverse energy, parameterizes the ALP quadratic coupling strength to nuclear spins and where we have introduced the effective field due to the quadratic axion coupling:

[TABLE]

Assuming the form of $a(r,t)$ given by Eq. (29):

[TABLE]

we obtain for the effective field:

[TABLE]

As for the linear coupling, the effect of $\boldsymbol{D}{{}_{\mbox{\scriptsize quad}}}$ acting on a single spin is equivalent to that of a magnetic field, whose amplitude in this case is given by:

[TABLE]

Specifically for formic acid, assuming equal coupling of the axion field to the proton and the neutron, we have

[TABLE]

We note that for the quadratic coupling case, the modulation frequency is twice the Compton frequency of the ALPs. Therefore, an experiment sensitive to a certain range of axion masses for the linear coupling is, at the same time, probing quadratic couplings of axions of half the mass.

S3.3 Dark-photon couplings to nuclear spins

There are two possible interactions of dark photons with nuclear spins that can be detected with CASPEr-ZULF: the coupling of the dark electric field to the dark electric dipole moment (dEDM) and the coupling of the dark magnetic field to the dark magnetic dipole moment (dMDM).

Assuming the dark photons make up the totality of the dark matter energy density, $\rho_{\rm{DM}}$, the energy stored in the dark electric field, $\boldsymbol{\boldsymbol{E}}{{}_{\mbox{\scriptsize dEDM}}}(t)$, is set equal to $\rho_{\rm{DM}}$, in a form analogous to that of the usual electromagnetism:

[TABLE]

where $\epsilon_{0}$ is permittivity of free space. The dark magnetic field is given by:

[TABLE]

The Hamiltonians describing the dark EDM and dark MDM are respectively given by:

[TABLE]

and:

[TABLE]

where $g{{}_{\mbox{\scriptsize dEDM}}}$ and $g{{}_{\mbox{\scriptsize dMDM}}}$ (of dimension of inverse energy) parametrize the strength of the dark photon couplings to the EDM and MDM, respectively.

Similarly to axions and ALPs, the above Hamiltonians are expressed by analogy to the Zeeman interaction by introducing the following magnetic fields:

[TABLE]

and

[TABLE]

where $\boldsymbol{\hat{\epsilon}}$ is the bosonic field’s polarization vector. Assuming equal coupling to the proton and neutron, yields for a two-spin (1H-13C) system:

[TABLE]

and

[TABLE]

Note that for ordinary photons, the coupling to the EDM is strongly suppressed since the electromagnetic interaction respects the CP symmetry. However, the dark photon need not respect the CP symmetry, in which case there would be no associated suppression [37].

S4 Signal processing: post-processing phase-cycling

The averaging method is illustrated in Fig. 2 in the publication main text. Subsequent transient acquisitions are separated by a time interval $\tau$. The $n^{\text{th}}$ transient acquisition yields a transient spectrum denoted FFTn. The operation is repeated $N$ times, yielding a set of $N$ transient spectra: \Big{\{}\text{FFT}_{n} where ${n}=1,2...$ {N}\Big{\}}. Once every $\text{FFT}_{n}$ has been acquired, a phase shift, $\phi$, incremented by the current acquisition number $n$, is applied to each $\text{FFT}_{n}$:

[TABLE]

The $N$ transient phase-shifted $\text{FFT}_{n}^{\phi}$ are then averaged:

[TABLE]

If the phase shift matches the accumulated phase of the dark-matter field between each transient acquisition,

[TABLE]

where $\omega{{}_{\mbox{\scriptsize DM}}}$ is the oscillation frequency of the dark-matter field, the sidebands will be averaged coherently. Generally, the “carrier” $J$-coupling peaks will be averaged to zero, with the exception of cases where $\omega_{DM}=2\pi/\tau$. Because the sidebands are coherently averaged, their SNR scales as SNR$(\text{T}_{\text{tot}})\propto\text{T}_{\text{tot}}^{1/2}$ while the $J$-coupling peaks are averaged away along with the uncorrelated noise.

The data shown in Fig. 2 of the main text were acquired in an experiment during which the AC-field frequency was set to $0.73$ Hz. The time interval between each measurement was set to $\tau=61$ s (including $30$ s of acquisition and $31$ s for polarization, shuttling and pulsing, combined). The optimal phase shift is computed using Eq. (52), yielding $\phi\approx 2.93$ rad. Note that as the frequency of the bosonic fields is unkown, this optimal phase shift cannot be computed during the actual dark-matter-search run. However, this optimal phase shift must necessarily lie within the [$-\pi,\pi$] interval. In practice we numerically try $2001$ phase increments within [$-\pi,\pi$] : \Big{\{}\phi_{k}=k\pi/1000 where k=-1000,-999,...1000\Big{\}}.. This method produces a set of $2001$ phase-shifted averaged FFTs: $\langle\text{FFT}^{\phi_{k}}\rangle$.

Figure 2 illustrates the result of the phase-shifting processing method for an experiment where the transient acquisition is repeated $850$ times. The $\langle\text{FFT}^{\phi_{k}}\rangle$ are ordered by phase increment $\phi_{k}$, and stacked in a two-dimensional plot. This plot is then examined and searched for sidebands. The phase increment of $\phi=0$ rad corresponds to averaging without applying any phase shift. In this case, the averaging is optimized for the $J$-resonances and sidebands cannot be seen (i.e., they are destructively averaged). However, a specific phase increment of $\phi=2.93$ rad (highlighted in orange), yields maximum sideband amplitude while the $J$-coupling peaks are averaged out. This phase value (and its negative counterpart) corresponds to the optimal value for which the benchmark-field coherence is restored and yields maximum SNR, matching the computed value using Eq. (52).

S5 Sideband amplitude determination

A benchmark sideband measurement was performed using a test field of amplitude $B_{0}=0.24$ nT, oscillating at the frequency $\omega_{0}/2\pi=0.73$ Hz. The transient acquisition data, consisting of $30$ s time-series sampled at $5$ kHz, are Fast-Fourier Transformed. The transient FFTs are multiplied by a calibration function expressing the magnetometer sensitivity within the [math] to $500$ Hz bandwidth. The calibration function is obtained by measuring the magnetometer response to an applied frequency-varying magnetic field of fixed amplitude.

Following the calibration, the baseline of the FFTs is fit and subtracted. The transient FFTs are then averaged using the phase-shifting scheme described in the previous section.

The next operation consists in extracting the sideband amplitudes from the averaged spectrum. To this end, six complex Lorentzians are fit to the averaged spectrum (two for the $J$-peaks, each possessing two sidebands):

[TABLE]

where $\alpha$, $\Gamma$, $\nu_{\text{c}}$ and $\phi$ are the amplitudes, linewidths, center frequencies and phases of the Lorentzians, respectively, and are real, freely fitted parameters. The benchmark sideband amplitude, centered at $\nu_{\text{c}}$ is given by $\text{A}_{s_{0}}=|L(\nu_{\text{c}})|$.

S6 Calibration: signal scaling versus bosonic-field amplitude and frequency

In order to experimentally verify that the sideband amplitude scales with the modulation index, $\text{A}_{s}(\text{B},\omega)\propto\text{B}/\omega$, two calibration experiments were performed.

A first calibration was performed by varying the amplitude of the AC-field with constant frequency. The results of this experiment are shown in Fig. S1 and show the linear dependence of the sideband amplitude on the amplitude of the AC-field. Each data point corresponds to the measured sideband amplitude (as defined is the previous section) for AC-calibration fields of amplitude varying from $0.095$ to $0.31$ nT at fixed frequency $0.73$ Hz.

The second calibration was performed by varying the frequency of the AC-calibration field from $0.45$ to $1.7$ Hz with constant amplitude of $0.16$ nT. The results of this calibration (Fig. S1) show the $1/\omega$ dependence of the sideband amplitude with the frequency of the AC-field.

Both calibration experiments were carried out with $2100$ seconds total integration time for each data point (corresponding to $70$ transient acquisitions of $30$ seconds).

S7 Detection threshold determination

The detection threshold, $\text{A}{{}_{\mbox{\scriptsize th}}}$, is defined as the required signal amplitude such that we can claim detection with a $90\%$ confidence in the case: $\text{A}_{s}(\text{B},\omega)>$ $\text{A}{{}_{\mbox{\scriptsize th}}}$, where $\text{A}_{s}(\text{B},\omega)$ is the measured signal as described in Eq. (16). To define the detection threshold, we use the standard $p$-value test described in Ref. [53]. To this aim we consider a complex spectrum containing noise only:

[TABLE]

where $\chi_{1}(\omega)$ and $\chi_{2}(\omega)$ are the real and imaginary part of the FFT, respectively, and $\omega/2\pi$ is the probed frequency. The signal detection threshold can be defined by studying the power spectral density (PSD) of this noisy spectrum:

[TABLE]

To this aim we use the result that if the Fourier coefficients $\chi_{1}(\omega)$ and $\chi_{2}(\omega)$ are normally distributed with zero mean, variance $\sigma^{2}$ and are frequency independent (consistent with white noise), then the PSD cumulative distribution function (CDF) follows an exponential distribution [53]. Hence, at any frequency $\omega/2\pi$, the probability that the measured power P is smaller than $z_{0}$ is:

[TABLE]

In practice, the value of $\sigma^{2}$ is obtained by fitting the evaluating the PSD CDF of a noise-only spectrum (in the present case via the MATLAB Kaplan-Meier estimator). We then fit the PSD CDF to the exponential distribution given in Eq. (56). An illustration of this procedure is shown in Fig. S2.

We now define a signal power $z_{0}$ such that, if at the probed frequency $\omega/2\pi$, the power of a spectrum containing both the signal and noise, P$(\omega)$, is greater that $z_{0}$, a detection can be claimed with a false-alarm probability $p_{0}$ (a false alarm corresponds to the case where the measured power P$(\omega)$ is due to noise fluctuations). This yields:

[TABLE]

Solving the above equation for $z_{0}$ yields:

[TABLE]

Setting now the false-alarm probability to $p_{0}=0.10$ (corresponding to a confidence level of $90\%$) we find the signal power threshold is terms of the noise fluctuation:

[TABLE]

As our analysis is done on the spectrum amplitude, defined as the square-root of the PSD, we define the detection threshold as:

[TABLE]

Any signal with amplitude $\text{A}{{}_{\mbox{\scriptsize s}}}(B,\omega)>\text{A}{{}_{\mbox{\scriptsize th}}}$ is inspected for detection. The case where no such signal is observed, yields the upper-bound exclusion line:

[TABLE]

Using Eq. (16) this limit can be expressed in term of magnetic field versus frequency:

[TABLE]

Such analysis is done in a frequency-dependent manner by computing the phase increment corresponding to the studied frequency. The averaged spectrum obtained using this phase increment is then studied as explained above to define the signal-detection threshold at $\omega$. The results of this analysis are given in Fig. S2.c, where Eq. (65) is evaluated for every frequency within the bandwidth.

The exclusion plots shown in Figs. 3, 4 and 5 in the publication main text, are obtained by converting $\omega$ to bosonic masses in eV and $B{{}_{\mbox{\scriptsize th}}}(\omega)$ to coupling constants in GeV*-1* using Eqs. (12), (13), (14) and (15) in the publication main text.

The exclusion lines show the 2.14$\sigma$ confidence level of the experiment ($90\%$ CL). Such frequency-dependent analysis yields higher detection threshold for frequencies for which the phase increment is close to [math] rad. Indeed, the spectra obtained from the [math] radian phase increment contains maximum amplitude $J$-couplings, which then raise the detection threshold to a higher value due to spectral leakage. This effect manifests itself in the exclusion lines by sporadic and sharp increases of the detection threshold.

We note that in the literature, similar analyses are done by applying a statistical penalty for searching multiple frequencies within a bandwidth of interest to account for the look-elsewhere effect (corresponding to raising Eq. (56) to the power of the number of frequencies probed). Our analysis does not require such a penalty. Indeed, the characteristic signature of the bosonic fields, a set of four sidebands around the $J$-coupling frequencies, allows one to differentiate false alarms from true detection events.

S8 Coherent averaging: signal scaling with integration time

Figure S3 shows the sidebands amplitude $\text{A}{{}_{\mbox{\scriptsize s}}}$, detection threshold $\text{A}{{}_{\mbox{\scriptsize th}}}$, and signal-to-noise ratio, SNR:= $\text{A}{{}_{\mbox{\scriptsize s}}}/\text{A}{{}_{\mbox{\scriptsize th}}}$, scaling with the total integration time using the phase-shifting procedure described in S4.

During this experiment, the AC-field is maintained at a frequency $\omega{{}_{\mbox{\scriptsize AC}}}/2\pi=0.73$ Hz and amplitude ${B}{{}_{\mbox{\scriptsize AC}}}=0.24$ nT. The experiment is repeated with increased number of transient acquisitions each of length $30$ s. All values are extracted using the averaged spectrum obtained from the optimal phase increment corresponding to $\omega{{}_{\mbox{\scriptsize AC}}}$.

Figure S3 shows that the sideband amplitude remains constant while the signal detection-thresholds decreases as $\text{A}{{}_{\mbox{\scriptsize th}}}($T${{}_{\mbox{\scriptsize tot}}})\propto$ T${{}_{\mbox{\scriptsize tot}}}^{-1/2}$ as expected while performing a coherent averaging. The SNR scales accordingly as SNR $\propto$ T${{}_{\mbox{\scriptsize tot}}}^{1/2}$.

S9 Bandwidth: accessible bosonic mass range

The bandwidth of the experiment is limited from below by the linewidth of the $J$-coupling peaks. For bosonic fields with frequencies lower than $45$ mHz (corresponding to bosons with mass $\lesssim 1.86\times 10^{-16}$ eV), the corresponding sidebands are located inside the $J$-resonance and cannot be resolved. As a result, we limit the lower end of the search bandwith to $45$ mHz, below which calibration could not be achieved (see Fig. S4). In principle, the bandwidth extends to frequencies up to $\sim 500$ Hz, corresponding to the maximum detectable frequency of the magnetometer (see Materials and Methods Sec. 4.1). However, the data are acquired during $\sim 14$ hours (corresponding to $850$ transient acquisitions each $30$ seconds long, with a duty cycle of $50$% imposed by the polarization step between each acquisition). The $850$ transient signals are coherently averaged for as long as the bosonic field remains phase-coherent. The coherence time of the bosonic fields is on the order of $10^{6}$ periods of oscillation, which for fields of $\sim 19$ Hz corresponds to $\sim 14$ hours. This limits the search bandwidth to fields of frequencies below $19$ Hz (corresponding to bosons with mass lower than $7.8\times 10^{-14}$ eV).

S10 Dark-matter-field directionality

S10.1 Daily and annual modulations

We note that the ALP field described in Eq. (29), seen in the Earth’s reference frame, has annual and daily modulations due to the motion of the Earth in the Solar System. These modulations are expected to 1- modulate the orientation of the experiment’s sensitive axis with respect to the dark-matter wind, thus modulating the overall experimental sensitivity and 2- change the lineshape of the experimental power spectra from pure Lorentzians to skewed Lorentzians. These effects are addressed and accounted for during the data processing by weighting transient spectra during the averaging sequence using the following method.

For axions and ALPs, the experiment is sensitive to the projection of the bosonic pseudo-magnetic field gradient, $\boldsymbol{\nabla}a(\boldsymbol{r},t)$, onto the direction of the leading magnetic field. We denote this axis in the laboratory frame by the unit vector $\boldsymbol{\hat{e}}{{}_{\mbox{\scriptsize z}}}$. As a result of the yearly revolution of the Earth around the Sun and the daily revolution of the Earth, the alignment of $\boldsymbol{\hat{e}}{{}_{\mbox{\scriptsize z}}}$ with $\boldsymbol{\nabla}a(\boldsymbol{r},t)$ varies over time. This effect can increase or totally suppress the expected signal depending on when transient acquisitions happen and must be included in the analysis to avoid overestimating the experiment’s sensitivity.

To this aim, we use the analysis proposed in Ref. [54], originally applied to nuclear recoil in WIMP-nucleus scattering but which directly translates to this sensitive-axis misalignment. Detailed and comprehensive calculations are given in Ref. [54], here we show how these results are applied to our current signal analysis. We recall that the ALP scalar-field gradient can be written as (ignoring its initial phase):

[TABLE]

where $\omega{{}_{\mbox{\scriptsize DB}}}=\sqrt{(m_{\text{DM}}c^{2})^{2}+(m_{\text{DM}}vc)^{2}}/\hbar$ is the bosonic particle de Broglie frequency, $\boldsymbol{k}=m_{\text{DM}}\boldsymbol{v}/\hbar$ is its wave-vector with $\boldsymbol{v}$, the average velocity of the particles in the laboratory frame and $m_{\text{DM}}$ is the mass of the particle. When shifting from the galactic-center to the laboratory frame, $\boldsymbol{v}$ can be written as the following vectorial sum of velocities:

[TABLE]

where $\boldsymbol{v}_{\text{DM}}$ is the velocity of the local bosonic field with respect to the center of the galaxy and is taken to be zero on average due to the isotropic structure of the dark-matter halo, $\boldsymbol{v}_{{}_{\text{SG}}}$ is the velocity of the Sun with respect to the center of the galaxy, i.e. $\boldsymbol{v}_{{}_{\text{SG}}}$ is the velocity of the Sun towards the Cygnus constellation and is approximatly constant, $\boldsymbol{v}_{{}_{\text{ES}}}(t)$ is the velocity of the center of the Earth with respect to the sun and $\boldsymbol{v}_{{}_{\text{LE}}}(t)$ is the velocity of the laboratory with respect to center of the Earth. The expected signal amplitude $S,$ is proportional to the projection of the field’s gradient onto the sensitive axis:

[TABLE]

During the calibration and benchmark acquisitions, an oscillating-magnetic field is continuously applied along the $\boldsymbol{\hat{e}}{{}_{\mbox{\scriptsize z}}}$ direction. After averaging the benchmark sideband amplitude $\text{A}_{\text{s}_{0}}$ is obtained from the fitting procedure. To avoid overestimating the expected signal amplitude and account for the previously discussed astrophysical motions, we weigh the transient spectra by the factor $\frac{|\boldsymbol{v}(t_{n})\cdot\boldsymbol{\hat{e}}{{}_{\mbox{\scriptsize z}}}|}{|v|}\approx\frac{|\boldsymbol{v}(t_{n})\cdot\boldsymbol{\hat{e}}{{}_{\mbox{\scriptsize z}}}|}{10^{-3}c}$ before including them in the average spectra, where $t_{n}$ is the starting time of the $n^{\text{th}}$ transient acquisition. This replicates the daily and yearly amplitude modulations of the signal induced by the motion of the Earth in the Solar System such as to avoid overestimating the expected signal amplitude from the bosonic fields.

During the search acquisitions, no oscillating-magnetic field is applied and this misalignment effect is naturally included in the data. Nonetheless we weight the transient spectra by the factor $\frac{\boldsymbol{v}(t_{n})\cdot\boldsymbol{\hat{e}}{{}_{\mbox{\scriptsize z}}}}{|v|}$ before averaging them. Indeed, a transient spectrum acquired when $\boldsymbol{v}(t_{n})\cdot\boldsymbol{\hat{e}}{{}_{\mbox{\scriptsize z}}}=0$, cannot possibly include any signal and thus contains only noise. The effect of weighting such a transient spectrum by this factor is to prevent the addition of noise in the averaged spectrum. On the contrary, spectra acquired when $\boldsymbol{v}(t_{n})\cdot\boldsymbol{\hat{e}}{{}_{\mbox{\scriptsize z}}}$ is maximal possess higher weight to account for the fact that the signal should be maximal during the acquisition.

For dark photons, the experiment is sensitive to the projection of the field’s polarization, denoted by the unit-vector $\hat{\boldsymbol{\epsilon}}$, onto the leading field direction $\boldsymbol{\hat{e}}{{}_{\mbox{\scriptsize z}}}$. We note that, at any point in space and time, $\hat{\boldsymbol{\epsilon}}$ can take any direction and is uniformly distributed along the unit-3D-sphere. Moreover, $\hat{\boldsymbol{\epsilon}}$ is assumed to be constant during one coherence time, which is always longer that the total measurement time. Therefore the experiment probes a unique value of $\hat{\boldsymbol{\epsilon}}$, modulated by the annual and daily motions of the sensitive axis of detection.

The analysis is performed assigning three equally probable orientations $\hat{\boldsymbol{\epsilon}}_{i}$; and by time propagating $\hat{\boldsymbol{\epsilon}}_{i}\cdot\boldsymbol{\hat{e}}{{}_{\mbox{\scriptsize z}}}(t_{0})$. This quantity is then used as a weight during the averaging of transient signal similar to the ALP case, thus yielding three corresponding limits.

To obtain limits comparable with other experiments, the three $\hat{\boldsymbol{\epsilon}}_{i}$, are given in the non-rotating Celestial Frame. The axis of this frame are practically fixed on any relevant time scale, since it takes the Solar System $\sim 250$ million years to complete one orbit about the galactic centre, so different experiments can be compared easily in this case (even if they are performed many years apart in time).

The three assigned directions are the usual orthonormal axis of the non-rotating Celestial Frame. The first assigned orientation, $\hat{\boldsymbol{\epsilon}_{1}}$, lies in the Celestial Equator Plane and points towards the Sun at the Vernal Equinox. The second orientation, $\hat{\boldsymbol{\epsilon}_{2}}$, points towards the North Celestial Pole (i.e. along the axis of rotation of the Earth). In this case the sensitivity of the experiment is null, as $\boldsymbol{\hat{e}}{{}_{\mbox{\scriptsize z}}}$ is pointing East in the laboratory frame and remains perpendicular to $\hat{\boldsymbol{\epsilon}_{2}}$ at any point in time. The third orientation is the following vector product $\hat{\boldsymbol{\epsilon}_{3}}=\hat{\boldsymbol{\epsilon}_{1}}\times\hat{\boldsymbol{\epsilon}_{2}}$, and lies in the Celestial Equator Plane. A complete mathematical description and schematic of this coordinate system can be found in Ref. [54].

S10.2 Search data time-stamp

The search data are composed of $850$ transient scans, each $30$ s long, spaced by a polarization time of $31$ s. The first transient scan was acquired on the $28^{\text{th}}$ of July $2018$ at $10$:$00$:$06$ (CET). Data acquired in Mainz, Germany.

S10.3 Expected sideband lineshape

The oscillating bosonic fields described in the Earth frame also exhibit annual and daily modulations due to the presence of the modulated wave-vector $\boldsymbol{k}$, and de Broglie frequency $\omega{{}_{\mbox{\scriptsize DB}}}$ in the arguments of the cosine function. The effect of such modulations is to modify the lineshape of the power spectra from the expected Lorentzian signal. The first modulation due to the motion of the Earth with respect to the galactic center enters in the de Broglie frequency $\omega{{}_{\mbox{\scriptsize DB}}}=\sqrt{(m_{\text{DM}}c^{2})^{2}+(m_{\text{DM}}vc)^{2}}/\hbar$, where $\boldsymbol{v}$ is modulated as in Eq. (67). The amplitude of the velocities in Eq. (67) are $|\boldsymbol{v}_{\text{DM}}|\approx 0$, $|\boldsymbol{v}_{{}_{\text{SG}}}|\approx 220$ km.s${}^{-1}\approx 10^{-3}c$, $|\boldsymbol{v}_{{}_{\text{ES}}}(t)|\approx 30$ km$/$s and $|\boldsymbol{v}_{{}_{\text{LE}}}(t)|\approx 0.5$ km$/$s [54]. Those components are therefore negligible compared to the rest mass-energy of the particle, so the usual approximation $\omega{{}_{\mbox{\scriptsize DB}}}\approx m_{\text{DM}}c^{2}/\hbar$ is used. The second modulation appears in the spatial argument of the cosine function, $\boldsymbol{k}\cdot\boldsymbol{r}$. This component is dominated by the velocity of the Sun towards the Cygnus constellation, $|\boldsymbol{v}_{{}_{\text{SG}}}|$, which is taken as constant both in amplitude and direction. Therefore, we neglect other components of the wave-vector. These two approximations yield a non-modulated cosine-form gradient, thus yielding a pure Lorentzian power-spectral line.