Muon g-2 and EDM experiments as muonic dark matter detectors

The detection of ultralight dark matter through interactions with nucleons, electrons, and photons has been explored in depth. In this work we propose to use precision muon experiments, specifically muon g-2 and electric dipole moment measurements, to detect ultralight dark matter that couples predominantly to muons. We set direct, terrestrial limits on DM-muon interactions using existing g-2 data, and show that a time-resolved reanalysis of ongoing and upcoming precession experiments will be sensitive to dark matter signals. Intriguingly, we also find that the current muon g-2 anomaly can be explained by a spin torque applied to muons from a pseudoscalar dark matter background that induces an oscillating electric dipole moment for the muon. This explanation may be verified by a time-resolved reanalysis.

Despite the presence of dark matter (DM) and its gravitational interactions being well established, its particle nature and non-gravitational interactions with the standard model (SM) are yet to be illuminated. While the elementary dark matter mass could span many orders of magnitude, the ultralight dark matter regime, 10 −22 eV ≤ m dm eV, has received much attention recently. These ultralight particles arise naturally in solutions to tuning problems, e.g. the axion [1] and the relaxion [2], as well as in the string landscape. Furthermore, they also have attractive production mechanisms -misalignment for scalars [3], inflationary production for vectors [4] and parametric resonance for both [5][6][7].
Traditional direct detection experiments targeting the WIMP scale are not sensitive to ultralight DM, so a plethora of experiments have been performed and proposed in recent years exploiting the wave-like properties of this mass regime. Yet these have exclusively tested dark matter couplings to photons [8], electrons [8], protons, and neutrons [9,10]. Meanwhile, the muon g-2 anomaly [11,12] has led to exploration of theories with dark forces that predominantly couple to muons and experimental proposals to find them [13][14][15][16]. Similarly, dark matter itself could dominantly couple to muons. In this work, we study such models and explore the possibility of precision muon experiments directly detecting such muophilic dark matter.
Muon g-2 and EDM experiments, such as the measurement done at BNL [11] in 2004, the ongoing work at Fermilab [17] and J-PARC [18], and the proposed frozen spin experiments [19,20] 1 , are precision efforts to track the time evolution of muon spins subject to an external magnetic field. The primary aim of the g-2 experiments [11,17,18] is the determination of the muon's magnetic dipole moment (MDM). However, they are sensitive to any new physics which sufficiently alters the precession dynamics of muon spins. For example, the existence of a muon electric dipole moment (EDM) has been constrained in this manner by the BNL experiment [22] and will be further tested at Fermilab and J-PARC. The frozen spin proposals are a more sensitive, dedicated search for this EDM signal. A coherent dark matter background may couple to muons in these experiments and alter their precession by applying a spin torque and by possibly altering their orbital trajectories. This results in a characteristic DM precession signal which is observable in these experimentswe thus propose to repurpose muon precession experiments as dark matter detectors.
DM perturbations to precession may yield a variety of signals in these experiments depending on the nature of the DM candidate. Some candidates would have noticeably altered the form of the precession signal in the existing analysis of BNL, allowing us to place immediate constraints. These limits will become more stringent with ongoing and future measurements. In addition, some candidates may leave the form of the signal unchanged while shifting the precession frequency or amplitude. This is intriguing, as it provides an effective contribution to the anomalous muon MDM or the muon EDM which is set by the local DM density. Such a DM MDM contribution may indeed explain the observed discrepancy between the BNL result and the SM prediction [11,12]. Finally, an ultralight DM perturbation is generally harmonic in time, resulting in a modulation of the precession signal on timescales set by the DM mass. The usual g-2 and EDM analysis is typically blind to this modulation as it averages over precession data spanning many DM modulation periods. However, the modulation may be revealed with a time-resolved reanalysis of precession data. This provides both a means of testing the background DM explanation of the muon g-2 anomaly, as well as a new opportunity for ultralight DM detection. 1 In the final stages of this work, [21] appeared which primarily considers frozen spin techniques with proton storage rings to constrain pseudoscalr DM-proton wind couplings, but also briefly considers the use of muons.
The rest of this paper is organized as follows. In Sec. II we provide an overview of muon precession experiments. In Sec. III we explore muon precession in the presence of a coherent dark matter field. In Sec. IV we describe the sensitivity of existing and upcoming muon precession experiments to characteristic DM signal shapes. In Sec. V we consider specific DM candidates and project limits. Concluding remarks are presented in Sec. VI.

II. OVERVIEW OF MUON SPIN PRECESSION EXPERIMENTS
This section will provide a criminally simplified description of the physics and techniques employed to measure the precession of muon spins. We discuss only what is necessary to reveal the implications of these measurements on DM-muon interactions. For more thorough reviews, see [23][24][25] A. Spin Tracking via Muon Decay The spin of a muon is imprinted on the angular and energy distribution of the positrons 2 produced by its decay. This is a consequence of the chiral structure of the Weak interaction. In the muon rest frame, the decay rate to positrons of energy E emitted into a solid angle dΩ alongn depends on the overlap ofn with the muon spin S: where the asymmetry factor A(E) is positive 3 at the relevant energies. The outgoing positron flux is emitted predominantly parallel to the muon spin, with the correlation becoming stronger for higher energy positrons [24]. The average spin of an ensemble of muons may thus be inferred by measuring the distribution of decay positrons. This technique is employed by the BNL, Fermilab, and J-PARC g-2 experiments. Two specific observables are measured in each experiment, a total count and a vertical count, each of which tracks a particular component of the muon spin. a. Total Count. In a lab frame the highest energy decay positrons are those emitted along the muon momentum p, so the lab frame energy may serve as a proxy for outgoing direction. As positrons are predominantly emitted parallel to the muon spin, it follows that more positrons will be produced at the highest possible energies if the muons' spin and momentum are anti-aligned than if they are aligned. The rate of positrons emitted over all directions with a lab frame energy E depends on the overlap ofŜ andp: The total count N T (t) is the number of positrons emitted above a carefully chosen energy threshold, which from Eqn. (2) has the form for an energy-dependent constant A and the dilated muon lifetime τ µ [23]. The time-evolution of N T (t) thus records the evolution of the projection of the muon spin along its momentum.
b. Vertical Count. The second observable is the difference in the number of positrons emitted with a velocity component parallel and anti-parallel to the vertical direction, defined as the direction of the experiment's large, static magnetic fieldB. From Eqn. (1), this is proportional toŜ ·B and thus probes the component of muon spin along the magnetic field. Instead of a differenatial count, an analgous quantity may be measured which is also proportional to the vertical component of the spin, such as the average vertical angle of outgoing positrons [22,26]. We will refer to this measurment generically as the 'vertical count' ∆N B (t), which has the form

