Testing Light Dark Matter Coannihilation With Fixed-Target Experiments

TL;DR

This paper proposes new fixed-target experimental methods to detect sub-GeV light dark matter that interacts via inelastic transitions, focusing on signals from scattering and decay, and analyzes existing and future experimental constraints.

Contribution

It introduces a novel fixed-target search strategy for inelastic light dark matter, including theoretical modeling, reanalysis of existing data, and projections for future experiments.

Findings

Current experiments already constrain significant parameter space.

Future experiments like Belle II and MiniBooNE can probe nearly all remaining viable regions.

Proposed methods can effectively detect inelastic dark matter via distinctive scattering and decay signals.

Abstract

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 heavier state, freeze-out proceeds via coannihilation and the unstable heavier state is depleted at later times. For sufficiently large mass splittings, direct detection is kinematically forbidden and indirect detection is impossible, so this scenario can only be tested with accelerators. Here we propose new searches at proton and electron beam fixed-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 model in which DM is a pseudo-Dirac fermion coupled to a hidden gauge field (dark photon), which kinetically mixes with the visible photon. We define theoretical targets in this framework and determine the existing bounds by reanalyzing results from previous experiments. We find that LSND, E137, and BaBar data already place strong constraints on the parameter space consistent with a thermal freeze-out origin, and that future searches at Belle II and MiniBooNE, as well as recently-proposed fixed-target experiments such as LDMX and BDX, can cover nearly all remaining gaps. We also briefly comment on the discovery potential for proposed beam dump and neutrino experiments which operate at much higher beam energies.

Equations

$$B_{\mu \nu },H^{†}H,LH,$$

$$U(1{)}_{B-L},U(1{)}_{{\ell }_i-{\ell }_j},U(1{)}_{3B-{\ell }_i},$$

$$\mathcal{L}=i\overline{\psi }/D\psi +M\bar{\psi }\psi +\lambda \varphi \overline{{\psi }^c}\psi +h.c.,$$

$$\mathcal{L}\supset g_D{A}_{\mu }^{\prime }{\overline{\chi }}_2{\gamma }^{\mu }{\chi }_1+h.c.,$$

$$\mathcal{L}\supset \frac{ϵ}{2\cos {\theta }_W}{F}_{\mu \nu }^{\prime }{B}^{\mu \nu }=ϵF_{\mu \nu }^{\prime }{F}^{\mu \nu }+ϵ\tan {\theta }_W{F}_{\mu \nu }^{\prime }{Z}^{\mu \nu },$$

$$e{A}_{\mu }{J}_{EM}^{\mu }\to e(A_{\mu }+ϵA_{\mu }^{\prime })J_{EM}^{\mu },$$

$$\sigma v\propto {ϵ}^2{\alpha }_D\frac{m_1^2}{m_{A^{\prime }}^4}=\left({ϵ}^2{\alpha }_D\frac{m_1^4}{m_{A^{\prime }}^4}\right)\frac{1}{m_1^2}\equiv \frac{y}{m_1^2},$$

$$y\equiv {ϵ}^2{\alpha }_D{\left(\frac{m_1}{m_{A^{\prime }}}\right)}^4,$$

$$y_{\exp }={ϵ}_{\exp }^2\times \left[{\alpha }_D{\left(\frac{m_1}{m_{A^{\prime }}}\right)}^4\right],$$

$$N_{{\chi }_1{\chi }_2}=N_{m^0}\times 2{ϵ}^2{\left(1-\frac{m_{A^{\prime }}^2}{m_{m^0}^2}\right)}^3\mathrm{BR}(m^0\to \gamma \gamma ),$$

$$\frac{d\sigma }{dxd{\Omega }_{A^{\prime }}}$$

$$U(x,{\theta }_{A^{\prime }})=E^{2}x{\theta }_{A^{\prime }}^2+m_{A^{\prime }}^2\frac{1-x}{x}+m_e^{2}x,$$

$$q_0=\frac{2U}{E(1-x)},q_f=m_{A^{\prime }}.$$

$$\frac{d\sigma }{dx}=4{\alpha }^3{ϵ}^2\overline{\Phi }(m_{A^{\prime }},E)\frac{x^2+3x(1-x)}{3m_{A^{\prime }}^2(1-x)},$$

$$\displaystyle\sigma(eN\to eNA^{\prime})\approx\frac{4\,\alpha^{3}\epsilon^{2}}{3\,m_{A^{\prime}}^{2}}\overline{\Phi}(m_{A^{\prime}},E)\bigl{[}\log\delta+{\cal O}(1)\bigr{]},~{}~{}~{}~{}~{}~{}~{}$$

$$\frac{d\sigma }{d{E}_e}=4\pi {ϵ}^2\alpha {\alpha }_D\frac{2m_e{E}_{\chi }^2-f(E_e)(E_e-m_e)}{(E_{\chi }^2-m_{\chi }^2)(m_{A^{\prime }}^2+2m_e{E}_e-2m_e^2{)}^2},$$

$$f(E_e)=2m_e{E}_{\chi }-m_e{E}_e+m_{\chi }^2+2m_e^2.$$

$${\sigma }_{\chi e}\approx \frac{8\pi \alpha {\alpha }_D{ϵ}^2{m}_e{E}_{\chi }}{m_{A^{\prime }}^4},$$

$$P_{\mathrm{scatter}}=n_e{\sigma }_{\chi e}\ell ,$$

$${\Gamma }_{{\chi }_2}\approx \frac{4{ϵ}^2\alpha {\alpha }_D{\Delta }^5}{15\pi m_{A^{\prime }}^4}$$

$$N_{\mathrm{signal}}\propto {ϵ}^2{P}_{\mathrm{survive}}{P}_{\mathrm{decay}},$$

$$P_{\mathrm{survive}}{P}_{\mathrm{decay}}=e^{-{\Gamma }_{{\chi }_2}d/\beta \gamma }\left(1-e^{-{\Gamma }_{{\chi }_2}\ell /\beta \gamma }\right)$$

$$R\equiv \frac{P_{\mathrm{decay}}}{P_{\mathrm{scatter}}}\simeq \frac{{\Delta }^5}{60{\pi }^2\gamma E_{\chi }{m}_e{n}_e},$$

$$R\simeq 3\left(\frac{m_1}{20}\right){\left(\frac{\Delta }{0.1m_1}\right)}^5{\left(\frac{10}{E_{\chi }}\right)}^2,$$

$${\Delta }^2\ge 2m_e^2\left(1+\cos {\theta }_{\min }\right)+2E_{\min }^2\left(1-\cos {\theta }_{\min }\right).$$

$$N_{\mathrm{signal}}\approx N_{e^{-}}{n}_W\sigma (eN\to eN{A}^{\prime })X_0{e}^{-{\Gamma }_{{\chi }_2}\ell /\beta \gamma },$$

$$\frac{d{Y}_{1,2}}{dx}$$

$${\lambda }_A$$

$${\lambda }_S$$

$${\lambda }_D$$

$${\Omega }_{\mathrm{DM}}=\frac{{\rho }_{{\chi }_1}}{{\rho }_{\mathrm{cr}}}=\frac{m_1{s}_0{Y}_1^{\infty }}{{\rho }_{\mathrm{cr}}},$$

$$\mathcal{A}=\frac{iϵe{g}_D}{s-m_{A^{\prime }}^2}\bar{u}(p_1){\gamma }^{\mu }v(p_2)\bar{u}(k_1){\gamma }_{\mu }v(k_2),$$

Full text

††preprint: FERMILAB-PUB-17-068-PPD, PUPT 2520

I Introduction

Understanding the particle nature of Dark Matter (DM) is among the highest priorities in all of physics. Perhaps the most popular DM candidate to date has been the Weakly Interacting Massive Particle (WIMP), which is charged under the electroweak force and naturally yields the observed cosmological abundance via thermal freeze-out (see Jungman et al. (1996) for a review). However, decades of null results from direct detection, indirect detection, and collider searches have cast doubt on this paradigm and motivated many alternative possibilities Alexander et al. (2016); Essig et al. (2013a).

Nonetheless, the thermal freeze-out mechanism remains compelling even if DM is not a WIMP. First and foremost, thermal DM is largely UV insensitive; its abundance is determined by the DM particle properties and is unaffected by the details of earlier, unknown cosmological epochs (e.g. reheating). Furthermore, unlike nonthermal production mechanisms, which can accommodate DM masses anywhere between $10^{-22}\,-10\,M_{\odot}$ (!), thermal DM is only viable between $\sim 5~{}-100~{}$, and is therefore more predictive. Dark matter masses outside this window are either too hot for acceptable structure formation Iršič et al. (2017) or violate perturbative unitarity Griest and Kamionkowski (1990). Finally, achieving the observed abundance requires a minimum interaction rate between DM and the Standard Model (SM), which provides a clear target for discovery or falsification. Thus, there is ample motivation to identify and study every viable realization of this mechanism.

One simple way to completely eliminate the tension between a thermal origin and experimental limits, in particular those from direct detection experiments, is for the DM to couple inelastically to SM particles Tucker-Smith and Weiner (2001). In this class of models, the halo DM species $\chi_{1}$ is the lightest stable particle in the dark sector and interacts with the SM only by transitioning to a slightly heavier state $\chi_{2}$. This class of models has several appealing features:

• •

Large Viable Couplings: If the inelastic interaction with the SM also determines the leading annihilation process, the relic abundance arises dominantly through $\chi_{1}\chi_{2}\to{\rm\small SM}$, a process dubbed coannihilation Griest and Seckel (1991). Since the heavier $\chi_{2}$ population is Boltzmann-suppressed during freeze-out, the requisite annihilation rate must compensate for this penalty.111For a general scenario where coannihilation proceeds without inelastic couplings, but an analogous enhanced thermal cross section appears, see Ref. D’Agnolo et al. . Thus, the coannihilation cross section always satisfies $\sigma v\geq 3\times 10^{-26}{\rm cm}^{3}{\rm s}^{-1}$.

