On sub-GeV Dark Matter Production at Fixed-Target Experiments

We analyze the sensitivity of fixed-target experiments to sub-GeV thermal relic dark matter models, accounting for variations in both mediator and dark matter mass, and including dark matter production through both on- and off-shell mediators. It is commonly thought that the sensitivity of such experiments is predicated on the existence of an on-shell mediator that is produced and then decays to dark matter. While accelerators do provide a unique opportunity to probe the mediator directly, our analysis demonstrates that their sensitivity extends beyond this commonly discussed regime. In particular, we provide sensitivity calculations that extend into both the effective field theory regime where the mediator is much heavier than the dark matter and the regime of an off-shell mediator lighter than a dark matter particle-antiparticle pair. Our calculations also elucidate the resonance regime, making it clear that all but a fine-tuned region of thermal freeze-out parameter space for a range of simple models is well covered.


INTRODUCTION
The evidence for a dark matter (DM) component in the universe, through its gravitational effects over many distance scales, provides the strongest indicator to-date for physics beyond the Standard Model (SM). Over the past decade, theoretical and experimental approaches to particle DM have evolved significantly, with a loosening of theoretical priors about how DM may be explained within particle physics. This has been driven on the one hand by the increasingly stringent LHC constraints on new TeV-scale degrees of freedom, and on the other by the recognition that the strong empirical evidence for DM motivates exploring all viable and testable scenarios for DM, not just those that are linked to other expectations for new physics (e.g., the naturalness problem motivating electroweak-scale WIMPs, or the strong CP problem motivating axions).
Within this broader perspective, there is strong theoretical motivation for considering thermal relic models where the DM was once in thermal equilibrium with ordinary matter. This assumption constitutes a "minimal" cosmology in two important ways: the cosmological evolution of DM has minimal dependence on initial cosmological conditions, and at the same time closely mimics the early-universe production processes known to accurately describe Big Bang nucleosynthesis (BBN). The full mass range for such thermal relics is wider than for conventional electroweak WIMPs, which must be heavier than a few GeV as a result of the Lee-Weinberg bound [1]. The light DM window from an MeV to a GeV is viable for DM that is neutral under SM gauge groups but possesses distinct interactions [4]. The lower boundary of the light DM mass range is motivated by the successful predictions of BBN and observations of the cosmic microwave background (CMB) (see Refs. [2,3] for an exception). The upper boundary is purely conventional, as heavier hidden-sector DM has increasingly WIMP-like phenomenology. The light DM paradigm has been explored in some depth over the past 5-10 years [5][6][7] in both direct detection [8][9][10][11][12][13][14] and accelerator-based experiments using electron [15][16][17][18][19][20][21][22] and proton [23][24][25][26][27][28][29][30][31][32][33][34][35] beams (see also Refs. [15, for studies of related dark sectors). Yet the multi-dimensional parameter space of light DM models poses a challenge for assessing both the status of the field and prospects for upcoming searches.
From a general effective field theory (EFT) perspective, the leading couplings of the SM to a fully neutral dark sector could arise from vector, Higgs, or neutrino portals, which rely on the presence of an associated vector (A µ ), scalar (S), or singlet fermion (N ) mediator, respectively. These portals are the unique marginal or relevant operators that can couple the two sectors, and the neutrino portal provides the simplest model of neutrino mass. Thermal freeze-out is generically characterized by the requirement of a weak scale annihilation cross section, σ ann v ∼ 1 pb. This ensures that the frozen out mass density matches the observed relic abundance of DM. For orientation, we can parametrize the interactions between DM fields, χ, and the SM through a set of higher-dimensional operators of the form where J µ... denotes a set of (pseudo-)scalar, (axial-)vector, or other currents associated with the mediation channels between DM (J D ) and SM fields (J SM ). In the minimal case that 2 → 2 processes mediate DM annihilation through dimension-six interactions, we have n = 2 in Eq. (1) and thermal freeze-out implies the parametric relation where c is a numerical factor that depends on spins and other quantum numbers of the currents. This scaling explains why WIMPs (with m χ ∼ Λ ∼ TeV) are archetypal DM candidates, but allows a broad range of DM masses when Λ is allowed to vary. The relation of Eq. (2) is a schematic form of the thermal target that models must satisfy in order to reproduce all of the observed DM and must at a minimum exceed in order not to overpopulate the matter content of the universe. An operator description, as in Eq. (1), is valid at energy scales below the mass of the force mediators responsible for coupling SM particles to DM. If the force mediator mass scale is experimentally accessible or relevant during thermal freeze-out, then an accurate description of DM phenomenology and/or cosmology requires a more complete theoryone that includes the mediator particle explicitly, rather than integrating it out as in Eq. (1). Indeed, the presence of light dark force mediators is a critical ingredient to bypassing the Lee-Weinberg bound. Thus, light DM must form part of a multi-component dark sector, and it is reasonable to expect that the mediator mass may be relatively near the DM mass.
Motivated by models of light DM, much theoretical and experimental effort over the past decade has focused on the vector portal, which at low energies involves kinetic mixing between the photon and the dark vector, ( /2) F µν F µν . This is in part because it is the least constrained scenario that allows for bilinear mixing. It also provides the most scope for model building, including what has become the benchmark model for sub-GeV DM [5,6]. We will consider the vector portal, and modest generalizations of it, as our primary examples below as we discuss the domain of validity of effective theory and the full parameter space of DM detection.
Under the assumptions of Eq. (2), the EFT is always valid for the kinematics of direct detection, which involves momentum transfers parametrically smaller than the DM particle mass. The most convenient description of the interaction is via a relativistic (or non-relativistic) EFT for DM interactions with nucleons or electrons, which is entirely determined by Eq. (1). For example, for fermionic DM χ coupled to SM vector currents, the leading interactions are typically of the form The scattering cross section on nucleons then depends parametrically on µ 2 N,χ /Λ 4 , with µ N,χ the reduced DM-nucleon mass, and so direct detection experiments simply constrain the cutoff scale Λ for a given DM mass.
By contrast, the usual analysis of accelerator-based searches starts from the assumption that the mediator is, in fact, light enough to be produced. We shall see that this assumption is often reasonable, but it is illustrative to instead begin from the EFT description (Eq. (3)) of the DM-SM interaction. The result of this analysis is illustrated in Fig. 1 for a particular operator (see Eq. (12) below) and a particular experiment (the proposed LDMX experiment, see Sec. 4 below). Other operators of the same dimension lead to results that are qualitatively very similar. For a given DM mass m χ on the horizontal axis of Fig. 1, smaller values of Λ −1 on the vertical axis correspond to greater sensitivity. This illustrates LDMX's sensitivity, within the EFT treatment, to the interaction strengths expected for thermal freeze-out for DM particle masses below ∼ 200 MeV.
The maximum momentum transfer in the kinematics explored by LDMX is √ s ∼ GeV, such that the EFT of Eq. (3) does not become strongly coupled anywhere within the range of Fig. 1. Nonetheless, there is reason to believe that the EFT has a smaller region of validity. Concretely, other constraints at GeV energy scales motivate the assumption that the interaction of Eq. (3) is mediated by tree-level exchange of a mediator (of mass m V ) that couples to SM fields with interaction strengths g SM and to DM with strength g D . In this case, we expect (other high-energy completions for the contact interaction of Eq. (3) would change the constant of proportionality, while maintaining Λ ∝ m V ). Constraints on g SM from electron self-interactions imply g e O(10 −3 ) − O(10 −2 ) depending on the precise mediator mass. The EFT treatment is valid only for m V > √ s -a constraint that (using Eq. (4) and the constraints on g e ) is inconsistent with large Λ −1 , and is also invalid in models where a small Λ −1 results from a light mediator with weaker coupling. A further complication to the EFT interpretation of accelerator-based searches is that different experiments (and even different production modes relevant within a single experiment) have different maximum energies √ s, so that the regions of consistency of the EFT vary from one experiment to another. We emphasize, however, that reaching beyond the regime where EFT is valid is not a deficit of accelerator-based 10 -3 10 -2 10 -1 1 10 -2 experiments -rather, it is a signal that they can potentially uncover the mediation mechanism behind the freeze-out process.
For these reasons, the experimental constraints should instead be evaluated in the larger parameter space: This same enlarged parameter space is also motivated by the calculation of the DM relic abundance resulting from thermal freeze-out. When m V m DM , the DM abundance reduces to a function of the cutoff Λ 2 ∝ m 2 V gSM g D , but for m V 2 m DM the thermal abundance has more complex dependence on the four parameters above.
Practical limitations motivate examining the sensitivity of different experimental probes in two-dimensional slices of this parameter space, of which several have been considered in the literature. The most commonly used is to fix m V = 3 m DM , which generally yields a conservative estimate of experimental sensitivity within the region where the EFT is approximately valid at freeze-out and mediator production is kinematically accessible. The resulting 3-dimensional parameter space is further reduced by either saturating g D at a large value near the perturbativity bound or fixing the product g SM · g D to the value predicted by thermal freeze-out for given m DM . However, fixing the ratio of m V to m DM obscures some important features of the parameter space and of experiments' sensitivities, particularly in the EFT limit and the resonance region (the latter was the focus of Ref. [61]).
With the experimental effort to fully explore light DM now becoming a reality, 1 it is timely to analyze the full thermal relic parameter space, in order to determine what capabilities are needed to fully test and either discover or exclude simple models of light thermal DM. In this paper, we initiate such an analysis and uncover some generic features of these models that were less visible in the conventional parameter slices considered thus far. In addition, the extended coverage allows the sensitivity to be tracked from the EFT limit, with heavy off-shell mediators, to the resonance region where these mediators may be produced on-shell in fixed-target experiments, and all the way to off-shell light mediators below the DM pair threshold. For concreteness, our analysis will focus on models that involve a dark photon mediator, but the features uncovered are generic for models within this dark sector paradigm.
In the next section, we summarize and motivate the thermal DM framework in more detail. In Sec. 3, we introduce a new parameter plane, which allows for a more comprehensive analysis of the reach of experiments. In Secs. 4 and 5, we summarize the fixed-target experimental facilities to be analyzed and the primary production modes. Our results, presenting a comprehensive picture of the reach of these facilities, appears in Sec. 6, and we finish with some concluding remarks in Sec. 7.

DARK SECTORS AND PREDICTIVE LIGHT DM MODELS
In this section, we will briefly review the assumptions made in constructing light DM models in the thermal relic window. Fermion models of DM are especially natural to consider for several reasons. First, all stable matter that we know of is comprised of fundamental fermions -it would be natural for this to be true of DM as well. Second, the stability of new fermions is rather easy to realize if the analog of lepton or baryon number, i.e., "dark lepton number," is preserved. This is a rather natural feature of many models, especially those where the fermions are coupled to a vector field A , as illustrated below. Third, fermion masses can be technically natural -radiatively insensitive to unknown UV physics.
The fermion model with the fewest number of degrees of freedom, and therefore the most "minimal" DM model in some sense, consists of a single two-component Weyl fermion ξ. In this case, the vector and scalar interactions are uniquely fixed by relativistic covariance to be, where we have grouped ξ into a 4-component Majorana spinor χ = (ξ ξ † ) T . Without additional structure, "ξ-number" is automatically conserved, so χ-particles will be stable. Majorana fermion DM arises naturally in supersymmetric dark sector models and more generally in scenarios involving additional U (1) D gauge groups, such that U (1) D -breaking mass terms are larger than U (1) D -preserving mass terms. A well-motivated hypothesis for the origin of DM is that it reached thermal equilibrium with ordinary matter in the hot early universe, then "froze out" of equilibrium as the universe cooled, leaving a residual DM abundance surviving to the present day. This thermal DM scenario is both simple, readily achieving the observed DM abundance, and predictive, requiring a precise cross section for DM annihilation into ordinary matter to match the observed density. In many models, this in turn implies a minimum interaction strength for DM, which serves as an important sensitivity benchmark for DM searches.
For DM lighter than a GeV, DM couplings to the lightest SM species (such as the SM fermions, f ) can easily explain DM's thermal origin. At low-energies, such couplings can be parametrized as four-Fermi contact operators suppressed by the mass-scale of the interaction where c v , c s are dimensionless coupling constants and J µ f , J f are SM currents. Examples of such currents include where f is a SM lepton ( ) or nucleon (N ). The DM annihilation rate near freeze-out is parametrically of size σv ∼ c 2 v,s m 2 χ /Λ 4 , implying that thermal freeze-out is consistent with the measured DM energy density for where T eq is the temperature at matter-radiation equality and m pl is the Planck mass. If DM annihilates efficiently to SM matter, as it must for the thermal freeze-out mechanism to produce the observed abundance, then we can look for evidence of this annihilation in the universe at late times. In the sub-GeV mass range, the tightest constraint comes from the annihilation of DM at eV-scale temperatures, shortly after recombination. In the early universe, the annihilation products could have reionized atomic hydrogen, altering the CMB power spectrum. The total power injected has been constrained by Planck, leading to the constraint [62,63] where σv T denotes the DM annihilation cross section to visible final states at a temperature T . It follows that, for MeV−GeV scale thermal DM, the annihilation rate to visible final states at recombination (T ∼ eV) must be considerably suppressed (by 2 − 5 orders of magnitude) relative to the annihilation near freeze-out (T ∼ m χ ). This is a revealing constraint, but not a problematic one. Many of the generic models for thermal DM involve a significant suppression of the annihilation rate at low temperatures. This is typically due to one of two effects: velocity-suppression or population-suppression. Annihilations of fermionic DM coupled to the SM through spin-0 currents or Majorana DM coupled to the SM through spin-1 currents, as in Eq. (7), are p-wave with σv ∝ v 2 and hence are velocity-suppressed at the time of recombination. This is also the case for scalar DM coupled to the SM through a spin-1 current, "Population-suppression" arises in scenarios involving U (1) D -breaking mass terms (δm) that are much smaller than U (1) D -preserving mass terms (m D ). Such models involve the vector coupling of nearly-degenerate Majorana massstates, χ 1,2 , If χ 2 is only slightly heavier than χ 1 , with general mass eigenvalues m 1,2 m D ∓ δm, the leading annihilation process during freeze-out involves the excited state, χ 2 , which is thermally depopulated at later times but well-before recombination.
Details of the DM model, such as the spin-type of the interaction and the structure of the DM mass matrix, can dramatically suppress signals relying on non-relativistic DM scattering, such as direct detection. Indeed, in many scenarios of DM lighter than a few GeV, the strong suppression of DM annihilation at low velocities, as required by the CMB, is accompanied by a suppression of low-velocity scattering. One virtue of accelerator-based searches is that they involve momentum transfers comparable to or greater than the DM mass. Their sensitivity is therefore independent of these details, allowing a robust and inclusive test of the sub-GeV thermal relic scenario.

A. Force mediators and UV-complete portals
The higher-dimensional operators discussed above are a valid description of a given process if the kinematic scale of that process is small compared to the mass-scale of the interaction, Λ. As summarized in the introduction, parametrically small couplings, as expected in dark sectors, can complicate the understanding of when the EFT is valid. Indeed, weakly-coupled light mediators will, when integrated out, lead to an EFT with a high characteristic cutoff scale Λ, e.g., as given in Eq. (4), even though the mediator mass scale is parametrically smaller.
Each of the EFT interaction channels above can be extended to a wider kinematic range by "integrating in" the mediator field at tree-level, so that schematically where V is the mediator field. In practice, while low-energy descriptions like this are possible for each of the EFT operators discussed above, not all can be rendered independent of unknown short-distance physics (namely ultraviolet or UV-complete) in a simple manner. In many cases, the interactions are secretly still higher-dimensional when accounting either for anomalies or simply the chiral structure induced by the SM electroweak sector. As is now well known, there are only three UV-complete interactions of this type that are consistent with the SM electroweak symmetry-breaking structure. They involve couplings to new dark photon (A ), dark scalar (S), or dark fermion (N ) degrees of freedom, and are characterized as follows: in terms of SM photon, Higgs and lepton doublet fields. On general EFT grounds, if these dark photon, scalar, or fermion degrees of freedom are present, we expect the three renormalizable portals to provide the leading sensitivity to them.
The simplest and most viable way in which light DM can interact with the SM is through the mixing of any force that DM is charged under ("dark forces") with the photon. The basic setup that realizes this possibility consists of a massive dark sector gauge boson A , with general Lagrangian containing, Here, is the kinetic mixing parameter, A is the massive "dark photon" of a broken U (1) D symmetry, A is the SM photon, F µν and F µν are the dark-photon and electromagnetic field-strength tensors, J µ em is the electromagnetic current, J µ D is the dark current, and m A is the dark photon's mass. If m A O(100) GeV, the dark photon predominantly mixes with the SM photon so that the visible sector acquires a millicharge under U (1) D , In this setup, the U (1) D symmetry is broken by the mass of the A . Thus, the effective theory analysis of different scenarios should include other sources of soft U (1) D breaking, which we include below. General discussions of hidden sector physics often include other possible portals to the dark sector: e.g., the UVcomplete Higgs or neutrino portals noted above, dark vectors coupled primarily to lepton or baryon numbers, or pseudoscalars coupled through the axion portal. Before proceeding, we pause briefly to explain why the focus on dark photons is appropriate here. The more general case of dark vector mediators is UV-sensitive, and thus theoretically distinct, but phenomenologically identical to dark photon physics except that the dark gauge boson couplings to SM particles are not proportional to charge. This changes the relative sensitivity of different DM production/scattering experiments in the expected way (leptophilic models favor DM-electron scattering and electron-production over DMnuclear scattering and hadron-production, while the opposite is true for leptophobic models) and introduces other constraints on the mediator-SM coupling strength that are beyond the scope of this work. The most minimal neutrino portal is not a viable portal for thermal DM production, though right-handed neutrinos are a viable DM candidate in their own right. Furthermore, predictive scenarios involving dark (pseudo-)scalars are strongly constrained by meson decay constraints [64]. Such considerations motivate the dark photon vector portal as the most viable, and we will focus on that case in the first instance, but will aim to illustrate the expectations for the other viable mediation channels as well.

INTRODUCING THE y − R PLANE
The DM reach of accelerator experiments is often presented in terms of the canonical parameter, . This is simply related to the EFT cutoff Λ (y = (16π 2 α em ) −1 (m χ /Λ) 4 0.9 (m χ /Λ) 4 ), and well-motivated from considerations of DM production in the early universe. For m A few×m χ , the DM annihilation rate near freeze-out is proportional to σv ∝ α em y/m 2 χ , implying that the cosmological χ abundance is in agreement with the observed DM energy density for where T eq is the temperature at matter-radiation equality and m pl is the Planck mass. Although this cosmologicallymotivated range for y is independent of m A /m χ for m A few × m χ , DM observables at low-energy accelerators often scale with different powers of m χ and m A . Furthermore, the particular scaling of such observables depends on, e.g., the hierarchy between the typical collisional center of mass energy, E cm , and m A . As a result, when presenting the sensitivities of terrestrial experiments in the y − m χ plane, it has become standard convention to fix the dark photon-to-DM mass ratio to one or several representative values, such as m A /m χ = 3, 5, etc. In this section, we further explore the range of DM parameter space spanned by y and R ≡ m A /m χ , including regimes where DM production proceeds through off-shell mediators, i.e., m A 2m χ and m A E cm . In Fig. 2, we schematically illustrate the characteristic sensitivity scaling of electron fixed-target (blue), proton fixed-target (orange), and high-energy collider (gray) experiments in the y − R plane for a fixed DM mass. Along the solid black line, the thermal freeze-out of χχ → A * → ff (where f is an electrically charged SM fermion) in the early universe leads to an abundance of χ that is in agreement with the observed DM energy density. For R 1, the cosmologically-favored value of y is given parametrically by Eq. (16). Note that this region is independent of the particular value of the DM-dark photon coupling, α D . Near R 2, χχ → A * → ff is resonantly enhanced, allowing for an adequate annihilation rate near freeze-out with much smaller couplings [61,65]. For R 2, the rate for -independent reactions, such as χχ → A A and χχA → A A , can significantly modify the process of thermal freeze-out, leading to an adequate abundance of χ for significantly smaller values of y (provided that α D is sufficiently large and that the -dependent decay, A → e + e − , occurs sufficiently rapidly) [66][67][68]. In this "secluded" regime, there is no sharp cosmological target and the cosmologically viable range of y is extended to much smaller couplings. The sector of parameter space where the secluded annihilation χχ → A A dominates over direct annihilations, χχ → ff , is shown as the solid green region of Fig. 2 (for a particular choice of α D ).
In the remainder of this section, we will briefly discuss the qualitative behavior of the fixed-target and collider sensitivities in the y − R plane, as shown in Fig. 2. More detailed discussion will be provided below in Secs. 4-6.
• For R 2, the on-shell decay A → χχ is kinematically forbidden, and hence DM production proceeds through reactions involving off-shell dark photons. Since in this case the typical accelerator center of mass energies are much greater than the dark photon mass, the DM production rate is independent of m A and, hence, the corresponding sensitivities scale as y reach ∝ R −4 .
• For R 2 and m A E cm , DM production proceeds through on-shell production and decay of dark photons. In this region of parameter space, the signal yield for fixed at high-energy collider experiments is independent of m A , so that y reach ∝ R −4 , while the yield at fixed-target experiments typically scales as m −2 A , in which case the sensitivity to y scales as y reach ∝ R −2 . In electron fixed-target experiments, the finite size of the nucleus suppresses high-m A production leading to a gradual flattening of the yield as m A approaches E cm . For proton beams, there is additional structure in the R-dependence that emerges from thresholds in production from hadronic decays, and also resonant enhancements. This will be reviewed in Sec. 5.
• For m A E cm , collisions of electron or proton beams are unable to produce on-shell dark photons and DM production again occurs through an off-shell dark photon. In this regime, the rate of DM production decouples as m −4 A . In this case, the dark photon can be integrated out from the low-energy theory and the sensitivity in the y − R plane asymptotes to an R-independent value, y reach ∝ constant.

SUMMARY OF EXPERIMENTAL SIGNALS AT FIXED-TARGET EXPERIMENTS
In this section, we summarize the existing and proposed electron and proton fixed-target experiments whose capabilities for light DM detection will be analyzed in subsequent sections. Here, we provide a brief introduction and summary of detection signatures and provide further details of production modeling in the next section.
• LDMX is a proposed experiment designed to search for signals of missing momentum in electron-nuclear fixedtarget collisions of a high repetition rate ( 50 MHz) and energetic ( GeV) electron beam [69]. The momentum of every beam electron is measured by silicon trackers before and after it scatters in a thin ∼ 10% radiationlength tungsten target. This allows for a direct measurement of the momentum transfer that occurs in each electron-nuclear collision. Electromagnetic and hadronic calorimeters are placed downstream of this target region, in order to detect visible activity, such as energy deposition from the recoiling beam electron and any other charged or neutral SM particles. DM that is produced in the electron-nuclear collision is characterized by a large energy loss and momentum exchange of the beam electron and the absence of measured activity in the downstream calorimeters, aside from the soft recoiling electron. The projected reach will be shown as blue lines/regions (see Figs. 3-6 in Sec. 6).
• NA64 is an existing experiment that searches for missing energy signals in electron-nuclear fixed-target collisions, utilizing the 100 GeV secondary electron beam at the CERN SPS on a thick lead target. Its setup is akin to that of the proposed LDMX configuration, without the downstream tracking sufficient to accurately measure missing momentum. Instead, DM that is produced in beam-target collisions registers as missing energy, which is inferred from the electromagnetic shower of the recoil electron in a downstream electromagnetic calorimeter.
In this work, we recast the recently published limits from the NA64 collaboration, corresponding to a total of 2.84 × 10 11 electrons on target (EOT) and background levels of 1 event [22]. The ensuing limits will be shown as solid gray (see Figs. 3-6 in Sec. 6).
• BDX is a proposed DM beam dump experiment at JLab [70]. It takes advantage of the high-current 11 GeV CEBAF electron beam on a thick aluminum target and is expected to acquire 10 22 EOT within a year of parasitic running downstream of JLab Hall-A. DM that is produced in the electron-nuclear collisions is detected upon scattering in a shielded electromagnetic calorimeter placed ∼ 25 m downstream of the target. This detector consists of ∼ 800 CsI(Tl) crystals over an equivalent volume of a ∼ cubic meter and is sensitive to both DMproton and DM-electron scattering. The projected reach will be shown as purple lines/regions (see Figs. 3-6 in Sec. 6).
• SBND is under development as the near-detector component of the Fermilab short baseline neutrino (SBN) facility. Fed by the 8.9 GeV Booster Neutrino Beam (BNB), with 10 21 protons on target (POT) expected for a 4 year run, SBND will be a 77 ton liquid argon (LAr) TPC, with a short baseline of 110 m. DM can be produced through proton-nucleon scattering in the beryllium target and absorber, via several channels including secondary π 0 and η decays. SBND is being designed to have low thresholds and will be sensitive to a variety of elastic and inelastic scattering signatures, both of neutrinos and potential light DM species. We focus on elastic electron scattering in this work, which allows for good neutrino background rejection via a tight forward angle cut, as utilized to good effect in the recent MiniBooNE off-target DM search run [32,33]. The projected reach will be shown as cyan lines/regions (see Figs. 3-6 in Sec. 6).
• NOνA is an operational long-baseline neutrino oscillation experiment, running from the 120 GeV NuMI beamline at Fermilab, and has now accumulated ∼ 10 21 POT in both neutrino and antineutrino mode. Of most interest here is the 300 ton near detector, located 990 m from the carbon target. DM can be produced in the target and absorber as in the BNB, but the higher energy beam allows more significant production from, e.g., proton bremsstrahlung, and leads to a forward peaked DM distribution. We again focus on electron scattering, following the analysis in Ref. [71], as the most effective means of limiting neutrino backgrounds. The projected reach is shown as yellow lines in Figs. 3 and 4 (see Sec. 6).
• COHERENT is an expanding series of detectors at the SNS, focused on measuring coherent neutrino-nucleus scattering. The 1 GeV proton beam produces 10 23 POT/year and a wide angle distribution of secondary mesons, and DM can be produced through π 0 decays or radiative pion capture. We utilize the recent analysis with the CENNS-10 liquid argon detector, a 27 kg shielded detector located at about 120 degrees to the beamline at a distance of approximately 28 m. The analysis placed a 1σ upper limit of 7.4 on the number of potential coherent elastic neutrino-nucleus scattering events [72], and we constrain light DM scenarios by assuming that all of these potential events were the product of DM scattering. This is shown as the gray region in the right panel of Fig. 5 (see Sec. 6). The COHERENT collaboration recently published projections for the DM sensitivity of a proposed next generation ton-scale liquid argon detector [73].
• CCM (Coherent CAPTAIN-Mills) is a newly commissioned experiment using the 800 MeV proton beamline at LANSCE-Lujan, Los Alamos which can produce 10 20 POT/yr. Similar to COHERENT at the SNS, it is designed to detect coherent neutrino scattering. However, the setup includes a 10 ton LAr detector placed 20 m from the tungsten target, which allows enhanced sensitivity, despite it being positioned at an angle of 120 degrees with respect to the beamline [74]. We again focus on DM produced through π 0 decays and coherently scattering off of nuclei in the detector. The projected reach will be shown as brown lines (see Figs. 3-5 in Sec. 6).
In addition to these fixed-target searches, we will also highlight in Figs. 3-6 the existing constraints from a monophoton search at BaBar [20,75] and electroweak precision measurements at LEP (shaded gray) [76,77], as well as the projected sensitivity of a monophoton search at Belle II with 20 fb −1 of integrated luminosity (red lines/regions) [78].

A. Off-shell production in electron fixed-targets
We model the primary production of DM in electron-nuclear collisions in a modified version of MadGraph5 [79]. This proceeds through on-or off-shell bremsstrahlung of an intermediate dark photon, i.e., eN → eN A * → eN χχ. In this section, we briefly describe the computational procedure that we adopt throughout this work.
Instead of generating independent signal event samples for each possible combination of the DM and dark photon mass, we generate one large statistical sample in the heavy-mediator/contact-operator limit for a discrete set of DM masses. As we will explain in more detail below, for each DM mass, the differential invariant mass spectrum can be rescaled accordingly for any dark photon mass. This allows for a more computationally efficient, but still accurate, calculation of the final signal yield.
As described above, the first step in our procedure is to generate large statistical samples of events for eN → eN χχ in the contact-operator limit for a given experimental configuration, DM mass, spin, and interaction type. For instance, pseudo-Dirac DM particles that interact through a spin-1 coupling to electrons are modeled through a Lagrangian of the form where g χ is a dimensionful coupling controlling the strength of the interaction, and similarly for other interaction types. We then bin, in the invariant mass of the DM pair m χχ , the subset of these events that pass the given experimental selection criteria. From this, we obtain the differential invariant mass spectrum, dσ EFT /dm χχ , in the heavy-dark-photon-limit. To obtain the invariant mass spectrum for any choice of dark photon mass or couplings, we rescale this distribution bin-by-bin by the squared propagator of an intermediate dark photon, where above, g χ is the value of the contact-operator coupling that was chosen in making the initial distribution and Γ A is the total dark photon width. In practice, we approximate the A to have an O(1) branching ratio to DM pairs, i.e., Γ A Γ(A → χχ), which is valid provided that α D α em 2 . The total signal rate (accounting for experimental selection efficiencies) is calculated by integrating the rescaled distributions over all relevant values of the DM invariant mass, When simulating the DM production in this manner for each electron beam fixed-target experiment, we assume negligible background and utilize the following experimental parameters: • LDMX -We adopt a beam energy of E beam = 16 GeV, a tungsten target thickness of 0.1 radiation-lengths, an integrated luminosity of 10 16 EOT, and a 50% signal efficiency. The signal region is defined by imposing that the energy of the recoil electron is E e 0.3 E beam .
• NA64 -In recasting the recently published limits from the NA64 collaboration, we adopt a beam energy of E beam = 100 GeV, a thick lead target, an integrated luminosity of 2.84 × 10 11 EOT, and a 50% signal efficiency. The signal region is defined by imposing that the energy of the recoil electron is E e 0.5 E beam .
• BDX -We assume a beam energy of E beam = 11 GeV, a thick aluminum target, 10 22 EOT, and a 20% signal detection efficiency. The target-detector distance is taken to be 25 m and we model the front-face size and length of the CSI(Tl) detector as 50 × 40 cm 2 and 250 cm, respectively. From our simulated sample of produced DM pairs, we analytically incorporate DM-electron scattering in the detector. The signal region is defined by imposing that the energy of the scattered electron is E e 500 MeV.

B. Off-shell Production in Proton fixed-targets
The DM production rate at proton beam facilities is determined by σ(pN → A * + · · · ). The primary production modes have been analyzed in the literature (see Refs. [23-26, 58, 80-83]) and include: • π 0 /η decay in flight -The channel π 0 , η −→ γ + A * (→ χχ) is often the dominant mode. Given the 3-body decay distribution in the contact operator limit, we can write the full differential distribution (with A on-or off-shell) as where q 2 = m 2 χχ is the invariant mass of a DM pair, i.e., the dark photon virtuality.
• Bremsstrahlung with resonant vector meson mixing -For this channel p + N → p + N + A * (→ χχ), standard approaches involving variants of the Weizsacker Williams approximation require that kinematic criteria such as E p , E A m V be satisfied [31], which are restrictive if A is highly off-shell. However, in practice this mode is significant only over the m A range for which the timelike form-factor of the radiated A is dominated by the vector resonances, e.g., m A ∼ m ρ,ω . Therefore, for the DM masses of interest here, the on-shell approximation is sufficient.
• Drell-Yan production from quark/gluon constituents -Once the dark photon virtuality q 2 is above the characteristic hadronic scale, parton-level Drell-Yan processes qq → A * (→ χχ) provide the relevant description. We have where L qq (q 2 ) is the parton luminosity.
We focus our attention on meson decays and bremsstrahlung as the two channels expected to dominate the low-mass regime to which proton beam dump experiments are most sensitive. We make use of the following parameters when simulating the signal at each experiment: • SBND -The Sanford Wang parameterization [84] was sampled to generate mesons produced using the BNB, with π 0 s simulated by the the mean of a π + and π − distribution and η's from the K 0 distribution. The overall production rates were estimated to be N π 0 = 0.9×POT and N η = N π 0 /30. For bremsstrahlung, we placed limits , adjusting the limits on z such that zE p , (1 − z)E p > 3m A and the transverse A momentum satisfies p T < 0.1 GeV. For the NCE electron scattering, we assume that all neutrino backgrounds can be rejected with a cos θ > 0.99 cut on the scattering angle of the recoil electron, and a sensitivity limit can be placed with 2.3 events. For the inelastic π 0 scattering case, we make no cuts on the recoil particles but assume that at least 60 events are required to overcome existing backgrounds. The event rate is calculated for 10 21 POT with a detection efficiency of 60% in both cases.
• NOνA -Both the π 0 and η production distributions from the NuMI beamline were simulated by sampling the BMPT distribution [85]. The π 0 production rate N π 0 was conservatively estimated at 1 per POT with N η = 0.078N π 0 . For bremsstrahlung, we put limits of z ∈ [0.1, 0.9] and p T < 1 GeV, once again limiting z such that zE p , (1 − z)E p > 3m A . Only NCE-like electron scattering was considered for NOνA. We follow the treatment of Ref. [71], selecting recoil electrons with E ∈ [0.5, 5] GeV and Eθ 2 < 0.005 GeV rad 2 . We place a limit on 41 scattering events with a detection efficiency of 50%.
• COHERENT and CCM -The Burman-Smith parametrization [86] was sampled to generate mesons produced along the SNS or LANSCE beams, with π 0 s simulated by a π + distribution due to the lack of a π − distribution at these energies. For COHERENT, we take N π 0 = 0.06 × POT, while for CCM we take it to be N π 0 = 0.0425 × POT. For COHERENT, we considered coherent DM scattering off of Liquid Argon with a minimum energy cut of 80 keVnr, 2 and placed a 1σ limit on 7.4 coherent DM scattering events over a 4.2 × 10 22 POT run. At CCM, we adopted a minimum energy cut of 10 keVnr and placed a limit on 10 coherent DM scattering events. In both cases, we assume 50% efficiency.
The projected sensitivity contours in Fig. 6 were generated by simulating the expected signal at each experiment for 2500 combinations of m A and m χ .

DM SENSITIVITY CONTOURS AT FIXED-TARGETS
We now combine the production modeling above with the relevant detection signatures at each experiment to produce sensitivity contours. The contours are chosen to reflect either known or estimated backgrounds, as outlined in each case. We will present two classes of figures, with the goal of exhibiting the full parameter space of light DM coupled to the SM through the vector portal.
A. The y − R plane In this section, we present the fixed target and collider sensitivities in the y − R plane, where R ≡ m A /m χ . For a fixed choice of DM mass m χ , this allows us to explore the extrapolation in sensitivity from the fully off-shell EFT regime at large masses (m A E cm ) to the opposite regime where on-shell decays may enhance the rate, or resonant effects become important. In Fig. 3, we show results for pseudo-Dirac DM and two choices of m χ (10 MeV and 100 MeV) and DM-dark photon coupling (α D = 0.1, 0.5). Many of the qualitative features were already discussed briefly in Sec. 3. We summarize some of the physical features below: • Thermal target -The parameters required to ensure the full DM relic abundance from freeze-out is shown as a solid/dotted black line. The annihilation process, χχ → A * → ff , is resonantly enhanced near R = 2. For R 2, "secluded" processes such as χχ → A A lead to an adequate abundance for much smaller couplings to the SM and there is effectively no sharp thermal target (shaded green).
• EFT region -The R 1 regions of the figures show that the sensitivity of fixed-target experiments asymptotes to a fixed value. This occurs beyond the value of R at which the production becomes predominantly off-shell and well-approximated by a contact operator. Note that the contours for LEP (and, in Fig. 3 (right), BaBar and Belle II) do not flatten within the range of the plot, as the A can still be produced on-shell for the parameter range shown.
• Resonant production -For proton fixed-target experiments, bremsstrahlung of the dark photon mediator (and thus DM production) is resonantly enhanced when m A is close to the mass of one of the SM vector resonances, such as the ρ meson.
• Light off-shell window -For R < 2, DM production necessarily involves off-shell dark photons, leading to the reduced sensitivity that is apparent in both panels of Fig. 3.
• R=3 -In much of the recent literature, the reach of various experiments is compared after fixing R = 3. We observe that this is relatively conservative, but does not fully illustrate the resonant freeze-out region near R = 2 or the off-shell freeze-out region of Eq. (16).
In the right-panel of Fig. 4, we show the results of a similar set of calculations for scalar DM, fixing m χ = 10 MeV. To aid the visual comparison between the two models, we show again the results for pseudo-Dirac DM in the leftpanel. We note that the accelerator reach for these two models is very comparable, aside from minor differences when the dark photon is highly off-shell, due to the distinct spin-structure of the DM-dark photon couplings. While the qualitative features of the thermal target are similar for pseudo-Dirac and scalar DM, DM annihilations to SM fermions is p-wave suppressed during freeze-out for scalar DM. As a result, the thermal target for scalar DM is shifted to slightly larger SM couplings.
Separately, in Fig. 5, we show the reach of experiments that depend only on leptonic or hadronic couplings, respectively. Compared to Figs. 3 and 4, in the right panel of Fig. 5 we also include existing constraints from a DM search at MiniBooNE [87] and monojet searches at high-energy proton colliders, such as the Tevatron and the LHC [88]. We choose to present the complementary sensitivity in this way, rather than working with fully UVcomplete models of this type, e.g., those involving mediators such as L µ − L τ or B − 3L τ [21], partly for simplicity, but also because such models entail a number of other more model-dependent constraints. It is important to note that in any specific model, there may be additional constraints on this parameter space that will be relevant. A precise specification of the model is also necessary to determine whether a thermal target lies in unconstrained parameter space.
The sensitivity contours of the electron beam experiments can be understood by analyzing the production channels for each experiment that are relevant as a function of R. The existing constraint from an electron-beam missing-energy search at NA64 is shown as shaded gray in Figs. 3 and 4. The projected sensitivity of a missing momentum search at LDMX is also shown as shaded blue. For m A m χ , the on-shell decay A → χχ is kinematically forbidden, and hence DM production proceeds through an off-shell dark photon, i.e., eN → eN A * → eN χχ. Therefore, the expected number of signal events scales as N signal ∝ α D 2 . Instead, for dark photons above the DM mass-threshold, m A 2m χ , but lighter than the typical center of mass energy of an electron-nucleus collision, m A E cm , DM production proceeds through standard on-shell processes. In this case, if the dark photon has an O(1) branching ratio to DM pairs, the signal rate scales as the on-shell A production rate, i.e., N signal ∝ 2 /m 2 A [15]. For dark photons much heavier than the typical center of mass energy, m A E cm , DM production once again proceeds through off-shell processes, but the rate is now further suppressed by the dark photon mass, N signal ∝ α D 2 /m 4 A . Therefore,  Fig. 4). Along the black lines, the abundance of χ agrees with the observed dark matter energy density. For R 2, the shaded green region corresponds to dark photon-to-dark matter mass ratios for which secluded annihilations dominate over direct annihilations to Standard Model particles. In this case, there is no sharp cosmological target in parameter space. for each of these various mass regimes, for a fixed DM mass the sensitivity of NA64 and LDMX in the y − R plane scales as For NA64 and LDMX, the typical center of mass energy in the electron-nucleus collisions is roughly 5 GeV and 2 GeV, respectively. The R and α D scaling of the different mass regions of Eq. (22) can be seen directly in the behavior of the NA64 and LDMX contours of Figs. 3 and 4. In particular, taking E cm ∼ 5 GeV and m χ = 10 MeV implies that NA64's sensitivity in y should asymptote to a constant value for m A /m χ 500.
The projected sensitivity of the BDX experiment is shown as shaded purple in Figs. 3 and 4. The α D and R scaling in the y − R plane can be understood in a manner very similar to the discussion regarding NA64 and LDMX above. However, the signal yield needs to be supplemented with the scattering rate for χe → χe in the downstream detector. The total cross section for the scattering process scales as σ(χe → χe) ∝ α D 2 (2m e E th +m 2 A ) −1 (2m e E beam +m 2 A ) −1 . E th ∼ 500 MeV is the minimum threshold energy of the detected scattered electron, so that m e E th ∼ (10 MeV) 2 . Similar to Eq. (22), we therefore have For the proton beam fixed-target experiments, CCM, SBND, and NOνA, the sensitivity can be understood in a similar manner to BDX, as the production rate needs to be convoluted with the scattering cross section on either electrons (SBND, NOνA) or nuclei (CCM); thus N signal ∝ y 2 . The dominant production modes for 2m χ < m π/η are via pseudoscalar meson decay, and so N signal ∼ 2 α D for m A < 2m χ , N signal ∼ 2 for m A > 2m χ provided that the dark photon is produced on-shell, or N signal ∼ α D 2 /m 4 A if m A > m π/η and the dark photon is off-shell. The scaling of the scattering cross section is as above for BDX (or replacing m e with m N for nucleon scattering). For SBND and NOνA, the best reach at low mass comes from electron scattering, although the contours on the right of Fig. 5 rely on nucleon scattering. We find that the reach in y is as follows: The primary exception to this scaling is near the m A ∼ m ρ/ω resonance region, where production via proton bremsstrahlung is resonantly enhanced through mixing with vector mesons. This resonant peak is apparent in Figs. 3 For each value of R, mχ, and αD, the kinetic mixing parameter, , is fixed such that the relic abundance of χ agrees with the observed DM energy density. Many of the important features can be inferred directly from Fig. 3. We note that if the pseudo-Dirac mass splitting is sufficiently small, direct detection constraints, e.g., from CRESST-III [13], would provide complementary sensitivity on the right hand side of these plots for mχ 300 MeV. and 4 for SBND and NOνA, but the beam energy is too low at CCM to access the resonant regime. The 1/R 2 scaling for CCM, SBND and NOνA is apparent for m A < m π , and similarly (with reduced reach) for SBND and NOνA for m π < m A < m η . Just above this scale, the reach is enhanced by resonant mixing with the ρ and ω mesons.
Compared to fixed-target experiments, higher energy accelerators, such as LEP, BaBar, and Belle II, are typically most sensitive to thermal relics for larger dark photon masses. This is due to the fact that the center of mass energy for BaBar and LEP is approximately E cm ∼ 10 GeV and E cm ∼ 200 GeV, respectively. Hence, for most of the parameter space that we consider, DM production proceeds through on-shell dark photons, and for O(1) A branching ratios to DM pairs the signal rate is independent of m A or m χ , i.e., N sig ∝ 2 . Therefore, for 2m χ m A E cm , the sensitivity in the y − R plane continues to scale as To further explore light thermal DM, we now turn to a second two-dimensional slice of parameter space, namely the R − m χ plane. In this case, for every value of R, m χ , and α D , we fix the kinetic mixing parameter such that the relic abundance of χ agrees with the observed DM energy density. This analysis requires extensive simulation and is shown for α D = 0.01, 0.1, and 0.5 in Fig. 6 for a representative subset of accelerator experiments. While we focus on pseudo-Dirac DM for concreteness, the results for other DM models such as those involving Majorana or scalar DM are qualitatively very similar (see Fig. 4). In fact, for these other models, the accelerator reach to thermal targets is somewhat enhanced due to the larger couplings needed for adequate freeze-out in the early universe [21].
Although the main features in Fig. 6 can be inferred directly from Fig. 3, we summarize some of the most relevant points below: • Resonant thermal relic gap -In each panel of Fig. 6, there is a gap in sensitivity for R ∼ 2. This is due to the fact that the resonantly enhanced annihilation near freeze-out in the early universe requires much smaller couplings, as shown in the y − R plane of Figs. 3-5. As a result, the sensitivity of any experiment to light thermal DM is inhibited in this parameter regime.
• EFT regime -As discussed in Sec. 6 A, for R 1 and fixed DM mass, both the fixed-target reach and thermal target scale as y ∝ constant. This saturation point occurs once DM annihilation/production is mediated through highly off-shell mediators. As a result, the sensitivity regions of Fig. 6 become insensitive to R for R 1. The figures are truncated (R 30) to focus on the most interesting parameter region. Extending to larger values of R simply pushes further into this EFT regime, which is already excluded by existing limits. Hence, Fig. 6, along with Fig. 1, demonstrates that for m χ few×100 MeV, a generic class of thermal DM models can be generically tested by various fixed-target experiments, independent of the dark sector coupling (α D ) and mediator-to-DM mass ratio (R).
• Resonant production -The projected sensitivity of proton fixed-target experiments approaches that of missingmomentum experiments, such as LDMX, only in special cases, where resonant mixing of the mediator with vector mesons enhances production through proton bremsstrahlung. In Fig. 6, this is apparent for SBND in the region m A ∼ m ρ .