B. Precession Signals
All the muon spin precession experiments we consider, observe decaying muons which are executing cyclotron orbits in a uniform, static magnetic field B. The muon spin precesses in B and any additional EM fields which are present. The experiments are designed to measure the intrinsic muon MDM and/or EDM, so we briefly describe here the expected precession signals in that case. This will elucidate the specific design and data analysis choices made in these experiments (see Sections II C and II D), as well as introduce the notions needed to derive the DM-induced precession signals in Section V.
In the lab frame, muons are held in circular orbits in a plane perpendicular to B. They orbit with the cyclotron frequency ω C , given by the vertical magnetic field B and possibly a radial electric field E [23]: Note that for radial E, ω C is parallel or anti-parallel to B. We ignore for the moment non-radial E and the possibility of muons having non-zero momentum along B, which would cause a deviation from circular orbits. It is useful to view the evolution of the muon spin in the rotating muon rest frame (RMRF). This is a noninertial frame in which the muon is at rest and the velocity of the lab always points in the same direction, which we take to be the y-direction. To reach this frame at a particular time t, we start with a Cartesian lab frame with B in the z-direction, rotate so the muon momentum is in the y-direction, and then boost alongŷ so the muon is at rest. For muons in circular, cyclotron orbits, the z-axis of the lab frame and RMRF coincide. The momentum and vertical components of S appearing in the decay counts Eqn. (3) and Eqn. (4) are respectively the y and z components of spin in the RMRF.
The muon spin S in the RMRF evolves according to a precession equation where we take t to be the lab time. The precession frequency ω a is given by three distinct contributions: ω τ is the result of the net torque on the muon spin in the RFMR, with the factor of γ due to taking the derivative with respect to lab time in Eqn. (6). In this case ω τ is due entirely to the EM fields E and B in that frame: where m µ , g µ are the muon mass and gyromagnetic, and d µ is the intrinsic muon EDM. ω T is the Thomas precession, arising from the accelerated motion of the muon. This may be computed in terms of the lab frame trajectory v(t) of the muon [27]: Finally, ω C is the cyclotron frequency Eqn. (5), which appears because the RMRF rotates at ω C relative to the lab. All of these contributions may be expressed in terms of the lab frame fields E and B, which yields where a µ = g µ /2 − 1. Note that the v · B term vanishes for circular orbits. The spin trajectory in the RMRF is uniform, circular precession with angular velocity ω a , since ω a is time-independent in that frame. We take the muon spin to be initially parallel or anti-parallel to the momentum, as is the case in the experiments considered. 4 The g-2 experiments are designed so that the first term in Eqn. (10) dominates. And with the simplifying assumption of vanishing EDM, ω a is in the z-direction so the spin precesses in the xy-plane. The vertical component is zero and the momentum component is harmonic: S z = 0 (12) where the oscillation frequency is the magnitude ω a = | ω a |. For a small but nonzero d µ , ω a is slightly tilted in the RMRF from the z-direction into the x-direction, by an angle proportional to d µ . The spin now precesses in a plane slightly tilted from the xy-plane and has a harmonic vertical component in addition the harmonic momentum component: We may therefore think of the total count Eqn. (3) as probing the precession magnitude | ω a | and the vertical count in Eqn. (4) as probing components of ω a which are perpendicular to B. Note that a nonzero EDM always increases the magnitude of ω a (see Eqn. (10)). However, from the total count alone this is indistinguishable from the muon having zero EDM and an anomalous gyromagnetic ratio instead [29]. Breaking this degeneracy is a key motivation for the vertical count [22]. The g-2 experiments allow a simultaneous measurement of a µ and d µ . However, better sensitivity to d µ can be achieved with a dedicated search. One example is the frozen spin technique, in which the experiment is designed so that all of the terms in Eqn. (10) cancel except for the d µ term. Precession is then entirely due to an EDM, and the expected trajectory is S z ≈ S sin (ω a t) .
Note that the amplitude of the vertical component is no longer suppressed by d e and now ω a ∝ d e . A measurement of the EDM can now be made by determining ω a from the vertical count.

C. Data Analysis
We consider first the analysis of the g-2 experiments. The anticipated uniform precession of Eqn. (13) and Eqn. (14) would yield total and vertical counts in the form of decaying harmonic oscillations, These signals are observed from a succession of muon bunches, with the number of bunches ranging from 10 6 to 10 8 and occurring over the course of years-long experimental run times (see Section II D). The time-series of positron counts for every individual bunch are recorded and timestamped with GPS timing [30].
The experiments seek to extract from the ensemble of single-bunch signals an estimate of ω a and d µ . Since these quantities are expected to be constant in time, a sensible technique is to align and sum the signals from each bunch, creating a stacked signal with a large signal-to-noise ratio (SNR). The alignment may be readily done with the total count, which has SNR > 1 even within each bunch [11]. But this cannot be done independently with the vertical counts, as the expected SM amplitude is much smaller than the noise. However, as the vertical count oscillation for an EDM has a fixed phase shift relative to the total count oscillation (see Eqns. (17) and (18)), the same alignment shifts used in the total count may be used to coherently sum the vertical count [22]. The two resulting stacked signals may then be fit to deduce ω a and d µ .
Stacking of the vertical counts may also be used in frozen spin experiments. In that case ω a is small, being proportional to d e , and only the leading-order behavior of Eqn. (16) is observed, S z ≈ S ω a t. Alignment is therefore not an issue, and the vertical counts may be summed and then fit for the slope ω a , which determines d e .