• •

Indirect Detection Shuts Off: Since the heavier state is unstable, its population is fully depleted at low temperatures, so there are no remaining coannihilation partners for the $\chi_{1}$. This effect turns off possible late-time indirect detection signals and alleviates the bound from cosmic microwave background (CMB) power injection, which otherwise naively rules out thermal DM for masses below $\sim 10$ GeV for $s$-wave annihilation Ade et al. (2015).222Another way to evade this bound is “forbidden DM” D’Agnolo and Ruderman (2015).

• •

Direct Detection Forbidden: For a nonrelativistic halo particle scattering off a stationary SM target, the energy available to upscatter into the heavier state is $\sim\mu v^{2}$, where $\mu$ is the reduced mass of the DM-target system. Thus, for typical halo velocities $v\sim 10^{-3}$, direct detection off of nuclei is forbidden if the mass splitting exceeds $\gtrsim 100~{}$ Tucker-Smith and Weiner (2001). There is the possibility of loop-level induced elastic scattering off of electrons, which may be very relevant for new proposals for electron direct detection Essig et al. (2012); Lee et al. (2015); Essig et al. (2016); Hochberg et al. (2016); Emken et al. (2017); Essig et al. (2017). However, we leave this question for a future study.

Since direct and indirect detection search strategies are not available, testing thermal coannihilation fundamentally requires accelerator-based techniques. For DM masses near the weak scale (few GeV – 100s of GeV), Refs. Bai and Tait (2012); Weiner and Yavin (2012); Izaguirre et al. (2016) proposed LHC and B-factory searches with sensitivity to thermal coannihilation over a wide range of masses and splittings. However, for DM masses below the GeV scale, these searches become ineffective and new tools are required to fully test the keV–GeV mass range over which thermal coannihilation remains viable, yet unexplored.

In this paper, we fill this large gap by recasting and proposing a series of fixed-target searches for both electron and proton beam facilities. In both cases, an incident beam impinges on a fixed target and produces a boosted $\chi_{1}\chi_{2}$ pair. Depending on the experimental setup, this system can give rise to a variety of possible signals:

• •

Beam Dumps: The $\chi_{1}\chi_{2}$ pair can be produced radiatively through the “dark bremsstrahlung” reactions $pN\to pN\chi_{1}\chi_{2}$ at proton beam dumps, or $eN\to eN\chi_{1}\chi_{2}$ at electron beam dumps, where $N$ is a nuclear target. At proton beam dumps, the DM system can also be produced from meson decays $\pi^{0}\to\gamma\chi_{1}\chi_{2}$ and $\eta\to\gamma\chi_{1}\chi_{2}$, which can be advantageous for low DM masses. Once produced, the $\chi_{1,2}$ can reach a downstream detector and scatter off electrons or nucleons to induce an observable recoil signature. Alternatively, the unstable $\chi_{2}$ can also decay in flight as it passes through the detector via $\chi_{2}\to\chi_{1}e^{+}e^{-}$, thereby depositing an observable signal. These processes are depicted schematically in Fig. 2.

• •

Electron Missing Energy/Momentum: As above, the $\chi_{1}\chi_{2}$ pair is produced in $e^{-}N\to e^{-}N\chi_{1}\chi_{2}$ “dark bremsstrahlung”, but the production takes place in a thin target embedded in a detector. The principal observable signature in this context is the recoiling final-state electron. If this electron emerges having lost most of its incident beam energy, with no additional activity deposited in a downstream detector, this process can be sensitive to DM production. As above, a $\chi_{2}$ decaying inside the detector provides an additional potential signature. This process is depicted schematically in Fig. 3.

We will show that existing data already rules out large portions of the direct coannihilation parameter space. Moreover, the dedicated searches we propose which exploit the unique signals from this class of models can significantly improve on the current levels of sensitivity. Crucially, as shown in Fig. 4, the $\chi_{2}$ has a macroscopic decay length over nearly all of the parameter space we are interested in, and detecting the electromagnetic energy deposited by a $\chi_{2}$ decay in the detector — a signal not present in elastic DM models — is sufficient to cover large portions of the thermal relic curve.333A similar search strategy was proposed for the decay of long-lived scalars in Ref. Schuster et al. (2010). Indeed, we show in Sec. III.3 that the sensitivity to the decay of the exited state $\chi_{2}$ typically dominates over scattering channels at experiments with beam energies below 10 GeV. Thus, the prospects are excellent for dedicated experiments sensitive to these signatures, and we show that current and planned experiments can confirm or rule out nearly the entire mass range allowed for thermal coannihilating DM. We note a related study from Ref. Morrissey and Spray (2014) that investigated previous limits from fixed target facilities in the context of a supersymmetric hidden sector. Recently, Ref. Kim et al. (2016) proposed the signal of $\chi_{1}\to\chi_{2}$ upscattering followed by $\chi_{2}$ decay for the case of (boosted) astrophysical DM, as opposed to DM produced in a beam dump.

This paper is organized as follows. In Sec. II we introduce a representative renormalizable model featuring a pseudo-Dirac fermion DM coupled to a kinetically-mixed dark photon, where the abundance of the former arises from thermal coannihilation. In Sec. III we make some general comments on the production modes and detection signals at proton- and electron-beam fixed-target experiments. In Sec. IV we describe our methods for reinterpreting existing data from LSND Auerbach et al. (2001) and E137 Bjorken et al. (1988), and we discuss projections for the future data at MiniBooNE Stancu et al. (2001) and the proposed BDX Battaglieri et al. (2014) and LDMX Izaguirre et al. (2015a) experiments. In Sec. V we discuss the bounds and reach projections from these experiments in the context of the thermal coannihilating inelastic DM. Finally, in Sec. VI we offer some concluding remarks. Additional details on the matrix elements and Monte Carlo methods used in determining the thermal parameter space and in our simulations can be found in the Appendices.

II Sub-GeV Thermal Coannihilation

In this section, we describe a class of models of coannihilating DM: DM that couples inelastically to the SM through a kinetically-mixed dark photon. We detail the early universe cosmology and freeze out of the model, as well as introduce a useful parametrization of the parameters of the model in which the thermal target is largely an invariant under variation of couplings and of mass hierarchies.

II.1 Mediator Model Building

Unlike weak-scale WIMPs, which realize successful freeze-out with only SM gauge interactions, sub-GeV DM is overproduced in the absence of light ($\ll m_{Z}$) new mediators to generate a sufficiently large annihilation rate Lee and Weinberg (1977); Boehm and Fayet (2004). To avoid detection thus far, such mediators must be neutral under the SM and couple non-negligibly to visible particles.

If SM particles are neutral under the new interaction, a renormalizable model (without additional fields) requires the mediator to interact with the SM through the hypercharge, Higgs, or lepton portals

[TABLE]

for vector, scalar, and fermionic mediators, respectively. However, coupling a fermionic mediator to the lepton portal requires additional model building444A fermionic mediator coupled to the lepton portal requires additional model building to simultaneously achieve a thermal contact through this interaction and yield viable neutrino textures; the coupling to the mediator must be suppressed by neutrino masses, so it is generically difficult for the interaction rate to exceed Hubble expansion. and scalar mediators, which mix with the Higgs are ruled out for predictive models in which DM annihilates directly to SM final states (see Sec. II.3 and Krnjaic (2015) for a discussion of this issue), so we restrict our attention to abelian vector mediators; a nonabelian field strength is not gauge invariant, so kinetic mixing is forbidden.

Alternatively, the mediator could couple directly to SM particles if both dark and visible matter are charged under the same gauge group. In the absence of additional fields, anomaly cancellation restricts the possible choices to be

[TABLE]

and linear combinations thereof. In most contexts, the relevant phenomenology in fixed-target searches is qualitatively similar to the vector portal scenario, so below we will ignore these possibilities without loss of essential generality. We note, however, that viable models for both protophobic Feng et al. (2016) and protophilic Tulin (2014) mediators exist, so the complementarity provided by both proton- and electron-beam experiments is highly advantageous.

II.2 Representative Model

Our representative dark sector contains a 4-component fermion $\psi$ that transforms under a hidden abelian $U(1)_{D}$ gauge group. The fermion couples to a vector mediator $A^{\prime}$ as

[TABLE]

where $\phi$ is a $U(1)_{D}$ symmetry breaking scalar whose vacuum expectation value $v_{D}$ gives $A^{\prime}$ a nonzero mass $m_{A^{\prime}}\sim g_{D}v_{D}$ and gives $\psi$ a Majorana mass $\sim\lambda v_{D}$. Diagonalizing the resulting Dirac and Majorana masses gives rise to fermion mass eigenstates $\chi_{1,2}$ with a small mass splitting $\Delta\equiv m_{2}-m_{1}$ and an off-diagonal coupling to $A^{\prime}$,

[TABLE]

where $g_{D}\equiv\sqrt{4\pi\alpha_{D}}$ is the dark coupling constant. Note that it is technically natural to have $\Delta\ll M$ since the Majorana mass breaks the global symmetry associated with $\psi$ number.555If, unlike the construction in Eq. (3), the Majorana masses for the two Weyl components in $\psi=(\xi,\eta^{\dagger})$ are different, there is also a subleading diagonal interaction of the form $(\delta/M_{D})\overline{\chi}_{i}\displaystyle{\not}{A^{\prime}}\chi_{i}$, where $\delta\equiv m_{\xi}-m_{\eta}$ is the difference of Majorana masses for the the interaction eigenstates. We neglect this interaction in our analysis, assuming the off-diagonal interaction dominates.

This sector can interact with the SM through a renormalizable and gauge-invariant kinetic mixing term between $U(1)_{D}$ and $U(1)_{Y}$ gauge fields,

[TABLE]

