Exploring Leptophilic Dark Matter with NA64-$\mu$

We investigate the prospects for detecting light leptophilic dark sectors with a missing-momentum experiment at NA64 running in muon mode. In particular, we consider models in which dark matter connects to the visible sector through a lepton- or muon-specific scalar mediator. These scalars can also account for the $\sim 3.5\sigma$ discrepancy between the measured and predicted values of $(g-2)_{\mu}$ . We emphasize the complementarity between NA64-$\mu$ and other terrestrial and astrophysical probes.


Introduction
The measured value of the anomalous magnetic moment of the muon, a µ ≡ (g − 2) µ /2, differs from the Standard Model (SM) prediction [1] by Here, the first and second error bars indicate the experimental and theoretical uncertainties, respectively. In terms of these uncertainties, the measured result represents a 3.5σ upward deviation from the SM prediction. This discrepancy first surfaced about 15 years ago [2] and currently remains unexplained. On-going efforts to measure (g − 2) µ more precisely at Fermilab [3,4] and J-PARC [5], along with improvements in the SM theoretical predictions from e.g. advancements in lattice quantum chromodynamics (see [6] for a recent review), may shed additional light on the discrepancy in the near future.
At the same time, astrophysical and cosmological observations have provided strong evidence for the existence of dark matter (DM). To date, its identity remains undetermined and only its gravitational interactions have been observed. Null results at dark matter direct and indirect detection experiments and collider searches targeting ∼ O(10) GeV − O(10) TeV DM masses may be pointing towards scenarios beyond the typical weakly interacting massive particle (WIMP) paradigm. Light dark matter (LDM), with masses at the GeV scale or below, has been recognized as a particularly compelling and well-motivated alternative, receiving considerable attention in the literature and motivating several dedicated experimental efforts (see [7,8] and references therein for an overview). Typically, dark matter at or below the GeV scale requires an additional light particle to mediate interactions with the SM and avoid overproduction through thermal freeze-out in the early universe [9][10][11][12][13].
Light mediators, built into thermal LDM models, can also explain the (g − 2) µ anomaly. Several of the simplest scenarios furnishing new light degrees of freedom, such as dark photons and scalars mixing with the SM Higgs, are already disfavored as explanations of the (g − 2) µ anomaly by existing measurements (see e.g. [8,14]). However, there remain several viable possibilities if the new mediator couples predominantly to leptons. These leptophilic mediators will be our focus here.
Arguably the most direct method for exploring new physics explanations of (g − 2) µ is through muon beams at accelerator facilities. Given the possible connection between dark matter and the (g − 2) µ anomaly, missing momentum searches at muon beam experiments are a particularly appealing possibility, since one might expect the mediator to decay invisibly to dark matter. It has been pointed out [15,16] that the NA64 experiment can run in muon mode (dubbed NA64-µ) and perform a muon missing momentum search with a µ + beam 1 supplied by the Super Proton Synchrotron (SPS) at CERN [17,18]. If approved, NA64-µ is expected to run after the CERN long shutdown (2021). More recently, M 3 , a compact muon missing momentum search experiment at Fermilab, has been proposed [19]. It aims at the same measurement with potentially more muons on target (MOT).
In this study, we consider leptophilic dark matter models with light scalar mediators as physics targets for these experiments, focusing in particular on the prospects for NA64-µ. We will be interested in models in which the mediator decays primarily to dark matter, and hence evades most searches for visibly-decaying light particles. In contrast to the gauged L µ − L τ scenarios often discussed in the literature, light scalar mediators need not couple to neutrinos and can feature substantial freedom in the couplings to different lepton flavors. In fact, the couplings to all but one flavor could be strongly suppressed [20]. These models thus represent a distinct class of viable but experimentally challenging explanations of both DM and the (g − 2) µ anomaly. We will show that the parameter space of these models resolving the (g −2) µ discrepancy can be conclusively probed with NA64-µ for mediator masses below ∼ 10 GeV. We will also emphasize the complementarity between NA64-µ and other existing and proposed probes of light leptophilic dark sectors with invisibly-decaying mediators. We believe our results strengthen the scientific case for experiments like NA64-µ and provide additional well-motivated physics targets that would be difficult to experimentally access otherwise.
Before proceeding, let us briefly comment on the relation of the present study to previous work appearing in the literature. The prospects for NA64-µ in exploring leptophilic vector bosons were studied in Refs. [15,16] in the context of gauged L µ − L τ models. We utilize and extend their NA64 sensitivity estimates to the models introduced below. Ref. [19] discussed M 3 projections and constraints primarily in terms of a gauged L µ − L τ dark matter model. The authors also consider detection prospects for invisibly-decaying scalars coupled to muons at M 3 , and we utilize their results when comparing NA64-µ against M 3 . Refs. [21,22] explored explanations of the (g − 2) µ discrepancy with light leptophilic scalars coupling to all lepton flavors and decaying visibly. The leptophilic model we consider resembles those considered in these studies, extended to include couplings to dark matter. Ref. [23] studied light dark matter models with scalar mediators coupled to electrons, and we incorporate some of their results in discussing the cosmological and astrophysical signatures of our models. Ref. [20] introduced muon-specific scalar models and discussed ultraviolet (UV) completions of these scenarios, again focusing on scalars that decay visibly. We make use of their results in motivating models with muon-specific mediators. Finally, Ref. [24] considered light leptophilic and muon-specific mediators and their detection at NA64-µ, focusing on the case with visible (but displaced) S decays. Our study can be viewed as an extension of this work to the case where the scalar mediates interactions of a dark sector with the SM and decays invisibly.
The remainder of this study is structured as follows. In Sec. 2, we introduce light dark matter models with leptophilic scalar mediators and discuss the associated cosmological and astrophysical consequences. Sec. 3 discusses the muon missing momentum search at NA64-µ and its prospects for exploring these models. In Sec. 4 we survey several other probes of these models and compare them to the sensitivities afforded by NA64-µ. We conclude in Sec. 5.

