Search for axion-like dark matter using solid-state nuclear magnetic resonance
Deniz Aybas, Janos Adam, Emmy Blumenthal, Alexander V. Gramolin, Dorian Johnson, Annalies Kleyheeg, Samer Afach, John W. Blanchard, Gary P. Centers, Antoine Garcon, Martin Engler, Nataniel L. Figueroa, Marina Gil Sendra, Arne Wickenbrock, Matthew Lawson, Tao Wang, Teng Wu, Haosu Luo, Hamdi Mani, Philip Mauskopf, Peter W. Graham, Surjeet Rajendran, Derek F. Jackson Kimball, Dmitry Budker, Alexander O. Sushkov
aa r X i v : . [ h e p - e x ] J a n Search for axion-like dark matter using solid-state nuclear magneticresonance
Deniz Aybas,
1, 2
Janos Adam, Emmy Blumenthal, Alexander V. Gramolin, Dorian Johnson, Annalies Kleyheeg, Samer Afach,
3, 4
John W. Blanchard, Gary P. Centers,
3, 4
Antoine Garcon,
3, 4
Martin Engler,
3, 4
Nataniel L. Figueroa,
3, 4
Marina Gil Sendra,
3, 4
Arne Wickenbrock,
3, 4
Matthew Lawson,
5, 6
Tao Wang, Teng Wu, Haosu Luo, Hamdi Mani, Philip Mauskopf, Peter W. Graham, Surjeet Rajendran, Derek F. Jackson Kimball, Dmitry Budker,
3, 4, 14 and Alexander O. Sushkov
1, 2, 15, ∗ Department of Physics, Boston University, Boston, MA 02215, USA Department of Electrical and Computer Engineering,Boston University, Boston, MA 02215, USA Helmholtz-Institut, GSI Helmholtzzentrum für Schwerionenforschung, 55128 Mainz, Germany Johannes Gutenberg-Universität Mainz, 55128 Mainz, Germany The Oskar Klein Centre for Cosmoparticle Physics, Department of Physics,Stockholm University, AlbaNova, 10691 Stockholm, Sweden Nordita, KTH Royal Institute of Technology and Stockholm University,Roslagstullsbacken 23, 10691 Stockholm, Sweden Department of Physics, Princeton University, Princeton, New Jersey, 08544, USA State Key Laboratory of Advanced Optical Communication Systems and Networks, Department of Electronics,and Center for Quantum Information Technology, Peking University, Beijing 100871, China Shanghai Institute of Ceramics, Chinese Academy of Sciences, China School of Earth and Space Exploration, Arizona State University, Tempe, AZ 85287, USA Stanford Institute for Theoretical Physics,Stanford University, Stanford, California 94305, USA Department of Physics & Astronomy, The Johns Hopkins University, Baltimore, Maryland 21218, USA Department of Physics, California State University - East Bay, Hayward, California 94542-3084, USA Department of Physics, University of California, Berkeley, California 94720-7300, USA Photonics Center, Boston University, Boston, MA 02215, USA (Dated: January 6, 2021)We report the results of an experimental search for ultralight axion-like dark matter in the massrange 162 neV to 166 neV. The detection scheme of our Cosmic Axion Spin Precession Experiment(CASPEr) is based on a precision measurement of
Pb solid-state nuclear magnetic resonance ina polarized ferroelectric crystal. Axion-like dark matter can exert an oscillating torque on
Pbnuclear spins via the electric-dipole moment coupling g d , or via the gradient coupling g aNN . Wecalibrated the detector and characterized the excitation spectrum and relaxation parameters of thenuclear spin ensemble with pulsed magnetic resonance measurements in a 4.4 T magnetic field. Weswept the magnetic field near this value and searched for axion-like dark matter with Comptonfrequency within a 1 MHz band centered at 39.65 MHz. Our measurements place the upper bounds | g d | < . × − GeV − and | g aNN | < . × − GeV − (95% confidence level) in this frequencyrange. The constraint on g d corresponds to an upper bound of 7 . × − e · cm on the amplitudeof oscillations of the neutron electric dipole moment, and 3 . × − on the amplitude of oscillationsof CP-violating θ parameter of quantum chromodynamics. Our results demonstrate the feasibilityof using solid-state nuclear magnetic resonance to search for axion-like dark matter in the nano-electronvolt mass range. The existence of dark matter is indicated byastronomical and cosmological evidence, but itsinteractions, aside from gravity, remain unde-tected [1, 2]. A number of theoretical models ofphysics at high energies, such as string theory,grand unified theories, and models with extra di-mensions, incorporate light pseudoscalar bosons(axion-like particles, ALPs), which are potential dark matter candidates [3–7]. Among these, theaxion is particularly compelling, because it also of-fers a solution to the strong CP problem of quan-tum chromodynamics (QCD) [7–11]. The axionor axion-like field a ( t ) = a cos ( ω a t ) oscillates atthe Compton frequency ν a = ω a / (2 π ) = m a c /h ,where c is the speed of light in vacuum, h is thePlanck constant, and m a is the unknown ALP FIG. 1. Experimental setup. (a) The sample was a cylindrical ferroelectric PMN-PT crystal with diameter 0 .
46 cmand thickness 0 .
50 cm. It was electrically polarized along the cylinder axis, indicated with the black arrow. Thepickup coil and the cancellation coil were coaxial with the crystal, and the axis of the Helmholtz excitation coilwas orthogonal. The vertical leading magnetic field B set the direction of the equilibrium spin polarization.Coils were supported by G-10 fiberglass cylinders shown in gray and pink. (b) Electrical schematic, showingthe excitation and pickup circuits. Excitation pulses generated with the digital-to-analog converter (DAC) wereamplified (A e ), and coupled to the excitation coil via a tuned tank circuit that included matching and tuningcapacitors, as well as a resistor to set the circuit quality factor. The pickup probe was also designed as a tunedtank circuit, coupling the voltage induced in the pickup coil to a low-noise cryogenic amplifier (A ), whose outputwas filtered, further amplified, and digitized with an analog-to-digital converter (ADC). (c) Pulsed NMR sequenceused for FID measurements. The spin-ensemble equilibrium magnetization, initially parallel to B , was tilted intothe transverse plane by the excitation pulse. The FID signal was recorded after the excitation pulse, as themagnetization precessed and its transverse component decayed. mass, which can be in a broad range, roughly be-tween 10 − eV and 10 − eV [12–14]. The fieldamplitude a is fixed by the assumption thatit dominates the dark matter energy density: ρ DM = m a a / ≈ . × − GeV [15, 16]. Ki-netic energy of the axion-like dark matter field in-troduces small corrections to its frequency spec-trum. The standard halo model predicts thespectral shape with linewidth ( v /c ) ν a ≈ − ν a ,where v ≈
220 km / s is the circular rotation speedof the Milky Way galaxy at the Sun’s location [17,18].Experimental searches for axion-like particlesrely on symmetry arguments about the natureof their interactions with Standard Model parti-cles [7, 16, 19, 20]. These interactions are sup-pressed by a large energy scale, set by the decayconstant f a , which could lie near the grand uni-fication, or the Planck scale [21]. Most exper-iments to date have focused on the electromag-netic interaction, which can mix photons with ax-ions and ALPs in the presence of a strong mag-netic field [22–32]. The Cosmic Axion Spin Pre-cession Experiments (CASPEr) search for differentinteractions: the electric dipole moment (EDM)interaction and the gradient interaction with nu- clear spin I [19, 33–37]. The gradient interac-tion Hamiltonian is H aNN = g aNN ∇ a · I , where g aNN is the coupling strength. The EDM inter-action arises from the defining coupling of the ax-ion to the gluon field [38]. Its Hamiltonian canbe written as H EDM = g d a E ∗ · I /I , where g d isthe coupling strength and E ∗ is an effective elec-tric field [19]. This interaction is equivalent tothat of a parity- and time-reversal-violating os-cillating EDM, given by d = g d a cos ( ω a t ). Thiscorresponds to an oscillating QCD θ parame-ter: θ ( t ) = ( a /f a ) cos ( ω a t ), with g d inversely pro-portional to f a [16, 39]. The EDM couplinggenerates axion mass, and for the QCD axion m a ≈ Λ QCD /f a , where Λ QCD ≈
200 MeV is theQCD confinement scale [16, 40].The sensitivity of static EDM experiments tothe oscillating EDM is suppressed, although datare-analysis has produced limits at low frequen-cies [41, 42]. Astrophysical constraints can be de-rived by analyzing the cooling dynamics of the su-pernova SN1987A [16, 43]. Constraints can also beextracted from analysis of He production duringBig Bang nucleosynthesis [44] and from analysis ofblack hole superradiance [45]. CASPEr-electric isa direct, model-independent search for the EDMand gradient interactions of axion-like dark mat-ter, with the potential to reach the sensitivity tothe QCD axion [19]. We search for the effects ofthese interactions on the dynamics of a spin ensem-ble in a solid with broken inversion symmetry [46–52]. The measurements focus on Pb ions,with nuclear spin I = 1 /
2, in a poled ferroelec-tric PMN-PT crystal with the chemical formula:(PbMg / Nb / O ) / − (PbTiO ) / [53]. Thenon-centrosymmetric position of the ions in thiscrystal gives rise to a large effective electric field,analogous to the effect in polar molecules [54–56].The EDM or gradient interaction with axion-likedark matter creates an oscillating torque on thenuclear spins. We quantify the magnitude of thistorque by the Rabi frequency Ω a , which is propor-tional to the corresponding interaction strength.For a spin ensemble polarized by an external biasmagnetic field, this torque tilts the spins, if it isresonant with their Larmor frequency. The ex-perimental observable is the oscillating transversemagnetization: M a = uM Ω a T cos ( ω a t ) , (1)where M is the equilibrium magnetization of the Pb nuclear spin ensemble, T is the nuclear spincoherence time, and u is a dimensionless spectralfactor that takes into account the inhomogeneousbroadening of the spin ensemble and the detuningbetween the ALP Compton frequency and the spinLarmor frequency [53].Our apparatus makes use of inductive detectionto measure the Pb spin precession, Fig. 1(a).We poled the cylindrical PMN-PT crystal alongits axis, aligned with the [1,1,1] crystal direction.This created the axial effective electric field E ∗ ,proportional to the remanent polarization P r . Wemounted the crystal inside a fiberglass tube, sothat E ∗ was perpendicular to the vertical biasmagnetic field B , created with a superconduct-ing solenoid. A pickup coil, wound around thetube, was coupled to a low-noise cryogenic pream-plifier with a tuned matching circuit, Fig. 1(b).We tuned the pickup probe to have its resonanceat 39.7 MHz with quality factor 26, and matchedits impedance to the 50 Ω input impedance of thepreamplifier [53]. A cylindrical copper shield at-tenuated external sources of RF interference. Weperformed all experiments with the apparatus sub-merged in a liquid helium bath at 4 . Pbpulsed nuclear magnetic resonance (NMR) mea-surements, Fig. 1(c). The spins were excited byresonant magnetic field pulses, created by deliv-ering current to the 2 × α = V / ( µ M ), where V is the recorded voltagereferred to the amplifier input, M is the trans-verse sample magnetization, and µ is the perme-ability of free space. Despite our efforts to min-imize the inductive and capacitive couplings be-tween the excitation and the pickup coils, we foundthat the cryogenic preamplifier saturated duringexcitation pulses, and its recovery time was toolong to observe the fast FID decay [53]. To addressthis problem, we placed a single-turn cancellationcoil near the pickup coil, Fig. 1(a), and deliveredto it a compensating current during the excitationpulses. The amplitude and phase of this compen-sating current were chosen to cancel the currentin the pickup probe during excitation, and preventpreamplifier saturation, without affecting spin ex-citation. This scheme is a substitute for the trans-mit/receive switch, often used in NMR detectors.We performed the NMR calibration measure-ments at the leading magnetic field B = 4 . M of the spin ensemble was µ M = 2 . . ± .
2) nT by saturating the spins, thenletting magnetization recover over approximatelyone population relaxation time [53]. We set theexcitation carrier frequency to 39 .
71 MHz, andrecorded the FID signals after excitation pulsesof variable width. The Fourier spectrum of oneof these FID signals is shown in Fig. 2(a). Wemodeled the FID lineshapes by numerically solv-ing the Bloch equations for a spin ensemble withan inhomogeneously-broadened excitation spec-trum [53]. By fitting the data, we extracted thetransverse coherence time of the nuclear spins: T = (16 . ± .
9) ms, and the pickup-circuit trans-fer coefficient α = (2 . ± . × V / T. We notethat there is a sharp central feature with linewidthon the order of the Rabi frequency, but the over-all FID spectral width is much greater than 1 /T , FIG. 2. Sensitivity calibration. (a) Measurements of
Pb FID following a spin excitation pulse of length t p = 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 Fouriertransform of the detected voltage, referred to the input of the pickup probe amplifier A . Data points were binnedand averaged, the error bars show one standard deviation for each bin. The lines show the best-fit simulation ofthe spin response, with the light-colored narrow bands indicating the range of simulation results if parametersare varied by one standard deviation away from their best-fit values. We performed the fitting simultaneouslyto three FID data sets, with excitation pulse lengths t p = 0 . , ,
20 ms, with free parameters includingthe spin coherence time T and pickup circuit transfer coefficient α [53]. (b) Measurement of the normalized Pb NMR excitation spectrum near Larmor frequency 39 .
71 MHz. Excitation pulses of length 1 . e = 0 .
88 rad / ms were delivered at the carrier frequencies shown on the x-axis. Data points show theamplitude of the spin FID response, normalized so that the integral of the spectrum is unity. The error barsindicate one standard deviation uncertainties of the FID spectrum fits. We model the excitation spectrum asa super-Gaussian of order 2 (red line) [53]. (c) Detector calibration for varying drive Rabi frequency. Datapoints show the amplitude of the spin FID response after an excitation pulse of length 20 ms, delivered at thecarrier frequency 39 .
71 MHz, with Rabi frequency Ω e plotted on the x-axis. The error bars indicate one standarddeviation 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 Blochequations with parameters extracted from data in panel (a). (d) Measurement of ferroelectric hysteresis in thePMN-PT single crystal. The remanent polarization P r persists after the applied voltage has been ramped downto zero. since the tilting pulse excites a broad frequencyband within the inhomogeneous spin distribution.The exact shape of the FID Fourier spectrum de-pends on the interplay between the excitation-pulse spectrum, the distribution of tipping anglesacross the spin ensemble, and the T coherencetime.We measured the inhomogeneous broaden-ing of the Pb nuclear spins in the sam- ple by sweeping the excitation pulse carrier fre-quency and recording the corresponding FIDspectra. The resulting NMR excitation spec-trum was centered at 39 .
71 MHz and had a fullwidth Γ / (2 π ) = (78 ±
2) kHz, Fig 2(b). Thisbroadening is consistent with the chemical shiftanisotropy (CSA) of
Pb observed in solid-stateNMR [57]. We measured the population relax-ation time T of the Pb nuclear spin ensemblewith a saturation-recovery measurement, obtain-ing T = (25 . ± .
6) min [53].The spin evolution in our pulsed NMR calibra-tion measurements was more complicated than theCW-like small spin-tip angle response to axion-like dark matter, described by Eq. (1). In or-der to confirm the validity of our NMR modelin the limit of small spin-tip angles, we recordedand analyzed FID data for a range of excitationRabi frequencies Ω e . For these measurements wekept the excitation pulse width at 20 ms – approx-imately the coherence time of axion-like dark mat-ter field with Compton frequency near 40 MHz.At small excitation amplitudes, the spin responsewas linear in Ω e , as described by Eq. (1) for thecase of the drive due to interaction with axion-like dark matter, Fig 2(c). The slope of the lin-ear response is proportional to the spectral fac-tor u = (3 . ± . × − , which is well approxi-mated by the ratio of the homogeneous linewidth π/T and the inhomogeneously-broadened excita-tion spectrum width Γ [53]. The deviation fromlinearity at larger Ω e is due to saturation of theresonant spins in the excitation spectrum, consis-tent with our Bloch-equation simulations.Prior to any measurements, the PMN-PT crys-tal was ferroelectrically poled at room tempera-ture by applying 3 . P r = (22 ± µ C / cm . We recorded hysteresisdata before and after the experiments searchingfor axion-like dark matter, and verified that thefractional degradation of polarization due to ther-mal cycling and fatigue was smaller than thequoted uncertainty. The effective electric field E ∗ is proportional to the ferroelectric polariza-tion [48, 54, 55]. In order to calculate the valueof E ∗ we considered the Schiff moment S of the Pb nucleus, induced by the oscillating QCD θ parameter [58, 59]. The dominant contribution tothe Schiff moment arises from the parity- and time-reversal-violating nuclear forces, resulting in thevalue S = 0 . θ e · fm [53, 60–64]. This corre-sponds to the magnitude of effective electric field E ∗ = 340 kV / cm. We estimate the theoretical un-certainty in E ∗ on the level of 50% [53].In order to search for axion-like dark matterwe swept the leading magnetic field B in 21 steps, corresponding to the search frequency range39 . . Pb nuclear spin excitation spec-trum, Fig. 2(b). The broad NMR excitation spec-trum reduced the necessary number of magneticfield steps for a given search frequency range. Ateach value of B we recorded 58 s of scan data sen-sitive to axion-like dark matter, followed by 58 sof re-scan data that were used in our analysis toidentify statistical fluctuations. In order to confirmthe experimental calibration, we performed pulsedNMR measurements at three values of the lead-ing field, corresponding to the extremes and themidpoint of the search frequency range [53].Data analysis consisted of several processing,correction, and signal-search steps. At each valueof the leading field B we divided the recorded scandata into 27 blocks, each of 2 .
15 s duration, chosento be much longer than the ≈
25 ms coherence timeof any potential ALP dark matter signal in our fre-quency range. We used the pickup-circuit transfercoefficient α to convert the recorded voltage valuesto magnetization, and performed a discrete Fouriertransform on each block, subsequently averagingthe power spectral densities (PSDs) of the blocks.Many of the spectra were contaminated with nar-rowband RF interference that penetrated our elec-tromagnetic shielding. We used Savitzky-Golaydigital filtering to identify and reject these narrow-band features, while preserving potential axion-likedark matter signals, whose spectral shape is pre-dicted by the standard halo model [25, 53, 65].We then processed the data to search for sig-nals due to the EDM and the gradient interactions.The first step was optimal filtering, performed byconvolving the PSD with the signal lineshape pre-dicted for the corresponding interaction [53]. Ateach value of B we retained the optimally-filtereddata points in a frequency bin, centered at thecorresponding Larmor frequency, with full width80 kHz, covering the excitation spectrum band-width. We modeled the histogram of these datapoints as the normal distribution with standarddeviation σ , inset of Fig. 3. We set the candi-date detection threshold to 3 . σ , equivalent to95% confidence interval for a 5 σ detection, andflagged all points above the threshold as candi-dates [32, 53, 65].There were 617 candidates for EDM coupling(636 for gradient coupling). In order to rejectresidual RF interference, we used the fact that RF FIG. 3. Results of the search for spin interactions with axion-like dark matter. The line shows the limits on theaxion-like dark matter EDM coupling (left y-axis) and nucleon gradient coupling (right y-axis) in the mass range162 neV −
166 neV. The region above the line is excluded at 95% confidence level. The inset shows the histogramof the optimally-filtered power spectral density of transverse sample magnetization within the frequency windowcentered at 39 .
16 MHz. The red line shows the Gaussian distribution model, and the vertical black dashed lineshows the 3 . σ candidate threshold at 17 fT . pickup is independent of the leading field B , whilean axion-like dark matter signal should only ap-pear when B is tuned to a value such that the spinexcitation spectrum overlaps with the ALP Comp-ton frequency. We compared the candidates fromdata sets taken at different values of B , rejecting569 candidates for EDM coupling (577 for gradientcoupling). The remaining 48 candidates for EDMcoupling (59 for gradient coupling) were shownto be statistical fluctuations, using a scan/re-scananalysis [53]. The search sensitivity was limited bythe ≈ .
05 nV / √ Hz input noise level of the ampli-fier, corresponding to a magnetic field sensitivityof ≈ / √ Hz.Our search did not yield a discovery of the EDMcoupling g d or the gradient coupling g aNN of axion-like dark matter. In the absence of a detection,in each frequency bin the 95% confidence intervallimit on magnitudes of these coupling constantscorresponds to the 5 σ value in the Gaussian distri-bution of the optimally-filtered PSD [32, 53, 65].The limits were normalized by the NMR excita-tion spectrum for each bin and concatenated toproduce constraints on g d and g aNN over the en-tire frequency search range, Fig. 3. Over the fre-quency range 39 . . | g d | is 7 . × − GeV − , corresponding to anupper bound of 7 . × − e · cm on the ampli-tude of oscillations of the neutron electric dipole moment, and 3 . × − on the amplitude of os-cillations of the QCD θ parameter. The con-straint on | g aNN | is 2 . × − GeV − . We arenot aware of any existing experimental limits onthese interactions in this ALP mass range. Anal-ysis of cooling dynamics of supernova SN1987Acan be used to estimate bounds g d . − GeV − and g aNN . − GeV − [19, 24, 43]. Howeverthese model-dependent bounds are subject to sig-nificant caveats and uncertainties, and may beevaded altogether, reinforcing the importance oflaboratory searches [66, 67]. Stringent experimen-tal limits on g d and g aNN exist at much lower ALPmasses [35, 36, 41, 42, 68–71].There are several ways to improve experimen-tal sensitivity to axion-like dark matter. Sincethe CSA-induced inhomogeneous broadening isproportional to the Larmor frequency, search-ing in a lower ALP mass range will reduce thelinewidth and therefore improve the search sen-sitivity. A search in the lower mass range willlikely also benefit from superconducting detectors,such as SQUIDs and quantum upconverters [72].Manipulation of light-induced transient paramag-netic centers may enable control over the nuclearspin population-relaxation time T , and nuclearspin hyperpolarization using dynamic polarizationtechniques. A dramatic sensitivity improvementcould be achieved by scaling up the sample volume.We estimate that with a sample size of ≈
80 cm, itmay be possible to reach the sensitivity necessaryto detect the QCD axion g d coupling strength inthe mass range between ≈ peV and ≈ ∗ [email protected][1] D. N. Spergel, Science , 1100 (2015).[2] G. Bertone and T. M. P. Tait,Nature , 51 (2018).[3] J. Preskill, M. B. Wise, and F. Wilczek, PhysicsLetters B , 127 (1983).[4] L. Abbott and P. Sikivie, Physics Letters B ,133 (1983).[5] M. Dine and W. Fischler, Physics Letters B ,137 (1983).[6] P. Svrcek and E. Witten, Journal of High EnergyPhysics , 051 (2006).[7] I. G. Irastorza and J. Redondo, Progress in Parti-cle and Nuclear Physics , 89 (2018).[8] R. D. Peccei and H. R. Quinn,Physical Review Letters , 1440 (1977).[9] S. Weinberg, Physical Review Letters , 223 (1978).[10] F. Wilczek, Physical Review Letters , 279 (1978).[11] D. DeMille, J. M. Doyle, and A. O. Sushkov,Science , 990 (2017).[12] P. W. Graham and A. Scherlis,Physical Review D , 035017 (2018).[13] A. Ernst, A. Ringwald, and C. Tamarit,Journal of High Energy Physics , 103 (2018).[14] K. Schutz, Physical Review D , 123026 (2020).[15] M. Tanabashi et al. (Particle Data Group), Phys.Rev. D , 030001 (2018).[16] P. W. Graham and S. Rajendran, Physical Review D , 035023 (2013).[17] M. S. Turner, Physical Review D , 3572 (1990).[18] N. W. Evans, C. A. O’Hare, and C. McCabe,Physical Review D , 023012 (2019).[19] D. Budker, P. W. Graham, M. Ledbetter, S. Ra-jendran, and A. O. Sushkov, Phys. Rev. X ,021030 (2014).[20] A. Arvanitaki and A. A. Geraci,Physical Review Letters , 161801 (2014).[21] P. W. Graham and A. Scherlis, (2018),arXiv:1805.07362.[22] P. Sikivie, Physical Review Letters , 1415(1983).[23] N. Du, N. Force, R. Khatiwada, E. Lentz, R. Ot-tens, L. J. Rosenberg, G. Rybka, G. Carosi,N. Woollett, D. Bowring, A. S. Chou, A. Son-nenschein, W. Wester, C. Boutan, N. S. Oblath,R. Bradley, E. J. Daw, A. V. Dixit, J. Clarke,S. R. O’Kelley, N. Crisosto, J. R. Gleason, S. Jois,P. Sikivie, I. Stern, N. S. Sullivan, D. B. Tanner,and G. C. Hilton, Physical Review Letters ,151301 (2018).[24] P. W. Graham, I. G. Irastorza, S. K. Lamoreaux,A. Lindner, and K. A. van Bibber, Annual Reviewof Nuclear and Particle Science , 485 (2015).[25] B. M. Brubaker, L. Zhong, Y. V. Gurevich, S. B.Cahn, S. K. Lamoreaux, M. Simanovskaia, J. R.Root, S. M. Lewis, S. Al Kenany, K. M. Backes,I. Urdinaran, N. M. Rapidis, T. M. Shokair, K. A.van Bibber, D. A. Palken, M. Malnou, W. F.Kindel, M. A. Anil, K. W. Lehnert, and G. Carosi,Physical Review Letters , 061302 (2017).[26] J. Choi, H. Themann, M. J. Lee, B. R. Ko, andY. K. Semertzidis, Physical Review D , 061102(2017).[27] P. Sikivie, N. Sullivan, and D. B. Tanner, PhysicalReview Letters , 131301 (2014).[28] S. Chaudhuri, P. W. Graham, K. Irwin, J. Mar-don, S. Rajendran, and Y. Zhao, Physical ReviewD , 075012 (2015).[29] Y. Kahn, B. R. Safdi, and J. Thaler, PhysicalReview Letters , 141801 (2016).[30] S. Chaudhuri, K. Irwin, P. W. Graham, andJ. Mardon, aXiV:1803.01627 (2018).[31] J. L. Ouellet, C. P. Salemi, J. W. Foster, R. Hen-ning, Z. Bogorad, J. M. Conrad, J. A. Formag-gio, Y. Kahn, J. Minervini, A. Radovinsky, N. L.Rodd, B. R. Safdi, J. Thaler, D. Winklehner, andL. Winslow, Physical Review Letters , 121802(2019).[32] A. V. Gramolin, D. Aybas, D. John-son, J. Adam, and A. O. Sushkov,Nature Physics (2020), 10.1038/s41567-020-1006-6.[33] A. Garcon, D. Aybas, J. W. Blanchard,G. Centers, N. L. Figueroa, P. W. Gra-ham, D. F. J. Kimball, S. Rajendran, M. G.Sendra, A. O. Sushkov, L. Trahms, T. Wang,A. Wickenbrock, T. Wu, and D. Budker,Quantum Science and Technology , 014008 (2018). [34] T. Wang, D. F. J. Kimball, A. O. Sushkov, D. Ay-bas, J. W. Blanchard, G. Centers, S. R. O. Kel-ley, A. Wickenbrock, J. Fang, and D. Budker,Physics of the Dark Universe , 27 (2018).[35] T. Wu, J. W. Blanchard, G. P. Centers, N. L.Figueroa, A. Garcon, P. W. Graham, D. F. J.Kimball, S. Rajendran, Y. V. Stadnik, A. O.Sushkov, A. Wickenbrock, and D. Budker,Physical Review Letters , 191302 (2019).[36] A. Garcon, J. W. Blanchard, G. P. Centers,N. L. Figueroa, P. W. Graham, D. F. J. Kim-ball, S. Rajendran, A. O. Sushkov, Y. V. Stad-nik, A. Wickenbrock, T. Wu, and D. Budker,Science Advances , eaax4539 (2019).[37] D. F. Jackson Kimball, S. Afach, D. Ay-bas, J. W. Blanchard, D. Budker, G. Cen-ters, M. Engler, N. L. Figueroa, A. Gar-con, P. W. Graham, H. Luo, S. Rajendran,M. G. Sendra, A. O. Sushkov, T. Wang,A. Wickenbrock, A. Wilzewski, and T. Wu,in Springer Proceedings in Physics , Vol. 245(Springer, 2020) pp. 105–121, arXiv:1711.08999.[38] P. Graham and S. Rajendran,Physical Review D , 055013 (2011).[39] M. Pospelov and A. Ritz,arXiv:hep-ph/9908508 , 1 (1999).[40] M. Baldicchi, A. V. Nesterenko, G. M.Prosperi, D. V. Shirkov, and C. Simolo,Physical Review Letters , 242001 (2007).[41] C. Abel, N. J. Ayres, G. Ban, G. Bison, K. Bodek,V. Bondar, M. Daum, M. Fairbairn, V. V. Flam-baum, P. Geltenbort, K. Green, W. C. Griffith,M. van der Grinten, Z. D. Grujić, P. G. Harris,N. Hild, P. Iaydjiev, S. N. Ivanov, M. Kasprzak,Y. Kermaidic, K. Kirch, H.-C. Koch, S. Kom-posch, P. A. Koss, A. Kozela, J. Krempel,B. Lauss, T. Lefort, Y. Lemière, D. J. E. Marsh,P. Mohanmurthy, A. Mtchedlishvili, M. Musgrave,F. M. Piegsa, G. Pignol, M. Rawlik, D. Rebreyend,D. Ries, S. Roccia, D. Rozp¸edzik, P. Schmidt-Wellenburg, N. Severijns, D. Shiers, Y. V. Stadnik,A. Weis, E. Wursten, J. Zejma, and G. Zsigmond,Physical Review X , 041034 (2017).[42] T. S. Roussy, D. A. Palken, W. B. Cairncross,B. M. Brubaker, D. N. Gresh, M. Grau, K. C. Cos-sel, K. B. Ng, Y. Shagam, Y. Zhou, V. V. Flam-baum, K. W. Lehnert, J. Ye, and E. A. Cornell,arXiv: 2006.15787 (2020).[43] G. G. Raffelt, in Axions: Theory, Cosmology, and Experimental Searches (Springer Berlin Heidelberg, Berlin, Heidelberg,2008) pp. 51–71.[44] K. Blum, R. T. D’Agnolo, M. Lisanti, and B. R.Safdi, Physics Letters B , 30 (2014).[45] A. Arvanitaki, S. Dimopoulos, S. Dubovsky,N. Kaloper, and J. March-Russell,Physical Review D , 123530 (2010).[46] A. J. Leggett, Physical Review Letters , 586 (1978).[47] W. Bialek, J. Moody, and F. Wilczek, Physical Review Letters , 1623 (1986).[48] T. N. Mukhamedjanov and O. P. Sushkov,Physical Review A , 34501 (2005),arXiv:0411226 [physics].[49] D. Budker, S. Lamoreaux,A. Sushkov, and O. Sushkov,Physical Review A , 022107 (2006).[50] A. O. Sushkov, S. Eckel, and S. K. Lamoreaux,Physical Review A , 022104 (2010).[51] K. Z. Rushchanskii, S. Kamba, V. Goian,P. Vanek, M. Savinov, J. Prokleska, D. Nuzhnyy,K. Knízek, F. Laufek, S. Eckel, S. K. Lamoreaux,A. O. Sushkov, M. Lezaić, and N. A. Spaldin,Nature materials , 649 (2010).[52] S. Eckel, A. Sushkov, and S. Lamoreaux,Physical Review Letters , 193003 (2012).[53] Supplementary Information .[54] J. A. Ludlow and O. P. Sushkov,Journal of Physics B: Atomic, Molecular and Optical Physics , 085001 (2013).[55] L. V. Skripnikov and A. V. Titov,Journal of Chemical Physics , 054115 (2016).[56] V. Andreev, D. G. Ang, D. DeMille, J. M.Doyle, G. Gabrielse, J. Haefner, N. R. Hutzler,Z. Lasner, C. Meisenhelder, B. R. O’Leary, C. D.Panda, A. D. West, E. P. West, and X. Wu,Nature , 355 (2018), arXiv:0901.2328.[57] L.-S. Bouchard, A. O. Sushkov,D. Budker, J. Ford, and A. Lipton,Physical Review A , 022102 (2008).[58] L. Schiff, Physical Review , 2194 (1963).[59] P. G. Sandars, Physical Review Letters , 1396 (1967).[60] O. P. Sushkov, V. V. Flambaum, and I. B.Khriplovich, Sov. Phys. JETP , 873 (1984).[61] V. V. Flambaum, I. B.Khriplovich, and O. P. Sushkov,Nuclear Physics, Section A , 750 (1986).[62] I. B. Khriplovich and S. K. Lamoreaux, CP Violation Without Strangeness (SpringerBerlin Heidelberg, Berlin, Heidelberg, 1997).[63] V. V. Flambaum and V. A. Dzuba,Physical Review A , 042504 (2020).[64] K. Yanase and N. Shimizu,arXiv: 2006.15142 (2020).[65] B. M. Brubaker, L. Zhong, S. K. Lamoreaux,K. W. Lehnert, and K. A. van Bibber, PhysicalReview D , 123008 (2017).[66] W. DeRocco, P. W. Graham, and S. Rajendran,Physical Review D , 075015 (2020).[67] N. Bar, K. Blum, and G. D’amico,Physical Review D , 123025 (2020).[68] G. Vasilakis, J. M. Brown, T. W.Kornack, and M. V. Romalis,Physical Review Letters , 261801 (2009).[69] E. G. Adelberger and W. A. Terrano,Physical Review Letters , 169001 (2019).[70] T. Wu, J. W. Blanchard, G. P. Centers, N. L.Figueroa, A. Garcon, P. W. Graham, D. F.Kimball, S. Rajendran, Y. V. Stadnik, A. O.Sushkov, A. Wickenbrock, and D. Budker, Physical Review Letters , 169002 (2019).[71] W. A. Terrano, E. G. Adelberger,C. A. Hagedorn, and B. R. Heckel,Physical Review Letters , 231301 (2019). [72] S. Chaudhuri,
The Dark Matter Radio: Aquantum-enhanced search for QCD axion darkmatter , Ph.D. thesis (2019). r X i v : . [ h e p - e x ] J a n Supplemental Material forSearch for axion-like dark matter using solid-state nuclear magnetic resonance
Deniz Aybas,
1, 2
Janos Adam, Emmy Blumenthal, Alexander V. Gramolin, Dorian Johnson, Annalies Kleyheeg, Samer Afach,
3, 4
John W. Blanchard, Gary P. Centers,
3, 4
Antoine Garcon,
3, 4
Martin Engler,
3, 4
Nataniel L. Figueroa,
3, 4
Marina Gil Sendra,
3, 4
Arne Wickenbrock,
3, 4
Matthew Lawson,
5, 6
Tao Wang, Teng Wu, Haosu Luo, Hamdi Mani, Philip Mauskopf, Peter W. Graham, Surjeet Rajendran, Derek F. Jackson Kimball, Dmitry Budker,
3, 4, 14 and Alexander O. Sushkov
1, 2, 15, ∗ Department of Physics, Boston University, Boston, MA 02215, USA Department of Electrical and Computer Engineering, Boston University, Boston, MA 02215, USA Helmholtz-Institut, GSI Helmholtzzentrum für Schwerionenforschung, 55128 Mainz, Germany Johannes Gutenberg-Universität Mainz, 55128 Mainz, Germany The Oskar Klein Centre for Cosmoparticle Physics, Department of Physics,Stockholm University, AlbaNova, 10691 Stockholm, Sweden Nordita, KTH Royal Institute of Technology and Stockholm University, Roslagstullsbacken 23, 10691 Stockholm, Sweden Department of Physics, Princeton University, Princeton, New Jersey, 08544, USA State Key Laboratory of Advanced Optical Communication Systems and Networks, Department of Electronics,and Center for Quantum Information Technology, Peking University, Beijing 100871, China Shanghai Institute of Ceramics, Chinese Academy of Sciences, China School of Earth and Space Exploration, Arizona State University, Tempe, AZ 85287, USA Stanford Institute for Theoretical Physics, Stanford University, Stanford, California 94305, USA Department of Physics & Astronomy, The Johns Hopkins University, Baltimore, Maryland 21218, USA Department of Physics, California State University - East Bay, Hayward, California 94542-3084, USA Department of Physics, University of California, Berkeley, California 94720-7300, USA Photonics Center, Boston University, Boston, MA 02215, USA (Dated: January 6, 2021)
I. EXPERIMENTAL SETUPA. Description of the apparatus
Our cryogenic nuclear magnetic resonance (NMR) setup is inside a liquid helium (LHe) bath cryostat with a solenoidalsuperconducting magnet (Cryomagnetics, Inc. Model 90-300-010), Fig. S1. The apparatus is built around a crystalthat 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). Theexperimental 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 fieldexerts a torque on the spins, whose magnitude is quantified by the excitation Rabi frequency Ω e . The excitation-probe transfer function is defined as κ = Ω e V e . (S1)The excitation pulse tilts Pb nuclear spins into the plane transverse to the leading field B , creating a crystalmagnetization M 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 atthe input of the low-noise preamplifier A (Fig. S2). The pickup probe transfer function is defined as α = V µ M , (S2) ∗ [email protected] FIG. S1. Schematic of the setup operating at 4.2 K. Cylindrical PMN-PT crystal is placed close to the center of the supercon-ducting magnet, and it is coupled to the mutually orthogonal excitation and pickup coils. The sample and the pickup probeare mounted inside a cylindrical electromagnetic shielding enclosure within the superconducting magnet bore. The low-noisepreamplifier is inside the liquid helium bath, above the superconducting magnet. where µ is the permeability of free space. The preamplifier A has gain of 40. Its output is connected to a low-passfilter LP and another amplifier stage A (gain = 15) mounted inside the cryostat near the top flange. After a thirdamplifier stage A (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 (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 couplinginto 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 inputof A e is connected to ground via a 50 Ω termination. FIG. S2. Full electrical schematic of the experimental setup incorporating the Spectrum Instrumentation m4i.6622-x8 digital-to-analog converter (DAC), RF Lambda coaxial reflective SP2T RFSP2TRDC06G switches ( S ), Stanford Research SystemsSIM954 inverting amplifiers ( A e ), coil inductances ( L c , L e and L p ), surface mount probe tuning capacitors ( C c , C e and C p ),surface mount impedance matching capacitors ( C c , C e and C p ), surface mount resistors ( R c , R e and R p ) that determinethe quality factor of the circuit resonances, ULF-LNA- A ), Mini-Circuits ZX75LP-50-S+ 50 MHz low-pass filter ( LP ), Mini-Circuits ZX60-P103LN+low-noise amplifier ( A ), Femto HVA-200M-40-B amplifier ( A ), and Spectrum Instrumentation m4i.4421-x8 analog-to-digitalconverter (ADC). B. The crosstalk minimization scheme
During experimental assembly we carefully adjust the orthogonal axes of the excitation and the pickup coils tominimize mutual inductance between them, to ≈ .
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 residualinductive and capacitive couplings between the coils. This crosstalk signal saturates the preamplifier, resulting in arecovery time of ≈ µ s, which complicates the detection of the FID signal. In order to prevent saturation, duringthe excitation pulse we apply a waveform to the cancellation coil that is optimized to compensate the crosstalk currentin the pickup probe with minimal effect on spin dynamics. The phase and amplitude of this waveform are optimizedby monitoring the current measured at the pickup probe and minimizing its magnitude. This is done at zero leadingmagnetic field to avoid spin excitation during optimization. We emphasize that only a small ( < . . 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, coilinductance L is in parallel with a tuning capacitor C and a resistor R , and this tank circuit is in turn in series witha matching capacitor C . The total probe impedance is Z = 1 iωL + R + iωC + 1 iωC = (cid:18) ( ωL ) RR (1 − ω LC ) + ( ωL ) (cid:19) + i (cid:18) ωLR (1 − ω LC ) R (1 − ω LC ) + ( ωL ) − ωC (cid:19) . (S3)In order to match the probe impedance to Z = R = 50 Ω at the resonance angular frequency ω r , we have to choosethe following values for the circuit elements: R = Qω r L,C = 1 ω r L − Q r R − R R ! , (S4) C = 1 ω r s R ( R − R ) , where Q is the resonance quality factor. We used fixed-value surface mount capacitors and resistors, so the probesare 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 pathcancellation. It has a radius r p = 3 . L p = 0 . µ H, and is tuned to the resonant frequency ω p / (2 π ) = 39 .
71 MHz with quality factor Q p = 26 at 4 . Q p to 26. The excitation coil has a Helmholtz geometry with N e = 2 × r e = 7 . L e = 0 . µ H, which is tuned to the resonant frequency ω e / (2 π ) = 42 .
01 MHzwith quality factor Q e = 1 . . r c = 4 . L c = 0 . µ H, which is tuned to resonant frequency ω c / (2 π ) = 40 .
31 MHz withquality factor Q c = 2 at 4 . 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 κ , definedin Eqs. (S1,S2). Using Faraday’s law we can estimate the voltage induced in the pickup coil by an oscillating transversemagnetization M : V p ≈
13 ( γB ) N p ( µ M )( πr s ) , (S5)where B = 4 . r s = 2 . γ is the Pbgyromagnetic ratio. We set the demagnetizing factor to 1 /
3, as for a sphere, as an approximation for a cylindricalsample with height ≈ diameter. For the pickup probe on resonance with the spin Larmor frequency, the resultingvoltage at the input of preamplifier A is calculated from circuit analysis [1]: V ≈ V p s Q p R ω p L p , (S6)while Q p ω p L p ≫ R . We can therefore estimate the pickup probe transfer function: α = V µ M ≈
13 ( γB ) N p ( πr s ) s Q p R ω p L p ! ≈ × V / T . (S7)For an excitation voltage V e , referred to the output of the DAC, the current through the excitation coil is calculatedfrom circuit analysis: I e ≈ A e V e r Q e R ω e L e , (S8)where | A e | = 4 is the gain of the SRS SIM954 amplifier. Note that the SIM954 has an output impedance 3 . ≪
50 Ω.The magnetic field produced by this current at the center of the Helmholtz excitation coil is B e = µ I e r e (cid:16)p ( r e / + r e (cid:17) N e . (S9)Assuming the excitation is resonant with the spin Larmor frequency, the Rabi frequency is Ω e = γ ( B e / / B e is resonant(rotating wave approximation). Therefore the excitation probe transfer function can be estimated as κ = Ω e V e ≈ γµ A e r Q e R ω e L e ! r e (cid:16)p ( r e / + r e (cid:17) N e ≈ . / (ms · V) . (S10)Section II C describes how we used pulsed NMR to measure the values of α and κ . The proximity of the measuredvalues 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 thelow-noise amplifier A is shielded with a 0.5-mm thick copper mesh sleeve. Another copper mesh sleeve shields thebundled RG316 coaxial cables that run up to the top flange of the cryostat. Shields are connected to the ground pinof the A amplifier used as a common ground. We estimate the RF interference noise reduction factor due to theshields to be on the order of 10 within the 1 MHz range centered at 39 .
71 MHz.
II. CHARACTERIZATION OF THE SETUP WITH NUCLEAR MAGNETIC RESONANCEA. Properties of the
Pb spin ensemble
The
Pb isotope has nuclear spin I = 1 / γ π = µhI = 9 .
03 MHz / T , (S11)where µ = 0 . µ N is the magnetic moment of Pb nucleus [2], and the nuclear magneton is µ N /h =7 . / T [3].The chemical formula of PMN-PT is (PbMg / Nb / O ) / − (PbTiO ) / . The number density of Pb nuclearspins in a PMN-PT crystal is given by: n = ρM · N A · .
221 = 3 . × m − , (S12)where ρ = 8 . / cm [4] is the mass density, M = 317 . / mole [2] is the molar mass, and N A is the Avogadro constant.The natural abundance of Pb is 22 .
1% [2].We perform our experiments in the leading magnetic field B = 4 . T = 4 . Pb nuclear spin ensemble is given by [5] µ M = µ nγ ~ I ( I + 1) B k B T = 2 . , (S13)where k B is the Boltzmann constant, µ 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 f ( ν ) = 6 .
33Γ exp − (cid:20) ν − ν )Γ / (2 π ) (cid:21) ln 2 ! , (S14)where Γ / (2 π ) is the full-width at half-maximum, ν is the excitation frequency, and ν 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 T We use the standard NMR saturation recovery scheme to measure the T relaxation time of the Pb nuclearspin ensemble in PMN-PT at 4 . .
66 MHz to 39 .
76 MHz, and whose Rabi frequencies are fixed at 0 .
88 rad / ms. Each pulse duration is 0 . . t , after which a pulsed NMR measurement is per-formed, with spin FID recorded after excitation pulses of 20 ms duration and 180 ms repetition time. The dependenceof the FID amplitude on recovery time t is modeled as an exponential 1 − e − t/T . The best-fit value for the populationrelaxation time is T = (25 . ± .
6) min.
FIG. S3. Saturation of the spin-ensemble excitation spectrum by the sequence described in section II B. Bloch equationsimulations are performed as described in section II C. The initial spin excitation spectrum is shown by the blue dashed line.The spectrum immediately after the saturation step is shown by the orange solid line.
C. Spin-dynamics simulations with Bloch equations
We use the Bloch equations to quantitatively describe the magnetic resonance dynamics of the
Pb nuclear spinensemble [5, 6]. We choose the direction of the z-axis to be along the static magnetic field B . The linearly-polarizedexcitation magnetic field B e = (2Ω e /γ ) cos ( ω t ) is applied in the x -direction. In the reference frame that rotates atthe angular frequency ω around the leading magnetic field, the Bloch equations readd ˜ M x d t = − ˜ M x T + ∆ ω ˜ M y , d ˜ M y d t = − ∆ ω ˜ M x − ˜ M y T − Ω e M z , (S15)d M z d t = Ω e ˜ M y − M z − M ′ T , where ∆ ω = ω − ω is the detuning of spin Larmor frequency ω from the rotating frame frequency, M ′ is the initialensemble magnetization, T is the transverse spin coherence time, and ˜ M x,y are the transverse spin magnetizationcomponents in the rotating frame. The transformation between magnetization in the laboratory and the rotatingframes is M x = ˜ M x cos( ω t ) − ˜ M y sin( ω t ) ,M y = ˜ M x sin( ω t ) + ˜ M y cos( ω t ) , (S16)where in the lab frame M = M x ˆ x + M y ˆ y + M z ˆ z .We numerically solve the Bloch equations using the Runge-Kutta method. The inhomogeneously-broadened spinensemble is represented by 3251 spins, with their Larmor frequencies uniformly distributed in an excitation bandwidthof 65 kHz with 0 .
02 kHz spacing. We simulate the dynamics of each spin independently, and add their contributionsto obtain the total magnetization.The simulation parameters are the spin coherence time T , 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 valuesof these parameters to achieve the minimum value of the goodness-of-fit parameter χ = χ + χ + χ , where thesubscript enumerates the measurements with different pulse duration t p = 0 . , ,
20 ms. For each measurement i = 1 , , χ i = X ν h Re ( F exp [ ν ] − F sim [ ν ]) + Im ( F exp [ ν ] − F sim [ ν ]) i , (S17)where F exp is the Fourier transform of the experimentally detected voltage and F sim is the Fourier transform ofsimulation results, converted into voltage using the transfer coefficient α , and the index ν labels discrete frequencypoints within the window shown in Fig. 2(a) of the main text and Fig. S4. The real part of the Fourier transformcorresponds to the in-phase quadrature, and the imaginary part corresponds to the out-of-phase quadrature of theFID, 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 datastarting at 5 µ s after the end of an excitation pulse. To improve the signal-to-noise ratio, we average the recordedFID response data for several consecutive excitation pulses: 10 data sets are averaged for t p = 0 . 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 perbin for t p = 0 . t p = 2 ms pulse duration, and 2 points per bin for t p = 20 mspulse duration. The error bars shown in Fig. 2(a) of the main text and Fig. S4 are the standard deviation of thepoints within each bin.The spin ensemble was saturated before every FID measurement, and the FID measurements started after await time ≈ T after saturation. Therefore the initial magnetization at the start of every FID measurement was µ M ′ = (0 . ± . µ M = (1 . ± .
2) nT, where M is the thermal equilibrium ensemble magnetization given byEq. (S13).Using the measurements shown in Fig. 2(a) in the main text and Fig. S4, we extract the best-fit parameter values: T = (16 . ± .
9) ms ,κ = (0 . ± . / (ms · V) , (S18) α = (2 . ± . × V / T . FIG. S4. Measurements of
Pb FID spectra following a spin excitation pulse of length t p , as indicated in the panels. Weperformed fitting simultaneously to in-phase (blue) and out-of-phase (orange) components of Fourier transforms of averagedFID from three data sets: with excitation pulse duration t p = 20 ms shown in Fig. 2(a) in the main text, (a) with excitationpulse duration t p = 0 . t p = 2 ms. Data points were binned, and the error bars showone standard deviation in each bin. The lines show the best-fit simulation of the spin response, with the light-colored narrowbands indicating the range of simulation results if parameters T , κ , and α are varied by one standard deviation away fromtheir best-fit values. The uncertainties are evaluated by bootstrapping: the frequency-domain data are down-sampled into 16 groups, andthe fit is performed independently on each data group; the uncertainty is given by the standard deviation of thebest-fit parameter values. The proximity of the best-fit values of transfer parameters α and κ to the estimates ineqns. (S7) and (S10) validates the analysis of the apparatus design in Sec. I D. 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 FIDdata for a range of excitation Rabi frequencies Ω e . For these measurements we keep the excitation pulse width at20 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 consecutiveexcitation pulses, spaced by 180 ms. After each pulse, we sample the FID voltage, starting 5 µ s after the end of thepulse, and lasting for 16 . µ s. We average the 100 FID data sets, and calculate the discrete Fourier transform F [ n ] ofthe averaged FID, where index n labels frequency points. Since we only sample the beginning of the FID, before it canstart to decay, we model it as a sinusoidal signal at the excitation carrier frequency. We extract the amplitude of thespin ensemble transverse magnetization by numerically integrating the power spectrum | F [ n ] | over a 400 kHz-widefrequency band centered at the excitation carrier frequency, and using the pickup probe transfer function α to convertthe voltage to magnetization. Uncertainties are calculated using bootstrapping: we group the 100 FID data sets into5 sets of 20 and perform analysis on these 5 sets independently. Error bars are set at the standard deviation of theresults for these 5 sets. To obtain the theory curve in Fig. 2(c) of main text, we use our Bloch equation model togenerate 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 ω , the steady-state transverse magne-tization of an unsaturated homogeneously-broadened spin ensemble is given by [5] M = L ( ω − ω ) M Ω e T cos ( ω t ) , (S19)where M is the longitudinal magnetization, T is the transverse coherence time, ω is the Larmor angular frequency,and L is the Lorentzian lineshape function: L ( ω − ω ) = 11 + ( ω − ω ) T . (S20)Let us describe the spin ensemble inhomogeneous broadening with the excitation lineshape h ( ω + ∆), normalizedsuch that Z ∞−∞ h ( ω + ∆) d ∆ = 1 . (S21)Under CW excitation, the steady-state transverse magnetization is then M = uM Ω e T cos ( ω t ) , (S22)where the spectral u factor is given by the integral over the lineshape: u = Z ∞−∞ L ( ω + ∆ − ω ) h ( ω + ∆) d ∆ . (S23)Let us estimate the value of u . Our NMR measurements indicate that the excitation spectrum is much broader than1 /T , therefore we can approximate the Lorentzian with the delta-function: L ( ω + ∆ − ω ) ≈ ( π/T ) δ ( ω + ∆ − ω ).Furthermore, we approximate the excitation spectrum as a rectangular function, centered at ω , with full width Γand height 1 / Γ. Then, provided | ω − ω | < Γ /
2, we can approximate u ≈ πh ( ω ) /T ≈ π/ (Γ T ) . (S24)In order to more accurately determine u , we solved the Bloch equations with the experimentally-determined values T = (16 . ± .
9) ms and excitation spectrum with Γ / (2 π ) = (78 ±
2) kHz (Fig. 2(b) in the main text). We obtained u = (3 . ± . × − , (S25)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 ensuregood electrical contact, we paint the faces with graphite paint, which is removed after polarization. We connect thecrystal to the Trek model 610E-G-CE high-voltage amplifier as shown in Fig. S5(a). The amplifier measures theapplied voltage and the current through the sample. In order to measure the ferroelectric hysteresis loop, we applytriangular voltage ramps with alternating polarities, Fig. S5(b). Current spikes are visible when the applied voltageis sufficient to reverse the ferroelectric polarization. In this experimental run the crystal started with a remanentpolarization corresponding to positive polarity, so there is no current spike during the first ramp. We obtain thesample polarization by integrating the current: P ( t ) = q ( t ) πr s = 1 πr s Z t I ( t ′ )d t ′ , (S26)where q ( t ) is the electric charge on the crystal surface and r s = 2 . P r persists after the voltage has been ramped down to zero. We verified that the remanentpolarization does not decay after thermal cycling of the sample. V. NUCLEAR SPIN DYNAMICS DUE TO THE EDM INTERACTION WITH AXION-LIKE DARKMATTERA. P,T-odd axion-like dark matter physics
Axion-like cold dark matter is a classical field: a ( t ) = a cos ( ω a t ), where ω a ≈ m a c / ~ . If the axion-field energydensity dominates dark matter, then ρ DM = m a a / ≈ . × − GeV [7]. In the QCD Lagrangian, this gives riseto an oscillating θ angle: θ ( t ) = af a = a f a cos ( ω a t ) . (S27) FIG. S5. (a) Ferroelectric polarization setup, showing the signal generator controlling the high-voltage amplifier (HVA) thatis connected to the electrodes in contact with the sample. TREK Model 610E high-voltage supply/amplifier/controller housesthe HVA, as well as an ammeter A , a unity gain buffer amplifier A , and a voltmeter V . (b) Voltage applied at the output ofHVA is measured at the voltmeter V and converted to the voltage across the sample (blue dashed line). The current measuredat the ammeter A is plotted as the orange full line. Let us consider the nucleon EDM induced by axion-like dark matter: d n = g d a = 2 . × − θ e · cm = 2 . × − θ e · fm , (S28)calculated with 40% accuracy [8, 9]. Here g d is the EDM coupling constant [9], introduced in the Lagrangian term: L EDM = − i g d a ¯Ψ N σ µν γ Ψ N F µν , (S29)where Ψ N is the nucleon wavefunction, F µν is the electromagnetic field tensor, and σ and γ are the standard Diracmatrices.From eqs. (S27,S28) we get the relationship between g d and f a : g d = 2 . × − e · cm f a = 3 . × − GeV − f a , (S30)where we used the natural unit conversions: 1 cm = 5 × GeV − and e = 0 . m a = 6 × − eV (cid:18) GeV f a (cid:19) , (S31)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–13] is defined as: S = e (cid:18) h r r i − Z h r ih r i (cid:19) , (S32)where e is the elementary electric charge, Z is the atomic number, and h r k i = R r k ρ ( r ) d r are the integrals overnuclear charge density ρ ( r ). The Schiff moment sources the P- and T-odd electrostatic potential ϕ ( r ) = 4 π ( S · ∇ ) δ ( r ) . (S33)Importantly, the definition of the Schiff moment in Ref. [14] differs from this one by a factor of 4 π . We adopt thedefinition 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]. Thecontribution of P,T-odd nuclear forces is larger than the contribution of nucleon EDM [12]. Let us consider the twocontributions separately, in the case of the Pb nucleus, whose ground state is I π = 1 / − , having a neutron 3 p / hole in a closed-shell magic nucleus.0
1. Nuclear Schiff moments induced by nucleon EDM
This contribution is due to non-coincident densities of nuclear charge and dipole moment. It can be estimated for
Pb using Eq. (8.76) in Ref. [14]: 4 πS EDM ≈ d n × π
25 ( K + 1) II ( I + 1) r , (S34)where K = ( ℓ − I + 1) = 1 and r = 1 . A / = 7 . S EDM ≈ d n × . More detailedcalculations [15, 16] give the result: S EDM = d n × . . (S35)In order to connect this to QCD axion physics, we use Eq. (S28): S EDM = g d a × . = 5 × − θ e · fm . (S36)
2. 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]: W = G F √ η m σ · ∇ ρ ( r ) , (S37)where G F ≈ − GeV − is the Fermi constant, m is the nucleon mass, σ is its spin, and ρ ( r ) is the density of corenucleons. A vacuum θ angle gives rise to this interaction via the P,T-odd pion-nucleon coupling constant [14, 17]: η = 1 . × θ. (S38)Next we need to calculate the nuclear Schiff moment induced by the interaction (S37). Reference [12] states thatthe 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
Hg is estimated as0 . × − η e · fm . 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 Hg is estimated as 2 . × − η e · fm , and the Schiff moment of Hg is estimated as − . × − η e · fm .The issue is complicated by nuclear many-body effects. These were numerically calculated for Hg in Refs. [19, 20],giving a factor ∼
10 reduction in the Schiff moment. However the physical origin of such a strong reduction is notclear. The only effect, not included in the shell model, that can change the value of the Schiff moment is the collectivenuclear octupole deformation, and, if anything, that should increase the Schiff moment. Reference [21] gives a resultfor
Hg that is ∼
10% away from the shell-model estimate. These authors attribute the Schiff moment suppressionin Ref. [20] to the mixing with the J π = 1 / − state, for which they get a small Schiff moment value. However thissmall 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 matrixelements are likely to be small unless a lot of excited states are carefully taken into account. This suggests that thecalculation may have large intrinsic uncertainties.Importantly,
Pb is close to a magic nucleus, which means that many-body effects should not play an importantrole here. Therefore, until the many-body effects can be better understood, for
Pb we retain the single-particleestimate of Ref. [18]: S η = 2 × − η e · fm = 0 . θ e · fm . (S39)Note that this is a factor of eight larger than the result (13) in Ref. [22], where the Pb Schiff moment was taken tobe the same as for the many-body suppressed
Hg. We can also see that this contribution is a factor of eight largerthan the EDM contribution in Eq. (S36). We therefore neglect the EDM contribution, and use the above estimate(S39).Similar estimates were performed for
Hg in Ref. [23].1
C. Nuclear Schiff moment-induced spin energy shift in ferroelectric PMN-PT
The energy shift of each nuclear spin sublevel of a Pb ion in ferroelectric PbTiO is estimated in Refs. [24, 25].The result of the full quantum chemistry calculation [26] is:∆ ǫ = 1 . × x .
58 Å
Sea B [eV] = 1 . × − x Å S e · fm [eV] , (S40)where x is the displacement of the Pb ion with respect to the center of the oxygen cage, S is the magnitude of theSchiff moment of the Pb nucleus, and a B = 0 .
53 Å is the Bohr radius. The nuclear spin is I = 1 /
2; each of the twonuclear 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:Ω a = 12 2∆ ǫ ~ = 1 . × x Å S e · fm [rad / s] , (S41)where we used ~ = 6 . × − eV · s. We note that the spin driving field is “linearly polarized”, and therefore theRabi frequency contains an extra factor of 1 /
2, which arises because only one of the two counter-rotating componentsof the linearly polarized drive is resonant (rotating wave approximation).Density functional theory calculations for PMN-PT give the Pb cation displacement from the center of the oxygencage: x = 0 .
39 Å, and the average polarization: P = 55 µ C / cm [27]. Our experiment was performed with thecrystal polarization P r = 22 µ C / cm , therefore we scale the average displacement to x = 0 .
16 Å.For
Pb in ferroelectric PMN-PT we can use Eq. (S39,S40) and x = 0 .
16 Å to get:∆ ǫ = 8 × − θ [eV] , (S42)Ω a = 1 . × θ [rad / s] . (S43)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∆ ǫ = d n [e · cm] × . × [V / cm] . (S44)We can extract the effective electric field (which includes the Schiff screening factor [28]): E ∗ = ∆ ǫ/d n = 340 [kV / cm] . (S45)For the drive Rabi frequency we obtain:Ω a = 1 . × g d a . × − [GeV − ] [rad / s] , (S46) ~ Ω a = 2 . × − ( g d a ) [GeV] , (S47)where g d is in GeV − and a = √ ρ DM /m a is in GeV. Let us introduce the sensitivity factor ξ , defined as ~ Ω a = ξg d a .Its estimated value is therefore ξ = 2 . × − [GeV ] . (S48)There are several contributions to the theoretical uncertainty in E ∗ and ξ . The uncertainty of the QCD calculationsis ≈
40% [8, 9]. The uncertainty of the solid-state calculation of the nuclear spin energy shift due to the Schiff momentis ≈
30% [24–26]. Therefore we estimate the total theoretical uncertainty in E ∗ and ξ at ≈ VI. NUCLEAR SPIN DYNAMICS DUE TO THE GRADIENT INTERACTION WITH AXION-LIKEDARK MATTER
The non-relativistic Hamiltonian for the gradient interaction of spin I with axion-like dark matter field a ( r , t ) is H aNN = g aNN ∇ a ( r , t ) · I , (S49)where g aNN is the coupling strength measured in units of GeV − , and we used natural units here ~ = c = 1 [9, 29]. Inthe first approximation we can write the axion-like dark matter field as: a ( r , t ) ≈ a cos ( ω a t − k · r ) , (S50)2where the field amplitude a is fixed by the assumption that it dominates the dark matter energy density: ρ DM = m a a / . × − GeV [7, 9]. We approximate the instantaneous value of the gradient ∇ a ≈ m a v a , where v isthe instantaneous value of the velocity of the ALP field in the laboratory frame. The Hamiltonian in natural unitsbecomes: H aNN = ( g aNN a ) m a v · I cos ( ω a t ) . (S51)The product g aNN a is dimensionless, so we can restore the values of fundamental constants by dimensional analysis: H aNN = ( g aNN a ) m a c v c · I cos ( ω a t ) . (S52)This interaction exerts a torque on nuclear spins, with the drive Rabi frequency given by ~ Ω a = 12 ( g aNN a ) m a c v ⊥ c , (S53)where v ⊥ is the component of the velocity perpendicular to the direction of the leading field B . As in the previoussection, the spin driving field is “linearly polarized”, and therefore the Rabi frequency contains an extra factor of1 /
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
Pb nuclearspins at carrier angular frequency ω a with Rabi frequency Ω a . The steady-state transverse spin magnetization thatdevelops under the action of this driving field is given by Eq. (1) of the main text. The resulting voltage recorded bythe ADC is: V a ( t ) = αµ M a = αuµ M Ω a T cos ( ω a t ) . (S54)The time-averaged power in this signal is h V a i = 12 ( αuµ M Ω a T ) . (S55)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-likedark matter halo that follows the standard halo model [30, 31]. In this model the ALP speeds v in the Galactic framefollow the Maxwell-Boltzmann distribution f gal ( v ) = 4 v √ πv e − v /v , (S56)where v ≈
220 km / s is the most probable speed [31]. The laboratory frame moves relative to the Galactic frame withthe average speed v lab ≈
232 km / s which has annual and daily modulations due to, respectively, Earth’s revolutionabout the Sun and Earth’s rotation around its axis [32]. The distribution of ALP speeds broadens the Fourier spectrumof the ALP field a ( t ), giving it a characteristic linewidth ≈ v ν a /c ≈ − ν a . The power spectrum of the ALP field a ( t ) is given by the function f ( ν ) = 2 c √ πv v lab ν a exp (cid:18) − c v ν − ν a ν a − v v (cid:19) sinh β, (S57)where β = 2 cv lab v s ν − ν a ) ν a . (S58)This spectral function is normalized so that Z ∞ ν a f ( ν ) dν = 1 . (S59)This is the spectral lineshape used in searches for ALP-photon interactions [33–36].3 A. EDM search
In our search for the ALP EDM interaction, the recorded voltage V EDM is directly proportional to the ALP field a ( t ).Therefore the voltage power spectral density will have the same spectral shape f ( ν ). We make use of Parseval’stheorem to ensure that the time-averaged power, Eq. (S55), matches the integral of the Fourier power spectrum, withthe lineshape normalized as in Eq. (S59). The result is the expression for the voltage power spectrum: V ( ν ) = 12 ( αuµ M Ω a T ) f ( ν ) = 12 ( αuµ M T ) (cid:18) ξg d a ~ (cid:19) f ( ν ) . (S60) B. Gradient search
In our search for the ALP gradient interaction, the recorded voltage V gr is proportional to the gradient of the ALPfield, which includes the velocity of the ALP field in the lab frame. Therefore the voltage power spectrum has adifferent form: f ( ν ) = 2 c v + v sin ζ ν − ν a ν a (cid:20) sin ζ + 1 β (cid:18) coth β − β (cid:19) (cid:0) − ζ (cid:1)(cid:21) f ( ν ) , (S61)where ζ is the angle between the vectors B and v lab . This spectral function is normalized so that Z ∞ ν a f ( ν ) dν = 1 . (S62)A detailed analysis of the ALP velocity distribution in the laboratory reference frame, resulting in the ALP gradientspectral line shape (S61), will be published elsewhere.Again, making use of Parseval’s theorem to ensure that the time-averaged power equals the integral of the Fourierpower spectrum, we write the expression for the voltage power spectrum: V ( ν ) = 12 ( αuµ M T ) (cid:18) g aNN a m a c ~ (cid:19) v + v sin ζc f ( ν ) . (S63) VIII. DATA ACQUISITION AND ANALYSIS FOR THE AXION-LIKE DARK MATTER SEARCH
The experimental search for axion-like dark matter took place on October 7, 2019. We varied the static magneticfield to sweep the spin Larmor frequency, starting at 40 .
16 MHz and ending at 39 .
16 MHz, in 21 steps with step size of50 kHz. This corresponds to magnetic fields between B = 4 .
45 T and 4 .
35 T. During the recording of data sensitiveto axion-like dark matter the superconducting magnet is in persistent mode with the power supply turned off and theexcitation probe and cancellation coil are terminated with a 50 Ω resistor. Data are recorded at the ADC samplingrate of 250 MS / s and saved to a hard drive using the first-in, first-out mode of the ADC. At each value of B we record58 s of “scan” data, immediately followed by 58 s of “re-scan” data, analyzed as described below. During the searchwe perform three pulsed NMR calibrations, at the first, the last, and the middle values of the magnetic field. Eachcalibration consists of FID data taken at five different excitation carrier frequencies near the corresponding Larmorfrequency, Fig. S6.The data sensitive to axion-like dark matter are analyzed using Matlab on the Shared Computing Cluster, whichis administered by Boston University’s Research Computing Services. The data-processing procedure consists of thefollowing steps.(1) Divide data at each value of the magnetic field B into 27 blocks. Each block contains 2 points and correspondsto 2 .
15 s of data. The block duration exceeds the axion-like dark matter coherence time of ≈
25 ms in the standardhalo model.(2) Perform the discrete Fourier transform on each data block without any windowing function, to obtain thespectral density F [ ν ]. Select the analysis frequency range between 39 . . | F [ ν ] | =Re( F [ ν ]) + Im( F [ ν ]) is a chi-square distribution with two degrees of freedom. Then, average | F [ ν ] | obtainedfrom 27 segments. Take a histogram and fit to a chi-square distribution and confirm that it has 54 degrees offreedom. Repeat this for the scan (and separately, re-scan) data at each resonant frequency.4 FIG. S6. NMR calibration at the three values of the bias field B . FID data are recorded after excitation pulses at Rabifrequency Ω e = 0 .
88 rad / ms and pulse length 20 ms. The excitation carrier frequency is plotted on the x-axis. Following theprocedure 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 onlyfree parameter is the central frequency. The best-fit values of the central frequency for the three calibration data sets are: ν = (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 rawpower spectral densities. The points where this difference is above a threshold are marked as narrow spectrallines 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 EDMinteraction f ( ν ) given in Eq. (S57). The separation between distinct ALP search frequencies is set to the ALPsignal linewidth 3( v /c ) ν /
4, where ν is the central Larmor frequency, determined by the value of the biasfield B [32, 37].(5) Model the histogram of the optimally-filtered power spectral density with 100 bins as the Gaussian distributionwith mean µ and standard deviation σ . Calculate the detection threshold at µ + 3 . σ , corresponding to a5 σ detection with 95% confidence level. Points above the threshold are ALP detection candidates. A detailedexplanation 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 in the scan.The spin response to an axion-like dark matter signal will only appear in the data set where B is such that theALP Compton frequency is within the magnetic resonance excitation spectrum. For each data set we use the 80 kHzfrequency band centered at Larmor frequency ν , corresponding to the excitation spectrum, to search for the ALPsignal, as described above. The rest of the spectral data within the 1 MHz scan range are used to reject residualbackground RF interference, which is not eliminated by the Savitzky-Golay filter. In addition, re-scan measurementsare analyzed to eliminate statistical fluctuations, which are expected, given the large bandwidth of our search (look-elsewhere effect). The analysis procedure is as follows.(a) At each value of bias magnetic field we consider ≈ . σ 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 excitationspectrum) to the candidate frequencies from the “background” data sets (for which the frequency is outside theexcitation spectrum). If the candidate frequency appears in one of the background data sets, it is rejected asRF interference. On average this eliminates ≈
28 candidates at each value of B .(c) We compare candidate frequencies from the scan and re-scan data sets. If a candidate frequency appears onlyin one of those data sets, it is rejected as a statistical fluctuation. On average this eliminates ≈ B .5This analysis procedure rejects all candidates above the 3 . σ threshold at all values of B . We do not detect anaxion-like dark matter signal. Therefore, for each value of B , we quote the g d coupling value that corresponds to the5 σ value of the power spectral density as the 95% confidence interval limit [36]. FIG. S7. Angle ζ between the bias magnetic field B and the laboratory velocity vector v lab as a function of time offset fromthe start of the experimental search for axion-like dark matter at 18:41 UTC on October 7, 2019 to 00:30 UTC on October 8,2019. Stars indicate the times at which data are recorded for different values of B . The magnitude of the laboratory velocityis v lab = 226 . ± . / s for the entire duration of data taking. The velocity v lab is calculated for the Physics Department atBoston University (42°20 ′ . , ′ . We search for the gradient coupling g aNN of axion-like dark matter using the same steps as described above, withthe standard halo model lineshape in step (4) replaced by the gradient coupling lineshape f ( ν ), given in Eq. (S61).We calculate the angle ζ at each value of B during the scan, based on the coordinates of our laboratory and the timeat which the data are recorded, Fig. S7.Our analysis for the gradient coupling g aNN rejects all candidates above the 3 . σ threshold at all values of B .Therefore, for each value of B , we quote the g aNN coupling value that corresponds to the 5 σ value of the powerspectral density as the 95% confidence interval limit. We note that the variation in ζ throughout the scan meansthat 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 thisdifference 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 signalswith the lineshape given by Eq. (S57). Figure S8(a) shows the spectrum with an injected signal at Compton frequency ν a = 39 . g d = 1 . × − GeV − . After optimal filtering, the injected signalshows up as a candidate with amplitude 101 fT , as shown in Fig. S8(b). The histogram of the optimally-filtereddata points shows that this injected signal is detected at 20 σ significance, Fig. S8(c). We test the recovery of thecoupling strength by injecting 10 simulated signals, whose coupling strength is varied between g d = 7 . × − GeV − and g d = 7 . × − GeV − and whose Compton frequencies are selected randomly between ν a = 39 . ν a = 39 . . ± . B. Projected sensitivity reach
Our experimental results demonstrate the feasibility of using solid-state nuclear magnetic resonance to search foraxion-like dark matter. There are several bounds on the relevant interactions of axion-like dark matter in this massrange, 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 evaded6
FIG. S8. Injecting simulated axion-like dark matter signals into experimental data. (a) A 400 Hz-wide band of experimentalpower spectrum, with an injected signal at Compton frequency ν a = 39 . g d = 1 . × − GeV − .Experimental data are shown as blue circles, the injected signal is shown as the orange line. (b) The optimally-filtered spectrumwithin the 50 kHz frequency band centered on ν a . (c) The histogram of optimally-filtered PSD data (blue circles) within the80 kHz band centered on the Larmor frequency ν = 39 .
16 MHz. Data are sorted into 100 bins, and the Gaussian fit is shownas the orange line. The 3 . σ detection threshold (vertical dashed black line) is at 17 fT . (d) Recovered coupling for injectedsignals with coupling strengths varying logarithmically from g d = 7 . × − GeV − to g d = 7 . × − GeV − at differentCompton frequencies, sampled randomly between ν a = 39 . ν a = 39 . . ± . altogether [44, 45]. Stringent experimental limits on g d and g aNN exist at much lower ALP masses [29, 46–50], butthe 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 Pb nuclear spin ensemble, and (2)the small size of our PMN-PT sample. We plan to circumvent the inhomogeneous broadening by concentrating ourfuture searches on the lower Compton frequencies ( ν a < T spincoherence time, rather than CSA. The long T relaxation time will allow us to pre-polarize the nuclear spins, retainingtheir polarization even at lower fields. We plan to use Superconducting Quantum Interference Devices (SQUIDs) todetect the transverse magnetization in this frequency range. The green dashed curves in Fig. S9 show the projectedexperimental sensitivity for the search with the same 4 . /T NMR linewidth, and the cutoff at high frequencies is set by the Larmor frequencyat 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 volumeof 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 ≈ ≈ FIG. S9. Existing bounds and sensitivity projections for the: (a) EDM and (b) gradient coupling of axion-like dark matter.The region shaded in red is the exclusion at 95% confidence level placed by this work (CASPEr-e). The purple line shows theQCD axion coupling band. The darker purple color shows the mass range motivated by theory [9]. The blue regions mark themass ranges where the ADMX and HAYSTAC experiments have probed the QCD axion-photon coupling [33, 34]. The greenregion is excluded by analysis of cooling in supernova SN1987A, with color gradient indicating theoretical uncertainty [9]. Thedashed green line marks the projected 5 σ sensitivity of our CASPEr-e search with a 4 . σ sensitivity of our CASPEr-e search with an 80 cm sample, operating at 100 mKtemperature. Implementing a resonant coupling circuit will enable operation with a smaller sample. The black dashed linemarks the sensitivity limited by the quantum spin projection noise [28]. This is sufficient to detect the EDM coupling of theQCD axion across the 6-decade mass range from ≈ . ≈
500 neV. The other bounds are as follows. (a) The pink regionis 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 theneutron 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 yellowregion is excluded by the new-force search with K- He comagnetometer [46]. The bounds are shown as published, althoughcorrections should be made to some of the low-mass limits, due to stochastic fluctuations of the axion-like dark matter field [51].[1] J. B. Miller, B. H. Suits, A. N. Garroway, and M. A. Hepp, Concepts in Magnetic Resonance , 125 (2000).[2] J. E. Sansonetti and W. C. Martin, NIST Standard Reference Database (2013).[3] P. J. Mohr, D. B. Newell, and B. N. Taylor, Journal of Physical and Chemical Reference Data (2016).[4] F. Kochary, M. D. Aggarwal, A. K. Batra, R. Hawrami, D. Lianos, and A. Burger,Journal of Materials Science: Materials in Electronics , 1058 (2007).[5] A. Abragam, The principles of nuclear magnetism (1961).[6] F. Bloch, Physical Review , 460 (1946).[7] M. Tanabashi et al. (Particle Data Group), Phys. Rev. D , 030001 (2018).[8] M. Pospelov and A. Ritz, arXiv:hep-ph/9908508 , 1 (1999).[9] P. W. Graham and S. Rajendran, Physical Review D , 035023 (2013).[10] L. Schiff, Physical Review , 2194 (1963).[11] P. G. Sandars, Physical Review Letters , 1396 (1967).[12] O. P. Sushkov, V. V. Flambaum, and I. B. Khriplovich, Sov. Phys. JETP , 873 (1984).[13] V. V. Flambaum and J. S. Ginges, Physical Review A , 032113 (2002).[14] I. B. Khriplovich and S. K. Lamoreaux, CP Violation Without Strangeness (Springer Berlin Heidelberg, Berlin, Heidelberg,1997).[15] V. A. Dzuba, V. V. Flambaum, J. S. Ginges, and M. G. Kozlov, Physical Review A , 7 (2002).[16] V. F. Dmitriev and R. A. Sen’kov, Physical Review Letters , 212303 (2003).[17] V. V. Flambaum, D. Demille, and M. G. Kozlov, Physical Review Letters , 103003 (2014), arXiv:1406.6479.[18] V. V. Flambaum, I. B. Khriplovich, and O. P. Sushkov, Nuclear Physics, Section A , 750 (1986).[19] V. F. Dmitriev and R. A. Sen’kov, Physics of Atomic Nuclei , 1940 (2003).[20] S. Ban, J. Dobaczewski, J. Engel, and A. Shukla, Physical Review C - Nuclear Physics , 015501 (2010).[21] K. Yanase and N. Shimizu, arXiv: 2006.15142 (2020).[22] V. V. Flambaum and V. A. Dzuba, Physical Review A , 042504 (2020). [23] Y. V. Stadnik and V. V. Flambaum, Physical Review D , 043522 (2014).[24] T. N. Mukhamedjanov and O. P. Sushkov, Physical Review A , 34501 (2005), arXiv:0411226 [physics].[25] J. A. Ludlow and O. P. Sushkov, Journal of Physics B: Atomic, Molecular and Optical Physics , 085001 (2013).[26] L. V. Skripnikov and A. V. Titov, Journal of Chemical Physics , 054115 (2016).[27] I. Grinberg and A. M. Rappe, Physical Review B , 220101 (2004).[28] D. Budker, P. W. Graham, M. Ledbetter, S. Rajendran, and A. O. Sushkov, Phys. Rev. X , 021030 (2014).[29] A. Garcon, J. W. Blanchard, G. P. Centers, N. L. Figueroa, P. W. Graham, D. F. J. Kimball, S. Rajendran, A. O. Sushkov,Y. V. Stadnik, A. Wickenbrock, T. Wu, and D. Budker, Science Advances , eaax4539 (2019).[30] M. S. Turner, Physical Review D , 3572 (1990).[31] N. W. Evans, C. A. O’Hare, and C. McCabe, Physical Review D , 023012 (2019).[32] J. W. Foster, N. L. Rodd, and B. R. Safdi, Physical Review D , 123006 (2018).[33] N. Du, N. Force, R. Khatiwada, E. Lentz, R. Ottens, L. J. Rosenberg, G. Rybka, G. Carosi, N. Woollett, D. Bowring,A. S. Chou, A. Sonnenschein, W. Wester, C. Boutan, N. S. Oblath, R. Bradley, E. J. Daw, A. V. Dixit, J. Clarke, S. R.O’Kelley, N. Crisosto, J. R. Gleason, S. Jois, P. Sikivie, I. Stern, N. S. Sullivan, D. B. Tanner, and G. C. Hilton, PhysicalReview Letters , 151301 (2018).[34] B. M. Brubaker, L. Zhong, Y. V. Gurevich, S. B. Cahn, S. K. Lamoreaux, M. Simanovskaia, J. R. Root, S. M. Lewis, S. AlKenany, K. M. Backes, I. Urdinaran, N. M. Rapidis, T. M. Shokair, K. A. van Bibber, D. A. Palken, M. Malnou, W. F.Kindel, M. A. Anil, K. W. Lehnert, and G. Carosi, Physical Review Letters , 061302 (2017).[35] J. L. Ouellet, C. P. Salemi, J. W. Foster, R. Henning, Z. Bogorad, J. M. Conrad, J. A. Formaggio, Y. Kahn, J. Minervini,A. Radovinsky, N. L. Rodd, B. R. Safdi, J. Thaler, D. Winklehner, and L. Winslow, Physical Review Letters , 121802(2019).[36] A. V. Gramolin, D. Aybas, D. Johnson, J. Adam, and A. O. Sushkov, Nature Physics (2020), 10.1038/s41567-020-1006-6.[37] B. M. Brubaker, L. Zhong, S. K. Lamoreaux, K. W. Lehnert, and K. A. van Bibber, Physical Review D , 123008 (2017).[38] B. Pelssers, M. Lawson, and N. L. Figueroa, “The axion search simulation for lots of experiments (tassle),” https://github.com/thevorpalblade/tassle (2020).[39] T. P. Robitaille et al. (Astropy Collaboration), Astronomy & Astrophysics , A33 (2013),arXiv:1307.6212 [astro-ph.IM].[40] A. M. Price-Whelan et al. (Astropy Collaboration), The Astronomical Journal , 123 (2018),arXiv:1801.02634 [astro-ph.IM].[41] G. G. Raffelt, in Axions: Theory, Cosmology, and Experimental Searches (Springer Berlin Heidelberg, Berlin, Heidelberg,2008) pp. 51–71.[42] P. W. Graham, I. G. Irastorza, S. K. Lamoreaux, A. Lindner, and K. A. van Bibber, Annual Review of Nuclear andParticle Science , 485 (2015).[43] K. Blum, R. T. D’Agnolo, M. Lisanti, and B. R. Safdi, Physics Letters B , 30 (2014).[44] W. DeRocco, P. W. Graham, and S. Rajendran, Physical Review D , 075015 (2020).[45] N. Bar, K. Blum, and G. D’amico, Physical Review D , 123025 (2020).[46] G. Vasilakis, J. M. Brown, T. W. Kornack, and M. V. Romalis, Physical Review Letters , 261801 (2009).[47] C. Abel, N. J. Ayres, G. Ban, G. Bison, K. Bodek, V. Bondar, M. Daum, M. Fairbairn, V. V. Flambaum, P. Geltenbort,K. Green, W. C. Griffith, M. van der Grinten, Z. D. Grujić, P. G. Harris, N. Hild, P. Iaydjiev, S. N. Ivanov, M. Kasprzak,Y. Kermaidic, K. Kirch, H.-C. Koch, S. Komposch, P. A. Koss, A. Kozela, J. Krempel, B. Lauss, T. Lefort, Y. Lemière,D. J. E. Marsh, P. Mohanmurthy, A. Mtchedlishvili, M. Musgrave, F. M. Piegsa, G. Pignol, M. Rawlik, D. Rebreyend,D. Ries, S. Roccia, D. Rozp¸edzik, P. Schmidt-Wellenburg, N. Severijns, D. Shiers, Y. V. Stadnik, A. Weis, E. Wursten,J. Zejma, and G. Zsigmond, Physical Review X , 041034 (2017).[48] T. Wu, J. W. Blanchard, G. P. Centers, N. L. Figueroa, A. Garcon, P. W. Graham, D. F. J. Kimball, S. Rajendran, Y. V.Stadnik, A. O. Sushkov, A. Wickenbrock, and D. Budker, Physical Review Letters , 191302 (2019).[49] W. A. Terrano, E. G. Adelberger, C. A. Hagedorn, and B. R. Heckel, Physical Review Letters , 231301 (2019).[50] T. S. Roussy, D. A. Palken, W. B. Cairncross, B. M. Brubaker, D. N. Gresh, M. Grau, K. C. Cossel, K. B. Ng, Y. Shagam,Y. Zhou, V. V. Flambaum, K. W. Lehnert, J. Ye, and E. A. Cornell, arXiv: 2006.15787 (2020).[51] G. P. Centers, J. W. Blanchard, J. Conrad, N. L. Figueroa, A. Garcon, A. V. Gramolin, D. F. J. Kimball, M. Lawson,B. Pelssers, J. A. Smiga, et al.et al.