where $\epsilon\ll 1$ is the kinetic mixing parameter and $F^{\prime}_{\mu\nu}$ and $B_{\mu\nu}$ are respectively the dark and hypercharge field strength tensors and the kinetic mixing interaction has been written in terms of the SM mass eigenstates $A$ and $Z$ after electroweak symmetry breaking. Diagonalizing the kinetic terms in Eq. (5) and rescaling the field strengths to restore canonical normalization induce a coupling between $A^{\prime}$ and the SM fermions Holdom (1986). To leading order in $\epsilon$, the $A^{\prime}$-SM interaction becomes

[TABLE]

where $J^{\mu}_{EM}$ is the SM electromagnetic current and all charged fermions acquire millicharges under $U(1)_{D}$. There is also an analogous $A^{\prime}$ interaction with the SM neutral current that arises from $A^{\prime}-Z$ mixing, but in our mass range of interest, $m_{A^{\prime}}\ll m_{Z}$, the mixing parameter is suppressed by an additional factor of $(m_{A^{\prime}}/m_{Z})^{2}$ Davoudiasl et al. (2012a, b, 2014); Kahn et al. (2016), so we neglect this interaction for the remainder of paper.

II.3 Direct Coannihilation vs. Secluded Annihilation

In the hot early universe ($T\gg m_{i},m_{A^{\prime}}$), all dark species are in chemical and kinetic equilibrium with the SM radiation bath; this initial condition is guaranteed as long as the DM-SM scattering rate exceeds the Hubble expansion rate at some point during cosmic history. If $m_{i}>m_{A^{\prime}}$, the freeze-out process is analogous to that of WIMP models. Below the freeze-out temperature $T_{f}\sim m_{1,2}/20$, the number densities of both species are depleted predominantly through $\chi_{i}\chi_{i}\to A^{\prime}A^{\prime}$ self-annihilation (which depends only on $\alpha_{D}$), not coannihilation, which depends on the combination $\epsilon^{2}\alpha_{D}$ and is greatly suppressed by comparison. Although both components undergo freeze-out separately, since $\chi_{2}$ is heavier and unstable, it will be depleted through downscattering and decays. Thus, up to order-one corrections, the requisite self-annihilation cross section satisfies the familiar WIMP-like requirement $\langle\sigma v\rangle\sim 3\times 10^{-26}{}^{3}\,{\rm s}^{-1}$ in order for $\chi_{1}$ to have the observed abundance at late times.

However, this secluded ($m_{i}>m_{A^{\prime}}$) regime has several drawbacks. Since the self-annihilation rate for fermions is $s$-wave, annihilation continues to occur out of equilibrium during recombination, which ionizes newly-formed hydrogen and thereby modifies the CMB power spectrum. For a thermal annihilation rate, this bound rules out DM below $\sim 10~{}$ Ade et al. (2015).666If instead, the DM is a scalar and annihilates directly to SM fermions through an $s$-channel vector mediator, its annihilation rate is $p$-wave suppressed, which can evade CMB bounds. Furthermore, since the secluded annihilation cross section scales as $\sigma v\sim\alpha_{D}^{2}/m_{i}^{2}$, the abundance is independent of the $A^{\prime}$ coupling to SM states, so there is no minimum interaction strength target for the DM search effort — at direct detection or accelerators, since these are sensitive to $\epsilon$.

By contrast, in the direct coannihilation regime $(m_{i}<m_{A^{\prime}})$, the two species annihilate each other via $\chi_{1}\chi_{2}\to f\overline{f}$, which has several compelling features:

• •

Predictive: The annihilation rate depends crucially on the mixing with the SM, $\sigma v\propto\epsilon^{2}\alpha_{D}$, so for dark couplings that satisfy perturbative unitarity, there is a minimum value of $\epsilon$ that is compatible with thermal freeze-out. Reaching experimental sensitivity to this minimum value for each viable DM mass suffices to discover or falsify this entire class of models.

• •

Large Cross Section: Unlike secluded annihilation, which involves the annihilation of equal mass particles, direct coannihilation requires both $\chi_{1}$ and $\chi_{2}$ in the initial state. However, $\chi_{2}$ typically becomes Boltzmann-suppressed before the nominal $\chi_{1}$ freeze-out temperature $\sim m_{1}/20$, so the coannihilation cross section must be larger to compensate for this reduction. Thus, we generically require $\langle\sigma v\rangle\gg 3\times 10^{-26}{}^{3}{\rm s}^{-1}$ to achieve the observed abundance, where the precise value grows sensitively with increasing $\Delta$ (compare the relic density curves in Fig. 6).

• •

Bounded Viable Range: Since the necessary coannihilation rate for freeze-out grows sharply as the mass splitting is increased, for sufficiently large $\Delta$, i.e., $\Delta\sim m_{1}$, the requisite $\epsilon$ is easily excluded independently of other model properties (e.g. by precision QED/electroweak measurements). Thus, as we will see below, this class of coannihilating models can be tested over the full viable parameter space.

Thus, for the remainder of this paper, our focus will be on the direct coannihilation scenario. For simplicity, without loss of essential generality, we will also demand that $m_{A^{\prime}}>m_{1}+m_{2}$ so that $A^{\prime}$ can decay to dark sector states — otherwise, the $A^{\prime}$ decays to SM final state, a signature for which there are abundant ideas to cover Alexander et al. (2016). Since $\epsilon\ll 1$, the $A^{\prime}$ branching ratio to SM states is correspondingly negligible in this regime.

Note that if we had chosen a neutral scalar mediator instead of the vector $A^{\prime}$, it could only couple to SM fermions by mixing with the Higgs after electroweak symmetry breaking. However, such a scenario cannot viably realize thermal coannihilation below the GeV scale because the requisite Higgs-mediator mixing angle must be $\sim{\cal O}(1)$ to overcome the Yukawa penalty in the annihilation diagram and yield a thermal annihilation rate. Such large couplings are ruled out by Higgs coupling measurements and rare meson decay searches Clarke et al. (2014); Krnjaic (2015). We note that $t$-channel annihilation into a light scalar mediator can still realize secluded annihilation, but this process is less predictive and beyond the scope of this work.

II.4 Covering the Thermal Target

The goal of this paper is to compute existing bounds on the parameter space for thermal freeze-out via coannihilation and present sensitivity projections for future experiments. One major challenge of such an effort is the large dimensionality of this parameter space: each viable model point is defined by a unique choice of the inputs $\{m_{1},m_{A^{\prime}},\Delta,\epsilon,\alpha_{D}\}$, constrained by the requirement that the coannihilation rate yields the observed DM density. Fortunately, in the nonrelativistic limit, for each choice of $\Delta$, the coannihilation rate scales as

[TABLE]

which is valid for $s\simeq(m_{1}+m_{2})^{2}$ and sufficiently away from the $s$-channel resonance at $m_{A^{\prime}}\sim m_{1}+m_{2}$. Here we have defined the dimensionless parameter

[TABLE]

which uniquely determines the freeze-out annihilation rate for a given choice of $m_{1}$ and $\Delta$.

The virtue of this parameterization is that the observed relic density is achieved only along a one-parameter curve $y(m_{1})$ which is insensitive to the relative sizes of $\epsilon,\alpha_{D}$, or $m_{1}/m_{A^{\prime}}$, reducing the dimensionality of the parameter space. However, the drawback is that some experiments only constrain a subset of model parameters and are not sensitive to the same combination of inputs that define the $y$ variable in Eq. (8). For example, lepton colliders can be used to constrain $\epsilon$ as a function of $m_{A^{\prime}}$ in searches for $e^{+}e^{-}\to\gamma{\displaystyle}{\not}{E}$, which can be interpreted as $e^{+}e^{-}\to\gamma A^{\prime}$ production followed by prompt $A^{\prime}$ decay to invisible particles Essig et al. (2013b); Jaegle (2015); Lees et al. (2017). Since ${\rm Br}(A^{\prime}\to{\rm invisible})\simeq 1$, the signal event yield depends on $\epsilon$, but is insensitive to $\alpha_{D}$ or the ratio $m_{1}/m_{A^{\prime}}$ for larger $\chi_{2}$ lifetimes.

Nonetheless, it is still possible to use this information to constrain the $y$ variable for a conservative comparison with the thermal target. The strategy is to construct the largest possible $y_{\rm exp}$ using the experimentally determined $\epsilon_{\rm exp}$,

[TABLE]

where the quantity in square brackets is chosen to be as large as possible while remaining consistent with both perturbative unitarity and the requirement that direct coannihilation remains the dominant annihilation process. Thus, the conservative prescription is to adopt order-one values for $\alpha_{D}$ and $m_{1}/m_{A^{\prime}}$, while taking $m_{A^{\prime}}>m_{1}+m_{2}$ and avoiding the $m_{A^{\prime}}\sim m_{1}+m_{2}$ resonance. As an illustration, Fig. 5 shows the reach in $\epsilon$ for $\chi_{2}$ decay at various experiments, demonstrating that the weakest reach with respect to the thermal target occurs at large $\alpha_{D}$. For our numerical results in Figs. 6 and 8, we chose representative benchmarks $\alpha_{D}=0.1,m_{1}/m_{A^{\prime}}=1/3$ where appropriate (e.g. for collider bounds). In Fig. 7 and in the discussion in Sec. V below, we show how the constrains on the parameter space change for different choices of these benchmark values, and demonstrate that our conclusions are qualitatively unchanged.

Note that for our main results, which are discussed in Sec. V and presented in Figs. 6, 7, and 8, we compute the $y(m_{1})$ curve by numerically solving the full Boltzmann system for this model as described in Appendix A.

III Production and Detection Basics

In this section we review the basic $A^{\prime}$ production formalism in proton- and electron-beam fixed-target collisions to develop some intuition for our later numerical results. We also present the formalism for various DM detection signatures that ensue from boosted $A^{\prime}\to\chi_{1}\chi_{2}$ decays.

III.1 Production at Proton Beam Dumps