Model setup and general considerations
Motivated by the long standing (g −2) µ anomaly, we consider dark sectors with leptophilic interactions with the Standard Model. For concreteness, we will take the dark matter candidate to be a Dirac fermion χ, and the scalar mediator, S, to only couple to SM leptons. In UV-complete models, one might also expect S couplings to quarks. Our setup should be understood as corresponding to models in which the scalar's couplings to quarks are small relative to its couplings to leptons, so that the latter dominate the phenomenology. Given these assumptions, the effective Lagrangian governing the interactions of the scalar with the SM leptons and a Dirac fermion dark matter candidate is taken to be (2.1) We will assume that S and χ are both light, with masses below O(10) GeV, and take S to be real. The couplings of S to leptons in Eq. (2.1) violates the SU (2) L × U (1) Y gauge invariance of the Standard Model, but can be understood as originating from the effective gauge-invariant dimension-5 operators Here Λ is the associated scale of new physics and c i is a Wilson coefficient for the flavor i. We will assume that the couplings are diagonal in the mass basis. While the relative sizes of the Wilson coefficients c i (and hence the effective couplings g i ) are undetermined a priori, a natural expectation might be that they are proportional to the corresponding Yukawa coupling y i , so that the effective g i are proportional to the corresponding lepton masses after electroweak symmetry breaking. This is the case in the framework of Minimal Flavor Violation (MFV), for example. It is also possible, however, to have new physics at the scale Λ with non-minimal flavor structure. As emphasized in [20], this could give rise to couplings of S predominantly to one flavor in a technically natural way that avoids dangerous flavor changing neutral currents. In our analysis below, we will consider both the MFV-motivated lepton-specific case, with mass-proportional couplings, and the muon-specific case, in which the couplings to electrons and taus are negligible, i.e., Scalar mediator (S) = Lepton-specific scalar: g e : g µ : g τ = m e : m µ : m τ Muon-specific scalar: g µ = 0, g e = g τ = 0.
Throughout our study we will remain agnostic about the particular UV completion of the effective operators in Eq. (2.1), focusing on the model-independent constraints and prospects for observation. Possible UV completions involving lepton-or muon-specific scalars have been proposed in the literature and include scenarios with new vector-like leptons and lepton-specific two-Higgs-doublet plus singlet models (see e.g. [20][21][22]). Adding a coupling of the scalar to dark matter in most of these models is trivial. From the EFT perspective, assuming Wilson coefficients proportional to the corresponding Yukawa couplings, c ∼ O(1) × y in Eq. (2.2), new physics scales Λ 1 TeV correspond to muon couplings g µ O(10 −4 − 10 −3 ). For low values of m S , we will see that the a µ -favored region falls in this regime, as can the thermal relic target for leptophilic dark matter. If the new physics responsible for generating the operators in Eq. (2.2) does not involve new colored states, the LHC is unlikely to constrain the corresponding UV completions for Λ near the TeV scale. Couplings g µ O(10 −3 ) correspond to O(100 GeV) scales of new physics (assuming c ∼ O(1)×y ), and so the EFT can break down at LHC energies and a UV completion should be specified when considering constraints at highenergy experiments. Nevertheless, since several UV completions have been shown to remain viable while generating couplings in this range [20][21][22], and the Wilson coefficients need not be proportional to the corresponding Yukawa coupling, we content ourselves with the UV-agnostic treatment in what follows. Note also that, as discussed in Ref. [20], the scenarios we consider can be consistent with technical naturalness provided new physics enters to cut off the quadratically-divergent contributions to m 2 S around the TeV scale. The S-lepton couplings in Eq. (2.1) also introduce effective scalar couplings to pairs of vector bosons (SV V ) at one-loop. For a light scalar with mass below the electroweak (EW) scale, the most relevant SV V coupling is the scalar-diphoton coupling (Sγγ). 2 This interaction is important when the decays of S → + − are kinematically forbidden and is accounted in our thermal relic density calculations below. Additionally, the Sγγ coupling can give rise to experimental signatures involving photons, also discussed below. We can parametrize the corresponding interactions with an effective Lagrangian, where α is the electromagnetic fine structure constant and F 1/2 is a form factor that depends on the four-momentum-square of one of the photons (q 2 ) and the scalar S (p 2 S ). In computing F 1/2 we will consider one of the photons to be on-shell but allow S and the other photon (with four-momentum q) to be off-shell. Expressions for F 1/2 are given in Appendix A. Combining Eq. (2.1) and Eq. (2.4), the total decay width for S is given by (2.5) 2 Couplings like g Zγ can be also relevant. We discuss the corresponding signatures in Sec. 4.
We will primarily focus on cases where g χ g , so that S decays primarily to χχ. For the decay of S to two photons, all three particles are on-shell, corresponding to q 2 = 0 and p 2 S = m 2 S in F 1/2 . As emphasized in Sec. 1, besides mediating interactions between the visible and hidden sectors, the scalar S can contribute to (g − 2) µ through its couplings to muons [25]: This contribution is positive, and can raise the predicted value of a µ so that it agrees with experiment, cf. Eq. (1.1). In our discussion of the model parameter space below, we will indicate the regions consistent with the central a µ value within 2σ, as well as regions expected to be favored by future a µ measurements, assuming that the central value of ∆a µ will remain unchanged while the experimental and theoretical uncertainties will be improved by a factor of 4 and 2, respectively [1,[3][4][5]26]. Turning our attention to the DM, there are three distinct possibilities for the relative sizes of m χ and m S that carry different phenomenological consequences: • m χ < m S /2: in this case, for g χ g e,µ,τ , the mediator S will primarily decay invisibly to χχ.
The thermal freeze-out relic abundance of DM is driven by s-channel annihilation into leptons, χχ → (or γγ) in the Early Universe. The annihilation rate is roughly given by and depends on both the dark sector coupling g χ and the visible sector coupling(s) g (a sum over lepton flavors is implicit above). v 2 rel is the thermal average of the relative DM velocity squared, and its presence above reflects the fact that annihilation is a p-wave process in this scenario. Since the annihilation rate depends on g , this scenario provides a well-defined thermal dark matter target that can be searched for in terrestrial experiments.
• m S /2 m χ m S : here again s-channel annihilation into leptons sets the relic abundance of χ, providing a thermal relic target. However, S → χχ decays are kinematically forbidden, and so S will decay visibly. This dramatically changes the constraints and prospects for detection at accelerator experiments, and has been discussed in detail elsewhere in the literature (see e.g. [21,22]). In particular, the NA64-µ projections of interest below can be taken from Ref. [24].
• m S < m χ : in this case, again S decays visibly. However, annihilation in the early Universe will primarily proceed through secluded annihilation, χχ → SS. The cross-section for this process only depends on the dark sector coupling g χ , and so there is no well-defined thermal relic target for terrestrial experiments. Nevertheless, this is a viable possibility, and again the prospects for discovery can be inferred from e.g. [21,24].
In what follows, we will focus on the first of these scenarios, with m χ < m S /2, since it provides concrete thermal targets and is generally the most difficult to test, given the invisible decays of S.

