Search for axion-like dark matter using solid-state nuclear magnetic resonance

We report the results of an experimental search for ultralight axion-like dark matter in the mass range 162 neV to 166 neV. The detection scheme of our Cosmic Axion Spin Precession Experiment (CASPEr) is based on a precision measurement of $^{207}$Pb solid-state nuclear magnetic resonance in a polarized ferroelectric crystal. Axion-like dark matter can exert an oscillating torque on $^{207}$Pb nuclear spins via the electric-dipole moment coupling $g_d$, or via the gradient coupling $g_{\text{aNN}}$. We calibrated the detector and characterized the excitation spectrum and relaxation parameters of the nuclear spin ensemble with pulsed magnetic resonance measurements in a 4.4~T magnetic field. We swept the magnetic field near this value and searched for axion-like dark matter with Compton frequency within a 1~MHz band centered at 39.65~MHz. Our measurements place the upper bounds ${|g_d|<7.0\times10^{-4}\,\text{GeV}^{-2}}$ and ${|g_{\text{aNN}}|<2.1\times10^{-1}\,\text{GeV}^{-1}}$ (95\% confidence level) in this frequency range. The constraint on $g_d$ corresponds to an upper bound of ${7.6\times 10^{-22}\,\text{e}\cdot\text{cm}}$ on the amplitude of oscillations of the neutron electric dipole moment, and ${3.2\times 10^{-6}}$ on the amplitude of oscillations of CP-violating $\theta$ parameter of quantum chromodynamics. Our results demonstrate the feasibility of using solid-state nuclear magnetic resonance to search for axion-like dark matter in the nano-electronvolt mass range.