D. Specific g-2 and EDM Experiments
While the BNL, Fermilab, J-PARC, and frozen spin experiments all follow the general strategy outlined in Sections II A and II B, they differ in their detailed implementation. We outline here the differences which are relevant to the detection of DM precession. Unless otherwise cited, the specific values used here are taken from the experimental documentation [11,17,18,20].
a. BNL. Muons were held on their cyclotron orbits with an additional electric field E, configured as a Penning trap. This field is radial in the plane of the orbit, as in Eqn. (5), and yields a vertical restoring force above and below the orbital plane. To minimize the need to carefully measure E, the muon momentum is chosen such that the v × E term in the precession frequency ω a in Eqn. (10) vanishes. The boost factor of these muons is known as the magic gamma, γ magic ≈ 29.3. This also removes any energy-dependence from ω a , which is now determined only by the magnetic field. A field B ≈ 1.7 T was used, which yields a SM precession period 2π/ω a ≈ 4 µs. Decay positrons were collected by 24 calorimeters stations located along the inner radius of the muon orbit.
Muon precession is observed in a succession of muon bunches. Each bunch produced an oscillatory decay signal of duration 660 µs, which is roughly ten muon lifetimes at γ magic and contained about 150 spin precession periods. Each data run lasted around 5 months, observing roughly 10 6 bunches and 10 9 decay positrons in total. There runs were completed in three consecutive years, from 1999 to 2001, which measured a µ to a precision of 0.5 ppm and found a 3.3σ discrepancy from the SM prediction [12,31]. Note that this experiment directly measured ω a in Eqn. (10), and a determination of a µ requires an independent measurement of the muon mass. This was taken from measurements of the hyperfine splitting of muonium performed a few years earlier at LAMPF [32].
Three different observables were used to obtain a vertical count [22]. The least systematically difficult of these was the average outgoing angle of decay positrons relative to the orbital plane, which was measured with a tracking detector placed in front of one calorimeter station. Fewer positrons were therefore detected in this count than in the total count. This allowed a limit to be set on the muon EDM; |d µ | < 1.9 · 10 −19 e · cm. Converting this into a relative precision for measuring the perpendicular, EDM-induced component of ω a , we have δω EDM /ω a ≈ 0.5 · 10 −3 .
b. Fermilab. The Fermilab measurement is very similar to that of BNL, seeking to improve primarily by increased statistics. It employs a Penning trap electric field and uses muons at γ magic . The static field is slightly smaller, B ≈ 1.45 T. Decay positrons are counted with 24 calorimeter stations along the inner orbit radius. A vertical count is made using the average positron decay angle, obtained with two tracking detectors that have significantly increased acceptance compared to that of BNL.
The bunch duration and the number of positrons detected per bunch is similar, however the average bunch cadence is increased, allowing about 10 8 bunches and 10 11 total positrons to be observed during a roughly 5 month run. This is expected to improve the precision on a µ to 0.1 ppm. a µ will be extracted from ω a using the same LAMPF muonium measurements as BNL [32]. The enhanced tracking detection will significantly improve the measurement of the EDM, with an expected limit of |d µ | 2 · 10 −21 e · cm or δω EDM /ω a ≈ 0.5 · 10 −5 .
c. J-PARC. The J-PARC experiment will take a difference approach than BNL and Fermilab, seeking a measurement of a µ and the muon EDM with qualitatively different systematics and experimental challenges. J-PARC employs no electric field, so ω a is again set only by the magnetic field, in this case B ≈ 3 T, while allowing the use of slower muons, γ ≈ 3. The muons will be held in orbit with a weak radial magnetic field, which vanishes in the orbital plane and varies along the vertical direction. Detection for both the total and vertical count will be done with tracking detectors that record the spiral trajectory of decay positrons in the static magnetic field.
The timescales involved in this approach are naturally shorted, as slower muon have a shorter dilated lifetime. Each bunch will last around 40 µs, which is roughly 6 muon lifetimes at γ ≈ 3 and contains about 20 spin precession periods. Each bunch is expected to result in about 10 3 detected positrons, with 10 8 bunches and 10 11 positrons observed in total. The final precision is expected to be similar to that of Fermilab and BNL, 0.5 ppm on a µ and |d µ | 2 · 10 −21 e · cm. In addition, J-PARC is planning to perform new measurements of muonium spectroscopy using their muon source [33] which may be used to deduce a µ from the g-2 data.

d. Frozen Spin EDM Experiments
The frozen spin technique is newer than the g-2 approach, and a muon EDM search using these methods is still conceptual. We follow [28], which studies the possibility of using slow muons of γ ≈ 1.5 in a compact storage ring of B ≈ 1 T. An applied, radial electric field is used to cancel the precession of B, so that ω a ∝ d e . With future, high-intensity muon sources, this search can reach a sensitivity of |d µ | 10 −25 e · cm or δω EDM /ω a ≈ 10 −9 . In order to estimate the sensitivity to an oscillating DM signal, we assume that such an experiment takes data over a 3 year timespan, with each muon bunch having a duration of about 50 µs.

III. DM PERTURBED PRECESSION
In this section, we consider the evolution of muon spins in a coherent, non-relativistic DM background. We follow the muon spin in the RMRF, defined in Section II B. The most general equation of motion for the spin is a precession equation with a possibly time-dependent precession frequency: In the g-2 experiments at BNL, Fermilab, and J-PARC the SM prediction for this frequency is constant in time and given by where B is the magnitude of the lab frame magnetic field, as described in Section II B and II D. In the frozen spin proposal the SM prediction is ω a = 0. We will refer to this prediction in either case as ω sm , the SM precession frequency. DM interactions may alter ω a (t) by either perturbing the muon's orbital trajectory or by effecting the torque on the muon spin in the RMRF. In either case, the small DM perturbations may be linearized and ω a (t) may be written as where ω dm (t) is the contribution from DM-muon interactions. The DM field value will oscillate at a frequency equal to the DM particle mass m dm , and so the frequency perturbation ω dm (t) will similarly contain oscillatory components. We review here the precession trajectories that result from a perturbation with a single harmonic component of frequency m. Note that for a particular DM candidate, the frequency m of the perturbation may not be m dm but rather a multiple of m dm . The direction of ω dm (t) plays a significant role, so we consider separately parallel perturbations for which ω dm (t) = ω dm (t)ẑ and perpendicular perturbations for which ω dm (t) ·ẑ = 0.

A. Parallel Perturbations
If ω dm (t) = ω dm (t)ẑ, the precession equatioṅ may be solved exactly. The spin precesses aboutẑ with an instantaneous angular speed ω sm + ω dm (t). A spin S which is initially parallel to the momentum and perpendicular to B precesses as This may be compared to the expected SM precession with d µ = 0, given in Eqns. (11) and (12). The parallel perturbation results in a pure frequency modulation of the total count, and does not produce a signal in the vertical count. For a harmonic perturbation ω dm (t) = ω dm cos (mt + α), this has the form

B. Perpendicular Perturbations
Next we consider a perturbation to the precession frequency which is perpendicular to ω sm . For concreteness we take this to lie in the x-direction of the RMRF, ω dm (t) = ω dm (t)x, which corresponds to a precession frequency perpendicular to both B and the muon momentum, as in the case of an EDM (see Eqn (10)). 5 We focus on a quasistatic perturbation, that is ω dm (t) which varies at a characteristic rate m ω sm . This is not true in the frozen spin setup, which we consider separately in Section III C. Then the spin executes circular precession locally in time with a slowly-evolving instantaneous frequency ω smẑ + ω dm (t)x. The WKB solution to Eqn. (19) at leading order in m/ω sm and ω dm /ω sm gives: for a spin initially parallel to the momentum. This may be compared to the expected precession with d µ = 0, given in Eqns. (11) and (12). The perpendicular perturbation produces a frequency modulation in the total count which scales as ω 2 dm . This is because the oscillation of the total count is sensitive only to the magnitude of ω a (t). The perturbation also yields a non-zero vertical count, which oscillates with a fixed phase shift relative to the total count and has an amplitude modulation which is linear in ω dm . This amplitude is independent of m as it is due to the tilting of ω a (t) away fromẑ, which is set by ω dm alone -taking ω dm (t) to be static in Eqn. (28) recovers the tilted precession signal of Eqn. (12).
For a harmonic perturbation ω dm (t) = ω dm cos (mt + α), the quadratic scaling of Eqn. (27) produces both a net frequency shift and a frequency modulation at frequency 2m. The resulting spin trajectory is where:ω = ω sm + 1 4