Relic abundance via thermal freeze-out
To begin our investigation of the parameter space of these models, we first determine the regions consistent with the observed relic abundance of dark matter. We assume the dark sector is in thermal equilibrium with the SM plasma in the Early Universe [23]. At later times, the annihilation rate for χχ → (or γγ) drops below the expansion rate of the Universe and the dark matter abundance freezes out. To accurately compute the resulting relic abundance, Ω χ h 2 , we use the MadDM 2.1 [27] incorporating a UFO-format model built in FeynRules 2.3.27 [28]. For concreteness, we compute Ω χ h 2 for two specific DM-mediator mass ratios, where we assume R 1, adopting the notation from [29]. The second mass ratio is chosen to illustrate the effects of the resonant enhancement present for the s-channel annihilation process when m χ ∼ m S /2. As emphasized in Ref. [29], this enhancement can dramatically affect the thermal target parameter space, allowing for smaller couplings consistent with the observed dark matter density, Ω DM h 2 ≈ 0.12 [30]. Results for the mass relations in Eq. (2.8) (Eq. (2.9)) in the lepton-specific and muon-specific models are shown in the left and right panels of Fig. 5 (Fig. 6), respectively. Regions below the red contours are excluded due to an over-abundance of DM for the indicated mass ratio. "Kinks" around m S = 2m µ (both left and right panels) and m S = 2m τ (right panel) occur when a new annihilation channel becomes kinematically accessible as indicated by Eq. (2.7). The contrast between the relic density-compatible regions in Fig. 5 and 6 illustrates the model-dependence of the thermal freeze-out constraint. Figs. 5 and 6 also show the corresponding a µ -favored regions in the m S − g µ parameter space. Comparing the relic density curves with the green bands suggests that, given a suitable choice for m S /m χ , a large portion of the parameter space favored by the (g − 2) µ measurement can also furnish a potentially viable dark matter candidate (which may or may not saturate the entire observed abundance).

Constraints on g χ from dark matter self-interactions and perturbativity
The presence of the scalar mediator also introduces a Yukawa-like attractive self-interaction between DM particles [31]. Such interactions can be constrained by astrophysical observations. Taking m χ < m S /2 and given that v rel O(10 −3 ) in most astrophysical systems, the resulting momentum-transfer self-interaction cross-section strength 3 is given by (2.10) if the interaction falls in the perturbative Born regime, α χ ≡ g 2 χ /4π m S /m χ (as well as m χ v rel /m S 1).
Given that Eq. (2.10) is velocity-independent, self-interaction constraints from halo systems at different scales should be taken into account simultaneously. On the dwarf galaxy scale, σ/m O(10 cm 2 /g) is allowed since gravothermal collapse of the halo is avoided [33,34]. For Milky Waysized galaxies, σ/m 1 cm 2 /g is allowed by halo morphism [32]. At the galaxy cluster scale, σ/m 0.7 − 7 cm 2 /g is allowed by halo mergers (see e.g. [35]). Given these considerations, we choose 1 cm 2 /g as a robust upper limit on σ T,χχ /m χ . The resulting constraint is shown in Fig. 1. As expected from Eq. (2.10), the constraint on g χ weakly depends on m χ (or the ratio of m χ /m S ). We see that g χ = 1 (0.1) is disfavored for m S < 20 MeV (1 MeV) and not constrained for larger scalar masses.
Large values of g χ can also lead to a breakdown of perturbation theory at low scales, rendering perturbative results invalid. In Fig. 1, we also indicate regions for which α χ ≥ 4π at the TeV scale

Resonant regime
Born regime E x c l u d e d b y from DM self-interactions and perturbativity requirements. The orange region is disfavored by astrophysical constraints such as halo morphism (see discussions in the text). The purple region indicates the resonant regime (α χ ≡ g 2 χ /4π m S /m χ and m χ v rel /m S 1) where Eq. (2.10) is not applicable. Above the gray dash-dotted line, α χ runs to a value ≥ 4π at the TeV scale or below. or below. These results make use of the one-loop beta function for a real scalar coupled to a Dirac fermion [36,37] where µ is the renormalization scale. Motivated by the self-interaction and perturbativity constraints in Fig. 1, we take g χ = 1 and restrict the scalar mass to be above 20 MeV in what follows.

Other cosmological and astrophysical constraints
There are other astrophysical and cosmological constraints on light leptophilic dark sectors. Observations of the cosmic microwave background (CMB) constrain the amount of energy injected into the Intergalactic Medium (IGM) through dark matter annihilation at late times. This energy injection can distort the CMB. However, this is only an issue if the DM annihilation cross-section is s-wave [30,38,39]. Scenarios with a light scalar mediator and Dirac fermion dark matter feature a p-wave annihilation cross-section and therefore the CMB constraint is not relevant for the scenarios we consider. Light dark matter can also affect the successful predictions of big bang nucleosynthesis (BBN). If dark matter is sufficiently light, it is expected to have been in equilibrium with the SM thermal bath, and hence relativistic, until after the onset of BBN. This affects the Hubble rate, and thus the abundances predicted by BBN, which are tightly constrained. For dark matter masses above about 1 MeV, this is not an issue, as freeze-out occurs sufficiently early, and thus this constraint does not significantly impact the scenarios we consider. Dark scalars with mass below O(100 MeV) can also be abundantly produced in a core-collapsed supernova (SN) through resonant or continuum production [20,23,40]. Once produced, S can provide a new cooling channel for the supernova as it streams out from the core, and could conflict with observations of SN 1987A. We will discuss the corresponding SN 1987A constraint in Sec. 4.2. Figure 2: Diagrams contributing to the signal in a muon missing momentum search for dark scalar bremsstrahlung (a) with the scalar decaying invisibly and (b) the scalar decaying visibly but outside the detector. Given the NA64-µ setup and the parameter space we interested in, the missing momentum signals in our models are dominated by (a).

Missing momentum measurements at NA64-µ
We have argued that light leptophilic dark sectors are compelling from the standpoint of explaining both the observed DM abundance and the (g −2) µ discrepancy. How might one explore these scenarios experimentally? The NA64-µ experiment was proposed in [15,41] and is currently planned to run after 2021. It was originally envisioned to search for new vector gauge bosons solving the (g − 2) µ puzzle. In this section, we study the possibility of using this experiment to search for light dark matter produced through the decay of a scalar mediator in the class of models motivated and described above. We show that NA64-µ will provide impressive experimental coverage of these scenarios that are otherwise difficult to probe.