We report the results of an experimental search for ultralight axion-like dark matter in the mass range 162 neV to 166 neV. The detection scheme of our Cosmic Axion Spin Precession Experiment (CASPEr) is based on a precision measurement of 207 Pb solid-state nuclear magnetic resonance in a polarized ferroelectric crystal. Axion-like dark matter can exert an oscillating torque on 207 Pb nuclear spins via the electric-dipole moment coupling g d , or via the gradient coupling g aNN . We calibrated the detector and characterized the excitation spectrum and relaxation parameters of the nuclear spin ensemble with pulsed magnetic resonance measurements in a 4.4 T magnetic field. We swept the magnetic field near this value and searched for axion-like dark matter with Compton frequency within a 1 MHz band centered at 39.65 MHz. Our measurements place the upper bounds |g d | < 7.0 × 10 −4 GeV −2 and |g aNN | < 2.1 × 10 −1 GeV −1 (95% confidence level) in this frequency range. The constraint on g d corresponds to an upper bound of 7.6 × 10 −22 e · cm on the amplitude of oscillations of the neutron electric dipole moment, and 3.2 × 10 −6 on the amplitude of oscillations of CP-violating θ parameter of quantum chromodynamics. Our results demonstrate the feasibility of using solid-state nuclear magnetic resonance to search for axion-like dark matter in the nanoelectronvolt mass range.
The existence of dark matter is indicated by astronomical and cosmological evidence, but its interactions, aside from gravity, remain undetected [1,2]. A number of theoretical models of physics at high energies, such as string theory, grand unified theories, and models with extra dimensions, incorporate light pseudoscalar bosons (axion-like particles, ALPs), which are potential dark matter candidates [3][4][5][6][7]. Among these, the axion is particularly compelling, because it also offers a solution to the strong CP problem of quantum chromodynamics (QCD) [7][8][9][10][11]. The axion or axion-like field a(t) = a 0 cos (ω a t) oscillates at the Compton frequency ν a = ω a /(2π) = m a c 2 /h, where c is the speed of light in vacuum, h is the Planck constant, and m a is the unknown ALP It was electrically polarized along the cylinder axis, indicated with the black arrow. The pickup coil and the cancellation coil were coaxial with the crystal, and the axis of the Helmholtz excitation coil was orthogonal. The vertical leading magnetic field B0 set the direction of the equilibrium spin polarization. Coils were supported by G-10 fiberglass cylinders shown in gray and pink. (b) Electrical schematic, showing the excitation and pickup circuits. Excitation pulses generated with the digital-to-analog converter (DAC) were amplified (Ae), and coupled to the excitation coil via a tuned tank circuit that included matching and tuning capacitors, as well as a resistor to set the circuit quality factor. The pickup probe was also designed as a tuned tank circuit, coupling the voltage induced in the pickup coil to a low-noise cryogenic amplifier (A1), whose output was filtered, further amplified, and digitized with an analog-to-digital converter (ADC). (c) Pulsed NMR sequence used for FID measurements. The spin-ensemble equilibrium magnetization, initially parallel to B0, was tilted into the transverse plane by the excitation pulse. The FID signal was recorded after the excitation pulse, as the magnetization precessed and its transverse component decayed. mass, which can be in a broad range, roughly between 10 −21 eV and 10 −3 eV [12][13][14]. The field amplitude a 0 is fixed by the assumption that it dominates the dark matter energy density: ρ DM = m 2 a a 2 0 /2 ≈ 3.6 × 10 −42 GeV 4 [15,16]. Kinetic energy of the axion-like dark matter field introduces small corrections to its frequency spectrum. The standard halo model predicts the spectral shape with linewidth (v 2 0 /c 2 )ν a ≈ 10 −6 ν a , where v 0 ≈ 220 km/s is the circular rotation speed of the Milky Way galaxy at the Sun's location [17,18].
Experimental searches for axion-like particles rely on symmetry arguments about the nature of their interactions with Standard Model particles [7,16,19,20]. These interactions are suppressed by a large energy scale, set by the decay constant f a , which could lie near the grand unification, or the Planck scale [21]. Most experiments to date have focused on the electromagnetic interaction, which can mix photons with axions and ALPs in the presence of a strong magnetic field [22][23][24][25][26][27][28][29][30][31][32]. The Cosmic Axion Spin Precession Experiments (CASPEr) search for different interactions: the electric dipole moment (EDM) interaction and the gradient interaction with nu-clear spin I [19,[33][34][35][36][37]. The gradient interaction Hamiltonian is H aNN = g aNN ∇a · I, where g aNN is the coupling strength. The EDM interaction arises from the defining coupling of the axion to the gluon field [38]. Its Hamiltonian can be written as H EDM = g d aE * · I/I, where g d is the coupling strength and E * is an effective electric field [19]. This interaction is equivalent to that of a parity-and time-reversal-violating oscillating EDM, given by d = g d a 0 cos (ω a t). This corresponds to an oscillating QCD θ parameter: θ(t) = (a 0 /f a ) cos (ω a t), with g d inversely proportional to f a [16,39]. The EDM coupling generates axion mass, and for the QCD axion m a ≈ Λ 2 QCD /f a , where Λ QCD ≈ 200 MeV is the QCD confinement scale [16,40].
The sensitivity of static EDM experiments to the oscillating EDM is suppressed, although data re-analysis has produced limits at low frequencies [41,42]. Astrophysical constraints can be derived by analyzing the cooling dynamics of the supernova SN1987A [16,43]. Constraints can also be extracted from analysis of 4 He production during Big Bang nucleosynthesis [44] and from analysis of black hole superradiance [45]. CASPEr-electric is a direct, model-independent search for the EDM and gradient interactions of axion-like dark matter, with the potential to reach the sensitivity to the QCD axion [19]. We search for the effects of these interactions on the dynamics of a spin ensemble in a solid with broken inversion symmetry [46][47][48][49][50][51][52]. The measurements focus on 207 Pb 2+ ions, with nuclear spin I = 1/2, in a poled ferroelectric PMN-PT crystal with the chemical formula: . The non-centrosymmetric position of the ions in this crystal gives rise to a large effective electric field, analogous to the effect in polar molecules [54][55][56]. The EDM or gradient interaction with axion-like dark matter creates an oscillating torque on the nuclear spins. We quantify the magnitude of this torque by the Rabi frequency Ω a , which is proportional to the corresponding interaction strength. For a spin ensemble polarized by an external bias magnetic field, this torque tilts the spins, if it is resonant with their Larmor frequency. The experimental observable is the oscillating transverse magnetization: where M 0 is the equilibrium magnetization of the 207 Pb nuclear spin ensemble, T 2 is the nuclear spin coherence time, and u is a dimensionless spectral factor that takes into account the inhomogeneous broadening of the spin ensemble and the detuning between the ALP Compton frequency and the spin Larmor frequency [53]. Our apparatus makes use of inductive detection to measure the 207 Pb spin precession, Fig. 1(a). We poled the cylindrical PMN-PT crystal along its axis, aligned with the [1,1,1] crystal direction. This created the axial effective electric field E * , proportional to the remanent polarization P r . We mounted the crystal inside a fiberglass tube, so that E * was perpendicular to the vertical bias magnetic field B 0 , created with a superconducting solenoid. A pickup coil, wound around the tube, was coupled to a low-noise cryogenic preamplifier with a tuned matching circuit, Fig. 1 We calibrated the pickup probe using 207 Pb pulsed nuclear magnetic resonance (NMR) measurements, Fig. 1(c). The spins were excited by resonant magnetic field pulses, created by delivering current to the 2 × 3-turn Helmholtz excitation coil, coupled to a matching circuit, tuned at 42 MHz with a quality factor 2. The axis of this coil was orthogonal to the pickup coil axis, Fig. 1(a). After each pulse, nuclear spin free induction decay (FID) was measured with the pickup probe, characterized by transfer coefficient where V 1 is the recorded voltage referred to the amplifier input, M 1 is the transverse sample magnetization, and µ 0 is the permeability of free space. Despite our efforts to minimize the inductive and capacitive couplings between the excitation and the pickup coils, we found that the cryogenic preamplifier saturated during excitation pulses, and its recovery time was too long to observe the fast FID decay [53]. To address this problem, we placed a single-turn cancellation coil near the pickup coil, Fig. 1(a), and delivered to it a compensating current during the excitation pulses. The amplitude and phase of this compensating current were chosen to cancel the current in the pickup probe during excitation, and prevent preamplifier saturation, without affecting spin excitation. This scheme is a substitute for the transmit/receive switch, often used in NMR detectors.
We performed the NMR calibration measurements at the leading magnetic field B 0 = 4.4 T, for which the value of the equilibrium thermal magnetization M 0 of the spin ensemble was µ 0 M 0 = 2.9 nT. Before every FID measurement the spin ensemble magnetization was initialized to (1.9 ± 0.2) nT by saturating the spins, then letting magnetization recover over approximately one population relaxation time [53]. We set the excitation carrier frequency to 39.71 MHz, and recorded the FID signals after excitation pulses of variable width. The Fourier spectrum of one of these FID signals is shown in Fig. 2(a). We modeled the FID lineshapes by numerically solving the Bloch equations for a spin ensemble with an inhomogeneously-broadened excitation spectrum [53]. By fitting the data, we extracted the transverse coherence time of the nuclear spins: T 2 = (16.7 ± 0.9) ms, and the pickup-circuit transfer coefficient α = (2.3 ± 0.2) × 10 4 V/T. We note that there is a sharp central feature with linewidth on the order of the Rabi frequency, but the overall FID spectral width is much greater than 1/T 2 , 4 FIG. 2. Sensitivity calibration. (a) Measurements of 207 Pb FID following a spin excitation pulse of length tp = 20 ms. The excitation carrier frequency was set to 39.71 MHz, and the Rabi frequency was Ωe = 0.88 rad/ms. The data points show the in-phase (blue circles) and the out-of-phase (orange squares) quadratures of the Fourier transform of the detected voltage, referred to the input of the pickup probe amplifier A1. Data points were binned and averaged, the error bars show one standard deviation for each bin. The lines show the best-fit simulation of the spin response, with the light-colored narrow bands indicating the range of simulation results if parameters are varied by one standard deviation away from their best-fit values. We performed the fitting simultaneously to three FID data sets, with excitation pulse lengths tp = 0.2 ms, 2 ms, 20 ms, with free parameters including the spin coherence time T2 and pickup circuit transfer coefficient α [53]. (b) Measurement of the normalized 207 Pb NMR excitation spectrum near Larmor frequency 39.71 MHz. Excitation pulses of length 1.6 ms and Rabi frequency Ωe = 0.88 rad/ms were delivered at the carrier frequencies shown on the x-axis. Data points show the amplitude of the spin FID response, normalized so that the integral of the spectrum is unity. The error bars indicate one standard deviation uncertainties of the FID spectrum fits. We model the excitation spectrum as a super-Gaussian of order 2 (red line) [53]. (c) Detector calibration for varying drive Rabi frequency. Data points show the amplitude of the spin FID response after an excitation pulse of length 20 ms, delivered at the carrier frequency 39.71 MHz, with Rabi frequency Ωe plotted on the x-axis. The error bars indicate one standard deviation uncertainties, obtained by grouping 100 consecutive FID measurements taken at each Ωe into 5 sets, and independently analyzing each set [53]. The orange line shows the spin response simulated using the Bloch equations with parameters extracted from data in panel (a). (d) Measurement of ferroelectric hysteresis in the PMN-PT single crystal. The remanent polarization Pr persists after the applied voltage has been ramped down to zero. since the tilting pulse excites a broad frequency band within the inhomogeneous spin distribution. The exact shape of the FID Fourier spectrum depends on the interplay between the excitationpulse spectrum, the distribution of tipping angles across the spin ensemble, and the T 2 coherence time.
We measured the inhomogeneous broadening of the 207 Pb nuclear spins in the sam-ple by sweeping the excitation pulse carrier frequency and recording the corresponding FID spectra. The resulting NMR excitation spectrum was centered at 39.71 MHz and had a full width Γ/(2π) = (78 ± 2) kHz, Fig 2(b). This broadening is consistent with the chemical shift anisotropy (CSA) of 207 Pb observed in solid-state NMR [57]. We measured the population relaxation time T 1 of the 207 Pb nuclear spin ensemble with a saturation-recovery measurement, obtaining T 1 = (25.8 ± 0.6) min [53].
The spin evolution in our pulsed NMR calibration measurements was more complicated than the CW-like small spin-tip angle response to axionlike dark matter, described by Eq. (1). In order to confirm the validity of our NMR model in the limit of small spin-tip angles, we recorded and analyzed FID data for a range of excitation Rabi frequencies Ω e . For these measurements we kept the excitation pulse width at 20 ms -approximately the coherence time of axion-like dark matter field with Compton frequency near 40 MHz. At small excitation amplitudes, the spin response was linear in Ω e , as described by Eq. (1) for the case of the drive due to interaction with axionlike dark matter, Fig 2(c). The slope of the linear response is proportional to the spectral factor u = (3.8 ± 0.3) × 10 −4 , which is well approximated by the ratio of the homogeneous linewidth π/T 2 and the inhomogeneously-broadened excitation spectrum width Γ [53]. The deviation from linearity at larger Ω e is due to saturation of the resonant spins in the excitation spectrum, consistent with our Bloch-equation simulations.
Prior to any measurements, the PMN-PT crystal was ferroelectrically poled at room temperature by applying 3.5 kV across the crystal faces. We measured the ferroelectric hysteresis loop by sweeping the applied voltage while recording the current flowing through the sample, and integrating it to find the polarization, Fig. 2(d). The resulting value of remanent polarization was P r = (22 ± 2) µC/cm 2 . We recorded hysteresis data before and after the experiments searching for axion-like dark matter, and verified that the fractional degradation of polarization due to thermal cycling and fatigue was smaller than the quoted uncertainty. The effective electric field E * is proportional to the ferroelectric polarization [48,54,55]. In order to calculate the value of E * we considered the Schiff moment S of the 207 Pb nucleus, induced by the oscillating QCD θ parameter [58,59]. The dominant contribution to the Schiff moment arises from the parity-and timereversal-violating nuclear forces, resulting in the value S = 0.04θ e · fm 3 [53, 60-64]. This corresponds to the magnitude of effective electric field E * = 340 kV/cm. We estimate the theoretical uncertainty in E * on the level of 50% [53].
In order to search for axion-like dark matter we swept the leading magnetic field B 0 in 21 steps, corresponding to the search frequency range 39.1 MHz to 40.2 MHz. The step size was chosen to correspond to 50 kHz, on the order of the width of the 207 Pb nuclear spin excitation spectrum, Fig. 2(b). The broad NMR excitation spectrum reduced the necessary number of magnetic field steps for a given search frequency range. At each value of B 0 we recorded 58 s of scan data sensitive to axion-like dark matter, followed by 58 s of re-scan data that were used in our analysis to identify statistical fluctuations. In order to confirm the experimental calibration, we performed pulsed NMR measurements at three values of the leading field, corresponding to the extremes and the midpoint of the search frequency range [53].
Data analysis consisted of several processing, correction, and signal-search steps. At each value of the leading field B 0 we divided the recorded scan data into 27 blocks, each of 2.15 s duration, chosen to be much longer than the ≈ 25 ms coherence time of any potential ALP dark matter signal in our frequency range. We used the pickup-circuit transfer coefficient α to convert the recorded voltage values to magnetization, and performed a discrete Fourier transform on each block, subsequently averaging the power spectral densities (PSDs) of the blocks. Many of the spectra were contaminated with narrowband RF interference that penetrated our electromagnetic shielding. We used Savitzky-Golay digital filtering to identify and reject these narrowband features, while preserving potential axion-like dark matter signals, whose spectral shape is predicted by the standard halo model [25,53,65].
We then processed the data to search for signals due to the EDM and the gradient interactions. The first step was optimal filtering, performed by convolving the PSD with the signal lineshape predicted for the corresponding interaction [53]. At each value of B 0 we retained the optimally-filtered data points in a frequency bin, centered at the corresponding Larmor frequency, with full width 80 kHz, covering the excitation spectrum bandwidth. We modeled the histogram of these data points as the normal distribution with standard deviation σ, inset of Fig. 3. We set the candidate detection threshold to 3.355σ, equivalent to 95% confidence interval for a 5σ detection, and flagged all points above the threshold as candidates [32,53,65].
There were 617 candidates for EDM coupling (636 for gradient coupling). In order to reject residual RF interference, we used the fact that RF pickup is independent of the leading field B 0 , while an axion-like dark matter signal should only appear when B 0 is tuned to a value such that the spin excitation spectrum overlaps with the ALP Compton frequency. We compared the candidates from data sets taken at different values of B 0 , rejecting 569 candidates for EDM coupling (577 for gradient coupling). The remaining 48 candidates for EDM coupling (59 for gradient coupling) were shown to be statistical fluctuations, using a scan/re-scan analysis [53]. The search sensitivity was limited by the ≈ 0.05 nV/ √ Hz input noise level of the amplifier, corresponding to a magnetic field sensitivity of ≈ 2 fT/ √ Hz.
Our search did not yield a discovery of the EDM coupling g d or the gradient coupling g aNN of axionlike dark matter. In the absence of a detection, in each frequency bin the 95% confidence interval limit on magnitudes of these coupling constants corresponds to the 5σ value in the Gaussian distribution of the optimally-filtered PSD [32,53,65]. The limits were normalized by the NMR excitation spectrum for each bin and concatenated to produce constraints on g d and g aNN over the entire frequency search range, Fig. 3. Over the frequency range 39.1 MHz to 40.2 MHz the constraint on |g d | is 7.0 × 10 −4 GeV −2 , corresponding to an upper bound of 7.6 × 10 −22 e · cm on the amplitude of oscillations of the neutron electric dipole moment, and 3.2 × 10 −6 on the amplitude of oscillations of the QCD θ parameter. The constraint on |g aNN | is 2.1 × 10 −1 GeV −1 . We are not aware of any existing experimental limits on these interactions in this ALP mass range. Analysis of cooling dynamics of supernova SN1987A can be used to estimate bounds g d 10 −8 GeV −2 and g aNN 10 −9 GeV −1 [19,24,43]. However these model-dependent bounds are subject to significant caveats and uncertainties, and may be evaded altogether, reinforcing the importance of laboratory searches [66, 67]. Stringent experimental limits on g d and g aNN exist at much lower ALP masses [35,36,41,42,[68][69][70][71].
There are several ways to improve experimental sensitivity to axion-like dark matter. Since the CSA-induced inhomogeneous broadening is proportional to the Larmor frequency, searching in a lower ALP mass range will reduce the linewidth and therefore improve the search sensitivity. A search in the lower mass range will likely also benefit from superconducting detectors, such as SQUIDs and quantum upconverters [72]. Manipulation of light-induced transient paramagnetic centers may enable control over the nuclear spin population-relaxation time T 1 , and nuclear spin hyperpolarization using dynamic polarization techniques. A dramatic sensitivity improvement could be achieved by scaling up the sample volume.
We estimate that with a sample size of ≈ 80 cm, it may be possible to reach the sensitivity necessary to detect the QCD axion g d coupling strength in the mass range between ≈ peV and ≈ 5 neV.
The authors thank Oyku Acican for her help with Fig. 1

A. Description of the apparatus
Our cryogenic nuclear magnetic resonance (NMR) setup is inside a liquid helium (LHe) bath cryostat with a solenoidal superconducting magnet (Cryomagnetics, Inc. Model 90-300-010), Fig. S1. The apparatus is built around a crystal that is inductively coupled to a pickup probe along one axis, and an excitation probe along an orthogonal axis, both in the plane transverse to the leading magnetic field created by the magnet (Fig. 1(a) in the main text). The experimental setup is used both when measuring pulsed NMR and when performing the axion search. During pulsed NMR calibration measurements, a digital-to-analog converter (DAC) generates a radio frequency (RF) voltage waveform V e , which is coupled into the excitation probe (Fig. S2). The resulting RF magnetic field exerts a torque on the spins, whose magnitude is quantified by the excitation Rabi frequency Ω e . The excitationprobe transfer function is defined as The excitation pulse tilts 207 Pb nuclear spins into the plane transverse to the leading field B 0 , creating a crystal magnetization M 1 that rotates at the Larmor frequency. After the excitation pulse ends, this magnetization decays (free induction decay, FID). The magnetization induces an oscillating current in the pickup coil, and voltage V 1 at the input of the low-noise preamplifier A 1 (Fig. S2). The pickup probe transfer function is defined as where µ 0 is the permeability of free space. The preamplifier A 1 has gain of 40. Its output is connected to a low-pass filter LP 1 and another amplifier stage A 2 (gain = 15) mounted inside the cryostat near the top flange. After a third amplifier stage A 3 (gain = 10) outside the cryostat, the signal is sent to an analog-to-digital converter (ADC).
The excitation signal is routed through a switch S 1 (Fig. S2) that is controlled with a transistor-transistor logic (TTL) pulse with the same duration as the excitation RF pulse. This prevents the DAC output noise from coupling into the pickup probe after the end of the excitation pulse, during FID detection. When the TTL state is high at 5 V, the DAC is connected to the excitation probe through amplifier A e , and when the TTL state is low at 0 V the input of A e is connected to ground via a 50 Ω termination.

B. The crosstalk minimization scheme
During experimental assembly we carefully adjust the orthogonal axes of the excitation and the pickup coils to minimize mutual inductance between them, to ≈ 1.5% of its maximum value for parallel axes. Despite these efforts, the excitation pulse induces a crosstalk current in the pickup coil, with amplitude and phase depending on the residual inductive and capacitive couplings between the coils. This crosstalk signal saturates the preamplifier, resulting in a recovery time of ≈ 200 µs, which complicates the detection of the FID signal. In order to prevent saturation, during the excitation pulse we apply a waveform to the cancellation coil that is optimized to compensate the crosstalk current in the pickup probe with minimal effect on spin dynamics. The phase and amplitude of this waveform are optimized by monitoring the current measured at the pickup probe and minimizing its magnitude. This is done at zero leading magnetic field to avoid spin excitation during optimization. We emphasize that only a small (< 1.5%) fraction of the excitation pulse RF field couples into the pickup probe and has to be compensated, therefore our compensation scheme has a correspondingly small effect on spin dynamics.
In many room-temperature NMR measurements, preamplifier saturation is prevented by using a transmit/receive (T/R) switch between the pickup probe and the preamplifier. Because our preamplifier is at 4.2 K temperature, we chose to use the compensation scheme discussed above, rather than designing and constructing a cryogenic T/R switch.

C. Tuning and matching of pickup probe, excitation probe, and cancellation coil
Our magnetic resonance probes are designed as series capacitance-tuned tank circuits, Fig. S2. In these circuits, coil inductance L is in parallel with a tuning capacitor C 1 and a resistor R, and this tank circuit is in turn in series with a matching capacitor C 2 . The total probe impedance is In order to match the probe impedance to Z = R 0 = 50 Ω at the resonance angular frequency ω r , we have to choose the following values for the circuit elements: where Q is the resonance quality factor. We used fixed-value surface mount capacitors and resistors, so the probes are not tunable after the setup is assembled. The pickup coil with N p = 9 turns of 26 AWG (American Wire Gauge) copper wire is a solenoid with return path cancellation. It has a radius r p = 3.2 mm, an inductance of L p = 0.5 µH, and is tuned to the resonant frequency ω p /(2π) = 39.71 MHz with quality factor Q p = 26 at 4.2 K. The width of the pickup probe resonance limits the frequency range over which we can search for axion-like dark matter without re-tuning the probe. This is why we limited Q p to 26. The excitation coil has a Helmholtz geometry with N e = 2 × 3 turns of 26 AWG copper wire with radius r e = 7.1 mm and inductance of L e = 0.3 µH, which is tuned to the resonant frequency ω e /(2π) = 42.01 MHz with quality factor Q e = 1.5 at 4.2 K. The cancellation coil is a single turn loop around r c = 4.8 mm radius with 26 AWG copper and inductance of L c = 0.01 µH, which is tuned to resonant frequency ω c /(2π) = 40.31 MHz with quality factor Q c = 2 at 4.2 K. All the probes are matched to 50 Ω.

D. Estimates of the pickup probe transfer function α and the excitation probe transfer function κ
Based on the electrical schematic described above, let us estimate the values of the transfer functions α and κ, defined in Eqs. (S1,S2). Using Faraday's law we can estimate the voltage induced in the pickup coil by an oscillating transverse magnetization M 1 : where B 0 = 4.4 T is the leading static magnetic field, r s = 2.3 mm is the radius of the sample, and γ is the 207 Pb gyromagnetic ratio. We set the demagnetizing factor to 1/3, as for a sphere, as an approximation for a cylindrical sample with height ≈ diameter. For the pickup probe on resonance with the spin Larmor frequency, the resulting voltage at the input of preamplifier A 1 is calculated from circuit analysis [1]: while Q p ω p L p ≫ R 0 . We can therefore estimate the pickup probe transfer function: For an excitation voltage V e , referred to the output of the DAC, the current through the excitation coil is calculated from circuit analysis: where |A e | = 4 is the gain of the SRS SIM954 amplifier. Note that the SIM954 has an output impedance 3.3 Ω ≪ 50 Ω. The magnetic field produced by this current at the center of the Helmholtz excitation coil is Assuming the excitation is resonant with the spin Larmor frequency, the Rabi frequency is Ω e = γ(B e /2). The factor of 1/2 arises because only one circular component of the linearly polarized excitation magnetic field B e is resonant (rotating wave approximation). Therefore the excitation probe transfer function can be estimated as Section II C describes how we used pulsed NMR to measure the values of α and κ. The proximity of the measured values to the estimates above validates the approximations used when analyzing the apparatus design shown in Fig. S2.

E. Shielding to reduce RF interference
The probes are mounted on a G-10 fiberglass cylinder, with a 0.2-mm thick copper sheet wrapped around the outside. The cylinder is positioned inside the magnet bore, Fig. S1. The copper sheet forms a closed shielding enclosure, together with aluminum end caps on top and bottom. The RG316 coaxial cable between the pickup probe and the low-noise amplifier A 1 is shielded with a 0.5-mm thick copper mesh sleeve. Another copper mesh sleeve shields the bundled RG316 coaxial cables that run up to the top flange of the cryostat. Shields are connected to the ground pin of the A 1 amplifier used as a common ground. We estimate the RF interference noise reduction factor due to the shields to be on the order of 10 within the 1 MHz range centered at 39.71 MHz.

A. Properties of the 207 Pb spin ensemble
The 207 Pb isotope has nuclear spin I = 1/2 and gyromagnetic ratio where µ = 0.5926µ N is the magnetic moment of 207 Pb nucleus [2], and the nuclear magneton is µ N /h = 7.6226 MHz/T [3]. The chemical formula of PMN-PT is (PbMg 1/3 Nb 2/3 O 3 ) 2/3 − (PbTiO 3 ) 1/3 . The number density of 207 Pb nuclear spins in a PMN-PT crystal is given by: where ρ = 8.2 g/cm 3 [4] is the mass density, M = 317.9 g/mole [2] is the molar mass, and N A is the Avogadro constant. The natural abundance of 207 Pb is 22.1% [2]. We perform our experiments in the leading magnetic field B 0 = 4.4 T and at temperature T = 4.2 K. The equilibrium magnetization of the 207 Pb nuclear spin ensemble is given by [5] where k B is the Boltzmann constant, µ 0 is the permeability of free space, and is the reduced Planck constant. We model the NMR excitation spectrum as a super-Gaussian distribution of order 2, given by where Γ/(2π) is the full-width at half-maximum, ν is the excitation frequency, and ν 0 is the center of the distribution. The scaling pre-factor is chosen to ensure that the area under the distribution is normalized to 1.

B. Saturation-recovery measurements of the relaxation time T1
We use the standard NMR saturation recovery scheme to measure the T 1 relaxation time of the 207 Pb nuclear spin ensemble in PMN-PT at 4.2 K. Each measurement begins with a saturation step, comprising 100 consecutive repetitions of a sequence with 101 pulses whose carrier frequencies vary across the width of the excitation spectrum from 39.66 MHz to 39.76 MHz, and whose Rabi frequencies are fixed at 0.88 rad/ms. Each pulse duration is 0.8 ms, and the pulse spacing is 1.4 ms. Bloch-equation simulations confirm that this step saturates the spin ensemble, Fig. S3. The saturation step is followed by a variable recovery wait time t, after which a pulsed NMR measurement is performed, with spin FID recorded after excitation pulses of 20 ms duration and 180 ms repetition time. The dependence of the FID amplitude on recovery time t is modeled as an exponential 1 − e −t/T1 . The best-fit value for the population relaxation time is T 1 = (25.8 ± 0.6) min.

C. Spin-dynamics simulations with Bloch equations
We use the Bloch equations to quantitatively describe the magnetic resonance dynamics of the 207 Pb nuclear spin ensemble [5,6]. We choose the direction of the z-axis to be along the static magnetic field B 0 . The linearly-polarized excitation magnetic field B e = (2Ω e /γ) cos (ω 1 t) is applied in the x-direction. In the reference frame that rotates at the angular frequency ω 1 around the leading magnetic field, the Bloch equations read where ∆ω = ω 1 − ω 0 is the detuning of spin Larmor frequency ω 0 from the rotating frame frequency, M ′ 0 is the initial ensemble magnetization, T 2 is the transverse spin coherence time, andM x,y are the transverse spin magnetization components in the rotating frame. The transformation between magnetization in the laboratory and the rotating frames is where in the lab frame M = M xx + M yŷ + M zẑ . We numerically solve the Bloch equations using the Runge-Kutta method. The inhomogeneously-broadened spin ensemble is represented by 3251 spins, with their Larmor frequencies uniformly distributed in an excitation bandwidth of 65 kHz with 0.02 kHz spacing. We simulate the dynamics of each spin independently, and add their contributions to obtain the total magnetization.
The simulation parameters are the spin coherence time T 2 , and the transfer functions α and κ, defined in Sec. I D. We perform fits to experimental FID spectra, shown in Fig. 2(a) of the main text and Fig. S4, by varying the values of these parameters to achieve the minimum value of the goodness-of-fit parameter χ 2 = χ 2 1 + χ 2 2 + χ 2 3 , where the subscript enumerates the measurements with different pulse duration t p = 0.2 ms, 2 ms, 20 ms. For each measurement i = 1, 2, 3 where F exp is the Fourier transform of the experimentally detected voltage and F sim is the Fourier transform of simulation results, converted into voltage using the transfer coefficient α, and the index ν labels discrete frequency points within the window shown in Fig. 2(a) of the main text and Fig. S4. The real part of the Fourier transform corresponds to the in-phase quadrature, and the imaginary part corresponds to the out-of-phase quadrature of the FID, relative to the carrier phase of the excitation pulse. The excitation pulses induce probe ringing with time constant ≈ 500 ns, therefore we use the FID response data starting at 5 µs after the end of an excitation pulse. To improve the signal-to-noise ratio, we average the recorded FID response data for several consecutive excitation pulses: 10 data sets are averaged for t p = 0.2 ms pulse duration, 4 data sets are averaged for t p = 2 ms pulse duration, and 4 data sets are averaged for t p = 20 ms pulse duration. After performing the discrete Fourier transform, data points are binned along the frequency axis, with 4 points per bin for t p = 0.2 ms pulse duration, 2 points per bin for t p = 2 ms pulse duration, and 2 points per bin for t p = 20 ms pulse duration. The error bars shown in Fig. 2(a) of the main text and Fig. S4 are the standard deviation of the points within each bin.
The spin ensemble was saturated before every FID measurement, and the FID measurements started after a wait time ≈ T 1 after saturation. Therefore the initial magnetization at the start of every FID measurement was µ 0 M ′ 0 = (0.67 ± 0.05)µ 0 M 0 = (1.9 ± 0.2) nT, where M 0 is the thermal equilibrium ensemble magnetization given by Eq. (S13).
Using the measurements shown in Fig. 2(a) in the main text and Fig. S4, we extract the best-fit parameter values: T 2 = (16.7 ± 0.9) ms, κ = (0.352 ± 0.007) rad/(ms · V), (S18) S4. Measurements of 207 Pb FID spectra following a spin excitation pulse of length tp, as indicated in the panels. We performed fitting simultaneously to in-phase (blue) and out-of-phase (orange) components of Fourier transforms of averaged FID from three data sets: with excitation pulse duration tp = 20 ms shown in Fig. 2(a)

D. NMR response as a function of the Rabi frequency Ωe
In order to confirm the validity of our NMR model in the limit of small spin-tip angles, we record and analyze FID data for a range of excitation Rabi frequencies Ω e . For these measurements we keep the excitation pulse width at 20 ms -approximately the coherence time of an axion-like dark matter field with Compton frequency near 40 MHz. We vary the Rabi frequency from 0.02 rad/ms to 0.88 rad/ms. At each Rabi frequency, we apply 100 consecutive excitation pulses, spaced by 180 ms. After each pulse, we sample the FID voltage, starting 5 µs after the end of the pulse, and lasting for 16.4 µs. We average the 100 FID data sets, and calculate the discrete Fourier transform F [n] of the averaged FID, where index n labels frequency points. Since we only sample the beginning of the FID, before it can start to decay, we model it as a sinusoidal signal at the excitation carrier frequency. We extract the amplitude of the spin ensemble transverse magnetization by numerically integrating the power spectrum |F [n]| 2 over a 400 kHz-wide frequency band centered at the excitation carrier frequency, and using the pickup probe transfer function α to convert the voltage to magnetization. Uncertainties are calculated using bootstrapping: we group the 100 FID data sets into 5 sets of 20 and perform analysis on these 5 sets independently. Error bars are set at the standard deviation of the results for these 5 sets. To obtain the theory curve in Fig. 2(c) of main text, we use our Bloch equation model to generate numerical time-domain FID data, which we analyze in the same way as we analyze experimental data.

III. SPECTRAL PROPERTIES OF THE CW NMR RESPONSE
Under CW excitation with Rabi frequency Ω e and carrier angular frequency ω 1 , the steady-state transverse magnetization of an unsaturated homogeneously-broadened spin ensemble is given by [5] where M 0 is the longitudinal magnetization, T 2 is the transverse coherence time, ω 0 is the Larmor angular frequency, and L is the Lorentzian lineshape function: (S20) Let us describe the spin ensemble inhomogeneous broadening with the excitation lineshape h(ω 0 + ∆), normalized such that Under CW excitation, the steady-state transverse magnetization is then where the spectral u factor is given by the integral over the lineshape: Let us estimate the value of u. Our NMR measurements indicate that the excitation spectrum is much broader than 1/T 2 , therefore we can approximate the Lorentzian with the delta-function: L(ω 0 + ∆ − ω 1 ) ≈ (π/T 2 )δ(ω 0 + ∆ − ω 1 ). Furthermore, we approximate the excitation spectrum as a rectangular function, centered at ω 0 , with full width Γ and height 1/Γ. Then, provided |ω 0 − ω 1 | < Γ/2, we can approximate In order to more accurately determine u, we solved the Bloch equations with the experimentally-determined values T 2 = (16.7 ± 0.9) ms and excitation spectrum with Γ/(2π) = (78 ± 2) kHz ( Fig. 2(b) in the main text). We obtained in agreement with the estimate in Eq. (S24).

IV. FERROELECTRIC POLARIZATION OF PMN-PT
We polarize the ferroelectric PMN-PT crystal by applying a voltage across its faces at room temperature. To ensure good electrical contact, we paint the faces with graphite paint, which is removed after polarization. We connect the crystal to the Trek model 610E-G-CE high-voltage amplifier as shown in Fig. S5(a). The amplifier measures the applied voltage and the current through the sample. In order to measure the ferroelectric hysteresis loop, we apply triangular voltage ramps with alternating polarities, Fig. S5(b). Current spikes are visible when the applied voltage is sufficient to reverse the ferroelectric polarization. In this experimental run the crystal started with a remanent polarization corresponding to positive polarity, so there is no current spike during the first ramp. We obtain the sample polarization by integrating the current: where q(t) is the electric charge on the crystal surface and r s = 2.3 mm is the base radius of the cylindrical sample. The hysteresis loop shown in Fig. 2(d) of the main text is the plot of polarization as a function of applied voltage. The remanent polarization P r persists after the voltage has been ramped down to zero. We verified that the remanent polarization does not decay after thermal cycling of the sample.

A. P,T-odd axion-like dark matter physics
Axion-like cold dark matter is a classical field: a(t) = a 0 cos (ω a t), where ω a ≈ m a c 2 / . If the axion-field energy density dominates dark matter, then ρ DM = m 2 a a 2 0 /2 ≈ 3.6 × 10 −42 GeV 4 [7]. In the QCD Lagrangian, this gives rise to an oscillating θ angle: Let us consider the nucleon EDM induced by axion-like dark matter: calculated with 40% accuracy [8,9]. Here g d is the EDM coupling constant [9], introduced in the Lagrangian term: where Ψ N is the nucleon wavefunction, F µν is the electromagnetic field tensor, and σ and γ are the standard Dirac matrices. From eqs. (S27,S28) we get the relationship between g d and f a : where we used the natural unit conversions: 1 cm = 5 × 10 13 GeV −1 and e = 0.303. For the QCD axion, the decay constant is related to its mass: m a = 6 × 10 −10 eV 10 16 GeV but for a generic ALP there is no such connection.

B. Nuclear Schiff moments induced by the EDM coupling of axion-like dark matter
The nuclear Schiff moment [10][11][12][13] is defined as: where e is the elementary electric charge, Z is the atomic number, and r k = r k ρ(r)d 3 r are the integrals over nuclear charge density ρ(r). The Schiff moment sources the P-and T-odd electrostatic potential ϕ(r) = 4π(S · ∇)δ(r).
Importantly, the definition of the Schiff moment in Ref. [14] differs from this one by a factor of 4π. We adopt the definition in Eq. (S32), noting the factor of 4π wherever we refer to Ref. [14]. The Schiff moment can be induced by a permanent EDM of a nucleon, or by P,T-odd nuclear forces [14]. The contribution of P,T-odd nuclear forces is larger than the contribution of nucleon EDM [12]. Let us consider the two contributions separately, in the case of the 207 Pb nucleus, whose ground state is I π = 1/2 − , having a neutron 3p 1/2 hole in a closed-shell magic nucleus.

Nuclear Schiff moments induced by P,T-odd nuclear forces
The P,T-odd nuclear interaction of a non-relativistic nucleon with nuclear core is parametrized by strength η [12]: where G F ≈ 10 −5 GeV −2 is the Fermi constant, m is the nucleon mass, σ is its spin, and ρ(r) is the density of core nucleons. A vacuum θ angle gives rise to this interaction via the P,T-odd pion-nucleon coupling constant [14,17]: Next we need to calculate the nuclear Schiff moment induced by the interaction (S37). Reference [12] states that the Schiff moment is suppressed by a factor ∼ 10 for nuclei with a valence neutron, compared to a valence proton, and only core polarization leads to a non-zero effect. For example, the Schiff moment of 201 Hg is estimated as 0.2 × 10 −8 η e · fm 3 . However in Ref. [18] it was realized that virtual excitations in the core eliminate this suppression, and, in fact, the results for a valence neutron and a valence proton should be comparable. Here the Schiff moment of 201 Hg is estimated as 2.4 × 10 −8 η e · fm 3 , and the Schiff moment of 199 Hg is estimated as −1.4 × 10 −8 η e · fm 3 . The issue is complicated by nuclear many-body effects. These were numerically calculated for 199 Hg in Refs. [19,20], giving a factor ∼ 10 reduction in the Schiff moment. However the physical origin of such a strong reduction is not clear. The only effect, not included in the shell model, that can change the value of the Schiff moment is the collective nuclear octupole deformation, and, if anything, that should increase the Schiff moment. Reference [21] gives a result for 199 Hg that is ∼ 10% away from the shell-model estimate. These authors attribute the Schiff moment suppression in Ref. [20] to the mixing with the J π = 1/2 − 2 state, for which they get a small Schiff moment value. However this small value itself is questionable. This state is an admixture of a soft quadrupole phonon (J = 2) to the ground state, resulting still in J = 1/2. The excited states do not have this quadrupole deformation, therefore the overlap matrix elements are likely to be small unless a lot of excited states are carefully taken into account. This suggests that the calculation may have large intrinsic uncertainties.
Importantly, 207 Pb is close to a magic nucleus, which means that many-body effects should not play an important role here. Therefore, until the many-body effects can be better understood, for 207 Pb we retain the single-particle estimate of Ref. [18]: Note that this is a factor of eight larger than the result (13) in Ref. [22], where the 207 Pb Schiff moment was taken to be the same as for the many-body suppressed 199 Hg. We can also see that this contribution is a factor of eight larger than the EDM contribution in Eq. (S36). We therefore neglect the EDM contribution, and use the above estimate (S39). Similar estimates were performed for 199 Hg in Ref. [23].

C. Nuclear Schiff moment-induced spin energy shift in ferroelectric PMN-PT
The energy shift of each nuclear spin sublevel of a 207 Pb 2+ ion in ferroelectric PbTiO 3 is estimated in Refs. [24,25]. The result of the full quantum chemistry calculation [26] is: where x is the displacement of the Pb 2+ ion with respect to the center of the oxygen cage, S is the magnitude of the Schiff moment of the 207 Pb nucleus, and a B = 0.53 Å is the Bohr radius. The nuclear spin is I = 1/2; each of the two nuclear spin states shifts by this amount, in opposite directions. Since θ and S exhibit sinusoidal time dependence, the experimentally relevant quantity is the Rabi angular frequency: where we used = 6.58 × 10 −16 eV · s. We note that the spin driving field is "linearly polarized", and therefore the Rabi frequency contains an extra factor of 1/2, which arises because only one of the two counter-rotating components of the linearly polarized drive is resonant (rotating wave approximation). Density functional theory calculations for PMN-PT give the Pb 2+ cation displacement from the center of the oxygen cage: x 0 = 0.39 Å, and the average polarization: P 0 = 55 µC/cm 2 [27]. Our experiment was performed with the crystal polarization P r = 22 µC/cm 2 , therefore we scale the average displacement to x = 0.16 Å.
For 207 Pb in ferroelectric PMN-PT we can use Eq. (S39,S40) and x = 0.16 Å to get: To connect with the EDM d n and the coupling constant g d , we use Eqs. (S27,S28,S30). For the energy shift we obtain We can extract the effective electric field (which includes the Schiff screening factor [28]): For the drive Rabi frequency we obtain: where g d is in GeV −2 and a 0 = √ 2ρ DM /m a is in GeV. Let us introduce the sensitivity factor ξ, defined as Ω a = ξg d a 0 . Its estimated value is therefore There are several contributions to the theoretical uncertainty in E * and ξ. The uncertainty of the QCD calculations is ≈ 40% [8,9]. The uncertainty of the solid-state calculation of the nuclear spin energy shift due to the Schiff moment is ≈ 30% [24][25][26]. Therefore we estimate the total theoretical uncertainty in E * and ξ at ≈ 50%.

VI. NUCLEAR SPIN DYNAMICS DUE TO THE GRADIENT INTERACTION WITH AXION-LIKE DARK MATTER
The non-relativistic Hamiltonian for the gradient interaction of spin I with axion-like dark matter field a(r, t) is where g aNN is the coupling strength measured in units of GeV −1 , and we used natural units here = c = 1 [9,29]. In the first approximation we can write the axion-like dark matter field as: where the field amplitude a 0 is fixed by the assumption that it dominates the dark matter energy density: ρ DM = m 2 a a 2 0 /2 = 3.6 × 10 −42 GeV 4 [7,9]. We approximate the instantaneous value of the gradient ∇a ≈ m a va, where v is the instantaneous value of the velocity of the ALP field in the laboratory frame. The Hamiltonian in natural units becomes: The product g aNN a 0 is dimensionless, so we can restore the values of fundamental constants by dimensional analysis: This interaction exerts a torque on nuclear spins, with the drive Rabi frequency given by where v ⊥ is the component of the velocity perpendicular to the direction of the leading field B 0 . As in the previous section, the spin driving field is "linearly polarized", and therefore the Rabi frequency contains an extra factor of 1/2, which arises because only one of the two counter-rotating components of the linearly polarized drive is resonant (rotating wave approximation).

VII. SPECTRAL PROPERTIES OF THE SPIN RESPONSE DUE TO AXION-LIKE DARK MATTER
In the first approximation we assume that the axion-like dark matter field is coherent, and drives the 207 Pb nuclear spins at carrier angular frequency ω a with Rabi frequency Ω a . The steady-state transverse spin magnetization that develops under the action of this driving field is given by Eq. (1) of the main text. The resulting voltage recorded by the ADC is: The time-averaged power in this signal is Note that we use the term "power" in the signal processing context, and this is proportional to the physical power. The Galactic axion-like dark matter halo field a(t) is not perfectly coherent. In this work we search for the axion-like dark matter halo that follows the standard halo model [30,31]. In this model the ALP speeds v in the Galactic frame follow the Maxwell-Boltzmann distribution where v 0 ≈ 220 km/s is the most probable speed [31]. The laboratory frame moves relative to the Galactic frame with the average speed v lab ≈ 232 km/s which has annual and daily modulations due to, respectively, Earth's revolution about the Sun and Earth's rotation around its axis [32]. The distribution of ALP speeds broadens the Fourier spectrum of the ALP field a(t), giving it a characteristic linewidth ≈ v 2 0 ν a /c 2 ≈ 10 −6 ν a . The power spectrum of the ALP field a(t) is given by the function where This spectral function is normalized so that This is the spectral lineshape used in searches for ALP-photon interactions [33][34][35][36].
FIG. S6. NMR calibration at the three values of the bias field B0. FID data are recorded after excitation pulses at Rabi frequency Ωe = 0.88 rad/ms and pulse length 20 ms. The excitation carrier frequency is plotted on the x-axis. Following the procedure used to obtain Fig. 2(b) in the main text, results are normalized so that the integral of the spectrum is unity.
The error bars show one standard deviation uncertainties of the FID spectrum fits, performed as described in section II C. Each spectrum is modeled as a super-Gaussian of order 2 (Eq. (S14)) and constant width 78 kHz (orange line). The only free parameter is the central frequency. The best-fit values of the central frequency for the three calibration data sets are: ν0 = (39159 ± 1) kHz, (39708 ± 1) kHz, (40160 ± 2) kHz.
(3) Search for narrow RF interference spectral lines using the Savitzky-Golay filter with order 2 and length 31 [37]. Spectral lines narrower than the ALP linewidth are distinguished by the difference between the filtered and raw power spectral densities. The points where this difference is above a threshold are marked as narrow spectral lines and are assigned the average value of their neighboring points.
(4) Optimally-filter data by convolving the power spectral density with the spectral lineshape for the ALP EDM interaction f 0 (ν) given in Eq. (S57). The separation between distinct ALP search frequencies is set to the ALP signal linewidth 3(v 2 0 /c 2 )ν 0 /4, where ν 0 is the central Larmor frequency, determined by the value of the bias field B 0 [32,37].
(5) Model the histogram of the optimally-filtered power spectral density with 100 bins as the Gaussian distribution with mean µ and standard deviation σ. Calculate the detection threshold at µ + 3.355σ, corresponding to a 5σ detection with 95% confidence level. Points above the threshold are ALP detection candidates. A detailed explanation of the choice of threshold value can be found in Refs. [36,37].
This analysis process is repeated for data taken at each of the 21 settings of bias magnetic field B 0 in the scan. The spin response to an axion-like dark matter signal will only appear in the data set where B 0 is such that the ALP Compton frequency is within the magnetic resonance excitation spectrum. For each data set we use the 80 kHz frequency band centered at Larmor frequency ν 0 , corresponding to the excitation spectrum, to search for the ALP signal, as described above. The rest of the spectral data within the 1 MHz scan range are used to reject residual background RF interference, which is not eliminated by the Savitzky-Golay filter. In addition, re-scan measurements are analyzed to eliminate statistical fluctuations, which are expected, given the large bandwidth of our search (lookelsewhere effect). The analysis procedure is as follows.
(a) At each value of bias magnetic field we consider ≈ 5000 frequency points (independent values of the ALP Compton frequency). For Gaussian-distributed data we expect two points to be above the 3.355σ threshold.
Typically we obtain ≈ 30 candidates above the threshold. The excess candidates are due to RF interference.
(b) We compare candidate frequencies from the "resonant" data set (for which the frequency is within the excitation spectrum) to the candidate frequencies from the "background" data sets (for which the frequency is outside the excitation spectrum). If the candidate frequency appears in one of the background data sets, it is rejected as RF interference. On average this eliminates ≈ 28 candidates at each value of B 0 .
(c) We compare candidate frequencies from the scan and re-scan data sets. If a candidate frequency appears only in one of those data sets, it is rejected as a statistical fluctuation. On average this eliminates ≈ 2 candidates at each value of B 0 .
This analysis procedure rejects all candidates above the 3.355σ threshold at all values of B 0 . We do not detect an axion-like dark matter signal. Therefore, for each value of B 0 , we quote the g d coupling value that corresponds to the 5σ value of the power spectral density as the 95% confidence interval limit [36]. We search for the gradient coupling g aNN of axion-like dark matter using the same steps as described above, with the standard halo model lineshape in step (4) replaced by the gradient coupling lineshape f 1 (ν), given in Eq. (S61). We calculate the angle ζ at each value of B 0 during the scan, based on the coordinates of our laboratory and the time at which the data are recorded, Fig. S7.
Our analysis for the gradient coupling g aNN rejects all candidates above the 3.355σ threshold at all values of B 0 . Therefore, for each value of B 0 , we quote the g aNN coupling value that corresponds to the 5σ value of the power spectral density as the 95% confidence interval limit. We note that the variation in ζ throughout the scan means that the shape of the limit curves for g d and for g aNN is slightly different in Fig. 4(b) of the main text, however this difference is smaller than the line thickness on the logarithmic plot.

A. Testing the data analysis procedure by injecting ALP signals
We test our data-analysis procedure by injecting into the experimental spectra synthetic axion-like dark matter signals with the lineshape given by Eq. (S57). Figure S8(a) shows the spectrum with an injected signal at Compton frequency ν a = 39.1586 MHz and with coupling strength g d = 1.4 × 10 −3 GeV −2 . After optimal filtering, the injected signal shows up as a candidate with amplitude 101 fT 2 , as shown in Fig. S8(b). The histogram of the optimally-filtered data points shows that this injected signal is detected at 20σ significance, Fig. S8(c). We test the recovery of the coupling strength by injecting 10 simulated signals, whose coupling strength is varied between g d = 7.0 × 10 −4 GeV −2 and g d = 7.0 × 10 −3 GeV −2 and whose Compton frequencies are selected randomly between ν a = 39.1185 MHz and ν a = 39.1985 MHz. The coupling strengths recovered from detected signals are shown in Fig. S8(d). We find that, on average, our analysis procedure results in a (2.7 ± 0.8)% suppression in the recovered coupling strength. This is due to the discrete sampling of the ALP search frequencies. If the injected ALP frequency falls between the search frequencies, there is a small mismatch in the lineshapes, which reduces the recovered coupling strength. The limits reported in the main text are corrected for this suppression.

B. Projected sensitivity reach
Our experimental results demonstrate the feasibility of using solid-state nuclear magnetic resonance to search for axion-like dark matter. There are several bounds on the relevant interactions of axion-like dark matter in this mass range, based on analysis of cooling dynamics of supernova SN1987A [28,41,42], and of Big-Bang nucleosynthesis [43]. However these model-dependent bounds are subject to significant caveats and uncertainties, and may be evaded altogether [44,45]. Stringent experimental limits on g d and g aNN exist at much lower ALP masses [29,[46][47][48][49][50], but the mass range probed in the current search has been, until now, experimentally unexplored.
The current sensitivity is not yet sufficient to reach the benchmark QCD axion level. The two main reasons are: (1) the CSA-induced inhomogeneous broadening of the NMR linewidth of the 207 Pb nuclear spin ensemble, and (2) the small size of our PMN-PT sample. We plan to circumvent the inhomogeneous broadening by concentrating our future searches on the lower Compton frequencies (ν a < 1 MHz), where the linewidth will be dominated by the T 2 spin coherence time, rather than CSA. The long T 1 relaxation time will allow us to pre-polarize the nuclear spins, retaining their polarization even at lower fields. We plan to use Superconducting Quantum Interference Devices (SQUIDs) to detect the transverse magnetization in this frequency range. The green dashed curves in Fig. S9 show the projected experimental sensitivity for the search with the same 4.6 mm sample as used in the current work. The cutoff at the low frequency end is set at the 1/T 2 NMR linewidth, and the cutoff at high frequencies is set by the Larmor frequency at the maximum magnetic field of 15 T.
In order to reach sufficient sensitivity to probe the QCD axion coupling strengths, we plan to scale up the volume of the ferroelectric sample. If the sample is coupled to the SQUID sensor with a broadband circuit, sample size of ≈ 80 cm and operation at ≈ 100 mK temperature are sufficient to reach the QCD axion line over ≈ 3 decades in mass, Fig. S9, blue dashed line. Implementing a resonant coupling circuit with a modest quality factor ≈ 1000 may allow us to reach this sensitivity level with a sample that is an order of magnitude smaller. The ultimate sensitivity limit is determined by the nuclear spin projection noise, Fig. S9, black dashed line. The region shaded in red is the exclusion at 95% confidence level placed by this work (CASPEr-e). The purple line shows the QCD axion coupling band. The darker purple color shows the mass range motivated by theory [9]. The blue regions mark the mass ranges where the ADMX and HAYSTAC experiments have probed the QCD axion-photon coupling [33,34]. The green region is excluded by analysis of cooling in supernova SN1987A, with color gradient indicating theoretical uncertainty [9]. The dashed green line marks the projected 5σ sensitivity of our CASPEr-e search with a 4.6 mm sample, as used in current work. The dashed blue line marks the projected 5σ sensitivity of our CASPEr-e search with an 80 cm sample, operating at 100 mK temperature. Implementing a resonant coupling circuit will enable operation with a smaller sample. The black dashed line marks the sensitivity limited by the quantum spin projection noise [28]. This is sufficient to detect the EDM coupling of the QCD axion across the 6-decade mass range from ≈ 0.3 peV to ≈ 500 neV. The other bounds are as follows. (a) The pink region is excluded by the neutron EDM (nEDM) experiment [47]. The blue region is excluded by the HfF + EDM experiment [50]. The yellow region is excluded by analysis of Big Bang nucleosynthesis (BBN) [43]. (b) The pink region is excluded by the neutron EDM (nEDM) experiment [47]. The blue region is excluded by the zero-to-ultralow field comagnetometer (ZULF CM) experiment [48]. The gray region is excluded by the zero-to-ultralow field sideband (ZULF SB) experiment [29]. The yellow region is excluded by the new-force search with K-3 He comagnetometer [46]. The bounds are shown as published, although corrections should be made to some of the low-mass limits, due to stochastic fluctuations of the axion-like dark matter field [51].