C. Resonance and Frozen Spin
The amplitude of the vertical count in the case of a perpendicular perturbation scales as ω dm /ω sm , as in Eqn. (30). The suppression by ω sm is due to the following mechanism. The action of a perpendicular ω dm in the RMRF is to rotate the spin out of the xy-plane, and this rotation is either towards the +ẑ direction or the −ẑ direction depending on the polar angle of the spin in the xy-plane. Specifically, the spin rotates towards the direction of ω dm × S. But the dominant motion of S is rotation in the xy-plane at frequency ω sm , and so the action of ω dm is not coherent -it raises S for half of the SM period T sm and then lowers it for the next half-period. The maximal vertical component S z that may develop is limited by the SM rotation to be S ω dm T sm ∼ S ω dm /ω sm . This suppression is not fundamental. It is the by-product of an experimental design optimized for the measurement of ω sm itself and can be removed by using a different approach. There are two natural possibilities for this: the frozen spin technique and resonance. We discuss the spin trajectory in each of these cases below, focusing only on the vertical component S z as the vertical count is the most sensitive in these setups. Both techniques can achieve maximal coherence in the vertical signal, i.e. an oscillation in S z with an amplitude ∼ S. Indeed, they are conceptually the same technique as they both involve matching the frequency ω sm to m, with the distinction being whether this results in ω sm ≈ 0 or ω sm = 0.
a. Frozen Spin The frozen spin technique was invented for measuring intrinsic, static EDMs [19], and is thus most sensitive to static perturbations. In our case, this means modulation frequencies m such that m t bunch 1, where t bunch is the duration of a single muon bunch. This method engineers ω sm = 0, i.e. it freezes the spin in the xy-plane (see Section II D). Eqns. (29) and (30) are no longer valid in this regime, however the trajectory may be readily found as the total precession frequency in the RMRF varies only in magnitude, analogous to the parallel perturbation of Section III A. The spin rotates aboutx with an instantaneous angular speed ω dm (t) = ω dm cos(mt + α). This yields: where we have chosen ω dm to be alongx and the spin initially alongŷ, as in Section III B.
In the static limit, this yields a vertical signal with no amplitude suppression. Note that this is a uniform rotation over a single bunch only. For a later bunch the value of α changes and the rotation frequency may have an opposite sign. For large m the oscillation of ω dm (t) introduces a new source of decoherence. In this case the vertical signal is where we have assumed m ω dm as well, which is true in the cases we consider. The amplitude is now suppressed by ω dm /m. This is due to the fact that the spin's rotation aboutx is oscillating between clockwise and counter-clockwise motion at the DM frequency m, and after integrating this angular speed the vertical displacement of the spin scales as m −1 . This effect is analogous but physically distinct from that which produces the ω dm /ω sm scaling of Eqn. (30). If m t bunch 1, the spin is again unable to develop a large vertical component.

b. Resonance
The decoherence due to m t bunch 1 may be removed by a resonance technique, that is by engineering ω sm = m. In this case, the rotation of the spin in the xy-plane occurs at the same frequency as the oscillation of ω dm (t), and as a consequence ω dm × S does not change sign over the course of a single muon bunch. The spin will steadily rotate out of the xy-plane. Near-resonance, ω sm ≈ m, the trajectory may be found by decomposing the harmonic perturbation ω dm = ω dm cos(mt)x into two counter-rotating perturbations, one clockwise and the other counter-clockwise in the xy-plane. One of the these circular components rotates with S and dominates the dynamics. Ignoring the other component and transforming to a frame rotating at m yields a frame in which the precession frequency is constant and the spin trajectory may be easily found. Transforming back to the RMRF, the vertical component is For m = ω sm , this recovers a form similar to the static, spin frozen case of Eqn. (35). Again the vertical oscillation on-resonance is uniform over one bunch, however its amplitude will vary and may change sign between bunches. This is because at the start of a new bunch the spin is initialized to lie alongŷ, which differs from the position that a spin from the prior bunch would have if it survived until the start of the new bunch.

IV. SENSITIVITY
In this section we determine the sensitivity of existing and upcoming muon precession experiments to the generic harmonic DM perturbations given in Section III. Such a DM signal may appear in muon precession data in three distinct ways: i) A time-resolved analysis of the ensemble of single-bunch signals may directly reveal temporal variation in the muon precession frequency ω a (t).
ii) Temporal variation of ω a (t) may cause the stacked data to noticeably deviate from the expected harmonic behavior described in Section II B.
iii) The stacked data may follow the harmonic forms of Section II B within current precession, but the observed frequency or precession tilt may receive a measurable contribution which depends on the local DM density.
The first of these is the most compelling and provides an opportunity for DM detection upon reanalysis of existing and future muon precession data. The second and third allow us to set limits on DM-muon interactions using published, stacked results, while the third may also provide an explanation of the g-2 anomaly observed at BNL. A DM-muon interaction may give rise to one or more of these three signals, depending on the form of the interaction and the timescale of the perturbation, i.e. the DM mass, relative to the various experimental timescales outlined in Section II D.
We begin with the signals and constraints resulting from the total electron count, which is applicable to g-2 experiments. We then consider the vertical count, which applies to both g-2 and future frozen spin experiments, and which admits a resonant enhancement. Many of the derivations for the vertical count follow closely an analogous total count derivation, in which case only the final result is given. These results are applied to specific DM candidates in Section V.

A. Total Count
Ultralight DM may generate a frequency modulation or a frequency shift in the total count, as in Eqns. (25) and (29). We may describe both cases as a DM-induced frequency modulation of amplitude δω and frequency m in the oscillation of the momentum-component of spin S y . A static frequency shift simply corresponds to m = 0. During the i th muon bunch this has the form where α i is the phase of the DM oscillation at the start of the i th bunch. The stacked signal is where N b ≈ 10 6 − 10 8 is the number of bunches observed per experimental run. Note that δω is distinct from the DM contribution to the vector precession frequency ω dm and m is distinct from the DM particle mass m dm . δω may scale either linearly or quadratically with the magnitude | ω dm |, and m may be equal to either m dm , a non-zero multiple of m dm , or it may vanish, depending on the form of the DM-muon interaction (see Section III).

Static Frequency Shift
A DM-induced shift in the precession frequency may be directly compared with the stacked results of muon precession experiments and the predicted SM value. The current discrepancy between theory and experiment makes this comparison more intriguing. The BNL experiment has measured ω a with a precision σ ωa ≈ 0.5 · 10 −6 ω a and found a discrepancy ∆ω a between their measurement and the SM prediction of ∆ω a = 3.3 σ ωa [11]. For a DM candidate which generates a frequency shift δω, we may immediately say the following: i) If δω > ∆ω a + σ ωa , this candidate is disfavored 6 by at least 1-sigma.
ii) If δω < σ ωa , the candidate is unconstrained by this observable.
iii) If δω lies within σ ωa of ∆ω a , it provides a 1-sigma explanation of the discrepancy. iv) In the window σ ωa < δω < ∆ω a −σ ωa , a candidate cannot be said to be disfavored nor would it explain the discrepancy. Such a candidate would provide a non-negligible contribution to ω a , but additional physics would be needed to fully explain the discrepancy.
These criteria are used for the constraints given in Section V. The Fermilab and J-PARC measurements anticipate a decrease in σ ωa by a factor of 4 (see Section II D), and of course may yield a change in ∆ω a , which will necessitate a slight update to those limits.

