Testing Light Dark Matter Coannihilation With Fixed-Target Experiments
Eder Izaguirre, Yonatan Kahn, Gordan Krnjaic, Matthew Moschella
FFERMILAB-PUB-17-068-PPD, PUPT 2520
Testing Light Dark Matter Coannihilation With Fixed-Target Experiments
Eder Izaguirre, ∗ Yonatan Kahn, † Gordan Krnjaic, ‡ and Matthew Moschella § Brookhaven National Laboratory, Upton, NY Princeton University, Princeton, NJ Fermi National Accelerator Laboratory, Batavia, IL
In this paper, we introduce a novel program of fixed-target searches for thermal-origin Dark Matter (DM),which couples inelastically to the Standard Model. Since the DM only interacts by transitioning to a heavierstate, freeze-out proceeds via coannihilation and the unstable heavier state is depleted at later times. For suffi-ciently large mass splittings, direct detection is kinematically forbidden and indirect detection is impossible, sothis scenario can only be tested with accelerators. Here we propose new searches at proton and electron beamfixed-target experiments to probe sub-GeV coannihilation, exploiting the distinctive signals of up- and down-scattering as well as decay of the excited state inside the detector volume. We focus on a representative modelin which DM is a pseudo-Dirac fermion coupled to a hidden gauge field (dark photon), which kinetically mixeswith the visible photon. We define theoretical targets in this framework and determine the existing bounds byreanalyzing results from previous experiments. We find that LSND, E137, and BaBar data already place strongconstraints on the parameter space consistent with a thermal freeze-out origin, and that future searches at BelleII and MiniBooNE, as well as recently-proposed fixed-target experiments such as LDMX and BDX, can covernearly all remaining gaps. We also briefly comment on the discovery potential for proposed beam dump andneutrino experiments which operate at much higher beam energies.
I. INTRODUCTION
Understanding the particle nature of Dark Matter (DM) isamong the highest priorities in all of physics. Perhaps themost popular DM candidate to date has been the Weakly In-teracting Massive Particle (WIMP), which is charged underthe electroweak force and naturally yields the observed cos-mological abundance via thermal freeze-out (see [1] for a re-view). However, decades of null results from direct detec-tion, indirect detection, and collider searches have cast doubton this paradigm and motivated many alternative possibilities[2, 3].Nonetheless, the thermal freeze-out mechanism remainscompelling even if DM is not a WIMP. First and foremost,thermal DM is largely UV insensitive; its abundance is de-termined by the DM particle properties and is unaffected bythe details of earlier, unknown cosmological epochs (e.g. re-heating). Furthermore, unlike nonthermal production mech-anisms, which can accommodate DM masses anywhere be-tween − eV − M (cid:12) (!), thermal DM is only viable be-tween ∼ keV − TeV, and is therefore more predictive.Dark matter masses outside this window are either too hot foracceptable structure formation [4] or violate perturbative uni-tarity [5]. Finally, achieving the observed abundance requiresa minimum interaction rate between DM and the StandardModel (SM), which provides a clear target for discovery orfalsification. Thus, there is ample motivation to identify andstudy every viable realization of this mechanism.One simple way to completely eliminate the tension be-tween a thermal origin and experimental limits, in particular ∗ [email protected] † [email protected] ‡ [email protected] § [email protected] those from direct detection experiments, is for the DM to cou-ple inelastically to SM particles [6]. In this class of models,the halo DM species χ is the lightest stable particle in thedark sector and interacts with the SM only by transitioning toa slightly heavier state χ . This class of models has severalappealing features:• Large Viable Couplings:
If the inelastic interac-tion with the SM also determines the leading annihi-lation process, the relic abundance arises dominantlythrough χ χ → SM , a process dubbed coannihila-tion [7]. Since the heavier χ population is Boltzmann-suppressed during freeze-out, the requisite annihilationrate must compensate for this penalty. Thus, the coan-nihilation cross section always satisfies σv ≥ × − cm s − .• Indirect Detection Shuts Off:
Since the heavier stateis unstable, its population is fully depleted at low tem-peratures, so there are no remaining coannihilation part-ners for the χ . This effect turns off possible late-time indirect detection signals and alleviates the boundfrom cosmic microwave background (CMB) power in-jection, which otherwise naively rules out thermal DMfor masses below ∼ GeV for s -wave annihilation[9]. • Direct Detection Forbidden:
For a nonrelativistic haloparticle scattering off a stationary SM target, the energyavailable to upscatter into the heavier state is ∼ µv ,where µ is the reduced mass of the DM-target system. For a general scenario where coannihilation proceeds without inelasticcouplings, but an analogous enhanced thermal cross section appears, seeRef. [8]. Another way to evade this bound is “forbidden DM” [10]. a r X i v : . [ h e p - ph ] A p r Thus, for typical halo velocities v ∼ − , direct detec-tion off of nuclei is forbidden if the mass splitting ex-ceeds (cid:38) keV [6]. There is the possibility of loop-level induced elastic scattering off of electrons, whichmay be very relevant for new proposals for electron di-rect detection [11–16]. However, we leave this questionfor a future study.Since direct and indirect detection search strategies are notavailable, testing thermal coannihilation fundamentally re-quires accelerator-based techniques. For DM masses near theweak scale (few GeV – 100s of GeV), Refs. [17–19] proposedLHC and B-factory searches with sensitivity to thermal coan-nihilation over a wide range of masses and splittings. How-ever, for DM masses below the GeV scale, these searches be-come ineffective and new tools are required to fully test thekeV–GeV mass range over which thermal coannihilation re-mains viable, yet unexplored.In this paper, we fill this large gap by recasting and propos-ing a series of fixed-target searches for both electron and pro-ton beam facilities. In both cases, an incident beam impingeson a fixed target and produces a boosted χ χ pair. Depend-ing on the experimental setup, this system can give rise to avariety of possible signals:• Beam Dumps:
The χ χ pair can be produced ra-diatively through the “dark bremsstrahlung” reactions pN → pN χ χ at proton beam dumps, or eN → eN χ χ at electron beam dumps, where N is a nu-clear target. At proton beam dumps, the DM system canalso be produced from meson decays π → γχ χ and η → γχ χ , which can be advantageous for low DMmasses. Once produced, the χ , can reach a down-stream detector and scatter off electrons or nucleons toinduce an observable recoil signature. Alternatively,the unstable χ can also decay in flight as it passesthrough the detector via χ → χ e + e − , thereby de-positing an observable signal. These processes are de-picted schematically in Fig. 2.• Electron Missing Energy/Momentum:
As above, the χ χ pair is produced in e − N → e − N χ χ “darkbremsstrahlung”, but the production takes place in athin target embedded in a detector. The principal ob-servable signature in this context is the recoiling final-state electron. If this electron emerges having lost mostof its incident beam energy, with no additional activitydeposited in a downstream detector, this process can besensitive to DM production. As above, a χ decayinginside the detector provides an additional potential sig-nature. This process is depicted schematically in Fig. 3.We will show that existing data already rules out large por-tions of the direct coannihilation parameter space. Moreover,the dedicated searches we propose which exploit the uniquesignals from this class of models can significantly improve onthe current levels of sensitivity. Crucially, as shown in Fig. 4,the χ has a macroscopic decay length over nearly all of theparameter space we are interested in, and detecting the elec-tromagnetic energy deposited by a χ decay in the detector a )FIG. 10: a) Scalar DM pair production in electron-nucleuscollisions. An on-shell A is radiated and decays o↵ diago-nally to ' h,` pairs. b) Inelastic up scattering of the lighter ' ` into the heavier state via A exchange inside the detector.For order-one (or larger) mass splittings, the metastable statepromptly de-excites inside the detector via ' h ! ' ` e + e .This process yields a target (nucleus, nucleon, or electron)recoil E R and two charged tracks, which is a instinctive, zerobackground signature, so nuclear recoil cuts need not be lim-iting. A f + f A T T f f + A FIG. 1. Leading order diagram for χ χ → f + f − coannihilation,which sets the DM relic abundance in the m A (cid:48) > m , regime. — a signal not present in elastic DM models — is sufficientto cover large portions of the thermal relic curve. Indeed, weshow in Sec. III C that the sensitivity to the decay of the exitedstate χ typically dominates over scattering channels at exper-iments with beam energies below 10 GeV. Thus, the prospectsare excellent for dedicated experiments sensitive to these sig-natures, and we show that current and planned experimentscan confirm or rule out nearly the entire mass range allowedfor thermal coannihilating DM. We note a related study fromRef. [21] that investigated previous limits from fixed target fa-cilities in the context of a supersymmetric hidden sector. Re-cently, Ref. [22] proposed the signal of χ → χ upscatteringfollowed by χ decay for the case of (boosted) astrophysicalDM, as opposed to DM produced in a beam dump.This paper is organized as follows. In Sec. II we introducea representative renormalizable model featuring a pseudo-Dirac fermion DM coupled to a kinetically-mixed dark pho-ton, where the abundance of the former arises from thermalcoannihilation. In Sec. III we make some general commentson the production modes and detection signals at proton-and electron-beam fixed-target experiments. In Sec. IV wedescribe our methods for reinterpreting existing data fromLSND [25] and E137 [26], and we discuss projections forthe future data at MiniBooNE [27] and the proposed BDX[28] and LDMX [23] experiments. In Sec. V we discuss thebounds and reach projections from these experiments in thecontext of the thermal coannihilating inelastic DM. Finally, inSec. VI we offer some concluding remarks. Additional detailson the matrix elements and Monte Carlo methods used in de-termining the thermal parameter space and in our simulationscan be found in the Appendices. II. SUB-GEV THERMAL COANNIHILATION
In this section, we describe a class of models of coannihi-lating DM: DM that couples inelastically to the SM througha kinetically-mixed dark photon. We detail the early universecosmology and freeze out of the model, as well as introducea useful parametrization of the parameters of the model in A similar search strategy was proposed for the decay of long-lived scalarsin Ref. [20]. Beam p ! Target Detector i A ⇡ ,⌘ A i j T T f f + A and/or Beam e/p ! Target/Dump Detector i A Ze/p e/p p Z A ⇡ ,⌘ A i j T T and/or f f + A FIG. 2. Inelastic DM production at electron and proton beam dump experiments via dark bremsstrahlung and meson decay. The resulting χ , χ pair can give rise to a number of possible signatures in the detector: χ can decay inside the fiducial volume to deposit electromagneticenergy; both χ and χ can scatter off detector targets T and impart visible recoil energies to these particles; or χ can upscatter into χ ,which can then decay promptly inside the detector to deposit a visible signal. e ! ECAL/HCALTarget Tracker e Invisible e ! Active Target (ECAL/HCAL) e Invisible A Ze e FIG. 3. Inelastic DM production at electron beam fixed-target missing energy/momentum experiments.
Left:
Setup for an LDMX stylemissing momentum experiment [2, 23] in which a ( ∼ few GeV) beam electron produces DM in a thin target ( (cid:28) radiation length) and therebyloses a large fraction of its incident energy. The emerging lower energy electron passes through tracker material and registers as a signal eventif there is no additional energy deposited in the ECAL/HCAL system downstream, which serves primarily to veto SM activity. Right:
Setupfor an NA64 style experiment in which the beam (typically at higher energies, ∼ GeV) produces the DM system by interacting with aninstrumented, active target volume [24]. As with LDMX, the instrumented region serves to verify that the beam electron has abruptly lost mostof its energy and that there is no additional SM activity downstream. which the thermal target is largely an invariant under varia-tion of couplings and of mass hierarchies.
A. Mediator Model Building
Unlike weak-scale WIMPs, which realize successfulfreeze-out with only SM gauge interactions, sub-GeV DM isoverproduced in the absence of light ( (cid:28) m Z ) new mediatorsto generate a sufficiently large annihilation rate [29, 30]. Toavoid detection thus far, such mediators must be neutral underthe SM and couple non-negligibly to visible particles.If SM particles are neutral under the new interaction, a renormalizable model (without additional fields) requires themediator to interact with the SM through the hypercharge,Higgs, or lepton portals B µν , H † H , LH, (1)for vector, scalar, and fermionic mediators, respectively.However, coupling a fermionic mediator to the lepton por-tal requires additional model building and scalar mediators,which mix with the Higgs are ruled out for predictive mod-els in which DM annihilates directly to SM final states (see A fermionic mediator coupled to the lepton portal requires additional Δ = . m Δ = . m Δ = . m
10 10 - - m [ MeV ] c τ [ y r e li c / y ][ c m ] Rest Frame Decay Length χ → χ f f , y = y relic FIG. 4. Rest frame χ decay length for various mass splittings ∆ ≡ m − m , as a function of the χ mass m . The vertical axis is nor-malized to values of the dimensionless coupling y ≡ (cid:15) α D ( m m A (cid:48) ) which achieve the observed relic density; see Sec. II for more details. Sec. II C and [31] for a discussion of this issue), so we restrictour attention to abelian vector mediators; a nonabelian fieldstrength is not gauge invariant, so kinetic mixing is forbidden.Alternatively, the mediator could couple directly to SMparticles if both dark and visible matter are charged underthe same gauge group. In the absence of additional fields,anomaly cancellation restricts the possible choices to be U (1) B − L , U (1) (cid:96) i − (cid:96) j , U (1) B − (cid:96) i , (2)and linear combinations thereof. In most contexts, the rele-vant phenomenology in fixed-target searches is qualitativelysimilar to the vector portal scenario, so below we will ignorethese possibilities without loss of essential generality. Wenote, however, that viable models for both protophobic [32]and protophilic [33] mediators exist, so the complementarityprovided by both proton- and electron-beam experiments ishighly advantageous. B. Representative Model
Our representative dark sector contains a 4-componentfermion ψ that transforms under a hidden abelian U (1) D gauge group. The fermion couples to a vector mediator A (cid:48) as L = iψ /Dψ + M ¯ ψψ + λφ ψ c ψ + h.c., (3) model building to simultaneously achieve a thermal contact through thisinteraction and yield viable neutrino textures; the coupling to the mediatormust be suppressed by neutrino masses, so it is generically difficult for theinteraction rate to exceed Hubble expansion. where φ is a U (1) D symmetry breaking scalar whose vacuumexpectation value v D gives A (cid:48) a nonzero mass m A (cid:48) ∼ g D v D and gives ψ a Majorana mass ∼ λv D . Diagonalizing the re-sulting Dirac and Majorana masses gives rise to fermion masseigenstates χ , with a small mass splitting ∆ ≡ m − m and an off-diagonal coupling to A (cid:48) , L ⊃ g D A (cid:48) µ χ γ µ χ + h.c., (4)where g D ≡ √ πα D is the dark coupling constant. Note thatit is technically natural to have ∆ (cid:28) M since the Majoranamass breaks the global symmetry associated with ψ number. This sector can interact with the SM through a renormaliz-able and gauge-invariant kinetic mixing term between U (1) D and U (1) Y gauge fields, L ⊃ (cid:15) θ W F (cid:48) µν B µν = (cid:15)F (cid:48) µν F µν + (cid:15) tan θ W F (cid:48) µν Z µν , (5)where (cid:15) (cid:28) is the kinetic mixing parameter and F (cid:48) µν and B µν are respectively the dark and hypercharge field strengthtensors and the kinetic mixing interaction has been written interms of the SM mass eigenstates A and Z after electroweaksymmetry breaking. Diagonalizing the kinetic terms in Eq. (5)and rescaling the field strengths to restore canonical normal-ization induce a coupling between A (cid:48) and the SM fermions[34]. To leading order in (cid:15) , the A (cid:48) -SM interaction becomes eA µ J µEM → e ( A µ + (cid:15)A (cid:48) µ ) J µEM , (6)where J µEM is the SM electromagnetic current and all chargedfermions acquire millicharges under U (1) D . There is alsoan analogous A (cid:48) interaction with the SM neutral current thatarises from A (cid:48) − Z mixing, but in our mass range of inter-est, m A (cid:48) (cid:28) m Z , the mixing parameter is suppressed by anadditional factor of ( m A (cid:48) /m Z ) [35–38], so we neglect thisinteraction for the remainder of paper. C. Direct Coannihilation vs. Secluded Annihilation
In the hot early universe ( T (cid:29) m i , m A (cid:48) ), all dark speciesare in chemical and kinetic equilibrium with the SM radia-tion bath; this initial condition is guaranteed as long as theDM-SM scattering rate exceeds the Hubble expansion rateat some point during cosmic history. If m i > m A (cid:48) , thefreeze-out process is analogous to that of WIMP models. Be-low the freeze-out temperature T f ∼ m , / , the numberdensities of both species are depleted predominantly through χ i χ i → A (cid:48) A (cid:48) self-annihilation (which depends only on α D ), not coannihilation, which depends on the combination (cid:15) α D and is greatly suppressed by comparison. Although both com-ponents undergo freeze-out separately, since χ is heavier and If, unlike the construction in Eq. (3), the Majorana masses for the two Weylcomponents in ψ = ( ξ, η † ) are different, there is also a subleading diago-nal interaction of the form ( δ/M D ) χ i (cid:54) A (cid:48) χ i , where δ ≡ m ξ − m η is thedifference of Majorana masses for the the interaction eigenstates. We ne-glect this interaction in our analysis, assuming the off-diagonal interactiondominates. unstable, it will be depleted through downscattering and de-cays. Thus, up to order-one corrections, the requisite self-annihilation cross section satisfies the familiar WIMP-like re-quirement (cid:104) σv (cid:105) ∼ × − cm s − in order for χ to havethe observed abundance at late times.However, this secluded ( m i > m A (cid:48) ) regime has severaldrawbacks. Since the self-annihilation rate for fermions is s -wave, annihilation continues to occur out of equilibrium dur-ing recombination, which ionizes newly-formed hydrogen andthereby modifies the CMB power spectrum. For a thermal an-nihilation rate, this bound rules out DM below ∼ GeV[9]. Furthermore, since the secluded annihilation cross sec-tion scales as σv ∼ α D /m i , the abundance is independentof the A (cid:48) coupling to SM states, so there is no minimum in-teraction strength target for the DM search effort — at directdetection or accelerators, since these are sensitive to (cid:15) .By contrast, in the direct coannihilation regime ( m i The annihilation rate depends crucially onthe mixing with the SM, σv ∝ (cid:15) α D , so for dark cou-plings that satisfy perturbative unitarity, there is a min-imum value of (cid:15) that is compatible with thermal freeze-out. Reaching experimental sensitivity to this minimumvalue for each viable DM mass suffices to discover orfalsify this entire class of models.• Large Cross Section: Unlike secluded annihilation,which involves the annihilation of equal mass particles,direct coannihilation requires both χ and χ in the ini-tial state. However, χ typically becomes Boltzmann-suppressed before the nominal χ freeze-out temper-ature ∼ m / , so the coannihilation cross sectionmust be larger to compensate for this reduction. Thus,we generically require (cid:104) σv (cid:105) (cid:29) × − cm s − to achieve the observed abundance, where the precisevalue grows sensitively with increasing ∆ (compare therelic density curves in Fig. 6).• Bounded Viable Range: Since the necessary coanni-hilation rate for freeze-out grows sharply as the masssplitting is increased, for sufficiently large ∆ , i.e., ∆ ∼ m , the requisite (cid:15) is easily excluded inde-pendently of other model properties (e.g. by precisionQED/electroweak measurements). Thus, as we will seebelow, this class of coannihilating models can be testedover the full viable parameter space.Thus, for the remainder of this paper, our focus will be on thedirect coannihilation scenario. For simplicity, without loss ofessential generality, we will also demand that m A (cid:48) > m + m so that A (cid:48) can decay to dark sector states — otherwise, the A (cid:48) decays to SM final state, a signature for which there are If instead, the DM is a scalar and annihilates directly to SM fermionsthrough an s -channel vector mediator, its annihilation rate is p -wave sup-pressed, which can evade CMB bounds. abundant ideas to cover [2]. Since (cid:15) (cid:28) , the A (cid:48) branchingratio to SM states is correspondingly negligible in this regime.Note that if we had chosen a neutral scalar mediator in-stead of the vector A (cid:48) , it could only couple to SM fermions bymixing with the Higgs after electroweak symmetry breaking.However, such a scenario cannot viably realize thermal coan-nihilation below the GeV scale because the requisite Higgs-mediator mixing angle must be ∼ O (1) to overcome theYukawa penalty in the annihilation diagram and yield a ther-mal annihilation rate. Such large couplings are ruled out byHiggs coupling measurements and rare meson decay searches[31, 39]. We note that t -channel annihilation into a light scalarmediator can still realize secluded annihilation, but this pro-cess is less predictive and beyond the scope of this work. D. Covering the Thermal Target The goal of this paper is to compute existing bounds on theparameter space for thermal freeze-out via coannihilation andpresent sensitivity projections for future experiments. Onemajor challenge of such an effort is the large dimensional-ity of this parameter space: each viable model point is definedby a unique choice of the inputs { m , m A (cid:48) , ∆ , (cid:15), α D } , con-strained by the requirement that the coannihilation rate yieldsthe observed DM density. Fortunately, in the nonrelativisticlimit, for each choice of ∆ , the coannihilation rate scales as σv ∝ (cid:15) α D m m A (cid:48) = (cid:18) (cid:15) α D m m A (cid:48) (cid:19) m ≡ ym , (7)which is valid for s (cid:39) ( m + m ) and sufficiently away fromthe s -channel resonance at m A (cid:48) ∼ m + m . Here we havedefined the dimensionless parameter y ≡ (cid:15) α D (cid:18) m m A (cid:48) (cid:19) , (8)which uniquely determines the freeze-out annihilation rate fora given choice of m and ∆ .The virtue of this parameterization is that the observed relicdensity is achieved only along a one-parameter curve y ( m ) which is insensitive to the relative sizes of (cid:15), α D , or m /m A (cid:48) ,reducing the dimensionality of the parameter space. However,the drawback is that some experiments only constrain a subsetof model parameters and are not sensitive to the same combi-nation of inputs that define the y variable in Eq. (8). For ex-ample, lepton colliders can be used to constrain (cid:15) as a functionof m A (cid:48) in searches for e + e − → γ (cid:54) E , which can be interpretedas e + e − → γA (cid:48) production followed by prompt A (cid:48) decay toinvisible particles [40–42]. Since Br( A (cid:48) → invisible) (cid:39) ,the signal event yield depends on (cid:15) , but is insensitive to α D orthe ratio m /m A (cid:48) for larger χ lifetimes.Nonetheless, it is still possible to use this information toconstrain the y variable for a conservative comparison with thethermal target. The strategy is to construct the largest possible y exp using the experimentally determined (cid:15) exp , y exp = (cid:15) × (cid:34) α D (cid:18) m m A (cid:48) (cid:19) (cid:35) , (9) E L S ND R e li c D e n s i t y B DX - - - - - - α D ϵ Decay Sensitivity, m = 15 MeV, m A' = m , Δ = m E L S ND R e li c D e n s i t y B DX M i n i B oo N E - - - - - α D ϵ Decay Sensitivity, m = 50 MeV, m A' = m , Δ = m FIG. 5. Decay reach in (cid:15) as a function of α D , with other model parameters held fixed, for various fixed-target experiments. The value of (cid:15) corresponding to the thermal target y for each α D is shown for comparison in black. The most conservative choice for α D relative to thethermal target is the largest α D consistent with perturbativity. Note that the different slope of the LSND curve at 50 MeV compared to the otherexperiments is a consequence of off-shell A (cid:48) decay, see Sec. III A. Also note that the MiniBooNE projection from the right panel is missingfrom the left panel because the mass splitting on the left is 1.5 MeV, which is below the threhsold to pass the experimental cut in Eq. (25). SeeSec. IV for further details on computing the reach. where the quantity in square brackets is chosen to be as largeas possible while remaining consistent with both perturba-tive unitarity and the requirement that direct coannihilationremains the dominant annihilation process. Thus, the con-servative prescription is to adopt order-one values for α D and m /m A (cid:48) , while taking m A (cid:48) > m + m and avoid-ing the m A (cid:48) ∼ m + m resonance. As an illustration,Fig. 5 shows the reach in (cid:15) for χ decay at various experi-ments, demonstrating that the weakest reach with respect tothe thermal target occurs at large α D . For our numerical re-sults in Figs. 6 and 8, we chose representative benchmarks α D = 0 . , m /m A (cid:48) = 1 / where appropriate (e.g. for col-lider bounds). In Fig. 7 and in the discussion in Sec. V below,we show how the constrains on the parameter space change fordifferent choices of these benchmark values, and demonstratethat our conclusions are qualitatively unchanged.Note that for our main results, which are discussed in Sec. Vand presented in Figs. 6, 7, and 8, we compute the y ( m ) curve by numerically solving the full Boltzmann system forthis model as described in Appendix A. III. PRODUCTION AND DETECTION BASICS In this section we review the basic A (cid:48) production formalismin proton- and electron-beam fixed-target collisions to developsome intuition for our later numerical results. We also presentthe formalism for various DM detection signatures that ensuefrom boosted A (cid:48) → χ χ decays. A. Production at Proton Beam Dumps At proton beam dumps, χ χ pairs can be produced eitherfrom dark bremsstrahlung pN → pN A (cid:48) , A (cid:48) → χ χ or fromthe neutral meson decays m → γA (cid:48) ( ∗ ) → γχ χ where m = π , η . The meson production rates and spectra dependon various experimental factors, but the strongest dependenceis on the proton beam energy, while the detailed material prop-erties of the beam dump can largely be ignored. At few-GeVbeam energies, N π ≈ N p and N η ≈ N p / [43] where N p is the total flux of protons. See Ref. [44] for a detailed discus-sion of production modes and distributions.If m A (cid:48) < m m , the A (cid:48) can be produced on-shell and thenumber of χ χ pairs via meson decays is given by N χ χ = N m × (cid:15) (cid:18) − m A (cid:48) m m (cid:19) BR( m → γγ ) , (10)where N m is the number of π or η produced. In either case,there is a kinematic threshold for production through mesons: m + m = 2 m + ∆ < m π ≈ 135 MeV for pions, and m + ∆ < m η ≈ 548 MeV for etas. At sufficiently high energies, direct A (cid:48) production from quarks and gluonsin the target becomes relevant, but this process is negligible for the beamenergies we consider in this work. If m A (cid:48) > m m , the meson decay proceeds through an off-shell A (cid:48) andthe rate scales as (cid:15) α D [45], in contrast to on-shell production which isindependent of α D . At low-energy experiments such as LSND with onlya single production channel, off-shell production may be important, but forhigher-energy experiments such as MiniBooNE with multiple productionchannels, on-shell production is dominant. The rate of dark bremsstrahlung production is proportionalto (cid:15) , but varies non-trivially with the dark photon mass m A (cid:48) .Using the simulation of Ref. [44], we find that the numberof A (cid:48) produced varies from ∼ − (cid:15) N p as m A (cid:48) → to ∼ − (cid:15) N p at m A (cid:48) = 3 GeV . A (cid:48) production can be enhancedthrough vector meson mixing near m A (cid:48) = m ρ , m ω , but sincethis enhancement is only in a limited mass range and dependson the precise choice of m /m A (cid:48) , we do not attempt to modelit precisely. The dark bremsstrahlung production mechanismhas no mass threshold and so heavy DM can still be producedup to the beam energy, m + ∆ < E beam . B. Production at Electron Beams At electron-beam experiments, mesons are no longer co-piously produced so the dominant process is A (cid:48) productionthrough dark bremsstrahlung followed by on-shell mediatordecay A (cid:48) → χ χ . The reaction eN → eN A (cid:48) , where N isa nucleus of atomic (mass) number Z ( A ) has been well stud-ied in Refs. [46–48]. The relevant features of the reaction canbe better illuminated by considering the calculation using theWeizsacker-Williams (WW) approximation, although we notethat all plotted results in this study employ a numerical calcu-lation; see Appendix D for more details of our simulations.In the WW approximation, the differential production crosssection of A (cid:48) is then given by dσdx d Ω A (cid:48) ≈ α (cid:15) Φ( q , q f ) E xπU (cid:20)(cid:18) − x + x (cid:19) − x (1 − x ) m A (cid:48) ( E θ A (cid:48) + m e ) U (cid:21) , (11)where d Ω A (cid:48) = 2 πd cos θ A (cid:48) is the A (cid:48) solid angle with respectto the beam axis in the lab frame and we have defined U ( x, θ A (cid:48) ) = E xθ A (cid:48) + m A (cid:48) − xx + m e x, (12) q = UE (1 − x ) , q f = m A (cid:48) . (13)Here E is the electron beam energy, x = E A (cid:48) /E the frac-tion of the beam energy carried by the A (cid:48) , and α Φ( q , q f ) /π is the WW photon flux for small-virtuality photons with mo-mentum − t bounded by q and q f [46]. For q values muchsmaller than the inverse nuclear size ≈ . GeV /A / , wehave Φ( q , q f ) ∼ Z up to an overall logarithmic factor.Ref. [47] found that for any given x , the angular integral isdominated by angles θ A (cid:48) such that Exθ A (cid:48) (cid:46) m A (cid:48) − xx + m e x .Using this approximation, we can derive a simpler expressionfor the differential cross section dσdx = 4 α (cid:15) Φ( m A (cid:48) , E ) x + 3 x (1 − x )3 m A (cid:48) (1 − x ) , (14)where Φ( m A (cid:48) , E ) ≡ Φ( q = m A (cid:48) / (2 E ) , q f = m A (cid:48) ) . Af-ter integrating, the total A (cid:48) production cross-section scalesroughly as σ ( eN → eN A (cid:48) ) ≈ α (cid:15) m A (cid:48) Φ( m A (cid:48) , E ) (cid:2) log δ + O (1) (cid:3) , (15)where δ ≡ min( m A (cid:48) /E , m e /m A (cid:48) , m e /E ) . C. Generic Detection Signals (Electron & Proton Beams) In our regime of interest, namely (cid:15) α < α D , A (cid:48) de-cays promptly primarily into the χ χ system, imparting anorder-one fraction of its energy to the daughter particles,which emerge from the target with large boosts. The boostedDM system gives rise to three classes of observable signa-tures: detector target scattering, χ decays, and missing en-ergy/momentum carried away by the DM system. 1. Detector Target Scattering The χ produced in the dump travels unimpeded, enters thedetector situated downstream of the target/dump, and scattersoff of target particles in the detector. It scatters via the re-action χ T → χ T , where T here could in principle be anelectron, a nucleon, or even a nucleus. Similarly, any popula-tion of surviving χ s produced in the target could downscatterthrough the reverse reaction χ T → χ T . The cleanest sig-nal occurs when T is an electron, with energy typically above10s of MeV. The production rate is proportional to (cid:15) and thescattering rate is proportional to (cid:15) α D , so the total yield inthis channel will be proportional to (cid:15) α D .Specializing to the case of electron targets in the E χ (cid:29) ∆ regime, the approximate differential cross section is [48, 49] dσdE e = 4 π(cid:15) αα D m e E χ − f ( E e )( E e − m e )( E χ − m χ )( m A (cid:48) + 2 m e E e − m e ) , (16)where E e is the electron recoil energy, E χ is the incident χ , energy, and we define f ( E e ) = 2 m e E χ − m e E e + m χ + 2 m e . (17)This cross section is valid up to corrections of order ∆ /E χ and ∆ /m A (cid:48) , both of which are very small compared for thebenchmarks we consider throughout (see Figs. 6, 8, and 7).However, our numerical results evaluate the exact expressionpresented following the derivation in App. B 2.Neglecting subleading corrections and integrating the recoilenergy up to E χ gives the approximate result σ χe ≈ παα D (cid:15) m e E χ m A (cid:48) , (18)so the corresponding scattering probability inside the detectorbecomes P scatter = n e σ χe (cid:96), (19)where (cid:96) is the DM path length through the detector and n e isthe electron number density. There are similar expressions forscattering off detector nucleons and nuclei, but, as we will see,most of the relevant bounds and projections below exploit theelectron channel. 2. Decay of excited state One of the most powerful channels at beam dump exper-iments is the direct decay of χ → χ e + e − , whose partialwidth satisfies Γ χ ≈ (cid:15) αα D ∆ πm A (cid:48) (20)in the m A (cid:48) (cid:29) m (cid:29) m e limit. The χ is produced in thebeam dump and decays in a downstream detector (depictedschematically in Fig. 2), so the signal yield scales as N signal ∝ (cid:15) P survive P decay , (21)where the product P survive P decay = e − Γ χ d/βγ (cid:16) − e − Γ χ (cid:96)/βγ (cid:17) (22)is the probability for χ to survive a distance d out to do thedetector and decay inside after traversing a path length (cid:96) withboost factor γ and velocity β = v/c . As with scattering, thetotal decay yield scales as (cid:15) α D , with (cid:15) coming from pro-duction and (cid:15) α D coming from expanding the exponentialsin Eq. (22) assuming the long lifetime limit (see Fig. 4) andusing Eq. (20).To estimate the relative reach of decay and scatteringsearches at beam dump experiments, it is useful to define R ≡ P decay P scatter (cid:39) ∆ π γE χ m e n e , (23)where n e is the number density of detector electrons. Herewe have used the approximate χ -electron scattering cross sec-tion from Eq. (18); note that this expression is independent ofdetector geometry or DM production rate.For most materials, n e ∼ cm − , so we have R (cid:39) (cid:16) m MeV (cid:17) (cid:18) ∆0 . m (cid:19) (cid:18) GeV E χ (cid:19) , (24)where E χ represents the typical energy of a χ within the de-tector acceptance. Note that the R ∼ ∆ scaling implies thatdecays dominate the reach for most of the splittings shown inFig. 6, and that for experiments with lower DM energies (e.g.at LSND where E χ ∼ few 100 MeV) the decay reach domi-nates the scattering reach by several orders of magnitude dueto the lower beam energy and smaller boost factors. 3. Upscatter followed by decay A third possible signal combines the phenomenology ofboth scattering and decay. χ can upscatter at the detector, For ∆ > m µ , χ can also decay to muons, but for the majority of theparameter space we consider in this paper, only the electron channel isallowed. producing a χ along with a recoiling target T . If the χ pro-duced in the upscatter decays inside the detector via the re-action χ → χ e + e − , and if the electron and positron finalstates are energetic and separated enough, the final state leadsto a signature that is not easily mimicked by environmentalor beam-related backgrounds [28, 50, 51]. However, we findthat the yields in this channel are subdominant to the othertwo and we will not consider it further. The consequencesof this signal for the case of boosted astrophysical DM, wherethe different kinematics and the absence of an (cid:15) productionpenalty make this channel an attractive background-free op-tion, were investigated in detail in Ref. [22]. 4. Missing Energy/Momentum A unique feature of A (cid:48) production in electron-beam fixed-target collisions is that the A (cid:48) is typically radiated with anorder-one fraction of the incident beam energy (see Sec. III Band more detailed discussions in [23, 24, 46, 47]). This en-ables a novel detection strategy where the target is now em-bedded in the detector, and which is based on comparing theenergy and momentum of the beam electron before and afterit undergoes dark bremsstrahlung. Fig. 3 illustrates a set-upfor this detection strategy. D. Kinematic Thresholds There are several kinematic thresholds which influencethe detectability of the various signals. Clearly, for ∆ < m e ≈ , the excited state χ can no longer decayvia χ → χ e + e − . The only possible decay mode remain-ing is to χ plus neutrinos or γ , both of which are suffi-ciently suppressed that χ is cosmologically long-lived [52].In this case the only signals are from the recoiling target in χ upscattering or χ downscattering, along with missing en-ergy/momentum. Otherwise, for detectors with finite energythresholds and angular resolution, we might require the de-cay signal χ → χ e + e − to have an energetic, well-separated e + e − pair satisfying E e > E min and θ e + e − > θ min . The re-quirement on ∆ to allow a visible χ decay is actually morestringent: ∆ ≥ m e (1 + cos θ min ) + 2 E (1 − cos θ min ) . (25)For example, requiring E min = 50 MeV and θ min = 2 ◦ as atMiniBooNE, we find that the minimum splitting we can probeis ∆ (cid:38) . IV. EXISTING BOUNDS AND PROJECTIONS In this section we discuss the key features of the variousrepresentative experiments we consider and list our kinematic If the recoiling target T is not visible, the signal in this channel is identicalto the decay signal, but with the added (cid:15) α D penalty of scattering. cuts used to compute reach curves. The bounds and projec-tions we derive are presented in Figs. 6, 7, and 8. A. Signals at Proton Beam Dump Experiments 1. LSND The Liquid Scintillator Neutrino Detector (LSND) [53]was a neutrino experiment at Los Alamos which ran from1993-1998. The extremely high-luminosity proton beam pro-duced the largest available fixed-target sample of neutral pi-ons, N π ∼ . The LSND proton beam had kineticenergy 800 MeV, small enough that η production and darkbremsstrahlung are negligible. Therefore, the DM is produceddominantly from π decays. We modeled the π production atLSND using the GEANT simulation from Ref. [45] and thendecayed the pions to DM using Monte Carlo. Due primarilyto the large luminosity, we expect DM signal event yields atLSND to dominate in the regions of parameters space wherepion decay is kinematically allowed.The LSND detector was an approximately cylindrical tankof scintillating mineral oil with length ∼ and radius ∼ , located ∼ from the beam stop. The close proxim-ity to the beam stop and its off-axis orientation gives a largegeometric acceptance of tracks originating at the beam stop.The detector used photomultiplier tubes to detect Cerenkovlight emitted by charged particles in the scintillator. For thepurposes of our simulations, LSND is only capable of identi-fying lepton tracks but is blind to nucleons, and electrons areindistinguishable from positrons. We take the electron detec-tion efficiency at LSND to be [25].We numerically computed the event yields at LSND usingDM events that were produced as outlined in Sec. III A, foreach of the available signal channels in Sec. III C. Using thetechniques outlined in App. D 2, we simulated the χ and χ scattering off electrons, nucleons, and nuclei as well as direct χ decays in the detector. This simulation accounted for thegeometric acceptance of the detector as well as kinematic cutsand detection efficiencies for the observable final state parti-cles.Previous dark matter limits at LSND have been derived us-ing the 55-event background-subtracted limit from the LSNDneutrino-electron elastic scattering analysis [25, 54]. To de-rive the LSND decay constraints for inelastic DM, we inter-preted the 55-event limit using χ decay events where the e + /e − opening angle was too small to distinguish the twoleptons, and thus could fake an elastic event passing the cutsdescribed in Ref. [25]. The angular resolution of LSND is ◦ for electrons above 20 MeV [53], so we require thatthe angular separation of the e + e − pair is θ e + e − < ◦ and that the total energy of the pair lies in the range 50 MeV with either the electron or positron sat-isfying cos θ < . . This is conservative, as we have ignoreda similar number of events where the electron and positronare well-separated; with access to the full LSND data set, onecould search for events where two leptons were seen simulta-neously inside the detector and potentially improve the reach. Indeed, we chose such a conservative prescription because, aswe will see, the limits derived from LSND are already domi-nant for ∆ > m e , m + ∆ < m π .For the LSND scattering constraints, we again use the 55-event limit applied to any event with at least one leptonpassing the elastic cuts in the final state. This includes χ down-scattering and χ upscattering with and without χ de-cay, though the dominant process for all but the largest DMmasses is χ downscattering, for which the signal is identicalto neutrino-electron elastic scattering. We require a recoilingelectron track to be forward-oriented with cos θ > . andhave energy between < E e < 50 MeV . Note that becauseLSND does not have particle ID, a potential background isneutral-current π production where the photons from the π decay fake electrons. 2. MiniBooNE MiniBooNE [55] is a proton beam neutrino experiment cur-rently operating at Fermilab with beam energy 8 GeV. It hasbeen previously noted that MiniBooNE is sensitive to lightDM [54, 56], and indeed, the first limits from a dedicated DMsearch have recently been published [57]. At MiniBooNE,the beam energy is large enough that η production and darkbremsstrahlung contribute to DM production, however, con-tributions from π decays dominate in regions of parameterspace where the π → γχ χ is allowed. Due to the smallerluminosity, the number of pions produced at MiniBooNE isonly N π ∼ , with the number of etas being further sup-pressed by a factor of ∼ . We used BdNMC [44] to gen-erate a distribution of π and η produced at MiniBoone andthen decayed the mesons to DM using Monte Carlo. To simu-late dark bremsstrahlung production at MiniBooNE, we used BdNMC to generate a distribution of on-shell A (cid:48) and then de-cayed the dark photons to DM using Monte Carlo. The rateof dark bremsstrahlung production decreases as m A (cid:48) is in-creased, but for m A (cid:48) (cid:38) 200 MeV , kinematic thresholds and α D suppression from off-shell meson decay are significantenough that the dark bremsstrahlung production of DM domi-nates. A more detailed discussion of the Monte Carlo methodsand matrix elements used in these simulations can be found inAppendices C and D.The MiniBooNE detector also uses scintillating mineral oil,and is approximately spherical with radius ∼ , located ∼ from the beam stop. The smaller detector and largerdistance from the beam stop means that the geometric accep-tance is smaller than LSND, but because the beam energy ofMiniBooNE is ten times that of LSND, the DM producedat MiniBooNE has boost factors that are roughly ten timesgreater, which to some extent compensates for the geomet-ric acceptance as the DM is more collimated. That said, thecombination of the smaller proton luminosity and the largernumber of background events for a DM search [57] results inthe MiniBooNE scattering reach being suppressed by severalorders of magnitude compared to LSND. We do not includethe MiniBooNE scattering curves in our reach plots as they do0not cover new parameter space. MiniBooNE is capable ofseeing both nucleon and lepton tracks, but electrons are indis-tinguishable from positrons and neutrons are indistinguishablefrom protons. We take the efficiency for lepton and nucleondetection to be [60].At MiniBooNE, we impose energy and angular cuts sim-ilar to those in Ref [57]. We require any lepton track to beforward-oriented with cos θ > . and have energy in therange < E e < 600 MeV . We also require the electronand positron to have an angular separation of at least ◦ [60]and an energy of at least 50 MeV to be visible. We deter-mined the decay reach assuming 95% one-sided Poisson c.l.,corresponding to 3 events. B. Signals at Electron Beam Dump Experiments 1. E137 The E137 experiment at SLAC [26] was designed to searchfor axion-like particles and with a 20 GeV beam delivering ∼ 30 C ( ∼ × electrons) of current onto an aluminum targetpositioned approximately 400 m upstream from an aluminumand plastic scintillator detector. Ref. [49] demonstrated thatthe existing data from this sample is sensitive to sub-GeV DMif the DM can be produced in the beam dump and scatter elas-tically off detector electrons.The null result of the E137 search can be used to placebounds on our inelastically coupled scenario. The DM yieldis computed by evaluating the χ flux produced in the beamdump via dark-brehmstrahlung (described in Sec. III B) andconsidering both DM scattering off electron targets in the de-tector (described in Sec. III C 1), as well as χ decays insidethe detector volume (Sec. III C 2). To comply with the searchcriteria from Ref. [26], we demand either signal process de-posit at least 1 GeV of energy inside the detector’s geometricacceptance, and we place a 3-event bound to account for thenull result of the E137 experiment. 2. BDX The proposed BDX experiment at Jefferson Laboratory[28, 48, 50, 51] is a dedicated effort to search for lightDM produced and detected in analogy with the procedure atE137. The setup involves placing a meter-scale CsI scintil-lator detector behind the beam dump at the upgraded 11 GeV SBND [58] uses the same beamline as MiniBooNE, but with a smallerdetector placed closer to the beam stop. We expect the decay reach tobe similar to MiniBooNE for all but the largest mass splittings, but thescattering reach should be enhanced due to the higher detector efficiency[59]. However, MiniBooNE can distinguish photons from leptons [60]. A similar idea involving lower-energy electron beams installed near largeunderground neutrino detectors was proposed in Ref. [61]. CEBAF beam. In comparison with E137, BDX has greaterluminosity, ( ∼ electrons on target), a shorter baseline( ∼ m distance from dump to detector), and a larger detec-tor volume ( × × cm ).In the BDX setup, the boosted χ χ system emerges fromthe beam dump and can either scatter (see Sec. III C 1) via χ i e → χ j e or the excited state can decay via χ → χ e + e − Sec. III C 2) as it passes through the ∼ meter long detector.We compute decay and scattering projections using 3 and 10event yield contours, respectively, for EM energy depositionsabove 300 MeV. C. Signals at Electron Missing Momentum/EnergyExperiments 1. NA64 The NA64 experiment at the CERN SPS [24], depictedschematically in Fig. 3 (right panel), is sensitive to DM pro-duction via dark bremsstrahlung as described in Sec. III B.In this setup, × electrons with 100 GeV energies arefired into an active ECAL target which measures the en-ergy/momentum and triggers on events with large missing en-ergy. In principle, this technique is sensitive to our scenario ofinterest because the excited state can decay outside the detec-tor, thereby giving rise to missing energy in the ECAL mea-surement. While we include this discussion for completeness,we find that the existing data sample does not currently con-strain new parameter space, so we do not discuss it further. 2. LDMX The proposed LDMX experiment at SLAC [2, 23] aims toproduce DM using the 4 or 8 GeV LCLS-2 electron beampassing through a thin target upstream of a dedicated trackerand ECAL/HCAL system designed to veto SM particles pro-duced in these collisions – the setup is depicted schematicallyin Fig. 3 (left panel). For our inelastically coupled DM sce-nario, an A (cid:48) is produced via dark bremsstrahlung (describedin Sec. III B) and decays via A (cid:48) → χ χ followed by a dis-placed χ → χ e + e − decay. If this decay occurs behind theECAL/HCAL system, it can mimic the missing energy signa-ture which LDMX is optimized to observe.A representative realization of this setup involves ∼ electrons with 4 GeV energies passing through a tungsten tar-get of thickness ∼ . X and emerging from the target withless than 1 GeV of energy. The signal yield scales as N signal ≈ N e − n W σ ( eN → eN A (cid:48) ) X e − Γ χ (cid:96)/βγ , (26)where N e = 10 is the number of electrons on target, X ∼ cm is the tungsten radiation length, n W is the tungsten num-ber density, σ ( eN → eN A (cid:48) ) is the A (cid:48) production cross sec-tion from Eq. (15) and β/γ are the χ velocity and boost fac-tors, respectively. Here the exponential factor is the χ sur-vival probability through a path length of (cid:96) through the LDMX1ECAL/HCAL system. The LDMX missing momentum pro-jections present the 3 event yield contours for various param-eter benchmarks.We note in passing that it may also be possible for bothNA64 and LDMX to be sensitive to a promptly decaying χ → χ e + e − inside their ECAL systems, but the back-grounds for this process are not known and understanding thissignal is beyond the scope of the present work. V. PLOTS & MAIN RESULTS We now have all the ingredients in place to assess the poten-tial sensitivity of currently proposed fixed-target experimentsto discover inelastic DM models. We begin with a brief dis-cussion of existing constraints.The parameter space of inelastically coupled DM for massscales beneath ∼ GeV is constrained by precision measure-ments, B-factories, and previous fixed-target experiments. Onthe precision front, the anomalous magnetic moment of themuon and electron constrains the interaction strength betweenthe dark photon A (cid:48) and the SM particles [62]. On the colliderfront, both LEP and BaBar set a bound for larger values of (cid:15) .The former arises from the shift in m Z induced from mixingwith A (cid:48) [63], and the latter from a monophoton and missingmass re-analysis [40, 42, 48].Some of the strongest constraints from elastic DM arisefrom E137 [26, 49], an electron beam dump experiment, andLSND [25, 54], a proton beam fixed-target neutrino produc-tion experiment. Here, we reinterpret the constraints in termsof coannihilating DM. As discussed in Sec. III C above, thereare two qualitatively different signals: a scattering signal,where χ , up- or downscatter at the detector and produce arecoiling target and possibly an e + e − pair, and a decay signal,where χ survives to the detector and decays inside, produc-ing an e + e − pair. The reach of these experiments dependson their ability to distinguish these multiple signals. WhileE137 is only sensitive to total energy deposits, the angularresolution of LSND means it is potentially sensitive to well-separated e + /e − pairs, which can be distinguished from thefake elastic events we used in estimating the sensitivity in thiswork and could enhance the sensitivity. However, this wouldrequire access to the LSND data as this signal of two chargedtracks in the detector is not present in any published analy-sis. These existing constraints are illustrated in Fig. 6 for α D = 0.1, m A (cid:48) = 3 m , and various values of ∆ . For all but thesmallest splittings, the combination of LSND and E137 coversa large portion of the thermal target in the 1-100 MeV range.However, for m + ∆ > m π , DM production through pionsis kinematically forbidden, so we see sharp kinematic cutoffsat the pion threshold.For comparison, Fig. 7 shows the effect of varying ourbenchmark parameters (each panel varies one detail relative tothe top-left panel of Fig. 6) to demonstrate that these bench-marks are conservative and representative of the viable param-eter space. In particular, the top row of Fig. 7 shows how theparameter space in the top-left plot of Fig. 6 changes as weincrease and decrease α D while holding all other parameters fixed. Although there is slightly more viable parameter spacefor the large value α D = 0 . , this choice is close to the pertur-bativity limit for abelian dark sectors [65], so we regard ourbenchmark choice as a representative and conservative value;choosing a smaller coupling excludes more parameter spaceon the y vs. m plane, as we see for α D = α in the samefigure. In the bottom row of Fig. 7, we vary our ∆ and m A (cid:48) benchmarks from Fig. 6. In the bottom-left plot, we show thenearly-elastic case of ∆ = 0 . m , where the decay signalshuts off and the constraints are dominated by scattering. Forcomparison, we also show the recent MiniBooNE elastic scat-tering results [57], for which the beam energy is sufficientlylarge that the small 1% mass splitting does not affect the reach.In the bottom-right plot, we show results for a larger hi-erarchy, m A (cid:48) /m = 10 . For a given m , ∆ , and α D ,the production rate is decreased as that event now arisesfrom a much heavier A (cid:48) . If we parameterize the produc-tion rates at m A (cid:48) /m = 3 and m A (cid:48) /m = 10 as N (cid:15) and N (cid:15) , respectively, the total decay or scattering yield scalesas N , (cid:15) /m A (cid:48) . Thus, for a fixed event yield, (cid:15) scales lin-early with m A (cid:48) but only as N / , . Far from any kinematicboundaries, the sensitivity in y ∝ (cid:15) /m A (cid:48) improves relativeto the thermal target since the scaling with m A (cid:48) dominatesthe scaling with ( N /N ) / . However, the reach at largemasses degrades as the A (cid:48) mass approaches the maximumavailable energy more rapidly and A (cid:48) production shuts off.We now turn to the potential of new proposals to largelycover the entire parameter space motivated by thermal inelas-tic DM. We focus on three experiments representative of thesetups we have previously discussed: MiniBooNE, BDX, andLDMX, which are proton beam dump, electron beam dump,and missing energy experiments, respectively. For compari-son, we have also estimated the Belle II [41] sensitivity onlyin the monophoton and missing mass channel by rescaling theBaBar result by the appropriate luminosity factor. As dis-cussed in Sec. III, the dominant signal at MiniBooNE is χ de-cay in the detector whenever it is kinematically allowed. SinceMiniBooNE has particle ID [60, 67], electrons can in princi-ple be distinguished from photons, and thus a well-separated e + /e − pair and no other activity in the detector is a signalwith few irreducible backgrounds. This stands in sharp con-trast to the case of elastic DM scattering at MiniBooNE [57],which must always contend with an irreducible neutrino back-ground. Note that the lower boundary of the decay curve isset by the energy threshold and angular resolution accordingto Eq. (25). We did not attempt a detailed simulation of the ρ resonance, and thus the reach in y of the region very close to m A (cid:48) ∼ m ρ may differ slightly from what we show. The miss-ing energy signal at LDMX dominates at low masses, whilethe decay reach for BDX is similar to that of MiniBooNE athigh masses. Indeed, while the A (cid:48) production rate at protonbeam experiments is a few orders of magnitude larger than atelectron beam experiments, the higher luminosity and beamenergy at BDX compensate to give a roughly similar reach.This is advantageous given that different models for the me-diator could enhance or suppress production at proton beamscompared to electron beams, so the combination of both ex-periments probes a wide range of models.2 LDMXmissing mom.LSNDdecayLSNDscatter MiniBooNEdecayE137decay BDXdecayBDXscatterE137scatter ( g - ) μ BelleII RelicDensity N e ff , ( m o d e l d e p . ) BaBarmono γ → LEP - - - - - - - - - - - m [ MeV ] y = ϵ α D ( m / m A ' ) Thermal iDM, Δ = m , m A ' = m , α D = LSNDscatter LSNDdecay LDMXmissing mom. ( g - ) μ MiniBooNEdecayE137decayE137scatter BDXdecayLEPBaBarmono γ BDXscatter N e ff , ( m o d e l d e p . ) Belle IImono - γ → R e li c D e n s i t y - - - - - - - - - - - m [ MeV ] y = ϵ α D ( m / m A ' ) Thermal iDM, Δ = m , m A ' = m , α D = LDMXmissing mom.LSNDdecayLEP ( g - ) μ MiniBooNEdecay L S N D s c a t t e r E137decay BDXdecayBDXscatterBaBarmono γ E137scatter N e ff , ( m o d e l d e p . ) Belle IImono - γ → R e li c D e n s i t y - - - - - - - - - - - m [ MeV ] y = ϵ α D ( m / m A ' ) Thermal iDM, Δ = m , m A ' = m , α D = LDMXmissing mom.LSNDdecay ( g - ) μ E137scatter E137decay L S N D s c a t t e r MiniBooNEdecay BDXdecayBDXscatterBaBarmono γ LEPBelle IImono - γ R e li c D e n s i t y → N e ff , ( m o d e l d e p . ) - - - - - - - - - - - m [ MeV ] y = ϵ α D ( m / m A ' ) Thermal iDM, Δ = m , m A ' = m , α D = FIG. 6. Parameter space compatible with thermal inelastic DM for different choices of ∆ with constraints and future projections presented.The black relic density curve is computed using the procedure described in Appendix A. For each choice of ∆ , the relic density curve isinsensitive to the relative values of (cid:15), α D , or m /m A (cid:48) , however, some constraints depend sensitively on these choices. Typical examples ofthis sensitivity are LEP and ( g − µ for which the curves shown here are based only on their limits on (cid:15) ; the observables in question donot depend on α D or the DM/mediator mass ratio. Thus, where appropriate, we have adopted the conservative prescription α D = 0 . and m A (cid:48) /m = 3 to place these constraints on this plot, thereby revealing the remaining viable parameter space; see text for a discussion. Thecolored curves in these plots represent new results computed in work: solid lines are existing constraints, and dashed lines are projections.The Belle II [41] projections are estimated by rescaling the BaBar luminosity for the process e + e − → γA (cid:48) → γχ χ in which the χ decaysoutside the detector. The gray shaded regions represent kinetic mixing constraints ( g − µ [62]; LEP [63]; and BaBar [19]. Finally, the verticaldashed line labeled N eff . is a model-dependent bound from DM freeze-out reheating photons preferentially over neutrinos during BBN [64],excluding parameter space to the left of the line; if there are other sources of dark radiation, this bound can be alleviated. Fig. 8 combines the results in this work with the results ofRef. [19] to show the thermal target over a wide range of DMmasses. We see that except for a few isolated masses, thethermal target for coannihilating DM could be well-coveredby all three planned experiments below ∼ ∼ ∆ = 2 m e , while the decay signals dominate when kinemati-cally allowed. VI. CONCLUSION In this paper we have studied the fixed-target phenomenol-ogy of thermal dark matter with inelastic couplings to the SMand proposed a series of new searches for these interactions.These models are an instance of the general case where therelic abundance arises from thermal coannihilation betweenthe halo DM candidate χ and an unstable excited state, χ .Since the heavier state decays away in the early universe, thereis no annihilation at later times, and therefore no indirect de-tection. Furthermore, if the mass difference between these3 LDMXmissing mom . LSNDdecay ( g - ) μ E137decayBDXdecay R e li c D e n s i t y B D X s c a t t e r L S N D s c a t t e r E137scatterLEP BelleIIBaBarmono γ ( g - ) μ N e ff , ( m o d e l d e p . ) → MiniBooNEdecay - - - - - - - - - - - m [ MeV ] y = ϵ α D ( m / m A ' ) Thermal iDM, Δ = m , m A' = m , α D = LDMXmissing mom.LSNDdecay BDXdecayBDXscatterE137scatter L S N D s c a t t e r LEPBaBarmono γ Belle II ( g - ) μ R e li c D e n s i t y N e ff , ( m o d e l d e p . ) E137decay → MiniBooNEdecay - - - - - - - - - - - m [ MeV ] y = ϵ α D ( m / m A ' ) Thermal iDM, Δ = m , m A' = m , α D = α LDMXmissing mom.LSNDscatter ( g - ) μ MiniBooNEscatter R e li c D e n s i t y LEPBaBarmono γ BDXscatter Belle II ( g - ) μ E137scatter N e ff , ( m o d e l d e p . ) → - - - - - - - - - - - m [ MeV ] y = ϵ α D ( m / m A ' ) Thermal iDM, Δ = m , m A ' = m , α D = LDMXmissing mom.LSNDdecay L S N D s c a t t e r MiniBooNEdecayBDXscatterE137scatter ( g - ) μ R e li c D e n s i t y BaBar mono γ BelleIIBDXdecay N e ff , ( m o d e l d e p . ) E137decay → LEP - - - - - - - - - - - m [ MeV ] y = ϵ α D ( m / m A ' ) Thermal iDM, Δ = m , m A' = m , α D = FIG. 7. Top row: Same as the top-left panel of Fig. 6, but with different choices for α D . For larger couplings near the perturbativity limit [65](left with α D = 0 . ) the viable parameter space increases slightly relative to Fig. 6. For smaller couplings (right with α D = α ) the thermaltarget is nearly closed. Bottom row: Same as the-top left panel in Fig. 6, but with different choices for the inelastic splitting ∆ = 0 . m (left) and the mediator mass m A (cid:48) = 10 m (right). Note that for very small mass splittings, the decay searches become ineffective and the bestlimits arise from scattering and collider searches, whose observables do not rely on a prompt χ decay. states exceeds ∼ keV, upscattering at direct detection ex-periments is kinematically forbidden and loop-induced elasticscattering is small, so this scenario can likely only be discov-ered or falsified using accelerators. We leave the possibility ofone-loop elastic scattering at recently proposed electron directdetection experiments for a future study.At fixed-target experiments, the inelastic interaction re-sponsible for setting the relic abundance yields a variety ofobservable signatures arising from the boosted χ χ system,which is produced in a proton or electron beam collision withtarget nuclei. Once produced, either state can scatter off parti-cles in a downstream detector, thereby generating an observ-able signal. In addition, the boosted χ can also survive outto the detector and decay semi-visibly via χ → χ e + e − todirectly deposit a visible signal as it passes through the active volume.Using these signatures, we have extracted existing con-straints on this scenario by reinterpreting old LSND and E137data. To this end, we have generalized the analyses in [54] (forLSND) and [49] (for E137), which focused on the scatteringsignatures of elastically coupled DM. In our analysis, we havedemonstrated that there are several new signatures to whichthese older experiments are sensitive if DM couples inelasti-cally. In particular, we find that E137 and LSND already placenontrivial bounds on the parameter space that yields sub-GeVthermal coannihilation for a variety of DM masses, mass split-tings, and coupling strengths.We have also studied the prospects for future decay andscattering searches at the existing MiniBooNE (proton beam)experiment and the proposed BDX and LDMX (electron4 LDMXmissing mom.LSNDdecay L S N D s c a t t e r MiniBooNEdecayBDXdecay LHCdisplacedBelle IImono γ BDXscatterLEP B a B a r d i s p l a ce d E137scatter ( g - ) μ N e ff , ( m o d e l d e p . ) → E137decay R e l i c D e n s i t y - - - - - - - - - - m [ MeV ] y = ϵ α D ( m / m A ' ) Thermal iDM, Δ = m , m A ' = m , α D = LDMXmissing mom. LSND decayE137scatter E137decayLSNDscatter R e l i c D e n s i t y MiniBooNEdecay BDXdecay LHC l + l - + METBaBardisplaced LHCdisplacedBDXscatter ( g - ) μ Belle IImono - γ N e ff , ( m o d e l d e p . ) LSNDscatter → LEP - - - - - - - - - - m [ MeV ] y = ϵ α D ( m / m A ' ) Thermal iDM, Δ = m , m A ' = m , α D = FIG. 8. Same parameter space as the top-left and bottom-right panels of Fig. 6, but with the mass range extended out to the electroweak scale.Here we combine the results of this paper with the LHC, BaBar, and Belle II constraints and projections presented in Ref. [19]. The combinedreach from the sum of these efforts suffices to cover nearly all remaining parameter space for thermal coannihilation; thermal DM models withmasses below the MeV scale suffer generic conflicts with N eff [66] and masses above ∼ 100 TeV generically violate perturbative unitarity[5]. The only gaps not covered by this program of searches occur at very small mass splittings ∆ (cid:28) . m , depicted in the lower left panelof Fig. 7. For such small splittings, the decay searches become weak on account of the Γ χ ∝ ∆ scaling, and are not even kinematicallyallowed at low masses since ∆ < m e . mon beam) experiments. We find that the combined reach ofall scattering and decay searches at these experiments cancomprehensively test nearly all remaining parameter space,thereby providing strong motivation for these efforts.This paper also extends earlier work [19], which studiedthe collider phenomenology of inelastic thermal coannihila-tion models over the GeV – TeV mass range. Our work com-plements this effort by working out the constraints and projec-tions for the MeV–GeV range, thereby providing a roadmapfor covering thermal coannihilation over nearly all masses forwhich a thermal origin is viable (lower masses are in conflictwith early universe cosmology and higher masses genericallyviolate perturbative unitarity in most models). This full cover- age, spanning the results of both papers, is presented in Fig. 8.Finally we note that other existing and future experimentsmay also have powerful reach to this class of models. In par-ticular, the proton beam experiments DUNE [68], SeaQuest[69] (see forthcoming work by [70]), SHiP [71], and T2K[72] all involve beam energies in excess of 100 GeV, whichcan produce far more boosted DM than the beams consid-ered in this work (all < 10 GeV and below). Higher energiesat these experiments can open up new production modes forthe DM candidates (e.g. deep inelastic scattering) and impartgreater boosts to the DM system, which can profoundly affectthe sensitivity projections for these setups. In addition, liquidargon detectors such as ICARUS [73] may be more optimized5for seeing the two tracks characteristic of the decay signal. However, working out the implications of these features is be-yond the scope of this paper.Collider experiments may also probe the low-mass sce-nario we consider in this paper. The superior estimated reachof LHCb to visibly-decaying dark photons [74, 75] suggeststhat LHCb will also have sensitivity to inelastic DM, espe-cially when χ undergoes displaced decays within the detec-tor. However, this requires a dedicated analysis, as the pairsof leptons in this case do not reconstruct a resonance due tothe invisible χ . We leave this interesting analysis to futurework. In addition, BaBar could be additionally sensitive to thistopology through dedicated searches for monophoton and dis-placed leptons, as proposed by Ref. [19]. Similarly, throughanalogous searches, the larger luminosities at Belle II couldprovide unprecedented sensitivity to both displaced and long-lived decays. Parts of the parameter space with larger masssplittings can also lead to displaced vertices from χ decay,but requires a dedicated analysis for Belle II using displaced leptons; we defer this for a future study.We look forward to the results of the numerous plannedexperiments on the horizon, and encourage them to pursue theinelastic DM implications of the signals we have discussed inthis work. ACKNOWLEDGMENTS Acknowledgments: We thank Brian Batell, Asher Berlin,Nikita Blinov, Maxim Pospelov, Philip Schuster, Brian Shuve,Tim Tait, Natalia Toro, Richard Van De Water, and YimingZhong for helpful conversations. We thank Jordan Smolin-sky for pointing out typos in an earlier version of this paper.Fermilab is operated by Fermi Research Alliance, LLC, underContract No. DE-AC02-07CH11359 with the US Departmentof Energy. EI is supported by the United States Departmentof Energy under Grant Contract desc0012704. [1] G. Jungman, M. Kamionkowski, and K. Griest, Phys. Rept. , 195 (1996), arXiv:hep-ph/9506380 [hep-ph].[2] J. Alexander et al. (2016) arXiv:1608.08632 [hep-ph].[3] R. Essig, J. A. Jaros, W. Wester, P. H. Adrian, S. Andreas, et al. ,(2013), arXiv:1311.0029 [hep-ph].[4] V. Irˇsiˇc et al. , (2017), arXiv:1702.01764 [astro-ph.CO].[5] K. Griest and M. Kamionkowski, Phys. Rev. Lett. , 615(1990).[6] D. Tucker-Smith and N. Weiner, Phys. Rev. D64 , 043502(2001), arXiv:hep-ph/0101138 [hep-ph].[7] K. Griest and D. Seckel, Phys. Rev. D43 , 3191 (1991).[8] R. T. D’Agnolo, C. Mondino, J. T. Ruderman, and P.-J.Wang, https://indico.cern.ch/event/394659/contributions/944014/attachments/790207/1083119/Ruderman.2015.pdf , to appear.[9] P. A. R. Ade et al. (Planck), (2015), arXiv:1502.01589 [astro-ph.CO].[10] R. T. D’Agnolo and J. T. Ruderman, Phys. Rev. Lett. ,061301 (2015), arXiv:1505.07107 [hep-ph].[11] R. Essig, A. Manalaysay, J. Mardon, P. Sorensen, and T. Volan-sky, Phys. Rev. Lett. , 021301 (2012), arXiv:1206.2644[astro-ph.CO].[12] S. K. Lee, M. Lisanti, S. Mishra-Sharma, and B. R. Safdi, Phys.Rev. D92 , 083517 (2015), arXiv:1508.07361 [hep-ph].[13] R. Essig, M. Fernandez-Serra, J. Mardon, A. Soto, T. Volansky,and T.-T. Yu, JHEP , 046 (2016), arXiv:1509.01598 [hep-ph].[14] Y. Hochberg, Y. Kahn, M. Lisanti, C. G. Tully, and K. M.Zurek, (2016), arXiv:1606.08849 [hep-ph].[15] T. Emken, C. Kouvaris, and I. M. Shoemaker, (2017),arXiv:1702.07750 [hep-ph].[16] R. Essig, T. Volansky, and T.-T. Yu, (2017), arXiv:1703.00910[hep-ph]. We thank Maxim Pospelov for pointing this out. We thank the anonymous referee for suggesting this analysis. [17] Y. Bai and T. M. P. Tait, Phys. Lett. B710 , 335 (2012),arXiv:1109.4144 [hep-ph].[18] N. Weiner and I. Yavin, Phys. Rev. D86 , 075021 (2012),arXiv:1206.2910 [hep-ph].[19] E. Izaguirre, G. Krnjaic, and B. Shuve, Phys. Rev. D93 , 063523(2016), arXiv:1508.03050 [hep-ph].[20] P. Schuster, N. Toro, and I. Yavin, Phys. Rev. D81 , 016002(2010), arXiv:0910.1602 [hep-ph].[21] D. E. Morrissey and A. P. Spray, JHEP , 083 (2014),arXiv:1402.4817 [hep-ph].[22] D. Kim, J.-C. Park, and S. Shin, (2016), arXiv:1612.06867[hep-ph].[23] E. Izaguirre, G. Krnjaic, P. Schuster, and N. Toro, Phys. Rev. D91 , 094026 (2015), arXiv:1411.1404 [hep-ph].[24] D. Banerjee et al. (NA64), Phys. Rev. Lett. , 011802 (2017),arXiv:1610.02988 [hep-ex].[25] L. Auerbach et al. (LSND Collaboration), Phys.Rev. D63 ,112001 (2001), arXiv:hep-ex/0101039 [hep-ex].[26] J. Bjorken, S. Ecklund, W. Nelson, A. Abashian, C. Church, et al. , Phys.Rev. D38 , 3375 (1988).[27] I. Stancu et al. (MiniBooNE collaboration), (2001).[28] M. Battaglieri et al. (BDX), (2014), arXiv:1406.3028[physics.ins-det].[29] B. W. Lee and S. Weinberg, Phys.Rev.Lett. , 165 (1977).[30] C. Boehm and P. Fayet, Nucl.Phys. B683 , 219 (2004),arXiv:hep-ph/0305261 [hep-ph].[31] G. Krnjaic, (2015), arXiv:1512.04119 [hep-ph].[32] J. L. Feng, B. Fornal, I. Galon, S. Gardner, J. Smolinsky,T. M. P. Tait, and P. Tanedo, Phys. Rev. Lett. , 071803(2016), arXiv:1604.07411 [hep-ph].[33] S. Tulin, Phys. Rev. D89 , 114008 (2014), arXiv:1404.4370[hep-ph].[34] B. Holdom, Phys.Lett. B166 , 196 (1986).[35] H. Davoudiasl, H.-S. Lee, and W. J. Marciano, Phys. Rev. D85 ,115019 (2012), arXiv:1203.2947 [hep-ph].[36] H. Davoudiasl, H.-S. Lee, and W. J. Marciano, Phys. Rev. Lett. , 031802 (2012), arXiv:1205.2709 [hep-ph].[37] H. Davoudiasl, H.-S. Lee, and W. J. Marciano, Phys. Rev. D89 , , 123 (2014),arXiv:1310.8042 [hep-ph].[40] R. Essig, J. Mardon, M. Papucci, T. Volansky, and Y.-M.Zhong, JHEP , 167 (2013), arXiv:1309.5084 [hep-ph].[41] I. Jaegle (Belle), Phys. Rev. Lett. , 211801 (2015),arXiv:1502.00084 [hep-ex].[42] J. P. Lees et al. (The BaBar Collaboration), (2017),arXiv:1702.03327 [hep-ex].[43] S. Teis, W. Cassing, M. Effenberger, A. Hombach, U. Mosel,and G. Wolf, Z. Phys. A356 , 421 (1997), arXiv:nucl-th/9609009 [nucl-th].[44] P. deNiverville, C.-Y. Chen, M. Pospelov, and A. Ritz, Phys.Rev. D95 , 035006 (2017), arXiv:1609.01770 [hep-ph].[45] Y. Kahn, G. Krnjaic, J. Thaler, and M. Toups, Phys. Rev. D91 ,055006 (2015), arXiv:1411.1055 [hep-ph].[46] K. J. Kim and Y.-S. Tsai, Phys. Rev. D8 , 3109 (1973).[47] J. D. Bjorken, R. Essig, P. Schuster, and N. Toro, Phys.Rev. D80 , 075018 (2009), arXiv:0906.0580 [hep-ph].[48] E. Izaguirre, G. Krnjaic, P. Schuster, and N. Toro, Phys.Rev. D88 , 114015 (2013), arXiv:1307.6554 [hep-ph].[49] B. Batell, R. Essig, and Z. Surujon, (2014), arXiv:1406.2698[hep-ph].[50] E. Izaguirre, G. Krnjaic, P. Schuster, and N. Toro, (2014),arXiv:1403.6826 [hep-ph].[51] M. Battaglieri et al. (BDX), (2016), arXiv:1607.01390 [hep-ex].[52] B. Batell, M. Pospelov, and A. Ritz, Phys. Rev. D79 , 115019(2009), arXiv:0903.3396 [hep-ph].[53] C. Athanassopoulos et al. (LSND Collaboration),Nucl.Instrum.Meth. A388 , 149 (1997), arXiv:nucl-ex/9605002[nucl-ex].[54] P. deNiverville, M. Pospelov, and A. Ritz, Phys.Rev. D84 ,075020 (2011), arXiv:1107.4580 [hep-ph].[55] A. A. Aguilar-Arevalo et al. (MiniBooNE), Nucl. Instrum.Meth. A599 , 28 (2009), arXiv:0806.4201 [hep-ex].[56] P. deNiverville, D. McKeen, and A. Ritz, Phys.Rev. D86 ,035022 (2012), arXiv:1205.3499 [hep-ph].[57] A. A. Aguilar-Arevalo et al. (MiniBooNE), Submitted to: Phys.Rev. Lett. (2017), arXiv:1702.02688 [hep-ex].[58] M. Antonello et al. (LAr1-ND, ICARUS-WA104, Micro-BooNE), (2015), arXiv:1503.01520 [physics.ins-det].[59] R. Van De Water, private communication.[60] R. B. Patterson, E. M. Laird, Y. Liu, P. D. Meyers, I. Stancu,and H. A. Tanaka, Nucl. Instrum. Meth. A608 , 206 (2009),arXiv:0902.2222 [hep-ex].[61] E. Izaguirre, G. Krnjaic, and M. Pospelov, Phys. Rev. D92 ,095014 (2015), arXiv:1507.02681 [hep-ph].[62] M. Pospelov, Phys.Rev. D80 , 095002 (2009), arXiv:0811.1030[hep-ph].[63] A. Hook, E. Izaguirre, and J. G. Wacker, Adv. High EnergyPhys. , 859762 (2011), arXiv:1006.0973 [hep-ph].[64] K. M. Nollett and G. Steigman, Phys. Rev. D89 , 083508 (2014),arXiv:1312.5725 [astro-ph.CO].[65] H. Davoudiasl and W. J. Marciano, Phys. Rev. D92 , 035008(2015), arXiv:1502.07383 [hep-ph].[66] D. Green and S. Rajendran, (2017), arXiv:1701.08750 [hep-ph].[67] H.-J. Yang, B. P. Roe, and J. Zhu, Nuclear Instrumentsand Methods in Physics Research A , 370 (2005),physics/0508045.[68] R. Acciarri et al. (DUNE), (2016), arXiv:1601.05471 [physics.ins-det].[69] K. Nakahara (Fermilab E906/SeaQuest), Proceedings, 7th In-ternational Workshop on Neutrino-nucleus interactions in thefew GeV region (NUINT 11): Dehradun, India, March 7-11,2011 , AIP Conf. Proc. , 179 (2011).[70] A. Berlin, S. Gori, P. Schuster, and N. Toro, to appear.[71] S. Alekhin et al. , Rept. Prog. Phys. , 124201 (2016),arXiv:1504.04855 [hep-ph].[72] Y. Hayato (T2K), Neutrino physics and astrophysics. Proceed-ings, 21st International Conference, Neutrino 2004, Paris,France, June 14-19, 2004 , Nucl. Phys. Proc. Suppl. , 269(2005).[73] M. Antonello et al. , (2013), arXiv:1312.7252 [physics.ins-det].[74] P. Ilten, J. Thaler, M. Williams, and W. Xue, Phys. Rev. D92 ,115017 (2015), arXiv:1509.06765 [hep-ph].[75] P. Ilten, Y. Soreq, J. Thaler, M. Williams, and W. Xue, Phys.Rev. Lett. , 251803 (2016), arXiv:1603.08926 [hep-ph].[76] P. Gondolo and G. Gelmini, Nucl. Phys. B360 , 145 (1991).[77] G. Duda, A. Kemper, and P. Gondolo, JCAP , 012 (2007),arXiv:hep-ph/0608035 [hep-ph].[78] J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni,O. Mattelaer, H. S. Shao, T. Stelzer, P. Torrielli, and M. Zaro,JHEP , 079 (2014), arXiv:1405.0301 [hep-ph].[79] C.-Y. Chen, M. Pospelov, and Y.-M. Zhong, (2017),arXiv:1701.07437 [hep-ph]. Appendix A: Relic Abundance The relic abundance of χ is governed by a Boltzmannequation whose collision terms involve thermally averagedcoannihilation, decay, and inelastic scattering processes.Defining the dimensionless comoving yield as Y i ≡ n i /s ,where s ( T ) = 2 π g s, ∗ T / is the entropy density and g ∗ ,s is the number of entropic degrees of freedom, the Boltzmannsystem can be written as dY , dx = − λ A x (cid:16) Y Y − Y (0)1 Y (0)2 (cid:17) ± xλ D (cid:16) Y − Y (0)2 Y (0)1 Y (cid:17) ± λ S x Y (0) f (cid:16) Y − Y (0)2 Y (0)1 Y (cid:17) , (A1)where x ≡ m /T is a dimensionless time variable and Y (0) i is the comoving equilibrium yield for species i . We define λ A,S,D , to be the dimensionless annihilation, scattering, anddecay rates λ A = s ( m ) H ( m ) (cid:104) σv ( χ χ → f ¯ f ) (cid:105) (A2) λ S = s ( m ) H ( m ) (cid:104) σv ( χ f → χ f ) (cid:105) (A3) λ D = (cid:104) Γ( χ → χ f ¯ f ) (cid:105) H ( m ) , (A4)and we use the Hubble rate H ( T ) (cid:39) . √ g ∗ T /m P l duringradiation domination, where g ∗ is the number of relativisticdegrees of freedom, and m P l = 1 . × GeV is the Planckmass. Solving this system yields the asymptotic value Y ∞ at7freeze-out near x ∼ , which determines the relic abundance Ω DM = ρ χ ρ cr = m s Y ∞ ρ cr , (A5)where ρ cr = 8 . h × − GeV is the critical density and s (cid:39) cm − is the present day CMB entropy. An exam-ple solution to this system for a representative model point ispresented in Figure 9. Appendix B: Scattering and Annihilation Rates1. Coannihilation Rate The amplitude for this χ ( p ) χ ( p ) → f ( k ) ¯ f ( k ) coan-nihilation is A = i(cid:15)eg D s − m A (cid:48) ¯ u ( p ) γ µ v ( p ) ¯ u ( k ) γ µ v ( k ) , (B1)so squaring and averaging initial state spins gives (cid:104)|A| (cid:105) = 2 (cid:15) e g D ( s − m A (cid:48) ) (cid:20) (cid:0) t − m − m f (cid:1) (cid:0) t − m − m f (cid:1) + (cid:0) u − m − m f (cid:1) (cid:0) u − m − m f (cid:1) +2 m m (cid:0) s − m f (cid:1) + 2 m f (cid:0) s − m − m (cid:1) +8 m m m f (cid:21) . (B2)The differential cross section for this process is dσd Ω ∗ = | (cid:126)k ∗ || (cid:126)p ∗ | (cid:104)|A| (cid:105) π s , (B3)where | (cid:126)p ∗ | and | (cid:126)k ∗ | are the initial and final momenta in theCM frame and | (cid:126)k ∗ || (cid:126)p ∗ | = (cid:115) ( s − m f ) − m f ( s − m − m ) − m m . (B4)Integrating this result, the total cross section becomes σ ( s ) = | (cid:126)k ∗ || (cid:126)p ∗ | π(cid:15) αα D s ( s − m A (cid:48) ) (cid:20) s + 23 | (cid:126)p ∗ | | (cid:126)k ∗ | + 12 m χ ( s − m f ) + 12 m f ( s − m χ ) + m χ m f (cid:21) , (B5)where we have taken the elastic limit m = m = m χ forillustrative purposes, but we retain the full inelastic mass de-pendence in our calculations.Generalizing the derivation in [76] for coannihilation, thethermally averaged cross section is (cid:104) σv (cid:105) = 1 N (cid:90) ∞ s ds σ ( s )( s − s ) √ sK (cid:18) √ sT (cid:19) , (B6)where s ≡ ( m + m ) , the thermal averaging factor is N = 8 m m T K (cid:16) m T (cid:17) K (cid:16) m T (cid:17) , (B7) Y ( ) Y ( ) Y Y 10 20 30 40 5010 - - - - - - - - - - - - x = m / T Y = n / s Thermal Freeze Out via χ χ → f f FIG. 9. Example solution to the Boltzmann system in Eq. (A1)for m = 50 MeV, ∆ = 0 . m , m A (cid:48) = 3 m , and (cid:104) σv (cid:105) = 1 . × − cm s − , corresponding to Ω χ h = 0 . at late times. and K , are modified Bessel functions of the first and secondkinds.Note that in our calculation of the relic density, we ac-count for annihilation to hadrons (e.g. χ χ → A (cid:48)∗ π + π − )by following the procedure described in [19] where the finalstate phase space is extracted from the measured distribution R ( s ) ≡ σ ( e + e − → hadrons) /σ ( e + e − → µ + µ − ) . 2. DM-SM Scattering Cross Section For pseudo-Dirac DM, with mass splitting ∆ ≡ m − m , the matrix element for the process χ ( p ) T ( p ) → χ ( k ) T ( k ) is given by A = (cid:15)eg D ( t − m A (cid:48) ) [¯ u ( k ) γ µ u ( p )][¯ u ( k ) γ µ u ( p )] , (B8)so the squared, spin-average matrix element is (cid:104)|A| (cid:105) = 128 π (cid:15) αα D ( t − m A (cid:48) ) (cid:20) ( k · k )( p · p ) + ( k · p )( p · k ) − m m ( k · p ) − m T ( p · k ) + 2 m m m T (cid:21) . (B9)The differential scattering cross section in the CM frame is dσd Ω ∗ = 12 π dσd cos θ ∗ = (cid:104)|A| (cid:105) π s | (cid:126)k ∗ || (cid:126)p ∗ | , (B10)where initial/final state momenta satisfy | (cid:126)p ∗ | = ( s − m T − m ) − m T m s , (B11) | (cid:126)k ∗ | = ( s − m T − m ) − m T m s . (B12)8In terms of lab frame quantities with a stationary target T , s = ( p + p ) = m + m T + 2 m T E p , (B13) k · p = − 12 (2 m T − m − m − m T E k )= E ∗ p E ∗ k − | (cid:126)p ∗ || (cid:126)k ∗ | cos θ ∗ , (B14)so we can change variables to obtain the differential recoildistribution d cos θ ∗ = m T | (cid:126)p ∗ || (cid:126)k ∗ | dE T , (B15)where E T ≡ E k is the targets recoil energy. Thus, we have dσdE T = m T (cid:104)|A| (cid:105) πs | (cid:126)p ∗ | , (B16)which serves as an input into all detector scattering calcula-tions for both proton and electron beam dump experiments. Appendix C: Matrix Elements for DM Production and Detection1. Meson Decay The matrix element for pseudoscalar meson decay m ( p ) → γ ( k ) χ ( k ) χ ( k ) ( m = π , η ) is given by A m → γχ χ = − i(cid:15)e g D π f m (cid:16) g νλ − q ν q λ m A (cid:48) (cid:17) s − m A (cid:48) + im A (cid:48) Γ A (cid:48) × (cid:15) µναβ k α q β (cid:15) ∗ µ ( k ) (cid:2) ¯ u ( k ) γ λ v ( k ) (cid:3) , (C1)where f m is the meson decay constant, Γ A (cid:48) is the total A (cid:48) width, q ≡ p − k = k + k , and s ≡ q . The spin-averagedsquare of the matrix element is (cid:104) (cid:12)(cid:12) A m → γχ χ (cid:12)(cid:12) (cid:105) = 4 (cid:15) α α D πf m s − m A (cid:48) ) + m A (cid:48) Γ A (cid:48) (cid:2) ( m m − s ) ( s + 2 m m ) − s ( k · k )( k · k )+2( m m − s )( m − m )( k · ( k − k )) (cid:3) . (C2)If m + m (cid:28) m A (cid:48) (cid:28) m m then we can reasonably makethe narrow width approximation [45] and take s − m A (cid:48) ) + m A (cid:48) Γ A (cid:48) → πm A (cid:48) Γ A (cid:48) δ ( s − m A (cid:48) ) (C3)The decay width is given by d Γ m → γχ χ = 14 πm m β ( m m , , s )32 π β ( s, m , m )32 π ×(cid:104) (cid:12)(cid:12) A m → γχ χ (cid:12)(cid:12) (cid:105) d Ω ∗ γ d Ω ∗ χ ds, (C4)where we have defined the function β ( s , s , s ) = (cid:115) − s + s ) s + ( s − s ) s . (C5)Here, d Ω ∗ γ refers to angles in the m rest frame and d Ω ∗ χ refersto angles in the χ χ CM frame. 2. Excited State Decay The matrix element for χ ( p ) → χ ( k ) f ( k ) ¯ f ( k ) isgiven by A χ → χ f ¯ f = (cid:15)eg D ( s − m A (cid:48) ) + im A (cid:48) Γ A (cid:48) (cid:18) g µν − q µ q ν m A (cid:48) (cid:19) × [¯ u ( k ) γ µ u ( p )] [¯ u ( k ) γ ν v ( k )] , (C6)where again Γ A (cid:48) is the total A (cid:48) decay rate, q ≡ p − k , and s ≡ q . In this paper we will consider f = e − exclusively;decay to muons is allowed only for the largest masses andsplittings but may provide a distinctive signal at higher-energyexperiments. The spin-averaged square is (cid:104) (cid:12)(cid:12) A χ → χ f ¯ f (cid:12)(cid:12) (cid:105) = 16 (cid:15) e g D ( s − m A (cid:48) ) + m A (cid:48) Γ A (cid:48) × [( k · k )( k · p ) + ( k · p )( k · k )+ m f ( k · p ) − m m ( k · k ) − m m m f (cid:3) . (C7)We note that because we only consider m A (cid:48) > m and ∆ 3. DM-SM Scattering The tree-level matrix elements for scattering χ i ( p ) T ( p ) → χ j ( k ) T ( k ) off a pointlike fermionictarget have already been computed above in App. B 2. Sincewe will also be interested in scattering off targets with sub-structure such as nucleons and nuclei, we consider the moregeneral scattering process where T is a fermionic target withboth a monopole and dipole coupling to electromagnetism.The matrix element for this process is given by A χ i T → χ j T = (cid:15)eg D t − m A (cid:48) (cid:18) g µν − q µ q ν m A (cid:48) (cid:19) × [¯ u ( k ) γ µ u ( p )] [¯ u ( k )Γ ν u ( p )] , (C9)where q ≡ p − k is the four-momentum carried by the virtualphoton and t ≡ q is the Mandelstam variable. The Lorentzstructure Γ µ at the target vertex is Γ µ = F ( q ) γ µ + F ( q ) iq ν σ µν m f . (C10)Here, σ µν ≡ i [ γ µ , γ ν ] , and F ( q ) and F ( q ) are the elec-tric monopole and dipole form factors which depend on thetarget T . For the purposes of this paper, we take F = − q /m p ) T = p T = nZ T = N (C11)9and F = κ p (1 − q /m p ) T = p κ n (1 − q /m n ) T = n T = N (C12)with κ p ≈ . and κ n ≈ − . [54].The spin-averaged square of the matrix element is (cid:104) (cid:12)(cid:12) A χ i T → χ j T (cid:12)(cid:12) (cid:105) = 16 π (cid:15) αα D ( t − m A (cid:48) ) × (C13) (cid:110) ( F + F ) (cid:2) t ( p + k ) + 3 t ( t − ∆ ) − ( m − m ) (cid:3) + (cid:0) F + τ F (cid:1) (cid:104) (( p + k ) · ( p + k )) + ( k + p ) ( t − ∆ ) (cid:105)(cid:111) , where τ ≡ − t m T .For an incoming χ i with the target at rest in the lab frame,the lab-frame differential cross section is then given by dσ χ i T → χ j T d Ω ∗ = 12 E χ m T | (cid:126)v χ | β ( s, m j , m T )32 π ×(cid:104) (cid:12)(cid:12) A χ i T → χ j T (cid:12)(cid:12) (cid:105) , (C14)where β is defined as in Eq. (C5).We note that for the case of coherent scattering off of nuclei( T = N ), we make the additional insertion of the Helm formfactor [50, 77]. Appendix D: Monte Carlo Techniques Our simulation performs two distinct operations: produc-tion of χ and χ pairs and the detection of χ χ pairs ina detector. Many of the generic techniques used in our sim-ulations such as numerical phase space integration, rejectionsampling of differential probability distributions, and compu-tations utilizing detector geometries were borrowed from orinfluenced by BdNMC [44]. 1. DM Production Simulation for Proton and Electron Beams We use two simulation pipelines, one for proton beams andone for electron beams. Our proton beam production simula-tion takes as input an unweighted list of four-momenta fromone of the DM progenitors that we considered: π , η , and A (cid:48) from dark bremsstrahlung. The output of the production sim-ulation is an unweighted list of χ and χ particles produced.Each χ χ four-momentum pair in the output list is randomlysampled from the differential decay rate of the progenitor viaa rejection sampling method similar to that used in BdNMC .For electron beam production, we used a modified versionof Madgraph@NLO [78] from Ref. [79], which creates a newphysics model which contains new particle containers for thenucleus, A (cid:48) and χ , system. It utilizes elastic and inelasticatomic form factors from Ref. [46]. The elastic component isgiven by G ,el ( t ) = (cid:18) a t a t (cid:19) (cid:18) 11 + t/d (cid:19) Z , (D1) where a = 111 Z − / /m e , and d = 0 . GeV A − / . Wealso include a quasi-elastic (inelastic) term: G ,in ( t ) = (cid:18) a (cid:48) t a (cid:48) t (cid:19) (cid:32) t m p ( µ p − t . 71 GeV ) (cid:33) Z, (D2)where a (cid:48) = 773 Z − / /m e , m p is the proton mass, and µ p =2 . . The general form factor is then Φ ≡ (cid:90) t max t min dt t − t min t ( G ,el ( t ) + G ,in ( t )) . (D3) 2. DM Detection Simulation for Neutrino Detectors Our detector simulation takes as input the list of χ χ four-momentum pairs from the production simulation. We assumethat the various production processes all take place at the beamstop, so that the trajectories of each DM particle is well de-fined. There are three ways that a χ χ pair produced at thebeam stop can be detected: the χ can reach the detector andeither decay or downscatter or the χ can reach the detectorand upscatter.For each χ in the input list, a decay length x is randomlyselected from the appropriate exponential distribution. If the χ trajectory intersects with the surface of the detector, wecalculate the path length l along the trajectory from the beamstop to the detector. If x > l , i.e. if the χ persists until itreaches the detector, then the χ can be detected via directdecay or downscattering. If x < l or if the χ does not inter-sect the detector, we decay the χ → χ e + e − via rejectionsampling of the differential decay rate. We then process theresulting χ four-momentum in exactly the same way as a χ produced from primary progenitors, but now accounting forthe fact that this χ trajectory starts at the point where the χ decayed rather than the beam stop.Each χ that reaches the detector can either decay or down-scatter. We process decay events by weighting by the totalprobability for the χ to decay inside the detector and per-forming the decay χ → χ e + e − via rejection samplingof the differential decay rate. We then process the finalstate e + e − pair by applying the kinematic cuts described inSec. IV. We process scattering events by summing over targetsand consider scattering off of each target T independently. Toavoid double counting with decay events, we make the con-servative assumption that the χ can only scatter if it does notdecay in the detector. We therefore weight by the probabilitythat the χ persists through the detector and the independentprobability that the χ scatters in the detector and performthe scattering χ T → χ T by sampling from a uniform dis-tribution in the center of mass variables cos θ and φ , whichuniquely specify the final state kinematics. Because we sam-ple our final state kinematic variables from a uniform distri-bution rather than the true distribution of final states, P (cos θ, φ ) = 1 σ T dσ T d Ω , (D4)which is proportional to the differential cross section, we mustalso weight each event by the factor π P ( x , ..., x n ) . This0weighting scheme enables a cancellation of the total cross sec-tion σ T , making a Monte Carlo computation of σ T unneces-sary and thereby saving significant computational complexity.The cancellation occurs because each event is also weightedby the total probability for scattering which is ∝ σ T whenTaylor-expanded. Once the scattering final state is sampledand the weights computed, we apply the cuts described in Sec. IV to the recoiling target.For each χ that reaches the detector, we process the scat-tering χ T → χ T using the same method as χ downscat-tering. Additionally, we weight by the probability that the up-scattered χ will decay in the detector and perform the decayvia rejection sampling. The e + e −−