Production of light dark matter
At muon-beam experiments, dark scalars, S, can be produced as initial and finial state radiation off of the incident muons when they scatter with the target nuclei (µ + N → µ + N S). Through this mechanism, and depending on the mass of S and its couplings to dark matter and the SM leptons, a missing momentum measurement at NA64-µ would be capable of probing both (a) invisible decays of S and (b) visible decays of S outside of the detector, illustrated in Fig. 2. In both cases, dark scalar bremsstrahlung can induce a significant amount of missing energy/missing momentum and leave a feebly scattered muon as a distinct signal. 4 In the scenarios of interest, the expected number of events corresponding to case (a) is much greater than that from case (b) for the following reasons: (1) large values of g χ are favored from the standpoint of obtaining the correct relic density without fine-tuning, while large values of g are significantly constrained for the parameter space we are interested in. Thus, we expect BR(S → χχ) > BR(S → ¯ ). (2) Case (b) surfers a exponential decay volume suppression in order for S to decay outside the detector. There is no such suppression for case (a). A detailed calculation of the dark scalar bremsstrahlung (µ + N → µ + N S) production rate is discussed in [24]. We will briefly review it here and adopt it to estimate the number of missing momentum signal events expected in our models. The calculation is based on the improved Weizsacker-Williams (IWW) approximation [42], which treats the exchanged virtual photon between the muon and nucleon as a real photon. It is a good approximation when the beam energy is much larger than the momentum transfer [43]. The total number of missing momentum events, N χ , can be approximated as assuming that dE µ /dy, the change of the muon energy E µ with respect to the penetration length y, is approximately constant. For the lead (Pb) target used in NA64-µ, this is a good approximation for the relevant energy range (∼ 100 GeV), with dE µ /dy ≈ 12.7 × 10 −3 GeV/cm [44]. In Eq. (3.1), N µ is the total number of incident muons and n atom is the atomic number density of the target (n atom = 3.3 × 10 22 / cm 3 for lead). The integration is over the penetrating muon energy E µ and the bremsstrahlung scalar energy E S with respect to E µ (x ≡ E S /E µ ). The lower integration limit for E µ , E µ,min = E µ,beam − L tg dE µ /dy , is set by the energy loss of a positive muon after passing through the entire target of length L tg (in the projection below, we use a thin target with L tg = 20 cm). E µ,beam , the initial muon beam energy, is 150 GeV. The lower integration limit on x, x min , is set by requirements for background rejection. Here, we choose x min = 1/3 as suggested by [15]. This amounts to requiring signal events to have missing energy larger than E µ /3 ∼ 50 GeV.
In most of the parameter space under consideration, BR(S → χχ) appearing in Eq. (3.1) is very close to one. Nevertheless, it can decrease substantially when visible decays are enhanced by either a large g or large phase space factors, and we account for this in our sensitivity estimates.
The differential signal production cross-section appearing in Eq. (3.1) is given by 5 : where the boost factors are β µ = 1 − m 2 µ /E 2 µ ≈ 1 and β S = 1 − m 2 S /(E S ) 2 for muons and S, respectively. The effective photon flux, χ, is given by where the virtuality t represents the momentum transfer squared and G 2 is the combined atomic and nuclear form factor. An explicit expression for G 2 can be found in Appendix A of Ref. [24]. Finally, Eq. (3.1) assumes ∼100% trigger and reconstruction efficiencies. The reader should bear this is mind in what follows.

Experimental setup
With the signal rate computed, we can consider the sensitivity achievable by the NA64-µ experimental setup. NA64-µ is equipped with a high energy and high intensity muon-beam from the CERN SPS [17,18]. The muon beam has a maximum momentum in the range of 100 − 225 GeV. A typical intensity of 2 × 10 8 µ + per spill can be achieved for beam energies between 160 and 190 GeV. The period of a SPS cycle is around 15 seconds, which includes the spill with a duration of 4.8 seconds. The full experimental setup of NA64-µ is detailed in [15]. We highlight relevant aspects of the detector segment in Fig. 3.
A calibrated muon beam with energy 150 GeV is injected into an active target with a length of around 20 cm (the length of the surrounding ECAL). The momentum and energy of the outgoing  Figure 3: The detector segment of the proposed NA64-µ experimental setup. See text for more details. The full configuration can be found in [15,45]. scattered muons are measured by a set of trackers. Ref. [15] proposed to use a set of eight straw-tube chambers with a momentum resolution of σ(p µ )/p µ = 3% (for muon momentum p µ = 100 GeV) and 1 mm length resolution. Subsequent studies of the experimental setup [45] investigated alternative options, such as incorporating micromegas chambers or silicon based trackers, to further improve the momentum resolution 6 . Photons and other secondary particles at large angles, generated during scattering events, are rejected by the electromagnetic calorimeter (ECAL) surrounding the target and two veto counters downstream. Other secondary particles at small angles are detected in the four hermetic hadronic calorimeter (HCAL) modules.
A missing momentum or energy signal consists of a single scattered muon with energy ≤ 100 GeV with no accompanying energy deposition in the vetoes and a small amount of energy deposition in the ECAL and HCAL (E ECAL + E HCAL ≤ 12 GeV). Ref. [15] performed a detector Monte Carlo (MC) simulation given the above signal criteria. It was found that the dominant background arises from muon trident events, µ + N → µ + N (µ + µ − ) (see Fig. 4) at an expected level of 10 −12 events per MOT. Other backgrounds are subdominant and are expected to be at the level of 10 −13 events per MOT.
The muon trident background can be challenging to eliminate when the momentum of the µ + of the µ + µ − pair is much larger than the momentum of the µ − in Fig. 4. It can fake a signal event if the soft µ − is missed in the detector and the hard µ + is so collinear with the scattered µ + that they only produce a single track along the central region of the HCAL. Recently, Ref. [19] proposed a "1 vs 2" method to further reject this background. The authors point out that the two collinear µ + s in the muon trident background are both minimum-ionizing particles (MIPs). Thus when passing through a layer of the HCAL, a background event would deposit roughly twice the energy and produce twice the number of photoelectrons as compared to a genuine signal event. Ref. [19] suggested that the fake rate can be suppressed by a factor of 10 −4 if the number of photoelectrons produced by a MIP is on the order of 100. While this "1 vs 2" method was originally proposed in the context of M 3 , we suspect that a similar strategy can be incorporated at NA64-µ, with a MIP producing 150 − 200 photoelectrons when passing through a single HCAL module [15]. Therefore, using this method, a significant suppression of the muon trident background may be achievable at NA64-µ.
Given the considerations above, we suggest two background-free scenarios for a muon missing momentum search at NA64-µ: (1) MOT = 10 12 . This number is based on the background analysis Figure 4: Feynman diagrams contributing to the SM muon trident processes.
in [15]. Given the high intensity of the muon beam, it can be achieved in a one-day run . (2) MOT = 10 13 . This would be a viable search option if one adopts the "1 vs 2" method and successfully reduces the muon trident background by at least an order of magnitude. This corresponds to about a nine-day run of the muon beam. For a given number of MOT, we estimate the 95% confidence level (C.L.) sensitivity by requiring ≥ 3 events, given our assumptions above. Longer run time with better background understanding can of course achieve even better sensitivity.