Stacked Envelope
To what extent is a modulation with m > 0 visible in the stacked signal? Averaging a collection of nearharmonic signals with similar frequencies will generically produce another near-harmonic signal whose frequency is an average of the individual frequencies and whose amplitude is given by an envelope that evolves at a rate given by the frequency spread of the individual signals. This is the phenomenon of beats. In our case, in the limit of a large number of bunches and m t run 1, the stacked signal S y is given by the average of Eqn. (38) over the DM phase α. Here t run is the duration of a full experimental run, spanning all of the bunches in the stack. This average may be done exactly, yielding 7 whereas the expected SM signal is S y = S cos (ω sm t).
The envelope in Eqn. (41) has the form of an additional decay of the signal. Such a decay would be noticed if sufficiently strong, however there is already present in the data a systematic effect which mimics this -muons escaping the orbital trap. These muon losses are found empirically at BNL to be f loss ≈ 10% [23]. We estimate that a stacked envelope will go unnoticed if it decays by no more than a fraction f loss over the span of the stacked bunch. This bounds the argument of the Bessel-envelope in Eqn. (41) to be f loss . For simplicity, we implement this constraint as yielding an allowable DM candidate if where t bunch is the bunch duration. If the modulation does not vary appreciably over a bunch duration, this bounds the modulation amplitude in the g-2 experiments to be smaller than ∼ 10 −3 ω sm . For larger m this weakens, as the envelope decay saturates due to the decoherence between the bunches.

Stacked Frequency Residual
Supposing that Eqn. (42) is satisfied, the stacked signal S y takes the form of a harmonic oscillation. The frequency of this oscillation is approximately ω sm , but only in so far as the discrete average of the bunches approximates the continuous, single-period average over DM phase of Eqn. (40). Given Eqn. (42), the discrete average is well-approximated by This follows from linearizing Eqn. (38) in the DM-induced phase shift. We will be primarily concerned with the case mt bunch 1, where the modulation is approximately static over a single bunch. Then we have, that is, the stacked frequency is simply the average of the frequencies of each bunch. Note that α i = α 0 + mt i , where t i is the starting time of the i th bunch. In most of our regime of interest, the average time between bunches t i+1 − t i is short compared to the modulation period m −1 , so the discrete average in Eqn. (44) may be approximated by an integral , where t run is the duration of the entire data-taking run, encompassing all bunches. If (t i+1 − t i ) m 1, the value of the discrete average of frequencies depends on the uniformity of the time interval between bunches. We assume that the duration of this interval may vary by O (1) between different pairs of bunches, in which case the discrete average becomes well-approximated by a random-walk, Taking the time interval between bunches to be given on average by t run /N b , the full result is This stacked frequency shift coincides with the static m = 0 case if mt run 1, for which the shift is simply ∼ δω as in Section IV A 1. For larger m this is suppressed as the DM oscillation averages out. The suppressed shift is still constrained in the same manner as described in Section IV A 1. A DM candidate is allowed if Note that frequency residual limit in Eqn. (48) is generally less constraining than the envelope limit considered above in Eqn. (42), as the DM averaging effects appear at a much smaller value of m for the frequency residual than they do for the envelope decay. Only for m 10 −20 eV does the frequency residual give the stronger limit.

Time-Resolved Frequency Tracking
A DM modulation with m > 0 may be directly revealed by a time-resolved analysis of muon precession using each unstacked bunch. There are many specific analysis techniques that one might use, and it is beyond the scope of this work to assess them in detail. We are concerned instead with understanding the general sensitivity of the g-2 data to a DM modulation signal. For simplicity we focus on the case m t bunch 1, corresponding to m 10 −12 eV for the BNL and Fermilab experiments, for which the modulated precession frequency is constant over the duration of one bunch. The opposite limit, m t bunch 1, may be probed as well with an analysis of modulation occurring within each bunch, however we leave that case to future work.
For m t bunch 1, one may determine a local precession frequency ω(t i ) for each bunch, where t i is the start time of the i th bunch. This may be done by fitting independently the oscillations observed in each bunch. The modulated precession frequencies ω(t i ) depend on the DM field, so this is a direct measurement of a possible DM background interacting with muons. Consider the Fourier spectrumω(Ω) of the time series ω(t i ). We denote the frequency of this spectrum as Ω, to avoid confusion with the precession frequency itself ω(t i ). The zero-mode of this spectrum is non-vanishing and corresponds to ω sm . We may normalizeω as so the zero-mode is indeedω(0) ≈ ω sm . A DM-induced modulation of the form of Eqn. (38) appears in the spectrum as a peak of height δω at Ω = m. This DM signal is detectable provided δω is sufficiently large relative to the noise inω. The fit which determines ω(t i ) differs from the fit done on the stacked data, described in see Section II C, only in the number of counts and thus the SNR of the individual bunch. The precision of such a fit scales inversely with the square root of the number of counts [23], so the noise in ω(t i ) is white and has an amplitude σ i ∼ σ ωa N 1/2 b , where σ ωa is the precision of the fit to the stacked signal and N b is the number of bunches. For the Fermilab and J-PARC measurement, σ ωa ≈ 10 −7 ω sm and σ i ∼ 10 −3 ω sm . The noise in each frequency bin ofω is thus σ ωa . This is sensible, as the stacked analysis corresponds to measuring the height of the peak in the spectrum at Ω = 0. The remaining modes Ω > 0 are currently unused, but may be utilized for a DM search.
The specific frequency modes Ω i to which g-2 data is sensitive is determined by the specific timing intervals of the bunches. This is complicated by the fact that the bunches are not uniformly spaced in time, and a full analysis requires knowledge of the intervals between each bunch. This is beyond the scope of the present work. We seek an estimate of the sensitivity of such an analysis, and for our purposes we simply take the bunches to be uniformly spaced by their average spacing, t run /N b . Thenω(Ω) probes modes spaced by t −1 run with a maximum frequency of N b t −1 run . These correspond to DM masses of 10 −23 eV and 10 −15 eV, respectively. The approximation of a uniform interval between bunches has little effect onω(Ω) at small Ω, but it sets the value of the maximal frequency N b t −1 run . In a full analysis, sensitivity will extend beyond N b t −1 run as some bunches are spaced much closer together than the average spacing.
The detection reach may then be estimated as follows. The DM modulation peak has a width δΩ ≈ mv 2 dm ≈ 10 −6 m, due to the finite width of the DM velocity distribution. If mv 2 dm < t −1 run then the DM oscillation is coherent over the course of an experimental run, or equivalently the DM peak inω lies entirely within a single frequency bin. The SNR of that bin is SNR = δω/σ ωa . If mv 2 dm > t −1 run then the phase of the DM oscillation will drift during the course of a run, and the resulting peak in the spectrum will span several frequency bins. The full SNR is now properly given by the quadrature-sum of the SNR of each of those bins, which is SNR = mv 2 dm t run −1/2 δω/σ ωa . The SNR covering both regimes is We take the detection reach to be given by SNR > 3. This is properly the reach only for a predetermined frequency m, which is of interest in the event that a candidate DM signal is found in other experiments. Accounting for the look-elsewhere effect in a search with no preferred modulation frequency requires taking SNR 15, with the exact threshold depending on the desired confidence. This amounts to a sensitivity which is about a factor of ∼ 5 worse than those shown in Section V.