At proton beam dumps, $\chi_{1}\chi_{2}$ pairs can be produced either from dark bremsstrahlung $pN\to pNA^{\prime}$, $A^{\prime}\to\chi_{1}\chi_{2}$ or from the neutral meson decays $m^{0}\to\gamma A^{\prime(*)}\to\gamma\chi_{1}\chi_{2}$ where $m^{0}=\pi^{0},\eta$.777At sufficiently high energies, direct $A^{\prime}$ production from quarks and gluons in the target becomes relevant, but this process is negligible for the beam energies we consider in this work. The meson production rates and spectra depend on various experimental factors, but the strongest dependence is on the proton beam energy, while the detailed material properties of the beam dump can largely be ignored. At few-GeV beam energies, $N_{\pi^{0}}\approx N_{p}$ and $N_{\eta}\approx N_{p}/30$ Teis et al. (1997) where $N_{p}$ is the total flux of protons. See Ref. deNiverville et al. (2017) for a detailed discussion of production modes and distributions.

If $m_{A^{\prime}}<m_{m^{0}}$, the $A^{\prime}$ can be produced on-shell and the number of $\chi_{1}\chi_{2}$ pairs via meson decays is given by

[TABLE]

where $N_{m^{0}}$ is the number of $\pi^{0}$ or $\eta$ produced.888If $m_{A^{\prime}}>m_{m^{0}}$, the meson decay proceeds through an off-shell $A^{\prime}$ and the rate scales as $\epsilon^{2}\alpha_{D}$ Kahn et al. (2015), in contrast to on-shell production which is independent of $\alpha_{D}$. At low-energy experiments such as LSND with only a single production channel, off-shell production may be important, but for higher-energy experiments such as MiniBooNE with multiple production channels, on-shell production is dominant. In either case, there is a kinematic threshold for production through mesons: $m_{1}+m_{2}=2m_{1}+\Delta<m_{\pi^{0}}\approx 135\ \mathrm{MeV}$ for pions, and $2m_{1}+\Delta<m_{\eta}\approx 548\ \mathrm{MeV}$ for etas.

The rate of dark bremsstrahlung production is proportional to $\epsilon^{2}$, but varies non-trivially with the dark photon mass $m_{A^{\prime}}$. Using the simulation of Ref. deNiverville et al. (2017), we find that the number of $A^{\prime}$ produced varies from $\sim 10^{-4}\epsilon^{2}N_{p}$ as $m_{A^{\prime}}\to 0$ to $\sim 10^{-9}\epsilon^{2}N_{p}$ at $m_{A^{\prime}}=3\ \mathrm{GeV}$. $A^{\prime}$ production can be enhanced through vector meson mixing near $m_{A^{\prime}}=m_{\rho},m_{\omega}$, but since this enhancement is only in a limited mass range and depends on the precise choice of $m_{1}/m_{A^{\prime}}$, we do not attempt to model it precisely. The dark bremsstrahlung production mechanism has no mass threshold and so heavy DM can still be produced up to the beam energy, $2m_{1}+\Delta<E_{\rm beam}$.

III.2 Production at Electron Beams

At electron-beam experiments, mesons are no longer copiously produced so the dominant process is $A^{\prime}$ production through dark bremsstrahlung followed by on-shell mediator decay $A^{\prime}\to\chi_{1}\chi_{2}$. The reaction $eN\rightarrow eNA^{\prime}$, where $N$ is a nucleus of atomic (mass) number $Z(A)$ has been well studied in Refs. Kim and Tsai (1973); Bjorken et al. (2009); Izaguirre et al. (2013). The relevant features of the reaction can be better illuminated by considering the calculation using the Weizsacker-Williams (WW) approximation, although we note that all plotted results in this study employ a numerical calculation; see Appendix D for more details of our simulations.

In the WW approximation, the differential production cross section of $A^{\prime}$ is then given by

[TABLE]

where $d\Omega_{A^{\prime}}=2\pi d\cos\theta_{A^{\prime}}$ is the $A^{\prime}$ solid angle with respect to the beam axis in the lab frame and we have defined

[TABLE]

Here $E$ is the electron beam energy, $x=E_{A^{\prime}}/E$ the fraction of the beam energy carried by the $A^{\prime}$, and $\alpha\Phi(q_{0},q_{f})/\pi$ is the WW photon flux for small-virtuality photons with momentum $-t$ bounded by $q_{0}^{2}$ and $q_{f}^{2}$ Kim and Tsai (1973). For $q_{0}$ values much smaller than the inverse nuclear size $\approx 0.4\,/A^{1/3}$, we have $\Phi(q_{0},q_{f})\sim Z^{2}$ up to an overall logarithmic factor.

Ref. Bjorken et al. (2009) found that for any given $x$, the angular integral is dominated by angles $\theta_{A^{\prime}}$ such that $Ex\theta_{A^{\prime}}^{2}\lesssim m_{A^{\prime}}^{2}\frac{1-x}{x}+m_{e}^{2}x$. Using this approximation, we can derive a simpler expression for the differential cross section

[TABLE]

where $\overline{\Phi}(m_{A^{\prime}},E)\equiv\Phi(q_{0}=m_{A^{\prime}}^{2}/(2E),q_{f}=m_{A^{\prime}})$. After integrating, the total $A^{\prime}$ production cross-section scales roughly as

[TABLE]

where $\delta\equiv\min(m_{A^{\prime}}^{2}/E^{2},m_{e}^{2}/m_{A^{\prime}}^{2},m_{e}^{2}/E^{2})$.

III.3 Generic Detection Signals (Electron & Proton Beams)

In our regime of interest, namely $\epsilon^{2}\alpha<\alpha_{D}$, $A^{\prime}$ decays promptly primarily into the $\chi_{1}\chi_{2}$ system, imparting an order-one fraction of its energy to the daughter particles, which emerge from the target with large boosts. The boosted DM system gives rise to three classes of observable signatures: detector target scattering, $\chi_{2}$ decays, and missing energy/momentum carried away by the DM system.

III.3.1 Detector Target Scattering

The $\chi_{1}$ produced in the dump travels unimpeded, enters the detector situated downstream of the target/dump, and scatters off of target particles in the detector. It scatters via the reaction $\chi_{1}T\rightarrow\chi_{2}T$, where $T$ here could in principle be an electron, a nucleon, or even a nucleus. Similarly, any population of surviving $\chi_{2}$s produced in the target could downscatter through the reverse reaction $\chi_{2}T\rightarrow\chi_{1}T$. The cleanest signal occurs when $T$ is an electron, with energy typically above 10s of MeV. The production rate is proportional to $\epsilon^{2}$ and the scattering rate is proportional to $\epsilon^{2}\alpha_{D}$, so the total yield in this channel will be proportional to $\epsilon^{4}\alpha_{D}$.

Specializing to the case of electron targets in the $E_{\chi}\gg\Delta$ regime, the approximate differential cross section is Izaguirre et al. (2013); Batell et al. (2014)

[TABLE]

where $E_{e}$ is the electron recoil energy, $E_{\chi}$ is the incident $\chi_{1,2}$ energy, and we define

[TABLE]

This cross section is valid up to corrections of order $\Delta/E_{\chi}$ and $\Delta/m_{A^{\prime}}$, both of which are very small compared for the benchmarks we consider throughout (see Figs. 6, 8, and 7). However, our numerical results evaluate the exact expression presented following the derivation in App. B.2.

Neglecting subleading corrections and integrating the recoil energy up to $E_{\chi}$ gives the approximate result

[TABLE]

so the corresponding scattering probability inside the detector becomes

[TABLE]

where $\ell$ is the DM path length through the detector and $n_{e}$ is the electron number density. There are similar expressions for scattering off detector nucleons and nuclei, but, as we will see, most of the relevant bounds and projections below exploit the electron channel.

III.3.2 Decay of excited state

One of the most powerful channels at beam dump experiments is the direct decay of $\chi_{2}\to\chi_{1}e^{+}e^{-}$, whose partial width satisfies

[TABLE]

in the $m_{A^{\prime}}\gg m_{1}\gg m_{e}$ limit.999For $\Delta>2m_{\mu}$, $\chi_{2}$ can also decay to muons, but for the majority of the parameter space we consider in this paper, only the electron channel is allowed. The $\chi_{2}$ is produced in the beam dump and decays in a downstream detector (depicted schematically in Fig. 2), so the signal yield scales as

[TABLE]

where the product

[TABLE]

is the probability for $\chi_{2}$ to survive a distance $d$ out to do the detector and decay inside after traversing a path length $\ell$ with boost factor $\gamma$ and velocity $\beta=v/c$. As with scattering, the total decay yield scales as $\epsilon^{4}\alpha_{D}$, with $\epsilon^{2}$ coming from production and $\epsilon^{2}\alpha_{D}$ coming from expanding the exponentials in Eq. (22) assuming the long lifetime limit (see Fig. 4) and using Eq. (20).

To estimate the relative reach of decay and scattering searches at beam dump experiments, it is useful to define

[TABLE]

where $n_{e}$ is the number density of detector electrons. Here we have used the approximate $\chi$-electron scattering cross section from Eq. (18); note that this expression is independent of detector geometry or DM production rate.

For most materials, $n_{e}\sim 10^{24}\ {}^{-3}$, so we have

[TABLE]

where $E_{\chi}$ represents the typical energy of a $\chi$ within the detector acceptance. Note that the $R\sim\Delta^{5}$ scaling implies that decays dominate the reach for most of the splittings shown in Fig. 6, and that for experiments with lower DM energies (e.g. at LSND where $E_{\chi}\sim$ few 100 MeV) the decay reach dominates the scattering reach by several orders of magnitude due to the lower beam energy and smaller boost factors.

III.3.3 Upscatter followed by decay