C. Comparison to Direct Detection
Although we began by contrasting the generic validity of the EFT for direct detection with the rich dependence of accelerator-based searches on the mediator properties, we have not explicitly shown constraints from direct detection, nor the prospects from the upcoming generation of electron scattering experiments (see, e.g., Refs. [8][9][10][11][12]14]). For the pseudo-Dirac and scalar models considered in this paper, a relative mass splitting as small as O(10 −6 ) would forbid tree-level DM scattering in direct detection experiments, such that their sensitivity is not phenomenologically relevant. Such mass splittings are not forbidden by an unbroken symmetry, and so there is no generic expectation that they should be small. Because of this model-dependence, we have not shown direct detection constraints on the plots in this section. Nonetheless, because their y-sensitivity is independent of R in the range of interest for simple thermal DM models, it is straightforward to summarize the sensitivity of direct detection experiments for the favorable case of splittings below this level (where elastic tree-level scattering is allowed). In this case, the strongest constraints at our m χ = 10 (100) MeV benchmark points come from XENON10 electron-scattering constraints [9,92,93], which translate to R-independent limit of roughly y > 2 × 10 −10 (10 −7 ) at the two benchmark DM masses. These constraints are not competitive with current accelerator-based limits, except in the region R < 2. Future experiments hope to improve on this sensitivity by four orders of magnitude or more, making them comparable in sensitivity to proposed accelerator-based experiments for models with small mass splittings (allowing elastic scattering) and velocityindependent interactions. Qualitatively similar conclusions hold for all m χ 300 MeV. For DM masses between 300 MeV and a GeV, constraints on DM-nucleon scattering from CRESST III [13] extend somewhat below the BaBar limits in Fig. 6. The strongest complementarity between accelerator-based and direct detection experiments lies in the fact that the former do not rely on the assumptions of velocity-independent interactions or near-zero DM mass splittings, while the latter have excellent sensitivity to models with parametrically light mediators (in our language R 1), which are not relevant for thermal freeze-out but can arise, for example, in DM freeze-in scenarios [94].

CONCLUDING REMARKS
In this paper, we have analyzed the sensitivity of fixed-target experiments to a broad class of sub-GeV thermal relic dark matter models, across a wide range of parameter space, including variations in both the mediator and dark matter mass. This allowed us to assess the reach of existing and planned electron and proton beam fixed-target experiments to light dark matter production through both on-and off-shell mediators. It has become conventional to work with the fixed slice in which dark matter can be produced by on-shell mediator decays, but our results indicate that the reach of fixed-target experiments extends well into the off-shell production regime, including the effective field theory regime where the mediator decouples. We investigated new slices of parameter space, such as the y − R and R − m χ planes (where R = m A /m χ ), in order to illustrate how the reach of different experiments scales as a function of mass, in relation to the thermal target. There is a resonance for m A ∼ 2m χ , which renders annihilation more efficient in the early universe, and our analysis indicates that there currently exists only a small band near this region that will remain untested for simple models coupled through the vector portal.