B. Vertical Count
A non-zero vertical count is generated only for perpendicular frequency perturbations. We consider here a harmonic DM signal of frequency m in the non-resonant case, which in the i th muon bunch is given by (see Section III B) where:ω = ω sm + 1 4 where α i is the phase of the DM oscillation at the start of the i th bunch and the stacked signal is The limits and detection reach in this case are analogous to those for the total count in Section IV A, with the distinction that in this case it is the amplitude, not the frequency, of the precession which is observed and the DM oscillation induces an amplitude modulation in the signal rather than a frequency modulation. In addition, as demonstrated in Section III, this signal is always accompanied by a static frequency shift in the total count of amplitude δω = ω 2 dm /8ω sm , which is subject to the constraints of Section IV A. That is, At its most stringent, this is ω dm 3 · 10 −3 ω sm for the g-2 experiments.

Stacked Amplitude Residual
For a perpendicular perturbation which satisfies Eqn. (55), the stacked vertical signal S z is well approximated by We have ignored the frequency modulation, as in this case it is subdominant to the amplitude modulation. The stacked amplitude is given by an average over samples of a sinusoid, analogous to the frequency residual in Eqn. (47). The typical stacked signal is thus Let σ ⊥ be the sensitivity of a static EDM search to the perpendicular component of precession frequency. For the existing BNL measurement, σ ⊥ ≈ 0.5 · 10 −3 ω sm (see Section II D). The sensitivity to the amplitude of a vertical oscillation is σ ⊥ S/ω sm and the null result of BNL implies that a DM candidate is allowed only if ω dm σ ⊥ Max (m t run , 1) (58)

Time-Resolved Amplitude Tracking
It is again possible to use a time-resolved analysis of the unstacked bunches to reveal the modulation induced by a DM background. As in Section IV A 4, we consider here the general sensitivity in the limit that m t bunch 1, where the precession is approximately uniform for the duration of each bunch.
We employ the same strategy outlined in Section IV A 4, fitting each bunch independently and then considering the Fourier spectrum of the outcome of those fits. In this case, the signal is expected to be of the form of Eqn. (51) in each bunch and the quantity of interest is the amplitude modulation. We may fit each bunch to the form for the amplitude A and construct a time series A(t i ), where t i is the start time of the i th bunch. The total count will oscillate at the same frequencyω and with a fixed phase shift relative to the vertical count (see Eqns. (29) and (30)). Thus the frequency and phase in Eqn. (59) may be determined by first fitting the higher-SNR total count, and the vertical count can be fit for only the amplitude A. Note that this is again the same procedure currently applied to the stacked signal, as described in Section II C, but now applied independently to each bunch.
We may consider the Fourier spectrumÃ(Ω) of A(t i ), normalized as: The DM modulation now appears as a peak of height ω dm /ω sm at frequency Ω = m. By an analogous argument to that given in Section IV A 4, the noise amplitude in each frequency bin ofÃ(Ω) is σ ⊥ /ω sm and the SNR of a DM modulation is For the upcoming Fermilab and J-PARC experiments, σ ⊥ ≈ 0.5 · 10 −5 ω sm . We set the threshold SNR for detection as in Section IV A 4.

Frozen Spin
For a frozen spin experiment, we consider an analogous time-resolved measurement to that of Section IV B 2. In the limit m t bunch 1, the signal has the form of Eqn. (35). ω dm is generally small, so that this is a signal which grows linearly in time, Simply averaging S z over each bunch yields a signalS z (t i ) which oscillates between bunches according to the DM phase α,S As in Section IV B 1, let σ ⊥ be the sensitivity of a spin frozen experiment to a static, perpendicular precession frequency. From the Fourier spectrum ofS z (t i )/S, the SNR of a DM modulation peak of frequency m is which follows from an analogous argument to that of Sections IV A 4 and IV B 2. For larger masses, m t bunch 1, the signal follows Eqn. (36) and the average over one bunch is suppressed: The SNR covering both regimes is and we set the threshold SNR for detection as in Section IV A 4.

Resonance
The amplitude of the vertical signal is enhanced if the DM modulation frequency m matches the SM rotation of the spin ω sm . For an experiment operating with fixed external fields and muon momentum, this results in an extended detection reach for perpendicular perturbations in a narrow frequency window around m = ω sm . In the previous and upcoming g-2 experiments, this corresponds to m ≈ 10 −10 eV. Following Eqn. (37), on resonance, the vertical spin component will grow linearly during each bunch, as the bunch duration is short compared to the on-resonance precession frequency of the spin. The angular spin velocity will vary between bunches according the DM phase, analogous to the frozen spin signal given in Eqn. (62). Following the time-resolved analysis procedure of Section IV B 3, the near-resonance SNR of this signal is This SNR is enhanced by a factor ω sm t bunch ≈ 100 relative to the non-resonant SNR of Eqn. (61). The reach is thus extended to ω dm /ω sm 10 −8 for the upcoming Fermilab and J-PARC measurements. From Eqn. (37), the frequency width of this enhancement is given by |m − ω sm | < 1/t bunch ≈ 10 −2 ω sm . This is very narrow compared to the range of m dm considered in Section V, and so we refrain from showing this peak in sensitivity in Figs. 3, 4, 5, and 6.
In addition to yielding a fixed sensitivity peak near m = ω sm in spin precession experiments, resonance may be used to extend the reach of a future DM search at a variety of frequencies by tuning ω sm to a desired search window. This would be useful for follow-up observations in the event that an ultralight DM signal is observed in other experiments. The most natural and sensitive setup for such a search is the proposed frozen spin EDM experiments, which plan to employ electric fields to tune ω sm and utilize future high-intensity muon sources (see Section II D). Then sensitivity of such a search matches that of a near-static frozen spin signal, given in Eqn. (64), as the resonant signal follows the same form as the non-resonant static signal. We show this reach in Section V for all m dm , indicating the peak reach of a narrow resonant search at the given m dm . In principle a future search may cover a wide range of m dm by systematically varying ω sm , in which case the sensitivity is as shown in Section V. There are important practical challenges to varying ω sm over a large range, which are beyond the scope of this work. The results of Section V represent the ideal limit of such an experiment.

V. CANDIDATES
In this section, we explore models of ultralight dark matter that would produce one of the signals enumerated in Sections III and IV. We consider models where the ultralight boson couples preferentially to muons so as to avoid strong tension with experiments and limits on couplings to electrons, photons, and nucleons. In the absence of a symmetry, the muon coupling will radiatively generate couplings to other SM particles. In this Section, we conservatively project only direct muon constraints and postpone a discussion of indirect constraints from radiatively generated couplings and fine-tuning, which are severe for models without a shift symmetry or gauge symmetry, to Appendix. A.
A. Scalars