A third possible signal combines the phenomenology of both scattering and decay. $\chi_{1}$ can upscatter at the detector, producing a $\chi_{2}$ along with a recoiling target $T$. If the $\chi_{2}$ produced in the upscatter decays inside the detector via the reaction $\chi_{2}\rightarrow\chi_{1}e^{+}e^{-}$, and if the electron and positron final states are energetic and separated enough, the final state leads to a signature that is not easily mimicked by environmental or beam-related backgrounds Izaguirre et al. (2014); Battaglieri et al. (2014, 2016). However, we find that the yields in this channel are subdominant to the other two and we will not consider it further.101010If the recoiling target $T$ is not visible, the signal in this channel is identical to the decay signal, but with the added $\epsilon^{2}\alpha_{D}$ penalty of scattering. The consequences of this signal for the case of boosted astrophysical DM, where the different kinematics and the absence of an $\epsilon^{2}$ production penalty make this channel an attractive background-free option, were investigated in detail in Ref. Kim et al. (2016).

III.3.4 Missing Energy/Momentum

A unique feature of $A^{\prime}$ production in electron-beam fixed-target collisions is that the $A^{\prime}$ is typically radiated with an order-one fraction of the incident beam energy (see Sec. III.2 and more detailed discussions in Kim and Tsai (1973); Bjorken et al. (2009); Izaguirre et al. (2015a); Banerjee et al. (2017)). This enables a novel detection strategy where the target is now embedded in the detector, and which is based on comparing the energy and momentum of the beam electron before and after it undergoes dark bremsstrahlung. Fig. 3 illustrates a set-up for this detection strategy.

III.4 Kinematic Thresholds

There are several kinematic thresholds which influence the detectability of the various signals. Clearly, for $\Delta<2m_{e}\approx 1\ \mathrm{MeV}$, the excited state $\chi_{2}$ can no longer decay via $\chi_{2}\to\chi_{1}e^{+}e^{-}$. The only possible decay mode remaining is to $\chi_{1}$ plus neutrinos or $3\gamma$, both of which are sufficiently suppressed that $\chi_{2}$ is cosmologically long-lived Batell et al. (2009). In this case the only signals are from the recoiling target in $\chi_{1}$ upscattering or $\chi_{2}$ downscattering, along with missing energy/momentum. Otherwise, for detectors with finite energy thresholds and angular resolution, we might require the decay signal $\chi_{2}\to\chi_{1}e^{+}e^{-}$ to have an energetic, well-separated $e^{+}e^{-}$ pair satisfying $E_{e}>E_{\rm min}$ and $\theta_{e^{+}e^{-}}>\theta_{\rm min}$. The requirement on $\Delta$ to allow a visible $\chi_{2}$ decay is actually more stringent:

[TABLE]

For example, requiring E_{\rm min}=50\ and $\theta_{\rm min}=2^{\circ}$ as at MiniBooNE, we find that the minimum splitting we can probe is $\Delta\gtrsim 2\ \mathrm{MeV}$.

IV Existing Bounds and Projections

In this section we discuss the key features of the various representative experiments we consider and list our kinematic cuts used to compute reach curves. The bounds and projections we derive are presented in Figs. 6, 7, and 8.

IV.1 Signals at Proton Beam Dump Experiments

IV.1.1 LSND

The Liquid Scintillator Neutrino Detector (LSND) Athanassopoulos et al. (1997) was a neutrino experiment at Los Alamos which ran from 1993-1998. The extremely high-luminosity proton beam produced the largest available fixed-target sample of neutral pions, $N_{\pi^{0}}\sim 10^{22}$. The LSND proton beam had kinetic energy 800 MeV, small enough that $\eta$ production and dark bremsstrahlung are negligible. Therefore, the DM is produced dominantly from $\pi^{0}$ decays. We modeled the $\pi^{0}$ production at LSND using the GEANT simulation from Ref. Kahn et al. (2015) and then decayed the pions to DM using Monte Carlo. Due primarily to the large luminosity, we expect DM signal event yields at LSND to dominate in the regions of parameters space where pion decay is kinematically allowed.

The LSND detector was an approximately cylindrical tank of scintillating mineral oil with length $\sim 10\mathrm{m}$ and radius $\sim 3\mathrm{m}$, located $\sim 30\mathrm{m}$ from the beam stop. The close proximity to the beam stop and its off-axis orientation gives a large geometric acceptance of tracks originating at the beam stop. The detector used photomultiplier tubes to detect Cerenkov light emitted by charged particles in the scintillator. For the purposes of our simulations, LSND is only capable of identifying lepton tracks but is blind to nucleons, and electrons are indistinguishable from positrons. We take the electron detection efficiency at LSND to be $19\%$ Auerbach et al. (2001).

We numerically computed the event yields at LSND using DM events that were produced as outlined in Sec. III.1, for each of the available signal channels in Sec. III.3. Using the techniques outlined in App. D.2, we simulated the $\chi_{1}$ and $\chi_{2}$ scattering off electrons, nucleons, and nuclei as well as direct $\chi_{2}$ decays in the detector. This simulation accounted for the geometric acceptance of the detector as well as kinematic cuts and detection efficiencies for the observable final state particles.

Previous dark matter limits at LSND have been derived using the 55-event background-subtracted limit from the LSND neutrino-electron elastic scattering analysis Auerbach et al. (2001); deNiverville et al. (2011). To derive the LSND decay constraints for inelastic DM, we interpreted the 55-event limit using $\chi_{2}$ decay events where the $e^{+}/e^{-}$ opening angle was too small to distinguish the two leptons, and thus could fake an elastic event passing the cuts described in Ref. Auerbach et al. (2001). The angular resolution of LSND is $12^{\circ}$ for electrons above 20 MeV Athanassopoulos et al. (1997), so we require that the angular separation of the $e^{+}e^{-}$ pair is $\theta_{e^{+}e^{-}}<12^{\circ}$ and that the total energy of the pair lies in the range $18<E_{e^{+}}+E_{e^{-}}<50\ \mathrm{MeV}$ with either the electron or positron satisfying $\cos\theta<0.9$. This is conservative, as we have ignored a similar number of events where the electron and positron are well-separated; with access to the full LSND data set, one could search for events where two leptons were seen simultaneously inside the detector and potentially improve the reach. Indeed, we chose such a conservative prescription because, as we will see, the limits derived from LSND are already dominant for $\Delta>2m_{e}$, $2m_{1}+\Delta<m_{\pi^{0}}$.

For the LSND scattering constraints, we again use the 55-event limit applied to any event with at least one lepton passing the elastic cuts in the final state. This includes $\chi_{2}$ down-scattering and $\chi_{1}$ upscattering with and without $\chi_{2}$ decay, though the dominant process for all but the largest DM masses is $\chi_{2}$ downscattering, for which the signal is identical to neutrino-electron elastic scattering. We require a recoiling electron track to be forward-oriented with $\cos\theta>0.9$ and have energy between $18<E_{e}<50\ \mathrm{MeV}$. Note that because LSND does not have particle ID, a potential background is neutral-current $\pi^{0}$ production where the photons from the $\pi^{0}$ decay fake electrons.

IV.1.2 MiniBooNE

MiniBooNE Aguilar-Arevalo et al. (2009) is a proton beam neutrino experiment currently operating at Fermilab with beam energy 8 GeV. It has been previously noted that MiniBooNE is sensitive to light DM deNiverville et al. (2011, 2012), and indeed, the first limits from a dedicated DM search have recently been published Aguilar-Arevalo et al. (2017). At MiniBooNE, the beam energy is large enough that $\eta$ production and dark bremsstrahlung contribute to DM production, however, contributions from $\pi^{0}$ decays dominate in regions of parameter space where the $\pi^{0}\to\gamma\chi_{1}\chi_{2}$ is allowed. Due to the smaller luminosity, the number of pions produced at MiniBooNE is only $N_{\pi^{0}}\sim 10^{20}$, with the number of etas being further suppressed by a factor of $\sim 30$. We used BdNMC deNiverville et al. (2017) to generate a distribution of $\pi^{0}$ and $\eta$ produced at MiniBoone and then decayed the mesons to DM using Monte Carlo. To simulate dark bremsstrahlung production at MiniBooNE, we used BdNMC to generate a distribution of on-shell $A^{\prime}$ and then decayed the dark photons to DM using Monte Carlo. The rate of dark bremsstrahlung production decreases as $m_{A^{\prime}}$ is increased, but for $m_{A^{\prime}}\gtrsim 200\ \mathrm{MeV}$, kinematic thresholds and $\alpha_{D}$ suppression from off-shell meson decay are significant enough that the dark bremsstrahlung production of DM dominates. A more detailed discussion of the Monte Carlo methods and matrix elements used in these simulations can be found in Appendices C and D.

The MiniBooNE detector also uses scintillating mineral oil, and is approximately spherical with radius $\sim 5\mathrm{m}$, located $\sim 500\mathrm{m}$ from the beam stop. The smaller detector and larger distance from the beam stop means that the geometric acceptance is smaller than LSND, but because the beam energy of MiniBooNE is ten times that of LSND, the DM produced at MiniBooNE has boost factors that are roughly ten times greater, which to some extent compensates for the geometric acceptance as the DM is more collimated. That said, the combination of the smaller proton luminosity and the larger number of background events for a DM search Aguilar-Arevalo et al. (2017) results in the MiniBooNE scattering reach being suppressed by several orders of magnitude compared to LSND. We do not include the MiniBooNE scattering curves in our reach plots as they do not cover new parameter space.111111SBND Antonello et al. (2015) uses the same beamline as MiniBooNE, but with a smaller detector placed closer to the beam stop. We expect the decay reach to be similar to MiniBooNE for all but the largest mass splittings, but the scattering reach should be enhanced due to the higher detector efficiency Van De Water . MiniBooNE is capable of seeing both nucleon and lepton tracks, but electrons are indistinguishable from positrons and neutrons are indistinguishable from protons.121212However, MiniBooNE can distinguish photons from leptons Patterson et al. (2009). We take the efficiency for lepton and nucleon detection to be $35\%$ Patterson et al. (2009).