Projections
The projected NA64-µ sensitivity to the models presented in Sec. 2 is shown in Figs. 5 and 6. The projections are shown both for 10 12 (solid) and 10 13 (dashed) MOT. NA64-µ has the potential to probe the a µ -favored region up to m S ∼ 10 GeV, due to the high beam energy. The sensitivity can be comparable to that of M 3 for m S < O(10 MeV), provided 10 13 MOT can be achieved while running background-free. Figs. 5-6 also show that NA64-µ can probe a considerable portion of the parameter space for which the relic abundance constraint for light leptophilic DM is satisfied via thermal freeze-out and without significant fine-tuning of m S /m χ .

Complementary probes
Given the sensitivities achievable by NA64-µ, it is natural to ask: what is the extent to which other current or planned experiments can explore the same models? We address this question here, showing that an NA64-µ-type experiment is expected to probe significantly more of the model parameter space explaining the (g − 2) µ discrepancy and consistent with thermal relic dark matter than any other experiment. In muon-specific models, other relevant experimental probes include measurements of (g − 2) µ itself and supernova cooling bounds. For lepton-specific mediators with couplings to electrons and taus, there are several additional accelerator probes that can be relevant, including future searches at B-factories like Belle-II and proposed lepton collider "Z-factories". We discuss these various complementary probes in turn.

Anomalous magnetic moments
Current measurements of a µ offer one of the most sensitive bounds on light scalars in muon-specific or lepton-specific models. In Figs. 5 and 6, the excluded parameter space for which |a NP µ −∆a central µ | > 5σ are shaded gray. Here, a NP µ is the new physics contribution to a µ from S and ∆a central µ represents the central value of ∆a µ in Eq. (1.1). The excluded regions correspond to g µ 10 −3 at m S ∼ O(10) MeV and g µ few × 10 −2 at m S ∼ 10 GeV.
Measurements of the anomalous magnetic moment of the electron, a e , can also constrain the lepton specific model in principle. However the resulting sensitivity is significantly weaker due to the mass- The left and right panels correspond to the muon-specific (g µ = 0, g e = g τ = 0) and lepton-specific (g e : g µ : g τ = m e : m µ : m τ ) models, respectively. For the two models, we assume m χ = m S /3 and g χ = 1. The solid and dashed purple lines represent the expected NA64-µ sensitivity with 10 12 and 10 13 MOT, respectively. The red curve indicates the parameters required to reach the correct thermal freeze-out DM abundance. The green and gray regions represent the 2σ-favored and 5σ-excluded regions based on current a µ measurements. The hashed regions indicate the projected sensitivity of future a µ measurements, assuming the current central value stays the same while the experimental and theoretical uncertainties will be improved by a factor of 4 and 2, respectively. The blue region with dotted boundaries represents approximate exclusions from the cooling of SN 1987A. In both panels, we include the sensitivity projections for M 3 [19] and a Belle-II mono-photon search (with 50 fb −1 data) [46]. For the lepton-specific model, we also include projections for LDMX (with 16 GeV electron beam) [47], BDX [8], mono-photon searches at BaBar [46] and Tera-Z [48], current LDM direct detection limits from XENON 10/100, DarkSide-50, and CDMS HVeV (collectively denoted as "current DMDD"), expected DMDD sensitivities at SENSEI, DAMIC-1K, and SuperCDMS-G2+, as well as constraints from exotic Z and τ decays. See text for further details.
proportional coupling hierarchy and is not competitive with other probes of the relevant parameter space in our models.

Cooling of SN 1987A
A potentially important constraint on light scalars coupled to leptons arises from supernova cooling. A core-collapsed supernova behaves like a proto-neutron star. Its core region consists of highly-degenerate and relativistic electrons, near-degenerate and non-relativistic nucleons, and perhaps some amount of muons [49]. It cools mainly through neutrino diffusion. The measured SN 1987A neutrino burst flux agrees with SN model predictions [50] and can be used to constrain dark scalars and leptophilic dark matter produced through its prompt decays 7 . Similar to the SN studies for dark photons, if the DMlepton coupling g is too small, no significant dark scalar population can be produced in the SN and hence S does not contribute significantly to cooling. On the other hand, if g is too large, the DM produced through S decays will be trapped inside the SN and, due to its frequent interactions with the SM plasma, again will not contribute to the cooling. This implies a window of couplings (and masses) that can be constrained by supernovae.
Dark scalars with mass less than the plasma frequency of the photon, ω p ( 20 MeV for SN 1987A), can be resonantly produced through mixing with the longitudinal mode of the photon [40]. The resulting energy loss per unit mass is constrained by the Raffelt bound, Raffelt = 10 19 erg g −1 s −1 . Requiring the energy loss to be smaller than this value yields an almost flat upper limit in the m S − g µ plane (the lower edge of the excluded region) around g e ≈ 10 −10 for m S 20 MeV [23]. The bound on g e can be translated into g µ = g e m µ /m e ≈ 10 −8 for the lepton-specific model. If the dark scalar mass is much greater than the plasma frequency (e.g. m S ∼ 100 MeV), resonant production is suppressed but the dark scalars can still be produced through continuum production. For example, the dark scalar can be produced through a Primakoff-like process γN → N S [20]. Ref. [46] considered this production mechanism for a pseudoscalar with photon coupling (−g γγ /4)aFF and obtained an almost flat upper limit around g γγ ≈ 6 × 10 −9 GeV −1 for m S 100 MeV. Ignoring the difference in CP quantum numbers and applying Eq. (2.4) 8 , the bound can be translated into g µ ≈ few × 10 −7 for the lepton-specific and the muon-specific models. Considering both the resonant and continuum production mechanisms, we conclude that an upper limit on g µ exists around O(10 −8 − 10 −7 ) for scalar masses up to at least 100 MeV for both models. Of course, a much more careful analysis is needed to combine all the possible production mechanisms and to properly account for plasma effects.
Once dark scalars are produced in the SN, they promptly decay into χχ given our parameter choice m S > 2m χ and g g χ = 1. Along their path of escape, DM interacts with the particles in the plasma, such as electrons, photons, and protons. Frequent interactions limit the DM outflow, and hence yield a lower limit (upper edge of the excluded region) on g for scalars light enough to be produced significantly in the SN. For concreteness, we adopt the "π/2-deflection criterion" proposed by [51]: trapping is assumed to be sufficiently efficient provided that the expected accumulative deflection angle of DM particles is ∼ π/2 along their path starting from the kinetic decoupling radius to the neutrino-gain radius ( 100 km). Here we consider the trapping induced by χe − → χe − interactions to set a preliminary constraint. The corresponding trapping bound can be evaluated via Eqs. (2.11-2.14) and Eqs. (C1-C3) of [51]. Note that Eqs. (C1-C3) need to be modified to account for the change from vector-mediated DM-nucleon interactions to scalar-mediated DM-electron interactions. The modification includes changing the mass, the number density, and the deflection angle per collision of nucleons to those of electrons. Also, one must include the correct electron distribution function in the phase-space integrals of Eqs. (C2-C3) to account for Fermi blocking.
For the muon-specific model, g γγ -induced processes such as χγ → χγ, χN → χN γ, and χe − → χe − γ are likely to provide efficient trapping. In this case, one needs to account for different polarizations of the plasma photons. A detailed analysis of these trapping mechanisms is beyond the scope of this work. However, given that g γγ is loop-and α-suppressed with respect to g e , we expect that the upper edge of the excluded region for the muon-specific model in the m S − g µ parameter space will arise at larger couplings than in the lepton-specific model. To give a rough estimate of the trapping bound in this case, we consider only the χγ → χγ process neglecting plasma effects, and adopt a simple mean-free-path criterion, requiring (n γ σ χγ→χγ ) −1 r core ≈ 1 km.
We show the resulting approximate SN 1987A bound for the lepton-specific model as a shaded blue region with dotted boundaries in Fig. 5 and Fig. 6. The lower edge extends down to g µ ∼ O(10 −8 − 10 −7 ) (below the plot range) and we cut the right edge off at 100 MeV for a conservative estimate. The upper edge is around g µ 10 −4 , increasing slightly for larger m S given that heavier dark matter (m χ ∝ m S ) is more difficult to trap. A similar exclusion region is shown for the muonspecific model. Here the upper limit is at O(10 −3 ). We do not include a dotted boundary line here, given that our estimate excludes potentially important plasma effects and other trapping processes. We again emphasize that our limits are preliminary, and defer a more comprehensive analysis for the two models to [52]. Nevertheless, the general message is clear: for light dark scalars with mass up to ∼ 100 MeV, SN probes are complementary to muon missing momentum searches.