φμµ
The scalar coupling we first consider is This operator has already been proposed to explain the muon g-2 anomaly (see for e.g. [14] and references therein), albeit through radiative corrections to muon g-2. This limits y 10 −3 for small enough m φ . Constraints could also be drawn from the anomalous cooling of SN1987A [34,35] owing to the presence of a non-trivial amount of muons inside. Finally, it may also result in 5th force constraints from neutron stars [36]. These, however, suffer from uncertainties in the muon abundance inside the neutron star and moreover can be avoided by introducing a quadratic coupling to nuclei, φ 2n n, which effectively screens the fifth force. There are also indirect constraints from couplings introduced at loop level which we discuss in Appendix. A.
If this scalar φ is DM, it induces an oscillating mass for the muon ω sm depends on m µ through Eqn. 10. Expanding in small y, we get, This is a parallel perturbation as discussed in Sec. III A.
Constraints and projections for this operator from different experiments are plotted in Fig. 1. The red shaded region corresponds to parameters that predict deviations not observed in the completed analysis at BNL and is ruled out at the 2σ level. At the smallest masses, the frequency shift is static as discussed in Sec. IV A 1. However, the limit is flat as it is only the change in the effective mass of the muon between the muonium experiments and the g-2 experiment which is observable here. The boundary of this region marked in green could explain the anomaly with 50% probability -it happens in the event that the scalar vev decreases in magnitude from the muonium measurement to the g-2 measurement, resulting in a lower muon mass. At scalar masses corresponding to frequencies larger than 1 year −1 , the red shaded region corresponds to deviations in muon g-2 larger in magnitude but in principle different in sign over the three different BNL runs. For this reason, the boundary is green-hatched to indicate the low probability that the three runs reported the same sign deviation. At masses larger than ∼ 10 −21 eV, there is noticeable change to the decay envelope (Sec. IV A 2). At even higher masses, coherence is lost over a bunch and only stacked frequency residuals set a limit (Sec. IV A 3). If time stamps of individual electron events are retained and used for a time-resolved analysis as described in detail in Sec. IV A 4, a projected detection reach shown by the blue line is obtained. Also shown are constraints from the virtual contribution to the g-2 measurement, cooling from SN, and 5th force constraints from NS mergers in gray. The red shaded region corresponds to deviations to the stacked analysis that would have already been seen in the g-2 analysis. The green (dashed) line corresponds to parameter space that can explain the observed g-2 anomaly with (12.5%) 50% probability. Shown in blue are projections for a time-resolved analysis. Shown in gray are constraints from virtual corrections to muon g-2 [14], SN cooling adapted from [37] and 5th force constraints from NS [36]. See Section V A 1 for details.

φ 2μ µ
In models where φ originally satisfies a Z 2 symmetry, we start with a Lagrangian, This operator is not as well constrained as the Yukawa case as the scalar appears with additional loops or in pairs and hence its effect is suppressed. Repeating the analysis above, we obtain, The constraints on this parameter space are derived similarly to the linear case and plotted in Fig. 2. Note that the constant term in Eqn. (72) does not contribute to the limits, as it is perfectly degenerate with the "intrinsic" muon mass m µ .  2: Limits and projections for a scalar DM candidate φ with coupling φ 2 Λμ µ using the same color coding discussed in Fig. 1.

∂αaμγ α γ5µ
We start with the axion-muon "wind" coupling, In a background axion field a, this interaction generates a spin torque described in the muon rest frame by the Hamiltonian term [9] H where S is the muon spin, and contributes an amount to the muon's rest-frame precession frequency. In its rest frame the muon spin precesses about the direction of the axion momentum p a , as ∇a ∼ a p a for a plane wave axion mode. In Eqn. (75), a is the axion field in the muon rest frame and the gradient is taken with respect to the rest frame coordinates. In the lab frame, the axion DM background is non-relativistic and has the form a ≈ a 0 cos (m a t) while the muon is relativistic. Thus in the muon rest frame the axion background is now relativistic and has momentum p a ≈ γm a v, where v and γ are the velocity and boost factor respectively of the muon in the lab frame. Then a ≈ a 0 cos (E t − p a · x ) in the muon rest frame, and where primes refer to muon rest frame coordinates and t is the lab frame time. This gives a perpendicular frequency perturbation via Eqn. (7), This perturbation is perfectly perpendicular as we have ignored the velocity of the axion DM in the lab frame. There is, in fact, also a parallel perturbation due to the DM velocity component along the vertical direction, however this is suppressed relative to Eqn. (77) by at least v dm ≈ 10 −3 and we may ignore it. Direct constraints on this coupling come from virtual corrections to the measured muon g-2 (this produces a wrong-sign contribution to muon g-2 and hence does not explain the anomaly), which gives Λ ≥ 1 TeV for small enough m a [35]. Constraints could also be drawn from the anomalous cooling of SN1987A [34,35,37] owing to the presence of a non-trivial amount of muons inside, yielding Λ ≥ 10 6 GeV. However there are sizable uncertainties in the muon abundance inside supernovae which translate to large uncertainties in these limits.
Constraints and projections for this operator are plotted in Fig. 3. As explained in Sec. III B, perpendicular perturbations are always accompanied by a static shift in the precession frequency which is positive definite. The green line corresponds to the parameter space that explains the anomaly and the region above marked in red would predict even larger g − 2 measurements which are disfavored. The perpendicular perturbations can also be seen in the vertical count, and the non-observation of a static EDM rules out the pink region (see Sec. IV B 1 for more detail). If a time-resolved analysis is carried out, as outlined in Sec. IV B 2, the BNL and Fermilab/J-PARC data could be used to constrain regions above the orange and blue lines respectively. Finally projections for the frozen spin method described in Sec. IV B 3 are shown in purple. Also shown are existing limits from virtual contribution to muon g-2, as well as SN cooling, that effectively rule out a DM explanation to the g-2 anomaly from this operator. However, the frozen spin method could be sensitive to new parameter space.

∂αa 2μ γ α γ5µ
We could instead consider a CP violating operator This produces a RMRF precession analogous to Eqn. (75) In the lab frame we still have a ≈ a 0 cos (m a t), so that Only the oscillatory term will contribute to Eqn. (79), as it gets a spacial gradient upon boosting to the RMRF. The frequency perturbation is Existing limits on Λ now are weaker than in the linear case. The pseudoscalar must be pair produced inside stars and it occurs in two loops in vertex corrections to muon g-2. The same set of constraints as discussed in Section V B 1 is applied to this operator and the results are plotted in Fig. 4.

aμσF γ5µ
Finally, let us consider a pseudoscalar coupling only to muons via the operator This generates a time dependent electric dipole moment for the muon given by, 10 -24 10 -20 10 -16 10 -12
The contribution to the time-dependent precession frequency can be obtained from Eqn. (10). Ignoring the electric field, which is subdominant to v × B [22], we have The DM perturbation is perpendicular to ω sm and is subject to the same limits and projections as considered in Section V B 1. These are shown in Fig. 5. The direct constraints on this operator from virtual contributions to muon g-2 are two-loop suppressed and are not shown. This model does not possess a shift symmetry and constraints from radiatively generated couplings are discussed in Appendix. A.