At MiniBooNE, we impose energy and angular cuts similar to those in Ref Aguilar-Arevalo et al. (2017). We require any lepton track to be forward-oriented with $\cos\theta>0.99$ and have energy in the range $50<E_{e}<600\ \mathrm{MeV}$. We also require the electron and positron to have an angular separation of at least $2^{\circ}$ Patterson et al. (2009) and an energy of at least $50\ \mathrm{MeV}$ to be visible. We determined the decay reach assuming 95% one-sided Poisson c.l., corresponding to 3 events.

IV.2 Signals at Electron Beam Dump Experiments

IV.2.1 E137

The E137 experiment at SLAC Bjorken et al. (1988) was designed to search for axion-like particles and with a 20 GeV beam delivering $\sim$ 30 C ($\sim 2\times 10^{20}$ electrons) of current onto an aluminum target positioned approximately 400 m upstream from an aluminum and plastic scintillator detector. Ref. Batell et al. (2014) demonstrated that the existing data from this sample is sensitive to sub-GeV DM if the DM can be produced in the beam dump and scatter elastically off detector electrons.

The null result of the E137 search can be used to place bounds on our inelastically coupled scenario. The DM yield is computed by evaluating the $\chi$ flux produced in the beam dump via dark-brehmstrahlung (described in Sec. III.2) and considering both DM scattering off electron targets in the detector (described in Sec. III.3.1), as well as $\chi_{2}$ decays inside the detector volume (Sec. III.3.2). To comply with the search criteria from Ref. Bjorken et al. (1988), we demand either signal process deposit at least 1 GeV of energy inside the detector’s geometric acceptance, and we place a 3-event bound to account for the null result of the E137 experiment.

IV.2.2 BDX

The proposed BDX experiment at Jefferson Laboratory Battaglieri et al. (2014, 2016); Izaguirre et al. (2013, 2014) is a dedicated effort to search for light DM produced and detected in analogy with the procedure at E137.131313A similar idea involving lower-energy electron beams installed near large underground neutrino detectors was proposed in Ref. Izaguirre et al. (2015b). The setup involves placing a meter-scale CsI scintillator detector behind the beam dump at the upgraded 11 GeV CEBAF beam. In comparison with E137, BDX has greater luminosity, ($\sim 10^{22}$ electrons on target), a shorter baseline ($\sim 20$ m distance from dump to detector), and a larger detector volume ($50\times 50\times 200\,^{3}$).

In the BDX setup, the boosted $\chi_{1}\chi_{2}$ system emerges from the beam dump and can either scatter (see Sec. III.3.1) via $\chi_{i}e\to\chi_{j}e$ or the excited state can decay via $\chi_{2}\to\chi_{1}e^{+}e^{-}$ Sec. III.3.2) as it passes through the $\sim 2$ meter long detector. We compute decay and scattering projections using 3 and 10 event yield contours, respectively, for EM energy depositions above 300 MeV.

IV.3 Signals at Electron Missing Momentum/Energy Experiments

IV.3.1 NA64

The NA64 experiment at the CERN SPS Banerjee et al. (2017), depicted schematically in Fig. 3 (right panel), is sensitive to DM production via dark bremsstrahlung as described in Sec. III.2. In this setup, $2\times 10^{9}$ electrons with 100 GeV energies are fired into an active ECAL target which measures the energy/momentum and triggers on events with large missing energy. In principle, this technique is sensitive to our scenario of interest because the excited state can decay outside the detector, thereby giving rise to missing energy in the ECAL measurement. While we include this discussion for completeness, we find that the existing data sample does not currently constrain new parameter space, so we do not discuss it further.

IV.3.2 LDMX

The proposed LDMX experiment at SLAC Izaguirre et al. (2015a); Alexander et al. (2016) aims to produce DM using the 4 or 8 GeV LCLS-2 electron beam passing through a thin target upstream of a dedicated tracker and ECAL/HCAL system designed to veto SM particles produced in these collisions – the setup is depicted schematically in Fig. 3 (left panel). For our inelastically coupled DM scenario, an $A^{\prime}$ is produced via dark bremsstrahlung (described in Sec. III.2) and decays via $A^{\prime}\to\chi_{1}\chi_{2}$ followed by a displaced $\chi_{2}\to\chi_{1}e^{+}e^{-}$ decay. If this decay occurs behind the ECAL/HCAL system, it can mimic the missing energy signature which LDMX is optimized to observe.

A representative realization of this setup involves $\sim 10^{16}$ electrons with 4 GeV energies passing through a tungsten target of thickness $\sim 0.1X_{0}$ and emerging from the target with less than 1 GeV of energy. The signal yield scales as

[TABLE]

where $N_{e}=10^{16}$ is the number of electrons on target, $X_{0}\sim 9$ cm is the tungsten radiation length, $n_{W}$ is the tungsten number density, $\sigma(eN\to eNA^{\prime})$ is the $A^{\prime}$ production cross section from Eq. (15) and $\beta/\gamma$ are the $\chi_{2}$ velocity and boost factors, respectively. Here the exponential factor is the $\chi_{2}$ survival probability through a path length of $\ell$ through the LDMX ECAL/HCAL system. The LDMX missing momentum projections present the 3 event yield contours for various parameter benchmarks.

We note in passing that it may also be possible for both NA64 and LDMX to be sensitive to a promptly decaying $\chi_{2}\to\chi_{1}e^{+}e^{-}$ inside their ECAL systems, but the backgrounds for this process are not known and understanding this signal is beyond the scope of the present work.

V Plots & Main Results

We now have all the ingredients in place to assess the potential sensitivity of currently proposed fixed-target experiments to discover inelastic DM models. We begin with a brief discussion of existing constraints.

The parameter space of inelastically coupled DM for mass scales beneath $\sim$ GeV is constrained by precision measurements, B-factories, and previous fixed-target experiments. On the precision front, the anomalous magnetic moment of the muon and electron constrains the interaction strength between the dark photon $A^{\prime}$ and the SM particles Pospelov (2009). On the collider front, both LEP and BaBar set a bound for larger values of $\epsilon$. The former arises from the shift in $m_{Z}$ induced from mixing with $A^{\prime}$ Hook et al. (2011), and the latter from a monophoton and missing mass re-analysis Izaguirre et al. (2013); Essig et al. (2013b); Lees et al. (2017).

Some of the strongest constraints from elastic DM arise from E137 Bjorken et al. (1988); Batell et al. (2014), an electron beam dump experiment, and LSND Auerbach et al. (2001); deNiverville et al. (2011), a proton beam fixed-target neutrino production experiment. Here, we reinterpret the constraints in terms of coannihilating DM. As discussed in Sec. III.3 above, there are two qualitatively different signals: a scattering signal, where $\chi_{1,2}$ up- or downscatter at the detector and produce a recoiling target and possibly an $e^{+}e^{-}$ pair, and a decay signal, where $\chi_{2}$ survives to the detector and decays inside, producing an $e^{+}e^{-}$ pair. The reach of these experiments depends on their ability to distinguish these multiple signals. While E137 is only sensitive to total energy deposits, the angular resolution of LSND means it is potentially sensitive to well-separated $e^{+}/e^{-}$ pairs, which can be distinguished from the fake elastic events we used in estimating the sensitivity in this work and could enhance the sensitivity. However, this would require access to the LSND data as this signal of two charged tracks in the detector is not present in any published analysis. These existing constraints are illustrated in Fig. 6 for $\alpha_{D}$ = 0.1, $m_{A^{\prime}}=3m_{1}$, and various values of $\Delta$. For all but the smallest splittings, the combination of LSND and E137 covers a large portion of the thermal target in the 1-100 MeV range. However, for $2m_{1}+\Delta>m_{\pi^{0}}$, DM production through pions is kinematically forbidden, so we see sharp kinematic cutoffs at the pion threshold.

For comparison, Fig. 7 shows the effect of varying our benchmark parameters (each panel varies one detail relative to the top-left panel of Fig. 6) to demonstrate that these benchmarks are conservative and representative of the viable parameter space. In particular, the top row of Fig. 7 shows how the parameter space in the top-left plot of Fig. 6 changes as we increase and decrease $\alpha_{D}$ while holding all other parameters fixed. Although there is slightly more viable parameter space for the large value $\alpha_{D}=0.5$, this choice is close to the perturbativity limit for abelian dark sectors Davoudiasl and Marciano (2015), so we regard our benchmark choice as a representative and conservative value; choosing a smaller coupling excludes more parameter space on the $y$ vs. $m_{1}$ plane, as we see for $\alpha_{D}=\alpha$ in the same figure. In the bottom row of Fig. 7, we vary our $\Delta$ and $m_{A^{\prime}}$ benchmarks from Fig. 6. In the bottom-left plot, we show the nearly-elastic case of $\Delta=0.01m_{1}$, where the decay signal shuts off and the constraints are dominated by scattering. For comparison, we also show the recent MiniBooNE elastic scattering results Aguilar-Arevalo et al. (2017), for which the beam energy is sufficiently large that the small 1% mass splitting does not affect the reach.

In the bottom-right plot, we show results for a larger hierarchy, $m_{A^{\prime}}/m_{1}=10$. For a given $m_{1}$, $\Delta$, and $\alpha_{D}$, the production rate is decreased as that event now arises from a much heavier $A^{\prime}$. If we parameterize the production rates at $m_{A^{\prime}}/m_{1}=3$ and $m_{A^{\prime}}/m_{1}=10$ as $N_{3}\epsilon^{2}$ and $N_{10}\epsilon^{2}$, respectively, the total decay or scattering yield scales as $N_{3,10}\epsilon^{4}/m_{A^{\prime}}^{4}$. Thus, for a fixed event yield, $\epsilon$ scales linearly with $m_{A^{\prime}}$ but only as $N_{3,10}^{1/4}$. Far from any kinematic boundaries, the sensitivity in $y\propto\epsilon^{2}/m_{A^{\prime}}^{4}$ improves relative to the thermal target since the scaling with $m_{A^{\prime}}$ dominates the scaling with $(N_{3}/N_{10})^{1/4}$. However, the reach at large masses degrades as the $A^{\prime}$ mass approaches the maximum available energy more rapidly and $A^{\prime}$ production shuts off.