Dark Matter Direct detection (DMDD)
In the leptophilic scenario described above, the scalar mediates scattering between χ and electrons. 9 This can lead to a signal in direct detection experiments sensitive to electron recoils. To determine the sensitivity of such experiments to the models of interest, we follow [8,53,54] and define a reference χ − e − scattering cross-section, σ e with momentum exchange set to q = αm e , and a form factor, F DM (q) to account for the q-dependence. For the mediator masses of interest, m S αm e (the where µ χe is the dark matter -electron reduced mass. In the parameter space of interest, µ χe ≈ m e and the direct detection limits are generally weakly sensitive to the DM-scalar mass ratio. Nevertheless, we show the predicted value of σ e for the leptophilic DM parameter space saturating the observed dark matter relic abundance with g χ = 1 and m χ = m S /3, m χ = m S /(2 √ 1 + R ) with R = 0.1 in the left-and right-hand panels of Fig. 7, respectively. We also show the corresponding σ e values for the a µ -favored regions.
The XENON10 [55][56][57], XENON100 [56,58], DarkSide-50 [59], and CDMS HVeV [60] experiments already set limits on σ e in the relevant mass range. Excluded values are shaded in Fig. 7. These experiments are not currently sensitive to the leptophilic parameter space saturating the relic abundance or explaining the (g − 2) µ discrepancy for our choices of parameters. There are, however, many proposed and planned direct detection experiments that are expected to substantially improve the reach in σ e for light dark matter (see [8] for an overview of these efforts). In Fig. 7, we show projections for three such experiments: SENSEI with an O(100 gram) detector, DAMIC-1K, and SuperCDMS; all projections are taken from [8]. There are other experimental proposals potentially sensitive to the relevant parameter space, but the projections shown demonstrate the range of cross-sections expected to be reached by next-generation experiments.
From Fig. 7, we see that next-generation direct detection experiments can probe the leptophilic DM parameter space with a lepton-specific scalar mediator consistent with the measured value of a µ for m S 100 MeV (assuming g e /g µ ∼ m e /m µ ). For the benchmark case m χ = m S /3, g χ = 1, this Figure 8: Diagrams sensitive to g µ or the g µ -induced g γγ coupling in the muon-specific model and contributing to (a) e + e − (qq) → γS, region is inconsistent with the observed relic abundance. However, this is not necessarily the case for m χ ∼ m S /2 (near the s−channel resonance), in which both future direct detection experiments and muon beam experiments such as NA64-µ can be sensitive to the a µ -favored regions, as illustrated by the results for the near-resonant benchmark point in Fig. 7. This opens up the exciting possibility of detecting both dark matter and its mediator to the SM within the next generation of experiments, as well as resolving the (g − 2) µ puzzle. Of course in the muon-specific case where the scalar does not couple to electrons, direct detection experiments will not be sensitive to dark matter interactions with the SM, and instead a muon missing momentum experiment will be crucial for the discovery of the dark sector.

Collider experiments
There are several collider probes of light leptophilic invisibly-decaying scalars. e + e − experiments can be sensitive to S through the process e + e − → γS(→ χχ) in mono-photon searches. In the muonspecific model, S couples to photons through a muon-loop (cf. Eq. (2.4)), giving rise to a signal in this channel through the diagram in Fig. 8a. In the lepton-specific scenario, the Sγγ coupling receives contributions from all three lepton flavors. The coupling of S to electrons results in additional t/u-channel diagram contributions to e + e − → γS(→ χχ). Across the parameter space of interest, however, we find that the Sγγ constribution represented by Fig. 8a dominates the production rate, even in the lepton-specific model. To place limits on the m S − g µ parameter space, we therefore utilize the 90% C.L. upper bounds and projections from BaBar and Belle-II in Ref. [46]. In particular, we compute the e + e − → γS(→ χχ) production cross-section, utilizing the form factor in Appendix A with q 2 = s, the center-of-mass energy for BaBar or Belle-II, and compare to the cross-section for axion-like pseudoscalars in Ref. [46]. We find that the bounds in our scalar case are weakened by a factor of ∼ √ 1.8 relative to the pseudoscalar bounds in Ref. [46] across the angular regions accessible by BaBar. Assuming a similar angular acceptance at Belle-II, the same factor of √ 1.8 applies 10 . We 10 If, for example, all angles were accepted, the factor of √ 1.8 would reduce to a factor of √ 1.5 and not significantly impact our results.
Decay process 95% C.L. upper limit on BR Table 1: Summary of rare decays that can probe g µ . The 95% C.L. upper limits on the decay branching fraction are calculated based on data in [1]. X represents weakly interacting particles.
therefore interpret the BaBar bounds and Belle-II projections on g aγγ from Ref. [46] as bounds on g γγ , including this factor of √ 1.8. In addition, we scale up the coupling by 1.96/1.64 to estimate 95% C.L. limit based on the 90% C.L. limit.
The resulting limits and projections are shown in Figs. 5-6. We find that mono-photon search at BaBar is not sensitive to any of the parameter space not already excluded by a µ measurements in either the muon-or lepton-specific models. In the lepton-specific case, mono-photon search at Belle-II with 50 ab −1 will be able to probe the a µ -favored regions of the model with m S ∼ 100 − 1000 MeV. Note that this measurment is an indirect probe of the couplings to leptons, and an experiment like NA64-µ would still be required to conclusively determine whether or not S couples to muons.
Other collider searches yield less sensitivity. Experiments such as BaBar and LHCb can probe the S-muon coupling through the process e + e − (qq) → µ + µ − S. BaBar has searched for e + e − → µ + µ − S(→ µ + µ − ) [61] and placed constraints on visibly-decaying mediators. However, we are not aware of analogous searches for µ + µ − + / E. This likely has to do with the fact that, in the visible case, the µ + µ − final state can be reconstructed as a narrow resonance that provides powerful discrimination from the background, whereas in the invisible case, no such discrimination is possible. The same reasoning likely applies to LHCb as well. KLOE [62] has performed a µ + µ − + / E search for e + e − → A S → (µ + µ − )(χχ) where a bump hunt in m µ + µ − can be performed. Given that S decays primarily invisibly in the parameter space we are interested in, we do not expect a relevant constraint can be inferred from the search.
S can also be radiated off of muons produced in kaon decays at accelerator experiments, providing another direct probe of the coupling g µ . The most relevant bound comes from searches for K + → µ + + / E, with the corresponding production mechanism shown in Fig. 8b. The E949 experiment at Brookhaven limits BR(K + → µ + + / E) (2 − 3) × 10 −6 , depending on the kinematics of the BSM contribution [1]. In Ref. [22], the branching ratio for K + → µ + ν µ S was calculated to be 10 −8 for couplings g µ large enough to explain the (g − 2) µ discrepancy, with smaller g µ resulting in even smaller branching ratios. We thus conclude that these kaon decay experiments are not sensitive to the parameter space of interest.

Z decays
At higher energies, precision measurements of Z boson properties can also place constraints on the g µ − m S plane. In the muon-specific model, the corresponding constraints are weak. The best sensitivity arises in regions of the parameter space where BR(S → µ + µ − ) is non-negligible, so that one can obtain constraints from Z → µ + µ − S, S → µ + µ − . This process contributes to the 4µ decay of the Z, which has a measured branching ratio of BR(Z → µ + µ − µ + µ − ) = (3.5 ± 0.4) × 10 −6 [1] . However, we find that deviations at the level of these uncertainties still only arise for large couplings and masses that are already excluded by (g − 2) µ measurements. There can also be a bound from BR(Z → γ + X), reflected in Tab. 1. The current limit is set by an L3 search at LEP [63], BR(Z → γ + X) < 1.1 × 10 −6 , where the energy of the photon is required to be greater than ∼ 30 GeV. The decay Z → γS occurs through the SZγ coupling induced through a muon-loop (analogous to the Sγγ coupling), followed by S → χχ. The expressions for the loop-induced Z → γS decay width can be found in Refs. [64,65].The corresponding constraint again turns out to be very weak due to the Yukawa and loop suppression.
In the lepton-specific case, the couplings to τ 's can yield stronger sensitivity. In particular, the process Z → τ + τ − S, S → χχ can contribute to the measured Z → τ + τ − width, assuming that the kinematics associated with radiating the scalar S still allow events to pass the signal acceptance criteria. The decay width Γ(Z → τ + τ − ) is measured to be Γ(Z → τ + τ − ) = 84.08 ± 0.22 MeV [1]. To obtain an approximate 95% C.L. limit, we simply require that the width Γ(Z → τ + τ − S, S → χχ) be less than 1.96 times the uncertainty of the Γ(Z → τ + τ − ) measurement. The resulting approximate limit is shown in Fig. 5. Despite the larger coupling to τ 's, this constraint is still not strong enough to exclude the a µ favored region. The bound cuts off around ∼ 10 GeV, since for larger masses the branching fraction of S → χχ becomes significantly less than one for g µ 0.1 due to the large g τ . Note that the corresponding constraints from Z → e + e − or Z → µ + µ − are much weaker given the coupling hierarchy. The lepton-specific case also induces Z → γS decays, which, when followed by S → χχ, can be bounded as in Tab. 1. Again, however, only large couplings are excluded, g µ 0.7.
Future lepton colliders such as FCC-ee or CEPC could place stronger bounds on Z → Sγ. Both proposals plan to run at the Z-pole and can potentially produce 10 9 -10 12 Z bosons. We infer a sensitivity projection for a Tera-Z factory (corresponding to 10 12 Z bosons) using the results of Fig.  9(C) in Ref. [48], which shows the reach for Z boson decays to a photon and pseudo-scalar, followed by the invisible decay of the pseudo-scalar. The resulting projection is shown in Figs. 5-6. Although the expected sensitivity mostly lies within the 5σ (g − 2) µ exclusion region, it is important to notice that a Tera-Z factory could potentially probe the a µ -favored region for m S greater than 10 GeV. The reach cuts off around m S ∼ 80 GeV due to the phase space suppression. Other Tera-Z searches, such as Z → τ + τ − / E, could also in principle be relevant. But the corresponding sensitivity can be largely restricted by experimental systematic uncertainties and we leave a more detailed study in future. Also, above the Z pole, other high-energy collider searches might apply, however, since we are focused on light dark sectors below the electroweak scale, we defer their consideration to future study.

µ and τ decays
Given non-zero g µ , the scalar S can be radiated off of a muon in µ − → e − ν µνe and τ − → µ − ν τνµ decays. The corresponding diagrams are shown in Figs. 8c and 8d-1. The coupling g µ can thus be constrained by measurements of the µ − → e − / E and τ − → µ − / E decay widths. We have computed the decay widths of Γ(µ − → e − ν µνe S) and Γ(τ − → µ − ν τνµ S) for m S ≤ m µ − m e and m τ − m µ , respectively, in Madgraph5 2.5.1 [66] and compared them with the upper limits on the branching fractions in Ref. [1] (listed in Tab. 1). We find that, at 95% C.L., the measured muon lifetime generally constrains couplings g µ larger than 1. The measured τ − → µ − ν τνµ branching ratio yields somewhat stronger bounds, constraining g µ 0.5 at 95% C.L. for m S ≈ 20 MeV. The bound weakens as m S increases. We conclude that these bounds do not constrain the parameter space of interest in the muon-specific model.
Similarly to the muon-specific case, S can be produced in µ and τ decays in the lepton-specific scenario. In addition to the diagrams in Fig. 8c and Fig. 8d-1, S can also be radiated off of taus in the initial state, as in Fig. 8d-2. 11 While µ − → e − ν µνe still only constrains g µ O(1), the decay τ − → µ − ν τνµ restricts g µ 0.02 for m S ≈ 20 MeV at 95% C.L.. The bound is stronger than in the muon-specific case due to the m τ /m µ enhancement of the S − τ coupling, but is still not competitive with the other probes discussed above. Nevertheless, we show the corresponding sensitivity in Figs. 5 and 6.

Other fixed target and beam dump experiments
Fixed target experiments with electron beams can be sensitive to the coupling g e through the process e − N → e − N S(→ χχ) in the lepton-specific model. This is analogous to the production mechanism discussed in Sec. 3, except with the scalar being radiated off of an electron beam. Currently, the most sensitive of these experiments is NA64 running in electron mode. We re-interpret NA64-e limits on the dark photon kinetic mixing parameter ε [8] as limits on an effective S-electron coupling ε S ≡ g e /e where e = √ 4πα. The resulting bound is weaker than e.g. bounds from current dark matter direct detection, and is not able to probe the region favored by the (g − 2) µ measurement. Of course, this conclusion could change if g e /g µ > m e /m µ .
Figs. 5-6 show the expected reach of the proposed LDMX experiment [47] with a 16 GeV electron beam. Here we translate the bound on y ≡ 2 ϕ α χ (m χ /m S ) 3 from [47] with ϕ ≡ g e /e into a bound on g µ . The expected reach of LDMX is impressive, but again, because of the coupling hierarchy, accessing the a µ -favored region for m S 300 MeV in these models will require a muon beam experiment such as NA64-µ or M 3 . Note that a related search strategy is that implemented in the proposed BDX experiment [67], which will be sensitive to the production of χχ through S and the subsequent scattering of the dark matter with a detector volume down-beam from the interaction point. The sensitivity of BDX to the lepton-specific model is translated from the corresponding dark photon projections [8] and shown Figs. 5-6. We use E χ ≈ E beam = 11 GeV and E e − ≈ E th = 350 MeV [67] in accounting for the amplitude difference between the scattering of the ultra-relativistic DM and electron via a scalar and vector-mediator (see also [47]).The projected BDX sensitivity falls within the (g − 2) µ excluded regions and does not exceed that of e.g. LDMX. Analogous constraints from past beam dump experiments, such as E137 [68,69], yield significantly weaker sensitivity and are not shown in Figs. 5-6.
Proton beam dump experiments copiously produce pions, which can decay to final states involving muons. These muons can again radiate a light scalar S. The SeaQuest experiment is expected to provide sensitivity to light scalars with couplings to muons that decay visibly, as discussed in Ref. [70]. However, if S decays primarily invisibly, as we have assumed across the parameter space of interest, then the relevant final state is µ + µ − + / E, and more difficult to constrain. We are not aware of any existing or planned searches for rare pion decays targeting this final state.

Summary
To summarize, in the muon-specific case and for couplings near the size required to explain the (g −2) µ discrepancy, none of above constraints are stronger than those from the a µ measurements. This is reflected in the left-hand panels of Figs. 5 and 6. Supernova cooling places an additional constraint on low masses and small couplings. Muon beam experiments, such as NA64-µ, will be thus critical in exploring the parameter space of these models motivated by the a µ discrepancy and/or thermal relic dark matter, absent significant fine-tuning. 11 There is also a similar diagram for S radiated off from τ or e in τ − → e − νµνe and τ − → e − ντνe. Given the similar branching ratio of BR(τ → e / E) and BR(τ → µ / E) and ge gµ in the lepton-specific model, the latter search yields a stronger bound.
The right-hand panels of Figs. 5 and 6 compare the impact of the most sensitive aforementioned probes to the reach afforded by NA64-µ in models where the mediator couples to electrons, muons, and taus. Of the former, the most sensitive are cooling bounds from SN 1987A, which already constrain low masses/small couplings, and future mono-photon searches at Belle-II or Tera-Z factories. Future direct detection experiments will be sensitive to masses below ∼ 100 MeV, and could, together with NA64-µ, provide conclusive evidence for light leptophilic dark matter. LDMX-type electron fixed target experiments can also be sensitive to low masses, but by far the most powerful probe of these models would be provided by muon missing momentum experiments at NA64-µ or M 3 .

Conclusion
We have argued that sub-GeV dark matter with light leptophilic scalar mediators is a compelling target for current and future experiments. In addition to providing a viable dark matter candidate, these models can explain the long-standing (g −2) µ discrepancy. Although not a UV complete scenario, light leptophilic dark matter can arise as a viable effective field theory with phenomenology that deviates in important ways from e.g. gauged L µ − L τ models. In particular, neutrino experimental constraints are absent, allowing for a larger range of masses consistent with (g − 2) µ and not currently ruled out by other experiments.
Light DM with leptophilic scalar mediators is generally difficult to test, but we have argued that missing momentum searches at muon fixed target experiments, such as NA64-µ, can provide valuable coverage of these scenarios. In particular, NA64-µ will be sensitive to the entire region of parameter space consistent with the measured value of a µ for mediator masses below ∼ 5 GeV. Even better sensitivity could be achieved utilizing the"1 vs 2" technique for background elimination to accommodate 10 13 muons on target in a background-free environment. In addition to the a µ -favored regions, NA64-µ would also be able to explore a significant portion of the parameter space consistent with the observed dark matter density without fine-tuning. As such, muon beam experiments may afford us a first glimpse at what lies beyond the Standard Model.
Note added: As this paper was being finalized, Ref. [47] appeared that also discusses searches for invisibly-decaying leptophilic scalar mediators at missing momentum experiments, and thus overlaps with some of our results and discussion above.

A Scalar couplings to photons
In this appendix, we discuss the effective coupling of S to photons relevant for computing the diphoton decay width, the s-channel contribution to e + e − → Sγ (represented by Fig. 8a), and other g γγ -induced processes. It can be implemented in MC generator as a form factor. The coupling g γγ is induced by muon loops in the muon-specific case, as well as electron and tau loops in the lepton-specific case. The computation proceeds similarly to the calculation of the Higgs-diphoton effective interaction in the SM, except allowing for general values of p 2 S and q 2 , the four-momentum-squared of S and one of the photons, respectively (the other photon is assumed to be on-shell). We compute the lepton loop diagram using Package-X 2.1 [71] and match onto the tree-level effective coupling presented in Eq. (2.4). We obtain the following result for the form factor F 1/2 defined in Eq. (2.4): The expression above is calculated using the (+, −, −, −) metric convention. For computing the contribution to e + e − → Sγ, p 2 S = m 2 S and q 2 = s, where s is the center-of-mass energy of the given collider, s (10.4 GeV) 2 , for BaBar, and s (10.58 GeV) 2 for Belle-II.
For computing the S → γγ decay width, p 2 S = m 2 S and q 2 = 0, which leads to the well-known result