Lµ − Lτ
We consider an L µ − L τ gauge boson as a vector DM candidate. With gauge coupling g µ−τ , this produces a local dark electric and magnetic field with magnitudes [38] E dm = 2ρ dm cos (m dm t + α) Limits and projections for a pseudoscalar DM candidate a with coupling −i a 2Λ 2μ σ αβ γ 5 µF αβ using the same color coding discussed in Fig. 3 These fields apply both a spin torque and a force to muons, and yield a contribution to the RMRF precession frequency which has the same form as Eqn. (10) where we have ignored any intrinsic muon EDM. It is helpful to decompose B dm and E dm into components along the vertical direction B dm,z , E dm,z , and components in the plane of the muon orbit B dm,⊥ , E dm,⊥ . We consider the effects of each of these four components in turn.
i) E dm,z contributes to ω dm through the v × E dm term of Eqn. (87). This term vanishes at BNL and Fermilab due to the use γ magic (see Section II D), but would otherwise yield This is a harmonic, perpendicular perturbation which may be detected as described in Section IV B. The projected detection reach of upcoming experiments is shown in Fig. 6 for the J-PARC experiment in blue and frozen spin experiments in purple.
ii) E dm,⊥ yields a parallel perturbation if γ = γ magic , in which case its amplitude is of the same order as Eqn. (88). Since the direction of E dm,⊥ is constant in the lab frame and ω dm ∼ | v × E dm,⊥ |, the DM precession frequency now contains a product of two oscillations, one at frequency m dm and the other at the cyclotron frequency ω C . This yields two harmonic components with frequencies ω C ± m dm . ω C is much faster than the frequency ω sm at which S itself rotates, and so generically both of these components ω C ±m dm are very rapid, which further suppressed the signal, as in Eqn. (25). This suppression is removed in a narrow frequency interval around m dm ≈ ω C in which case one of the components is nearly static. This signal is not presently observable in frozen spin experiments and may only be seen in the J-PARC total count, however even in this case the signal is to weak to be observed at the projected sensitivity.
iii) B dm,z is a harmonic, parallel perturbation with which is too small to be observed by current sensitivity.. This is considerably weaker than the E dm,z effect, as it is suppressed by both v dm and a µ . iv) B dm,⊥ produces a perpendicular perturbation with an amplitude of the same order as that of B dm,z in Eqn. (89). Similar to case of E dm,⊥ , this produces perturbations which oscillate at frequencies ω C ± m dm . In this case, the two components of ω dm rotate in the RMRF. By an argument analogous to Section III C, the vertical precession amplitude is then generally suppressed by an additional factor ω sm /ω C ≈ 10 −3 which renders these perturbations unobservable with current sensitivity. This may be avoided in one of two narrow mass windows, either |ω C − m dm |/ω C 10 −3 in which case one of the components is slower then ω sm and the signal follows Eqn. (30), or |ω C − m dm − ω sm |/ω dm 1 which is the resonance regime discussed in Sections III C and IV B. We do not plot these cases as they are extremely narrow.

D. Other dark relics
The results presented thus far assume all of DM to be composed of the ultralight candidate under consideration. However, subcomponent dark matter may be easily tested as well -the limits and projections presented here may be simply rescaled in the coupling plotted on the y-axis, either linearly or as the square-root of the DM fraction, depending on the candidate. For this reason we allow the mass range in our results to extend below the existing limit on fuzzy DM from dwarf galaxies [40]. In principle, these experiments are also sensitive to background fields that redshift differently than cold DM, such as dark radiation and dark energy. We leave a careful study of these candidates for future work.

VI. CONCLUSION
We have shown that experiments designed to measure the muon g-2 and EDM are uniquely sensitive to DM models that interact predominantly with muons. DM-induced variations in the properties of muons and DMapplied spin torques and forces on muons leads to time-dependent variations in the muon precession frequencies which are measured in these experiments. While an ultralight boson making up O(1) DM was the focus of this work, subcomponent DM, dark radiation, or even dark energy could in principle be observed through these precession experiments.
Existing data from the muon g-2 experiments can be readily used to draw constraints on DM models that provide a perpendicular perturbation to the precession frequency, as these result in a net positive shift of the observed g-2 frequency. These models include the pseudoscalar wind couplings as well as pseudoscalar EDMlike couplings. Interestingly, a part of this parameter space also provides a unique explanation for the observed Frozen spin experiments have a projected detection reach shown in dark (light) purple for a static (resonant) measurement. Shown in gray are constraints from virtual corrections to muon g-2 [13] and SN cooling [39]. Existing constraints from BNL and projections for Fermilab are suppressed due to their use of γ = γ magic and are not shown. See Section V C 1 for details. muon g-2 anomaly, which is distinct from solutions that invoke radiative corrections and which typically involve larger couplings between BSM and SM. This proposition could be tested by studying timing data of electron counts in existing EDM measurements at BNL or at the currently running Fermilab experiment. Dark matter models that contribute parallel perturbations are unlikely to explain the muon g-2 anomaly, but could also be tested using timing data. Lastly, vector DM produces an electric field whose effects are suppressed at BNL and Fermilab, which employ muons at the magic momentum. This effect could instead be discerned at the J-PARC experiment or with a frozen spin measurement, which uses slower muons. The most powerful detection opportunity available in the near future is the use of a time-resolved analysis in the frozen spin experiments proposed to measure the muon EDM, either in their intended static mode or repurposed as a resonant search. Such an experiment can detect ultralight DM-muon interactions with unheralded sensitivity.

φμµ
The operators induced at 1-loop by the Yukawa operator are: L ⊃ yφ 2α 3m µ F F + y e y µ 4πē e + y n y µ 4πn n + y 2 m µ φ 2 y e y µ 4πē e + y n y µ 4πn n (A1) Here y e is the SM electron Yukawa and y N is the effective Yukawa of the nucleon. The Yukawa type couplings, to a pair of photons, electrons and nucleons induced above have limits from stellar cooling, EP tests and also from atomic clocks if φ makes up all of dark matter. These are shown in Fig. 7. The φ 2n n and φ 2ē e couplings induce a mass for the scalar in the presence of large SM number densities and can prevent the scalar from percolating into the earth. The estimate for this is, and is labeled in Fig. 7 as "shielded from φ 2n n". The Coleman Weinberg potential generates L ⊃ y 2 4π 2 Λ 2 UV φ 2 + y 3 24π 2 m µ φ 3 + y 4 24π 2 φ 4 (A3) The mass term in the CW potential tells us how tuned the scalar is and in general depends on the UV scale Λ UV . The quartic coupling generated needs to be small enough in order for φ to redshift like dark matter [8]. (A4) This is plotted as the "Quartic" line in Fig. 7. These curves together show that the new Yukawa parameter space that can be probed by muon g-2 experiments is finely tuned and clever model building has to be performed in order to explain the absence of additional operators that are severely constraining.

φ 2μ µ
This radiatively generates, L ⊃ φ 2 Λ 2α 3m µ F F + y e y µ 4πē e + y n y µ 4πn n +  Fig. 1. EP tests (in black) and DD limits (in brown) from [41] for photon and electron couplings are shown. Induced φnn can lead to shielding on earth above the relevant brown line. φ redshifts as DM only below the brown "Quartic" line.
Just like the Yukawa case, φ 2n n and φ 2ē e can prevent scalar from percolating into the earth. This is given by, Finally, depending on the details of UV physics, the EFT is safe only for field values well below the cutoff scale, i.e. φ DM ≤ Λ. These constraints are plotted in Fig. 8.
Finally, a 2n n and a 2ē e can prevent percolation into the earth. These radiatively generated couplings are: y e y µ 4πē e + y n y µ 4πn n (A12) which give a correction, Also shown is the brown line above which a 2 ≥ Λ where the EFT might not be well-defined. The region below the black tuning line corresponds to natural parameter space where the coupling is weak enough to accommodate light masses naturally.