We now turn to the potential of new proposals to largely cover the entire parameter space motivated by thermal inelastic DM. We focus on three experiments representative of the setups we have previously discussed: MiniBooNE, BDX, and LDMX, which are proton beam dump, electron beam dump, and missing energy experiments, respectively. For comparison, we have also estimated the Belle II Jaegle (2015) sensitivity only in the monophoton and missing mass channel by rescaling the BaBar result by the appropriate luminosity factor. As discussed in Sec. III, the dominant signal at MiniBooNE is $\chi_{2}$ decay in the detector whenever it is kinematically allowed. Since MiniBooNE has particle ID Yang et al. (2005); Patterson et al. (2009), electrons can in principle be distinguished from photons, and thus a well-separated $e^{+}/e^{-}$ pair and no other activity in the detector is a signal with few irreducible backgrounds. This stands in sharp contrast to the case of elastic DM scattering at MiniBooNE Aguilar-Arevalo et al. (2017), which must always contend with an irreducible neutrino background. Note that the lower boundary of the decay curve is set by the energy threshold and angular resolution according to Eq. (25). We did not attempt a detailed simulation of the $\rho$ resonance, and thus the reach in $y$ of the region very close to $m_{A^{\prime}}\sim m_{\rho}$ may differ slightly from what we show. The missing energy signal at LDMX dominates at low masses, while the decay reach for BDX is similar to that of MiniBooNE at high masses. Indeed, while the $A^{\prime}$ production rate at proton beam experiments is a few orders of magnitude larger than at electron beam experiments, the higher luminosity and beam energy at BDX compensate to give a roughly similar reach. This is advantageous given that different models for the mediator could enhance or suppress production at proton beams compared to electron beams, so the combination of both experiments probes a wide range of models.

Fig. 8 combines the results in this work with the results of Ref. Izaguirre et al. (2016) to show the thermal target over a wide range of DM masses. We see that except for a few isolated masses, the thermal target for coannihilating DM could be well-covered by all three planned experiments below $\sim$1 GeV, and by collider experiments from $\sim$1 GeV to 1 TeV. The scattering signals dominate at low mass below the kinematic threshold $\Delta=2m_{e}$, while the decay signals dominate when kinematically allowed.

VI Conclusion

In this paper we have studied the fixed-target phenomenology of thermal dark matter with inelastic couplings to the SM and proposed a series of new searches for these interactions. These models are an instance of the general case where the relic abundance arises from thermal coannihilation between the halo DM candidate $\chi_{1}$ and an unstable excited state, $\chi_{2}$. Since the heavier state decays away in the early universe, there is no annihilation at later times, and therefore no indirect detection. Furthermore, if the mass difference between these states exceeds $\sim 100$ keV, upscattering at direct detection experiments is kinematically forbidden and loop-induced elastic scattering is small, so this scenario can likely only be discovered or falsified using accelerators. We leave the possibility of one-loop elastic scattering at recently proposed electron direct detection experiments for a future study.

At fixed-target experiments, the inelastic interaction responsible for setting the relic abundance yields a variety of observable signatures arising from the boosted $\chi_{1}\chi_{2}$ system, which is produced in a proton or electron beam collision with target nuclei. Once produced, either state can scatter off particles in a downstream detector, thereby generating an observable signal. In addition, the boosted $\chi_{2}$ can also survive out to the detector and decay semi-visibly via $\chi_{2}\to\chi_{1}e^{+}e^{-}$ to directly deposit a visible signal as it passes through the active volume.

Using these signatures, we have extracted existing constraints on this scenario by reinterpreting old LSND and E137 data. To this end, we have generalized the analyses in deNiverville et al. (2011) (for LSND) and Batell et al. (2014) (for E137), which focused on the scattering signatures of elastically coupled DM. In our analysis, we have demonstrated that there are several new signatures to which these older experiments are sensitive if DM couples inelastically. In particular, we find that E137 and LSND already place nontrivial bounds on the parameter space that yields sub-GeV thermal coannihilation for a variety of DM masses, mass splittings, and coupling strengths.

We have also studied the prospects for future decay and scattering searches at the existing MiniBooNE (proton beam) experiment and the proposed BDX and LDMX (electron beam) experiments. We find that the combined reach of all scattering and decay searches at these experiments can comprehensively test nearly all remaining parameter space, thereby providing strong motivation for these efforts.

This paper also extends earlier work Izaguirre et al. (2016), which studied the collider phenomenology of inelastic thermal coannihilation models over the GeV – TeV mass range. Our work complements this effort by working out the constraints and projections for the MeV–GeV range, thereby providing a roadmap for covering thermal coannihilation over nearly all masses for which a thermal origin is viable (lower masses are in conflict with early universe cosmology and higher masses generically violate perturbative unitarity in most models). This full coverage, spanning the results of both papers, is presented in Fig. 8.

Finally we note that other existing and future experiments may also have powerful reach to this class of models. In particular, the proton beam experiments DUNE Acciarri et al. (2016), SeaQuest Nakahara (2011) (see forthcoming work by Berlin et al. ), SHiP Alekhin et al. (2016), and T2K Hayato (2005) all involve beam energies in excess of 100 GeV, which can produce far more boosted DM than the beams considered in this work (all $<$ 10 GeV and below). Higher energies at these experiments can open up new production modes for the DM candidates (e.g. deep inelastic scattering) and impart greater boosts to the DM system, which can profoundly affect the sensitivity projections for these setups. In addition, liquid argon detectors such as ICARUS Antonello et al. (2013) may be more optimized for seeing the two tracks characteristic of the decay signal.141414We thank Maxim Pospelov for pointing this out. However, working out the implications of these features is beyond the scope of this paper.

Collider experiments may also probe the low-mass scenario we consider in this paper. The superior estimated reach of LHCb to visibly-decaying dark photons Ilten et al. (2015, 2016) suggests that LHCb will also have sensitivity to inelastic DM, especially when $\chi_{2}$ undergoes displaced decays within the detector. However, this requires a dedicated analysis, as the pairs of leptons in this case do not reconstruct a resonance due to the invisible $\chi_{1}$. We leave this interesting analysis to future work.151515We thank the anonymous referee for suggesting this analysis.

In addition, BaBar could be additionally sensitive to this topology through dedicated searches for monophoton and displaced leptons, as proposed by Ref. Izaguirre et al. (2016). Similarly, through analogous searches, the larger luminosities at Belle II could provide unprecedented sensitivity to both displaced and long-lived decays. Parts of the parameter space with larger mass splittings can also lead to displaced vertices from $\chi_{2}$ 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 planned experiments on the horizon, and encourage them to pursue the inelastic DM implications of the signals we have discussed in this work.

Acknowledgements.

Acknowledgments: We thank Brian Batell, Asher Berlin, Nikita Blinov, Maxim Pospelov, Philip Schuster, Brian Shuve, Tim Tait, Natalia Toro, Richard Van De Water, and Yiming Zhong for helpful conversations. We thank Jordan Smolinsky for pointing out typos in an earlier version of this paper. Fermilab is operated by Fermi Research Alliance, LLC, under Contract No. DE-AC02-07CH11359 with the US Department of Energy. EI is supported by the United States Department of Energy under Grant Contract desc0012704.

Appendix A Relic Abundance

The relic abundance of $\chi_{1}$ is governed by a Boltzmann equation whose collision terms involve thermally averaged coannihilation, decay, and inelastic scattering processes. Defining the dimensionless comoving yield as $Y_{i}\equiv n_{i}/s$, where $s(T)=2\pi^{2}g_{s,*}T^{3}/45$ is the entropy density and $g_{*,s}$ is the number of entropic degrees of freedom, the Boltzmann system can be written as

[TABLE]

where $x\equiv m_{2}/T$ is a dimensionless time variable and $Y_{i}^{(0)}$ is the comoving equilibrium yield for species $i$. We define $\lambda_{A,S,D}$, to be the dimensionless annihilation, scattering, and decay rates

[TABLE]

and we use the Hubble rate $H(T)\simeq 1.66\sqrt{g_{*}}T^{2}/m_{Pl}$ during radiation domination, where $g_{*}$ is the number of relativistic degrees of freedom, and $m_{Pl}=1.2\times 10^{19}\,$ is the Planck mass. Solving this system yields the asymptotic value $Y_{1}^{\infty}$ at freeze-out near $x\sim 20$, which determines the relic abundance

[TABLE]

where $\rho_{\rm cr}=8.1h^{2}\times 10^{-47}\,{}^{4}$ is the critical density and $s_{0}\simeq 2969\,^{-3}$ is the present day CMB entropy. An example solution to this system for a representative model point is presented in Figure 9.

Appendix B Scattering and Annihilation Rates

B.1 Coannihilation Rate

The amplitude for this $\chi_{1}(p_{1})\chi_{2}(p_{2})\to f(k_{1})\bar{f}(k_{2})$ coannihilation is

[TABLE]

so squaring and averaging initial state spins gives

[TABLE]

The differential cross section for this process is

[TABLE]

where $|\vec{p}_{*}|$ and $|\vec{k}_{*}|$ are the initial and final momenta in the CM frame and

[TABLE]

Integrating this result, the total cross section becomes

[TABLE]

where we have taken the elastic limit $m_{1}=m_{2}=m_{\chi}$ for illustrative purposes, but we retain the full inelastic mass dependence in our calculations.

Generalizing the derivation in Gondolo and Gelmini (1991) for coannihilation, the thermally averaged cross section is

[TABLE]

where $s_{0}\equiv(m_{1}+m_{2})^{2}$, the thermal averaging factor is

[TABLE]

and $K_{1,2}$ are modified Bessel functions of the first and second kinds.

Note that in our calculation of the relic density, we account for annihilation to hadrons (e.g. $\chi_{1}\chi_{2}\to{A^{\prime}}^{*}\pi^{+}\pi^{-}$) by following the procedure described in Izaguirre et al. (2016) where the final state phase space is extracted from the measured distribution $R(s)\equiv\sigma(e^{+}e^{-}\to{\rm hadrons})/\sigma(e^{+}e^{-}\to{\rm\mu^{+}\mu^{-}})$.

B.2 DM-SM Scattering Cross Section

For pseudo-Dirac DM, with mass splitting $\Delta\equiv m_{2}-m_{1}$, the matrix element for the process $\chi_{1}(p_{1})T(p_{2})\to\chi_{2}(k_{1})T(k_{2})$ is given by

[TABLE]

so the squared, spin-average matrix element is

[TABLE]

The differential scattering cross section in the CM frame is

[TABLE]

where initial/final state momenta satisfy

[TABLE]

In terms of lab frame quantities with a stationary target $T$,

[TABLE]

so we can change variables to obtain the differential recoil distribution

[TABLE]

where $E_{T}\equiv E_{k_{2}}$ is the targets recoil energy. Thus, we have

[TABLE]

which serves as an input into all detector scattering calculations for both proton and electron beam dump experiments.

Appendix C Matrix Elements for DM Production and Detection

C.1 Meson Decay

The matrix element for pseudoscalar meson decay $m^{0}(p_{1})\to\gamma(k_{1})\chi_{1}(k_{2})\chi_{2}(k_{3})$ ($m^{0}=\pi^{0},\eta$) is given by

[TABLE]

where $f_{m^{0}}$ is the meson decay constant, $\Gamma_{A^{\prime}}$ is the total $A^{\prime}$ width, $q\equiv p_{1}-k_{1}=k_{2}+k_{3}$, and $s\equiv q^{2}$. The spin-averaged square of the matrix element is

[TABLE]

If $m_{1}+m_{2}\ll m_{A^{\prime}}\ll m_{m^{0}}$ then we can reasonably make the narrow width approximation Kahn et al. (2015) and take

[TABLE]

The decay width is given by

[TABLE]

where we have defined the function

[TABLE]

Here, $d\Omega_{\gamma}^{*}$ refers to angles in the $m^{0}$ rest frame and $d\Omega_{\chi}^{*}$ refers to angles in the $\chi_{1}\chi_{2}$ CM frame.

C.2 Excited State Decay

The matrix element for $\chi_{2}(p_{1})\to\chi_{1}(k_{1})f(k_{2})\bar{f}(k_{3})$ is given by

[TABLE]

where again $\Gamma_{A^{\prime}}$ is the total $A^{\prime}$ decay rate, $q\equiv p_{1}-k_{1}$, and $s\equiv q^{2}$. In this paper we will consider $f=e^{-}$ exclusively; decay to muons is allowed only for the largest masses and splittings but may provide a distinctive signal at higher-energy experiments. The spin-averaged square is

[TABLE]

We note that because we only consider $m_{A^{\prime}}>m_{1}$ and $\Delta<m_{1}$ in this paper, the $A^{\prime}$ is always off-shell and we never make the narrow width approximation for $\chi_{2}$ decays.

The decay width is given by

[TABLE]

where $\beta$ is defined as in Eq. (52), $d\Omega_{1}^{*}$ refers to angles in the $\chi_{1}$ rest frame, and $d\Omega_{f}^{*}$ refers to angles in the $f\bar{f}$ CM frame.

C.3 DM-SM Scattering

The tree-level matrix elements for scattering $\chi_{i}(p_{1})T(p_{2})\to\chi_{j}(k_{1})T(k_{2})$ off a pointlike fermionic target have already been computed above in App. B.2. Since we will also be interested in scattering off targets with substructure such as nucleons and nuclei, we consider the more general scattering process where $T$ is a fermionic target with both a monopole and dipole coupling to electromagnetism. The matrix element for this process is given by

[TABLE]

where $q\equiv p_{1}-k_{1}$ is the four-momentum carried by the virtual photon and $t\equiv q^{2}$ is the Mandelstam variable. The Lorentz structure $\Gamma^{\mu}$ at the target vertex is

[TABLE]

Here, $\sigma^{\mu\nu}\equiv\frac{i}{2}\left[\gamma^{\mu},\gamma^{\nu}\right]$, and $F_{1}(q^{2})$ and $F_{2}(q^{2})$ are the electric monopole and dipole form factors which depend on the target $T$. For the purposes of this paper, we take

[TABLE]

and

[TABLE]

with $\kappa_{p}\approx 1.79$ and $\kappa_{n}\approx-1.9$ deNiverville et al. (2011).

The spin-averaged square of the matrix element is

[TABLE]

where $\tau\equiv-\frac{t}{4m_{T}^{2}}$.

For an incoming $\chi_{i}$ with the target at rest in the lab frame, the lab-frame differential cross section is then given by

[TABLE]

where $\beta$ is defined as in Eq. (52).

We note that for the case of coherent scattering off of nuclei ($T=N$), we make the additional insertion of the Helm form factor Izaguirre et al. (2014); Duda et al. (2007).

Appendix D Monte Carlo Techniques

Our simulation performs two distinct operations: production of $\chi_{1}$ and $\chi_{2}$ pairs and the detection of $\chi_{1}$ $\chi_{2}$ pairs in a detector. Many of the generic techniques used in our simulations such as numerical phase space integration, rejection sampling of differential probability distributions, and computations utilizing detector geometries were borrowed from or influenced by BdNMC deNiverville et al. (2017).

D.1 DM Production Simulation for Proton and Electron Beams

We use two simulation pipelines, one for proton beams and one for electron beams. Our proton beam production simulation takes as input an unweighted list of four-momenta from one of the DM progenitors that we considered: $\pi^{0}$, $\eta$, and $A^{\prime}$ from dark bremsstrahlung. The output of the production simulation is an unweighted list of $\chi_{1}$ and $\chi_{2}$ particles produced. Each $\chi_{1}\chi_{2}$ four-momentum pair in the output list is randomly sampled from the differential decay rate of the progenitor via a rejection sampling method similar to that used in BdNMC.

For electron beam production, we used a modified version of Madgraph@NLO Alwall et al. (2014) from Ref. Chen et al. (2017), which creates a new physics model which contains new particle containers for the nucleus, $A^{\prime}$ and $\chi_{1,2}$ system. It utilizes elastic and inelastic atomic form factors from Ref. Kim and Tsai (1973). The elastic component is given by

[TABLE]

where $a=111\,Z^{-1/3}/m_{e}$, and $d=0.164\mbox{ GeV}^{2}A^{-2/3}$. We also include a quasi-elastic (inelastic) term:

[TABLE]

where $a^{\prime}=773\,Z^{-2/3}/m_{e}$, $m_{p}$ is the proton mass, and $\mu_{p}=2.79$. The general form factor is then

[TABLE]

D.2 DM Detection Simulation for Neutrino Detectors

Our detector simulation takes as input the list of $\chi_{1}\chi_{2}$ four-momentum pairs from the production simulation. We assume that the various production processes all take place at the beam stop, so that the trajectories of each DM particle is well defined. There are three ways that a $\chi_{1}\chi_{2}$ pair produced at the beam stop can be detected: the $\chi_{2}$ can reach the detector and either decay or downscatter or the $\chi_{1}$ can reach the detector and upscatter.

For each $\chi_{2}$ in the input list, a decay length $x$ is randomly selected from the appropriate exponential distribution. If the $\chi_{2}$ trajectory intersects with the surface of the detector, we calculate the path length $l$ along the trajectory from the beam stop to the detector. If $x>l$, i.e. if the $\chi_{2}$ persists until it reaches the detector, then the $\chi_{2}$ can be detected via direct decay or downscattering. If $x<l$ or if the $\chi_{2}$ does not intersect the detector, we decay the $\chi_{2}\to\chi_{1}e^{+}e^{-}$ via rejection sampling of the differential decay rate. We then process the resulting $\chi_{1}$ four-momentum in exactly the same way as a $\chi_{1}$ produced from primary progenitors, but now accounting for the fact that this $\chi_{1}$ trajectory starts at the point where the $\chi_{2}$ decayed rather than the beam stop.

Each $\chi_{2}$ that reaches the detector can either decay or downscatter. We process decay events by weighting by the total probability for the $\chi_{2}$ to decay inside the detector and performing the decay $\chi_{2}\to\chi_{1}e^{+}e^{-}$ via rejection sampling of the differential decay rate. We then process the final state $e^{+}e^{-}$ pair by applying the kinematic cuts described in Sec. IV. We process scattering events by summing over targets and consider scattering off of each target $T$ independently. To avoid double counting with decay events, we make the conservative assumption that the $\chi_{2}$ can only scatter if it does not decay in the detector. We therefore weight by the probability that the $\chi_{2}$ persists through the detector and the independent probability that the $\chi_{2}$ scatters in the detector and perform the scattering $\chi_{2}T\to\chi_{1}T$ by sampling from a uniform distribution in the center of mass variables $\cos\theta$ and $\phi$, which uniquely specify the final state kinematics. Because we sample our final state kinematic variables from a uniform distribution rather than the true distribution of final states,

[TABLE]

which is proportional to the differential cross section, we must also weight each event by the factor $4\pi\mathcal{P}(x_{1},...,x_{n})$. This weighting scheme enables a cancellation of the total cross section $\sigma_{T}$, making a Monte Carlo computation of $\sigma_{T}$ unnecessary and thereby saving significant computational complexity. The cancellation occurs because each event is also weighted by the total probability for scattering which is $\propto\sigma_{T}$ when Taylor-expanded. Once the scattering final state is sampled and the weights computed, we apply the cuts described in Sec. IV to the recoiling target.

For each $\chi_{1}$ that reaches the detector, we process the scattering $\chi_{1}T\to\chi_{2}T$ using the same method as $\chi_{2}$ downscattering. Additionally, we weight by the probability that the upscattered $\chi_{2}$ will decay in the detector and perform the decay via rejection sampling. The $e^{+}e^{-}$ pair from the de-excitation and the recoiling target are both included when we apply the kinematic cuts described in Sec. IV.