capbtabboxtable[][]

# Spin-1 Thermal Targets for Dark Matter Searches at Beam Dump and Fixed
Target Experiments

Riccardo Catena  and Taylor R. Gray

###### Abstract

The current framework for dark matter (DM) searches at beam dump and fixed
target experiments primarily relies on four benchmark models, the so-called
complex scalar, inelastic scalar, pseudo-Dirac and finally, Majorana DM
models. While this approach has so far been successful in the interpretation
of the available data, it a priori excludes the possibility that DM is made of
spin-1 particles – a restriction which is neither theoretically nor
experimentally justified. In this work we extend the current landscape of sub-
GeV DM models to a set of models for spin-1 DM, including a family of
simplified models (involving one DM candidate and one mediator – the dark
photon) and an ultraviolet complete model based on a non-abelian gauge group
where DM is a spin-1 Strongly Interacting Massive Particle (SIMP). For each of
these models, we calculate the DM relic density, the expected number of signal
events at beam dump experiments such as LSND and MiniBooNE, the rate of energy
injection in the early universe thermal bath and in the Intergalactic Medium
(IGM), as well as the helicity amplitudes for forward processes subject to the
unitary bound. We then compare these predictions with experimental results
from Planck, CMB surveys, IGM temperature observations, LSND, MiniBooNE, NA64,
and BaBar and with available projections from LDMX and Belle II. Through this
comparison, we identify the regions in the parameter space of the models
considered in this work where DM is simultaneously thermally produced,
compatible with present observations, and within reach at Belle II and, in
particular, at LDMX. We find that the simplified models considered here are
strongly constrained by current beam dump experiments and the unitarity bound,
and will thus be conclusively probed (i.e. discovered or ruled out) in the
first stages of LDMX data taking. We also find that the vector SIMP model
explored in this work predicts the observed DM relic abundance, is compatible
with current observations and within reach at LDMX in a wide region of the
parameter space of the theory.

##  1 Introduction

The lack of discovery of Weakly Interacting Massive Particles (WIMPs) at dark
matter (DM) direct detection experiments has motivated the exploration of a
variety of alternative theoretical and experimental paradigms over the past
decade [1, 2]. In this exploration, emphasis has been placed on probing DM
candidates lying outside the canonical WIMP mass window, with most of these
efforts focusing on the MeV – GeV mass range [3]. This choice is supported by
at least three reasons [4]. First, a DM candidate lighter than a nucleon would
not carry enough kinetic energy to induce an observable nuclear recoil in a
direct detection experiment, thereby explaining in a simple and economical way
the lack of discovery of WIMPs. Second, the present cosmological density of
particles in this mass range can match the one observed for DM by the Planck
satellite [5]. This can occur via the chemical decoupling mechanism if the new
sub-GeV states have interactions that involve new particle mediators in the
same mass range, thus evading the Lee-Weinberg bound [6]. Finally, the sub-GeV
DM hypothesis can be tested experimentally using existing methods, including
direct detection experiments sensitive to DM-induced electronic excitations in
materials, as well as beam dump and fixed target experiments. Especially
important for this work are the operating beam dump experiments LSND [7] and
MiniBooNE [8], and the fixed target experiment LDMX [9], as they can
effectively probe sub-GeV DM models with DM-electron or -nucleon scattering
cross sections that are suppressed by small DM velocities or momentum
transfers.

The theoretical framework currently used in the analysis of operating beam
dump experiments, as well as in assessing the prospects of next-generation
beam dump and fixed target experiments consists of four benchmark models [10,
11, 12, 13, 14, 15, 16], often referred to as complex scalar DM, scalar
inelastic DM, pseudo-Dirac DM, and, finally, Majorana DM. While reviewing
these four models goes beyond the scope of this introduction, it is apparent
that this theoretical framework a priori excludes the possibility that DM
consists of spin-1 particles. However, as pointed out by different groups in a
series of recent works focusing on the direct detection of vector DM [17, 18,
19, 20, 21], there is no theoretical or experimental argument supporting this
restriction.

The main purpose of this work is to extend the current framework for DM
searches at beam dump and fixed target experiments to the case of spin-1 DM.
As a first step towards this extension, we study the phenomenology of a set of
simplified models for vector DM featuring a kinetic mixing between an ordinary
and a “dark photon” [22, 23]. The latter is responsible for mediating the
interactions between DM and the known electrically charged particles. These
simplified models conceptually extend the Standard Model (SM) of particle
physics in the same minimal way as the four benchmark models listed above.
Next, we will focus on the phenomenology of an ultraviolet complete model
where DM is made of Strongly Interacting Massive Particles (SIMPs) [24]. In
both cases, we identify the regions in the parameter space of the theory where
DM can be thermally produced, is not excluded by current experiments and,
finally, is within reach at LDMX [25]. One usually refers to these regions as
thermal targets.

We find that the simplified models for spin-1 DM that we consider here are
subject to strong constraints from existing beam dump experiments, as well as
from the unitarity of the S-matrix. We also find that the regions in the
parameter space of these models that are not already ruled out by existing
theoretical and experimental constraints will soon be probed at LDMX, which
will conclusively discover or exclude this family of models in the early stage
of data taking. In contrast, in the case of vector SIMP DM, we find that the
observed DM cosmological density can be reproduced in a broader region of
parameter space. A significant fraction of this is not excluded by existing
beam dump experiments, and is within reach at LDMX.

This article is organised as follows. We start by introducing the spin-1 DM
models explored in our work in Sec. 2. We then review the experimental and
theoretical constraints these models have to fulfill in Sec. 3. The main
results of our analysis are reported in Sec. 4, where we identify the regions
of the parameter space of our models in which DM is not ruled out by current
experiments, is thermally produced, and is within reach at LDMX. Finally, we
summarise and conclude in Sec. 5. Useful scattering cross section formulae are
listed in the appendices.

##  2 Models for vector dark matter

In this section, we provide a brief review of the spin-1 DM models that we
consider in this work. We focus on so-called “simplified models” featuring one
DM candidate and a single new particle mediator in the mass spectrum, as well
as on a renormalisable, ultraviolet complete model where the DM candidate is a
Strongly Interacting Massive Particle (SIMP). The first framework enables us
to extend the study of DM at fixed target and beam dump experiments in a way
that can directly be compared with the existing literature on scalar and
fermionic DM. As we will see, this first approach is constrained by existing
data and subject to strong bounds from the unitarity of the S-matrix, and will
therefore be conclusively probed in the early stages of LDMX data taking. The
second framework is by construction compatible with the unitarity of the
S-matrix, and it complements our simplified model analysis by focusing on a
different mechanism to explain the present DM cosmological density.
Specifically, the relic abundance of SIMPs is set by so-called $3\rightarrow
2$ processes and $2\rightarrow 2$ forbidden annihilations, whereas in the
simplified models we consider here the relic abundance is set by pair
annihilation into visible particles.

###  2.1 Simplified models

We start our exploration of spin-1 DM by considering a general set of
simplified models [26]. These models extend the Standard Model (SM) of
particle physics by a complex vector field, $X^{\mu}$, playing the role of DM,
and a mediator particle described by the real vector field $A^{\prime\mu}$.
The following Lagrangian specifies the interactions between the complex and
real vector fields $X^{\mu}$ and $A^{\prime\mu}$, respectively, and the Dirac
spinors $f$ associated with the electrically charged SM fermions,

$\displaystyle\mathscr{L}=$
$\displaystyle-\left[ib_{5}X_{\nu}^{\dagger}\partial_{\mu}X^{\nu}A^{\prime\mu}+b_{6}X_{\mu}^{\dagger}\partial^{\mu}X_{\nu}A^{\prime\nu}+h.c.\right]$
$\displaystyle-\left[b_{7}\varepsilon_{\mu\nu\rho\sigma}\left(X^{\dagger\mu}\partial^{\nu}X^{\rho}\right)A^{\prime\sigma}+h.c.\right]$
$\displaystyle-h_{3}A^{\prime}_{\mu}\overline{f}\gamma^{\mu}f,$ (2.1)

where, as anticipated, $f$ includes electrons, muons, taus, and quarks
(neutrinos are not included). The first line in Eq. (2.1) describes
interactions that can be generated in models for non-abelian spin-1 DM, as
reviewed in Sec. 2.2 and shown in detail in [24]. The strength of these
interactions is parametrized by the coupling constants $b_{5}$ and $b_{6}$.
Without loss of generality $b_{5}$ can be taken to be real, while the coupling
constant $b_{6}$ is in general complex. The second line in Eq. (2.1) describes
interactions that can arise from abelian spin-1 DM models [27], and is
characterised by the in general complex coupling constant $b_{7}$. Finally,
the last line in Eq. (2.1) corresponds to the coupling between the
electrically charged SM fermions and a “dark photon”, here associated with the
vector field $A^{\prime\mu}$. In order to make the analogy with the dark
photon model explicit, one could identify $h_{3}$ with $h_{3}=e\varepsilon$,
where, $\varepsilon$ is the so-called kinetic mixing parameter, which enters
the Lagrangian of the dark photon model via the term
$-(\varepsilon/2)F_{\mu\nu}F^{\prime\mu\nu}$, $F_{\mu\nu}$ and
$F^{\prime\mu\nu}$ being the field strength tensors of the ordinary and dark
photon, respectively. In our numerical applications, we consider the following
cases in which only one $b$-coupling at a time and $h_{3}$ are non-zero:
$(h_{3},b_{5})\neq 0$, $(h_{3},\mathrm{Re}[b_{6}])\neq 0$,
$(h_{3},\mathrm{Im}[b_{6}])\neq 0$, $(h_{3},\mathrm{Re}[b_{7}])\neq 0$, and
$(h_{3},\mathrm{Im}[b_{7}])\neq 0$. To further facilitate the comparison
between our results and the existing literature on the dark photon model, we
later write the non-zero DM couplings in terms of
$\alpha_{D}=g_{D}^{2}/(4\pi)$, where $g_{D}$ is one of $b_{5}$,
$\mathrm{Re}[b_{6}]$, $\mathrm{Im}[b_{6}]$, $\mathrm{Re}[b_{7}]$ or
$\mathrm{Im}[b_{7}]$. As far as the DM and mediator particle mass are
concerned, we denote them by $m_{X}$ and $m_{A^{\prime}}$, respectively. In
the simplified models of this section, they are free, independent parameters.

In this analysis, we restrict ourselves to the case where $m_{A^{\prime}}\geq
2m_{X}$. Under this mass hierarchy, and since $g_{D}\gg e\varepsilon$, dark
photons are produced on shell and decay dominantly to DM at fixed target
experiments. This choice allows for easy comparison with previous literature
on invisible signatures at fixed target experiments [25, 28]. If
$m_{A^{\prime}}\leq 2m_{X}$, the predicted signal yield for missing
energy/momentum signatures is suppressed, resulting in less sensitive
predicted exclusion bounds [10, 29] and visible dark photon decay signatures
becoming more relevant [10].

###  2.2 Non-abelian SIMPs

As a second framework for vector DM, we consider a model where the DM
candidate is a Strongly Interacting Massive Particle, or SIMP [24]. In this
model, the SM gauge group is extended by a local $SU_{X}(2)\times
U(1)_{Z^{\prime}}$ symmetry group under which none of the SM particles is
charged. The corresponding gauge couplings are $g_{X}$ and $g_{Z^{\prime}}$.
The model also features an extended Higgs sector including a scalar singlet,
$S$, and a second scalar $H_{X}$ transforming non-trivially under
$SU_{X}(2)\times U(1)_{Z^{\prime}}$. The gauge bosons associated with the new
symmetry group are denoted by $X_{i,\mu}$, $i=1,2,3$ and $Z^{\prime}_{\mu}$,
or, equivalently, by $X_{\mu}\equiv(X_{1,\mu}+iX_{2,\mu})/\sqrt{2}$,
$X^{\dagger}_{\mu}\equiv(X_{1,\mu}-iX_{2,\mu})/\sqrt{2}$, $X_{3,\mu}$ and
$Z^{\prime}_{\mu}$, respectively. The $SU(2)\times U(1)_{Z^{\prime}}$ gauge
group is spontaneously broken by the vacuum expectation values of $S$, i.e.
$v_{S}$, and $H_{X}$, i.e. $v_{X}$. This generates one complex and two real
mass eigenstates, corresponding to $X_{\mu}$ and two linear combinations of
$Z^{\prime}_{\mu}$ and $X_{3,\mu}$, denoted here by $\tilde{Z}^{\prime}_{\mu}$
and $\tilde{X}_{3,\mu}$, respectively. Their masses are [24],

$\displaystyle m^{2}_{X}$ $\displaystyle=\frac{1}{2}g_{X}^{2}Iv_{X}^{2}\,,$
$\displaystyle m^{2}_{{\tilde{Z}}^{\prime}}$
$\displaystyle=g^{2}_{X}I^{2}v_{X}^{2}\left(1-\cot\theta^{\prime}_{X}\,g_{Z^{\prime}}/g_{X}\right)\,,$
$\displaystyle m^{2}_{{\tilde{X}}_{3}}$
$\displaystyle=g^{2}_{X}I^{2}v_{X}^{2}\left(1+\tan\theta^{\prime}_{X}\,g_{Z^{\prime}}/g_{X}\right)\,,$
(2.2)

where

$\displaystyle\tan(2\theta^{\prime}_{X})=\frac{2c_{X}s_{X}}{c^{2}_{X}-\alpha
s^{2}_{X}}\,,$ (2.3)

with $s_{X}=g_{Z^{\prime}}/\sqrt{g_{X}^{2}+g_{Z^{\prime}}^{2}}$,
$c_{X}=g_{X}/\sqrt{g_{X}^{2}+g_{Z^{\prime}}^{2}}$ and, finally $\alpha\equiv
1+q^{2}_{S}v^{2}_{S}/(I^{2}v^{2}_{X})$. Here $q_{S}$ is the
$U_{Z^{\prime}}(1)$ charge of $S$, while $I$ labels the representation of
$SU_{X}(2)$ under which $H_{X}$ transforms, e.g. $I=1/2$ for a doublet, $I=1$
for a triplet, and $I=3/2$ for a quadruplet. The interaction Lagrangian
containing the cubic self-interactions between these mass eigenstates is [24],

$\displaystyle\mathscr{L}_{3}=$ $\displaystyle-
ig_{X}\cos\theta^{\prime}_{X}\bigg{[}\left(\partial^{\mu}X^{\nu}-\partial^{\nu}X^{\mu}\right)X^{\dagger}_{\mu}{\tilde{X}}_{3,\nu}-\left(\partial^{\mu}X^{\nu\dagger}-\partial^{\nu}X^{\mu\dagger}\right)X_{\mu}{\tilde{X}}_{3,\nu}$
$\displaystyle+X_{\mu}X^{\dagger}_{\nu}\left(\partial^{\mu}{\tilde{X}}^{\nu}_{3}-\partial^{\nu}{\tilde{X}}^{\mu}_{3}\right)\bigg{]}$
$\displaystyle-
ig_{X}\sin\theta^{\prime}_{X}\bigg{[}\left(\partial^{\mu}X^{\nu}-\partial^{\nu}X^{\mu}\right)X^{\dagger}_{\mu}{\tilde{Z}}^{\prime}_{\nu}-\left(\partial^{\mu}X^{\nu\dagger}-\partial^{\nu}X^{\mu\dagger}\right)X_{\mu}{\tilde{Z}}^{\prime}_{\nu}$
$\displaystyle+X_{\mu}X^{\dagger}_{\nu}\left(\partial^{\mu}{\tilde{Z}}^{\prime\nu}-\partial^{\nu}{\tilde{Z}}^{\prime\mu}\right)\bigg{]}.$
(2.4)

The model also predicts quartic interactions between mass eigenstates. These
are relevant in relic density calculations, and have explicitly been
calculated in [24]. Finally, “neutral current” interactions between
$\tilde{Z}^{\prime}_{\mu}$, $\tilde{X}_{3\mu}$ and the charged SM fermions
arise from a kinetic mixing term added to the Lagrangian of this
$SU_{X}(2)\times U(1)_{Z^{\prime}}$ gauge model. Specifically, they are given
by [24]

$\displaystyle\mathscr{L}_{\rm mix}$
$\displaystyle=-e\varepsilon\cos(\theta_{X}^{\prime})\,\tilde{Z}^{\prime}_{\mu}\,\overline{f}\gamma^{\mu}f+e\varepsilon\sin(\theta_{X}^{\prime})\,\tilde{X}_{3\mu}\,\overline{f}\gamma^{\mu}f\,.$
(2.5)

Eq. (2.2) shows that for $\tan\theta^{\prime}_{X}<0$ and
$g_{Z^{\prime}}/g_{X}|\tan\theta^{\prime}_{X}|<1/2$, the predicted mass
hierarchy is

$\displaystyle m_{X}^{2}<m^{2}_{\tilde{X}_{3}}<m_{\tilde{Z}^{\prime}}^{2}\,.$
(2.6)

For example, for $\sin(2\theta_{X}^{\prime})=-0.1$ and $\alpha_{D}\equiv
g_{X}^{2}/(4\pi)=0.5$, Eq. (2.6) is always satisfied for perturbative values
of $g_{Z^{\prime}}$. Consequently, when the mass hierarchy in Eq. (2.2) is
realised, $X_{\mu}$ is a stable DM candidate, while $\tilde{X}_{3}$ and
$\tilde{Z}^{\prime}$ mediate the interactions between $X_{\mu}$ and the SM
fermions via Eq. (2.5). Finally, let us note that for
$\cos(\theta^{\prime}_{X})=1$, Eq. (2.4) reduces to the first line in Eq.
(2.1) with $b_{5}=-2\mathrm{Im}[b_{6}]=g_{X}$, if one integrates by part the
second line in Eq. (2.4) and uses the equation of motion for $X_{\mu}$ (i.e.
Proca equation), which implies $\partial_{\mu}X^{\mu}=0$.

###  2.3 Vector DM phenomenology: general considerations

In this subsection, we highlight the main differences between the models
defined in Sec. 2.1 and in Sec. 2.2 focusing on the predicted DM relic
abundance, kinetic equilibrium, mediator production and decay rates, and,
finally, relativistic DM-electron and -nucleon scattering cross sections.

####  2.3.1 Relic density

Let us denote the DM particle associated with $X_{\mu}$ by $X^{+}$ and the
corresponding DM antiparticle by $X^{-}$. Also, let $n_{X^{+}}$ be the
cosmological number density of DM particles and $n_{X^{-}}$ be the
corresponding cosmological density of DM antiparticles. In the case of the CP-
preserving simplified or SIMP DM models that we consider here,
$n_{X^{+}}=n_{X^{-}}$, and the total DM number density,
$n_{X}=n_{X^{+}}+n_{X^{-}}$, evolves with time according to the following
Boltzmann equation [24],

$\displaystyle\dot{n}_{X}+3Hn_{X}=$ $\displaystyle-\frac{1}{2}\langle\sigma
v_{\rm rel}\rangle_{X^{+}X^{-}\rightarrow
f{\overline{f}}}\,(n_{X}^{2}-n_{X,{\rm eq}}^{2})$
$\displaystyle-\frac{1}{2}\langle\sigma v^{2}_{\rm
rel}\rangle_{X^{+}X^{+}X^{-}\rightarrow
X^{+}\tilde{X}_{3}}\,\left(n_{X}^{3}-n_{X}n_{X,{\rm eq}}^{2}\right)$
$\displaystyle+2\,\langle\sigma v_{\rm
rel}\rangle_{\tilde{X}_{3}\tilde{X}_{3}\rightarrow
X^{+}X^{-}}\,n_{\tilde{X}_{3},{\rm
eq}}^{2}\left(1-\frac{n^{2}_{X}}{n^{2}_{X,{\rm eq}}}\right)\,.$ (2.7)

In general, Eq. (2.7) shows that $n_{X}$ evolves from an equilibrium
configuration to a constant co-moving value as a result of the DM chemical
decoupling from the thermal bath, i.e. the so-called freeze-out or chemical
decoupling mechanism. The first term in the right-hand-side of Eq. (2.7)
describes the time evolution of the DM number density due to DM pair
annihilation into SM fermions. For $\varepsilon\sim 10^{-4}$ ($10^{-6}$) and
$m_{X}\sim~{}100$ MeV ($1$ MeV), this term alone can account for the entire DM
relic density. We will explore this DM production channel within the framework
of simplified models for vector DM introduced above. The second and third line
in Eq. (2.7) describes the time evolution of the DM number density due to
$3\rightarrow 2$ processes and $2\rightarrow 2$ forbidden annihilations that
are specific to the vector SIMP model introduced above (i.e. they are zero in
the simplified model framework). The DM to $\tilde{X}_{3}$ mediator mass ratio
determines whether the DM relic abundance is set by $3\rightarrow 2$ processes
or $2\rightarrow 2$ forbidden annihilations, as shown in Fig. 3 of [24]. We
will explore the interplay of these two DM production mechanisms in the case
of SIMP DM.

####  2.3.2 Kinetic equilibrium

In Eq. (2.7) we implicitly assumed the mass hierarchy of Eq. (2.6) for the
SIMP DM model we introduced in Sec. 2.2. We also assumed that the associated
mediator ${\tilde{X}}_{3}$ is in kinetic equilibrium with the SM thermal bath
during the DM freeze-out. Interestingly, the in-equilibrium decay of
${\tilde{X}}_{3}$ into SM particles, combined with the effective scattering of
DM with the ${\tilde{X}}_{3}$ mediator induced by the cubic and quartic
interactions introduced above, serve as mechanisms to keep the DM particles in
kinetic equilibrium during freeze-out, and thereby satisfy the strong
constraints from structure formation on the DM kinetic decoupling temperature
[24]. For a given $g_{X}$, $\sin(2\theta^{\prime}_{X})$, $m_{X}$ and
$m_{\tilde{X}_{3}}$, the value of $\varepsilon$ required for ${\tilde{X}}_{3}$
and, consequently, for the DM particles to be in kinetic equilibrium at the
freeze-out temperature $T_{f}$ can be estimated from

$\displaystyle n_{\tilde{X}_{3},\rm
eq}(T_{f})\Gamma_{\tilde{X}_{3}}>H(T_{f})n_{X,\rm eq}(T_{f})\,,$ (2.8)

where $H(T_{f})$ is the Hubble rate at $T_{f}$, and the equilibrium densities
for $X$ and ${\tilde{X}}_{3}$ are given by

$\displaystyle n_{X,{\rm eq}}$
$\displaystyle=\frac{45x^{2}}{2g_{*s}(T)\pi^{4}}\,s\,K_{2}(x)\,,$
$\displaystyle n_{{\tilde{X}}_{3},{\rm eq}}$
$\displaystyle=\frac{45m_{\tilde{X}_{3}}^{2}x^{2}}{4g_{*s}(T)\pi^{4}m_{X}^{2}}\,s\,K_{2}\Big{(}\frac{m_{\tilde{X}_{3}}x}{m_{X}}\Big{)}\,,$
(2.9)

where $s$ is the entropy density, $g_{*s}(T)$ is the effective number of
entropic relativistic degrees of freedom at the temperature $T$, $K_{2}$ is a
modified Bessel function of the second kind and $x\equiv m_{X}/T$. In Eq.
(2.8), $\Gamma_{\tilde{X}_{3}}$ is the total decay rate of $\tilde{X}_{3}$,
for which the expression is given in [24].

####  2.3.3 Mediator production and decay

Key to our exploration of spin-1 DM is the production and subsequent decay of
mediator particles in fixed target and beam dump experiments. In the case of
the simplified models of Sec. 2.1, we generically expect that $A^{\prime}$
particles are produced through the kinetic mixing term in Eq. (2.1), either
via dark bremsstrahlung or via meson decays. These $A^{\prime}$s are then
expected to decay into a DM particle/anti-particle pair for values of
$\varepsilon$ that are consistent with the observed relic density set via the
first term in the right-hand side of Eq. (2.7). In contrast, in the case of
the vector SIMP DM model introduced in Sec. 2.2, both the $\tilde{X}_{3}$ and
the $\tilde{Z}^{\prime}$ mediator can in principle be produced in fixed target
and beam dump experiments. In our calculations we will focus on
$\sin(2\theta^{\prime}_{X})=-0.1$, which implies that only
$\tilde{Z}^{\prime}$ particles can significantly be produced via the
interaction in Eq. (2.5), and that, once produced, these $\tilde{Z}^{\prime}$
particles will dominantly decay invisibly via gauge bosons self-interactions,
as long as $\varepsilon$ is small compared to $g_{X}$. Consequently, for
$\sin(2\theta^{\prime}_{X})=-0.1$, the simplified and SIMP models behave
similarly from the point of view of “dark vector boson” production at fixed
target and beam dump experiments.

####  2.3.4 Scattering by electrons and nuclei

Finally, we also need the (relativistic) cross sections for DM-electron and
-nucleon scattering in order to compare the simplified models of Sec. 2.1 and
the SIMP DM model of Sec. 2.2 with the experimental constraints reviewed in
Sec. 3. This in particular applies to the analysis of DM direct detection
experiments and of beam dump experiments, where DM particles produced by the
decay of $A^{\prime}$ mediators (in the case of simplified models) or
$\tilde{Z}^{\prime}$ mediators (in the case of SIMPs) are searched for in
electron or nuclear recoil events in a downstream detector.

We calculate these cross sections by implementing the models of Sec. 2 in
FeynRules [30] and then using CalcHEP [31] to generate analytic expressions
for the squared modulus of the spin-averaged scattering amplitudes. We finally
validate the outcome of this symbolic calculation through direct analytical
calculations of a subset of selected cross sections. In appendices B and A, we
list the relativistic DM-electron and -nucleon scattering cross sections that
we find for the simplified models of Sec. 2.1. In appendix C, we report the
scattering cross sections that we obtain for the SIMP DM model of Sec. 2.2 as
explained above. In the calculation of DM-nucleon scattering cross sections,
we use the mediator-nucleon interaction Lagrangian,

$\displaystyle\mathscr{L}_{N}=e\varepsilon F_{1}\bar{\psi}\gamma^{\mu}\psi
A^{\prime}_{\mu}+e\varepsilon\frac{F_{2}}{2m_{N}}\left[\bar{\psi}\sigma^{\mu\nu}\left(\partial_{\nu}\psi\right)+\left(\partial_{\nu}\bar{\psi}\right)\sigma^{\mu\nu}\psi\right]A^{\prime}_{\mu}\,,$
(2.10)

in the case of simplified models. For the case of SIMP DM, we use Eq. (2.10)
but with $F_{1}(F_{2})\rightarrow\cos(\theta^{\prime}_{X})F_{1}(F_{2})$ and
$A^{\prime}\rightarrow\tilde{Z^{\prime}}$. Here, $F_{1}$ and $F_{2}$ are
nuclear form factors, and we list them in appendix A.

##  3 Experimental and theoretical constraints

In this section, we introduce a selection of constraints and projections from,
respectively, operating and future DM search experiments that apply to the
models of Sec. 2. By complementing these constraints and projections with
bounds from the unitarity of the S-matrix, in Sec. 4 we will identify the
regions of the parameters space of the vector DM models we consider in this
work where DM is simultaneously: 1) thermally produced, 2) experimentally
allowed and, finally 3) detectable. We refer to these regions in parameter
space as thermal targets.

Figure 1: Contours consistent with the observed DM abundance for each of the
simplified models for spin-1 DM introduced in Sec. 2.1 (coloured solid
contours) and for scalar, pseudo-Dirac, and Majorana DM (dashed, dotted, and
dash-dot) with $m_{A^{\prime}}=3m_{X}$, and $\alpha_{D}=g_{D}^{2}/(4\pi)=0.5$,
where $g_{X}$ is the corresponding $b$-coupling in the legends. Results are
presented in the plane spanned by the DM particle mass and the standard
“coupling constant” $y$.

###  3.1 Relic density

Accurate measurements of the Cosmic Microwave Background (CMB) angular power
spectrum by the Planck collaboration set strong constraints on the present DM
cosmological density [5]. The spin-1 DM models considered in this work,
introduced in Sec. 2, are capable of producing the observed DM relic abundance
consistent with Planck by the freeze-out mechanism. In the case of simplified
models (Sec. 2.1), the freeze-out of DM pair annihilations into visible SM
particles sets the DM relic density, whereas in the case of vector SIMP DM
(Sec. 2.2) the present DM cosmological density arises from the freeze-out of
$3\rightarrow 2$ and forbidden annihilations, as one can see from Eq. (2.7).
Below, we discuss the two scenarios separately.

In the case of simplified models, we are interested in the region of parameter
space where $m_{A^{\prime}}>2m_{X}$. In this region, the relic abundance is
set dominantly by direct DM annihilation into SM fermions through an
s-channel. Comparing our theoretical predictions based on Eq. (2.1) with CMB
data, we calculate the thermally averaged cross sections for direct
annihilation using the MicrOMEGAS software [32]. We then compute the DM relic
density by using our own Boltzmann solver, which relies on the freeze-out
approximation from [33]. Our software results agree with MicrOMEGAS, although
we include the contributions from DM annihilation into hadronic final states
[34, 35, 36]. Fig. 1 shows the contour lines consistent with the Planck DM
abundance, or relic targets, for each of the simplified models introduced in
Sec. 2.1 in addition to three benchmark models from [10] (and mentioned in the
introduction, namely complex scalar, pseudo-Dirac, and Mjorana DM). Each of
the simplified models for spin-1 DM in Fig. 1 is defined as having all
couplings in Eq. (2.1) set to 0 except for $h_{3}$ and the $b$-coupling
specified in the legends. There are kinks in the curves occurring at
$m_{X}\approx m_{\mu}$, where DM annihilation into muons become kinematically
accessible. The resonance features appear because of resonances in the cross
sections for DM annihilation into hadrons.

Figure 2: DM relic abundance, $\Omega h^{2}$, as a function of $m_{X}$ and
$m_{\tilde{X}_{3}}/m_{X}$ for SIMP DM thermally produced by the freeze-out of
$3\rightarrow 2$ processes. Here we take $\alpha_{D}=0.5$. The horizontal
dashed line corresponds to the observed value of $\Omega h^{2}$. Coloured
lines have been obtained by evaluating Eq. (40) from [24].

In the case of vector SIMP DM, we rely on previus results from [24]. For each
given $\alpha_{D}=g_{X}^{2}/(4\pi)$ ($g_{X}$ is one of the gauge couplings
introduced in Sec. (2.2)) and $\sin(2\theta^{\prime}_{X})$ we extract the
value of the $m_{\tilde{X}_{3}}/m_{X}$ ratio that gives the correct DM relic
density from Fig. 3 in [24]. While for $\alpha_{D}=0.1$ the correct relic
abundance can be obtained for DM masses above approximately 10 MeV, for
$\alpha_{D}=0.5$ vector SIMP DM can account for the whole cosmological
abundance of DM only for masses above about 80 MeV. This effect is illustrated
in Fig. 2, which shows the DM relic abundance as a function of the DM mass and
of the $m_{\tilde{X}_{3}}/m_{X}$ ratio, focusing on the case in which
$3\rightarrow 2$ are the dominant production mechanism. As one can see from
this figure, for $m_{\tilde{X}_{3}}/m_{X}=2$ the DM relic abundance drops to
zero, independently of the DM particle mass. This is due to the fact that for
$m_{\tilde{X}_{3}}/m_{X}>2$, the $3\rightarrow 2$ process is kinematically not
allowed.

###  3.2 Direct Detection

In principle, DM direct searches via electronic transitions in detector
materials located deep underground also place important constraints on sub-GeV
DM models, e.g. [37]. For spin-1 DM, the current most competitive constraints
are from Xenon1T and Xenon10, which only appear in our plots for the $b_{5}$
and $\mathrm{Re}[b_{7}]$ models. We take the 90% C.L. exclusion limits from
the work of [37].

###  3.3 Energy Injection

DM annihilations in the early universe can inject energetic particles into the
photon-baryon plasma or into the intergalactic medium (IGM), altering the CMB
or the IGM temperature, respectively. However, CMB limits are only competitive
for s-wave annihilating DM [38, 39]. Among the models introduced in Sec. 2,
only the model with $\mathrm{Im}[b_{7}]\neq 0$ gives rise to s-wave dominant
annihilation cross sections, thus we include 95% C.L. exclusion bounds from
the CMB for this model. For all remaining spin-1 models in Sec. 2, including
the SIMP DM model, the predicted annihilation cross section is p-wave. CMB
limits on p-wave annihilation cross sections have been computed here [40]. We
find that they are weaker than the IGM limits on the spin-1 models we consider
in this work. Consequently, we set 95% C.L. exclusion limits from measurements
of the IGM temperature extracted from Lyman-$\alpha$ observations on our
p-wave annihilating DM models. Specifically, we require that the predicted
thermally averaged annihilation cross section is smaller than the upper bounds
reported in [41] as a function of the DM mass. In this analysis, we calculate
the relevant annihilation cross sections analytically, and then compare them
with the output of MicrOMEGAS. These IGM upper limits only appear in the top
left corner of our parameter space for certain mass ratios
($m_{A^{\prime}}/m_{X}$) and models.

###  3.4 Beam Dumps and Fixed Target Experiments

A beam of protons or electrons incident on a fixed target creates cascades of
interactions at beam dump experiments. The goal is to produce DM in these
cascades, which can then be detected in a downstream detector. Calculating the
experimental reach of these experiments involves modeling the processes that
give rise to DM production and detection, which is done through Monte Carlo
(MC) simulations. These simulations include the details of the detector, the
beam, the interactions which produce DM, and the DM model. With this
objective, we use a modified version of the numerical tool BdNMC, a beam dump
Monte Carlo software package [12], in which we implement our spin-1 DM models
to simulate DM production and interactions. BdNMC has the benefit of providing
a simple and rigorous framework for simulating the relevant experiments Mini-
Boone and LSND. In the following section, 90% C.L. limits on the model
parameter space calculated from these simulations are presented, showing the
reach of current experiments including beam dump and missing energy/momentum
experiments which give rise to the most competitive constraints on spin-1 DM.
Below, each experiment considered is introduced.

#### LSND and MiniBooNE

At proton beam dump experiments LSND [16, 7] and MiniBooNE [8], a proton beam
is incident on a target and the following chain of interactions occur
producing DM: $pp\to X\pi^{0};\pi^{0}\to\gamma A^{\prime};A^{\prime}\to DMDM$,
where in the inclusive process, $pp\to X\pi^{0}$, $X$ denotes an
unspecified/unmeasured set of particles, whereas $\pi^{0}$ can also be $\eta$
in the case of MiniBooNE where the energy is sufficient. In addition, dark
proton bremsstrahlung also produces DM in the case of MiniBooNE. The produced
DM can then be detected in the downstream detector through DM-electron and for
MiniBooNE also DM-nucleon scattering. We compute the expected number of signal
events using a modified version of the software BdNMC [13], a MC simulation
tool for beam dumps. We added to BdNMC the model dependent DM-electron and
-nucleon scattering cross section reported in the App. A and B as well as the
relevant branching ratios for dark photon decay.

At LSND, 55 non standard events were observed at 90% confidence level. A
factor of 2 is included to account for the uncertainty in the pion production
rate, so we take the 90% C.L. at 110 events. We perform MC simulations of DM
events at LSND, taking into account our model dependent cross sections, and we
draw the contours in $y$ vs $m_{X}$ space that corresponds to 110 signal
events. In this work, we adopt the standard notation
$y=\varepsilon^{2}\alpha_{D}(m_{X}/m_{A^{\prime}})^{4}$.

Similarly, we perform MC simulations of MiniBooNE to calculate the expected
number of signal events. No events were observed at MiniBooNE, thus we take
the contour in $y$ vs $m_{X}$ space at 2.3 events which corresponds to a 90%
C.L. exclusion limit.

#### E137

Limits from the electron beam dump experiment E137 are not competitive with
MiniBooNE and LSND, thus we do not include them in our analysis.

#### NA64

The missing energy experiment NA64 [42], where an electron beam is incident on
a target, aims to produce a DM flux from dark bremsstrahlung and detect these
signals by their missing energy. No signal events have been observed at this
experiment. We report 90% confidence limits extracted from [42], which we
project onto the $y$ vs $m_{X}$ plane.

#### LDMX

LDMX is a future missing momentum experiment [25, 28] with an 8 GeV electron
beam incident on a tungsten target with $10^{16}$ EOT. The ultimate reach at
90% C.L. is calculated using the expected number of DM signal events from
simulations and we project this expected limit onto our parameter space. In
addition, we include the projected 90 % C.L. exclusion limits from the
analysis of [43], where bremsstrahlung photons are converted into vector
mesons and then decay to invisible states, giving an extended sensitivity
reach in the larger DM mass regime.

###  3.5 Monophoton searches

We also include searches for single photon events at $e^{+}e^{-}$ colliders,
where DM is produced through the process $e^{+}e^{-}\to\gamma
A^{\prime},A^{\prime}\to DMDM$.

#### BaBar

The BaBar detector at the PEP-II B-factory searches for a narrow peak in the
missing mass distribution in the events with one high energy photon. BaBar has
observed no signal events, and we take the 90% C.L. limits on $\varepsilon$ vs
$m_{X}$ from [44] and project them onto our parameter space.

#### Belle-II

For the future experiment Belle-II, we use the expected 90% C.L. limits on
$\varepsilon$ vs $m_{X}$ for phase 3 of Belle-II reported in [45], to obtain
the expected exclusion limits on our spin-1 DM model.

###  3.6 Unitarity bound

In general, simplified DM models can violate perturbative unitarity in some
regions of parameter space [46]. Due to the energy dependence introduced in
the cross section by the longitudinal component of the vector DM polarization
vectors, as well as from the underlying derivative couplings, the simplified
models for spin-1 DM introduced in Sec. 2.1 predict large and un-physical
scattering amplitudes. Therefore, their parameter space is subject to a
“unitarity bound”, or in other words a bound of theoretical validity.
Violations of this bound indicate either that the theory is non perturbative
(all terms in the perturbative expansion of the S-matrix are equally
important), or that it is not complete, and thus additional fields have to be
included to cancel out energy dependent, un-physical contributions to
scattering cross sections. Unitarity violation from DM self-scattering was
investigated in [47], while here we calculate the DM-$e^{-}$ scattering
amplitude to determine at which parameters unitarity is violated. Following
[46], we implement the unitary bound by requiring that in the $(y,m_{X})$
plane,

$\displaystyle\lvert\mathrm{Im}(M_{i\rightarrow i}^{J})\rvert\leq 1\,,$
$\displaystyle 2\lvert\mathrm{Re}(M_{i\rightarrow i}^{J})\rvert\leq 1\,,$
(3.1)

where $i\rightarrow i=Xe^{-}\to Xe^{-}$, $J=0$, and

$\mathcal{M}_{Xe^{-}\to
Xe^{-}}^{0}(s)=\frac{\beta}{32\pi}\int^{1}_{-1}d\cos\theta\mathcal{M}_{Xe^{-}\to
Xe^{-}}(s,\cos\theta)\,.$ (3.2)

(a)

(b)

(c)

(d)

Figure 3: Panels (3(a)) and (3(c)). Exclusion limits from MiniBooNE, LSND,
NA64, BaBar, IGM temperature observations and direct detection experiments
(coloured shaded regions) and projected reach for Belle II and LDMX (dashed
line) on the simplified model with $b_{5}\neq 0$ (see Sec. 2.1) in the
parameter space, $y$ vs $m_{X}$ for $\alpha_{D}=0.5$. The contours consistent
with the observed relic abundance ($\Omega h^{2}\approx 0.12$) are drawn in
solid black. The top panel is for $m_{A^{\prime}}=3m_{X}$, whereas the bottom
panel assumes $m_{A^{\prime}}=2.5m_{X}$. The model is strongly constrained by
MiniBooNE. Panels (3(b)) and (3(d)). Same as the left panels, now for complex
scalar DM.

(a)

(b)

(c)

(d)

Figure 4: Same as Fig. 3, now for the $\mathrm{Re}[{b_{6}}]$,
$\mathrm{Im}[{b_{6}}]$, $\mathrm{Re}[{b_{7}}]$, and $\mathrm{Im}[{b_{7}}]$
models with $\alpha_{D}=0.5$ and $m_{A^{\prime}}=3m_{\chi}$. The models with
$\mathrm{Re}[{b_{7}}]\neq 0$ and $\mathrm{Im}[b_{6}]\neq 0$ are strongly
constrained by beam dump experiments, while the model with
$\mathrm{Im}[{b_{7}}]\neq 0$ is ruled out by CMB observations. The
$\mathrm{Re}[b_{6}]$ model predicts the current DM relic abundance while also
being compatible with observations and within reach at LDMX in the mass range
between about 40 MeV and 200 MeV.

(a)

(b)

(c)

(d)

Figure 5: Same as Fig. 4, now for the mass ratio
$m_{A^{\prime}}=2.5~{}m_{\chi}$. For this choice of parameters, the
$\mathrm{Re}[b_{6}]$ model predicts the current DM relic abundance, is
compatible with observations, and is within reach at LDMX in the mass range
between about 10 MeV and 300 MeV.  Figure 6: Current and expected exclusion
limits on vector SIMP DM in the ($m_{X},y$) plane from NA64, MiniBooNE, LSND,
and BaBar, as well as from projections for LDMX and Belle II. The colour code
is the one used in previous figures. The strongest bound on $y$ arises from
N64, rather than from MiniBooNE or LSND. This is due to the fact that the N64
bound on $y$ is model-independent, while the cross section for DM-nucleon and
-electron scattering for SIMP DM is suppressed by cancellations between
contributions from the $m_{\tilde{Z}^{\prime}}$ and
$m_{\tilde{X}^{\prime}_{3}}$ mediators. In this figure, we assume
$\alpha_{D}=0.5$, $m_{A^{\prime}}/m_{X}=3$, and $m_{\tilde{X}_{3}}/m_{X}\simeq
2$. Above the grey area, SIMPs are in kinetic equilibrium at chemical
decoupling, and can constitute the entire cosmological density of DM on the
right of the vertical brown line. Figure 7: Same as Fig. 6, now for
$\alpha_{D}=0.1$.

##  4 Thermal targets identification

We now present our spin-1 thermal targets for DM searches at beam dump and
fixed target experiments. As anticipated, in these regions of parameter space
DM is simultaneously thermally produced, not excluded by existing experimental
results, and within reach at Belle II or LDMX. We discuss the simplified
models of Sec. 2.1 and the vector SIMP model of Sec. 2.2 separately.

###  4.1 Simplified Models

Figs. 3, 4, and 5 summarize the constraints and projections on the $(m_{X},y)$
plane that we obtain as explained in Sec. 3 for the simplified models
introduced in Sec. 2.1. Here,
$y=\varepsilon^{2}\alpha_{D}(m_{X}/m_{A^{\prime}})^{4}$ and the black curves
on these figures show the contours that are consistent with the observed
abundance of DM measured by Planck [5]. The scientific potential of the
upcoming experiment LDMX, the area above the red dashed curves of the figures,
is significant for sub-GeV DM, since it is projected to probe down to much
smaller couplings than previous experiments.

Specifically, the left panels in Fig. 3 show current and expected exclusion
limits on the spin-1 model with $b_{5}\neq 0$ for both
$m_{A^{\prime}}/m_{X}=3$ and $m_{A^{\prime}}/m_{X}=2.5$. The excluded areas in
orange, green, blue and purple correspond to the MiniBooNE, LSND, NA64, and
BaBar experiments, respectively. For comparison, we also include the exclusion
limits on the $y$ coupling obtained in [37] from the null result reported by
the Xenon10 and Xenon1T experiments. Finally, the expected exclusion limits
for LDMX are compared with the expected reach of Belle II. Independently of
the $(m_{A^{\prime}}/m_{X})$ ratio, this model is ruled about by current beam
dump experiments. The right panels in Fig. 3 shows the same set of constraints
and projections now for the familiar complex scalar DM model, which we
recompute as described in Sec. 3 to calibrate our codes, and for comparison.
Indeed, the complex scalar DM model exhibits the same derivative coupling
between DM and the dark photon as in the model with $b_{5}\neq 0$. The only
difference between the two models is in the Feynman rules for incoming and
outgoing DM particles, which in the latter case involves momentum dependent
polarisation vectors. This difference produces a cross section for DM-nucleon
scattering that is enhanced by a factor of $(E_{\boldsymbol{p}}/m_{X})^{2}\gg
1$ relative to the case of the complex scalar DM model, where
$E_{\boldsymbol{p}}$ is the relativistic energy of the incoming DM particle in
the rest frame of the downstream detector. As a result, the MiniBooNE
constraints on the $b_{5}$ model are much stronger than in the case of complex
scalar DM.

Let us now focus on Fig. 4. This figure shows current and expected exclusion
limits on the spin-1 models with $\mathrm{Re}[b_{6}]\neq 0$ and
$\mathrm{Im}[b_{6}]\neq 0$ (top left and top right panels) as well as on the
models with $\mathrm{Re}[b_{7}]\neq 0$ and $\mathrm{Im}[b_{7}]\neq 0$ (bottom
left and bottom right panels). In all panels we assume
$m_{A^{\prime}}/m_{X}=3$, while the colour code is the same as in Fig. 3.
Current beam dump experiments, LSND and MiniBooNE (green and orange shaded
regions respectively), and BaBar (pink shaded regions) are able to rule out
large parts of the $(m_{X},y)$ plane for the $\mathrm{Re}[b_{7}]$,
$\mathrm{Im}[b_{6}]$ and $\mathrm{Im}[b_{7}]$ models, while leaving the
$\mathrm{Re}[b_{6}]\neq 0$ spin-1 model still compatible with observations in
the mass range between approximately 40 MeV and 200 MeV. Remarkably, the
constraints from MiniBooNE on $y$ are much stronger in the case of the
$\mathrm{Im}[b_{6}]$ model than for the $\mathrm{Re}[b_{6}]$ model. This is
due to the fact that for the $\mathrm{Im}[b_{6}]$ model the cross section for
DM-nucleon scattering is enhanced by a $(E_{\boldsymbol{p}}/m_{X})^{2}\gg 1$
factor relative to the analogous cross section for the $\mathrm{Re}[b_{6}]$
model. Here $E_{\boldsymbol{p}}$ is the energy of the DM particle in the
downstream detector rest frame. As anticipated above, a similar enhancement is
also present in the case of the $b_{5}\neq 0$ spin-1 model. Finally, in the
case of the $\mathrm{Re}[b_{6}]$ model, top left panel in Fig. 4, we also
report the unitary bound on $y$ arising from the helicity amplitude for DM-
electron scattering, which is one of the processes directly entering the
calculation of the constraints in Fig. 4.

It should also be noticed that for the $\mathrm{Im}[b_{7}]$ spin-1 model,
Figs. 4(d) and 5(d), the DM annihilation cross section is s-wave dominant in
contrast to being p-wave dominant as in the other models we consider, leading
to a strong constraint from CMB measurements on this scenario.

Fig. 5 shows current and expected exclusion limits on the same models as in
Fig. 4, now assuming the different mass ratio, $m_{A^{\prime}}/m_{X}=2.5$.
While this leaves are conclusions qualitatively unchanged, the range of DM
masses that are compatible with all observations in the case of the
$\mathrm{Re}[b_{6}]$ model is now significantly broader, varying from about 10
MeV to 300 MeV.

###  4.2 Non-abelian SIMPs

We now turn our attention to the case of vector SIMP DM. We start by reviewing
some of the features of the model, and its differences with the simplified
models. In the SIMP model, the DM relic abundance is set by the freeze-out of
$3\rightarrow 2$ processes and forbidden annihilations, and is thus
independent of the coupling $y$. At the same time, for the thermal production
of SIMP DM to work, DM has to be in kinetic equilibrium in the early universe.
As anticipated, this occurs via effective DM-$\tilde{X}_{3}$ mediator
scattering processes combined with in-equilibrium $\tilde{X}_{3}$ decays into
SM particles. Since the decay rate of $\tilde{X}_{3}$ depends on $\varepsilon$
[24], only above a certain critical value for $y$ vector SIMP DM can
effectively be thermally produced in the right amount. As described in Sec. 2,
we calculate this lower bound by using Eq. (2.8). The result of this
calculation defines a region in the parameter space of the model above which
SIMPs are in kinetic equilibrium at their chemical decoupling, and the freeze-
out mechanics can in principle successfully predict the present DM
cosmological density. Whether or not the whole cosmological DM abundance is in
the form of vector SIMPs depends on the choice of $\alpha_{D}$. For example,
for $\alpha_{D}=0.1$ vector SIMPs can constitute the entire DM in our universe
for masses above about 10 MeV, whereas for $\alpha_{D}=0.5$ this is only
possible for masses above about 80 MeV. Here and in the figures below, we
assume the benchmark value for the mixing angle
$\sin(2\theta^{\prime}_{X})=-0.1$.

Keeping these general considerations on SIMP DM in mind, we are now ready to
compare the predictions of the vector SIMP DM model with the experimental and
theoretical constraints of Sec. 3. Fig. 6 shows the current and expected
exclusion limits on the $y$ coupling as a function of the DM mass that we
obtain from a reanalysis of data collected at NA64, MiniBooNE, LSND, and
BaBar, as well as from projections for LDMX and Belle II. We obtain these
exclusion limits by comparing theory and observations as described in Sec. 3,
while the colour code in the figure is the one of the figures in the previous
sections. Remarkably, in the case of SIMP DM the strongest bound on $y$ arises
from NA64, rather than from MiniBooNE or LSND. This is due to the fact that
the NA64 bound on $y$ is model-independent, while the cross section for DM-
nucleon and -electron scattering for SIMP DM is suppressed by cancellations
between contributions from the $m_{\tilde{Z}^{\prime}}$ and
$m_{\tilde{X}^{\prime}_{3}}$ mediators. In Fig. 6, we assume $\alpha_{D}=0.5$,
$m_{A^{\prime}}/m_{X}=3$, and $m_{\tilde{X}_{3}}/m_{X}\simeq 2$. Fig. 7
reports the results of an analogous analysis where we assume a different
combination of parameters, namely $\alpha_{D}=0.1$, $m_{A^{\prime}}/m_{X}=3$,
and $m_{\tilde{X}_{3}}/m_{X}\simeq 2$. The main difference between the
analyses is the range of masses for which the whole cosmological DM abundance
is in the form of vector SIMPs. As explained above, this difference is due to
the different assumptions made for $\alpha_{D}$ in the two figures.

##  5 Conclusion

In this analysis, we extended the current landscape of sub-GeV DM models
considered in the context of various experiments such as MiniBooNE, LSND and
LDMX to a set of models for spin-1 DM, including a general family of
simplified models (involving one DM particle and one mediator – the dark
photon) and an ultraviolet complete model based on a non-abelian gauge group
(now including two mediators and an extended Higgs sector) where DM is a
vector SIMP. For each of these models, we calculated the DM relic density, the
expected number of signal events in beam dump experiments such as LSND and
MiniBooNE, the rate of energy injection in the early universe thermal bath and
in the IGM, as well as the helicity amplitudes for forward processes subject
to the unitary bound. We then compared these predictions with a number of
different experimental results from Planck, CMB observations, direct detection
experiments (Xenon10 and Xenon1T), data on the IGM temperature from Lyman
alpha observations, LSND, MiniBooNE, NA64, and BaBar and with available
projections from LDMX and Belle II. Through this comparison, we identified the
regions in the parameter space of the models considered in this work where DM
is simultaneously thermally produced, compatible with present observations,
and within reach at Belle II and, in particular, at LDMX.

We found that the simplified models for spin-1 DM investigated in our analysis
are strongly constrained by LSND and MiniBooNE, as well as from bounds on the
unitarity of the S-matrix. The only model not already excluded by these
experimental and theoretical constraints is the one characterised by the
coupling constants $h_{3}$ and $\mathrm{Re}[b_{6}]$. For a dark photon to DM
mass ratio of 3 (2.5), this model is compatible with current observations,
within reach at LDMX and admits a thermal DM candidate for DM masses in a
window between about 40 (10) MeV and 200 (300) MeV. In this mass range, the
model has a relic density contour lying very close to current 90% C.L.
exclusion limits from beam dump experiments in the ($m_{X},y$) plane, and will
thus be conclusively probed (i.e. excluded or discovered) in the first LDMX
run.

At the same time, we found that the vector SIMP model explored in this work
admits thermal DM candidates that are not ruled out by beam dump experiments
and within reach at LDMX in a wide region of the underlying parameter space.
The model features a DM production mechanism that is complementary to the
freeze-out of DM pair annihilations into SM particles of simplified models,
and is based on the interplay of $3\rightarrow 2$ processes and forbidden
annihilations. It also exhibits a lower bound on the DM particle mass arising
from the relic density constraint. The larger $\alpha_{D}$, the larger the
minimum admissible DM particle mass.

Ultimately, our investigation bridges a gap in the current knowledge of sub-
GeV DM by providing the DM community with new sub-GeV spin-1 thermal targets
lying in the experimentally accessible region of next-generation beam dump and
fixed target experiments such as LDMX.

## Acknowledgments

We would like to thank Felix Kahlhoefer and Chris Chang for pointing out the
importance of the unitarity bound on simplified models for spin-1 DM. We would
also like to thank Avik Banerjee and Gabriele Ferretti for drawing our
attention to the work by Choi et al. [24] on vector SIMP DM. In addition, we
greatly appreciate the important discussions we had with Patrick deNiverville
on beam dump experiment simulations, in particular his software BdNMC. The
research contained in this article was performed within the Knut and Alice
Wallenberg project grant Light Dark Matter (Dnr. KAW 2019.0080). We would like
to thank all participants in the project for valuable discussions on sub-GeV
DM and the physics reach of LDMX during our weekly collaboration meetings.
R.C. acknowledges support from individual research grants from the Swedish
Research Council, Dnr. 2018-05029 and Dnr. 2022-04299.

##  Appendix A Cross sections for relativistic dark matter-nucleon scattering

In this appendix, we provide analytic expressions for the differential cross
sections for relativistic DM-nucleon scattering for the spin-1 DM models of
Sec. 2. As anticipated, we obtain these cross sections by implementing the
models of Sec. 2 in FeynRules [30] and then using CalcHEP [31] to generate
analytic expressions for the squared modulus of the spin-averaged scattering
amplitude. We finally validate the outcome of this symbolic calculation
through direct analytical calculations of a subset of selected cross sections.
In the laboratory frame, we find

$\displaystyle\frac{{\rm
d}\sigma_{NX}(E_{\boldsymbol{p}},E_{\boldsymbol{p}^{\prime}})}{{\rm
d}E_{\boldsymbol{p}^{\prime}}}=\frac{1}{32\pi
m_{N}(E_{\boldsymbol{p}}^{2}-m_{X}^{2})}\overline{|\mathcal{M}_{NX}(E_{\boldsymbol{p}},E_{\boldsymbol{p}^{\prime}})|^{2}}\,,$
(A.1)

where $E_{\boldsymbol{p}}$ ($E_{\boldsymbol{p}^{\prime}}$) is the initial
(final) DM particle energy and the squared scattering amplitude is given by

$\displaystyle\overline{|\mathcal{M}_{NX}(E_{\boldsymbol{p}},E_{\boldsymbol{p}^{\prime}})|^{2}}$
$\displaystyle=\frac{\eta^{2}\left[F_{1}^{2}\mathcal{A}_{\eta}(E_{\boldsymbol{p}},E_{\boldsymbol{p}^{\prime}})+F_{2}^{2}\mathcal{B}_{\eta}(E_{\boldsymbol{p}},E_{\boldsymbol{p}^{\prime}})+F_{1}F_{2}\mathcal{C}_{\eta}(E_{\boldsymbol{p}},E_{\boldsymbol{p}^{\prime}})\right]}{3m_{X}^{4}\left[2m_{N}(E_{\boldsymbol{p}}-E_{\boldsymbol{p}^{\prime}})+m_{A^{\prime}}^{2}\right]^{2}}\,.$
(A.2)

Here, a dependence of the nucleon form factors $F_{1}(q^{2})$ and
$F_{2}(q^{2})$ on the momentum transfer $q=p-p^{\prime}$ is understood. The
three functions,
$\mathcal{A}_{\eta}(E_{\boldsymbol{p}},E_{\boldsymbol{p}^{\prime}})$,
$\mathcal{B}_{\eta}(E_{\boldsymbol{p}},E_{\boldsymbol{p}^{\prime}})$ and
$\mathcal{C}_{\eta}(E_{\boldsymbol{p}},E_{\boldsymbol{p}^{\prime}})$ in Eq.
(A.2) are model dependent and have dimension of mass to the eighth power, so
that
$\overline{|\mathcal{M}_{NX}(E_{\boldsymbol{p}},E_{\boldsymbol{p}^{\prime}})|^{2}}$
is dimensionless. $\eta$ is the coupling constant characterising the
underlying DM model. Below, we specify
$\mathcal{A}_{\eta}(E_{\boldsymbol{p}},E_{\boldsymbol{p}^{\prime}})$,
$\mathcal{B}_{\eta}(E_{\boldsymbol{p}},E_{\boldsymbol{p}^{\prime}})$ and
$\mathcal{C}_{\eta}(E_{\boldsymbol{p}},E_{\boldsymbol{p}^{\prime}})$ for
different choices of $\eta$.

  * •

For $\eta=b_{5}$, we find

$\displaystyle\mathcal{A}_{b_{5}}(E_{\boldsymbol{p}},E_{\boldsymbol{p}^{\prime}})$
$\displaystyle=8m_{N}\left[E_{\boldsymbol{p}^{\prime}}\left(2E_{\boldsymbol{p}}m_{N}+m_{X}^{2}\right)-E_{\boldsymbol{p}}m_{X}^{2}\right]\sum_{s}\varepsilon^{s}_{\mu}\varepsilon^{s\mu*}\,,$
$\displaystyle\mathcal{B}_{b_{5}}(E_{\boldsymbol{p}},E_{\boldsymbol{p}^{\prime}})$
$\displaystyle=2m_{N}\left(E_{\boldsymbol{p}}-E_{\boldsymbol{p}^{\prime}}\right)\left[E_{\boldsymbol{p}}^{2}+2(E_{\boldsymbol{p}}+m_{N})E_{\boldsymbol{p}^{\prime}}+E_{\boldsymbol{p}^{\prime}}^{2}-2E_{\boldsymbol{p}}m_{N}-4m_{X}^{2}\right]$
$\displaystyle\times\sum_{s}\varepsilon^{s}_{\mu}\varepsilon^{s\mu*}\,,$
$\displaystyle\mathcal{C}_{b_{5}}(E_{\boldsymbol{p}},E_{\boldsymbol{p}^{\prime}})$
$\displaystyle=-8m_{N}\left(E_{\boldsymbol{p}}-E_{\boldsymbol{p}^{\prime}}\right)\left(-m_{N}E_{\boldsymbol{p}^{\prime}}+E_{\boldsymbol{p}}m_{N}+2m_{X}^{2}\right)\sum_{s}\varepsilon^{s}_{\mu}\varepsilon^{s\mu*}\,,$
(A.3)

where the sum over spin configurations of the product of DM polarisation
vectors is given by

$\displaystyle\sum_{s}\varepsilon^{s}_{\mu}\varepsilon^{s\mu*}=\left\\{E_{\boldsymbol{p}}^{2}m_{N}^{2}+m_{N}E_{\boldsymbol{p}^{\prime}}\left[m_{N}E_{\boldsymbol{p}^{\prime}}-2\left(E_{\boldsymbol{p}}m_{N}+m_{X}^{2}\right)\right]+2E_{\boldsymbol{p}}m_{N}m_{X}^{2}+3m_{X}^{4}\right\\}\,.$
(A.4)

  * •

For $\eta=\mathrm{Re}(b_{6})$, we obtain

$\displaystyle\mathcal{A}_{\mathrm{Re}(b_{6})}(E_{\boldsymbol{p}},E_{\boldsymbol{p}^{\prime}})$
$\displaystyle=8m_{N}^{2}m_{X}^{2}\left(E_{\boldsymbol{p}}-E_{\boldsymbol{p}^{\prime}}\right)\left[E_{\boldsymbol{p}^{\prime}}\left(m_{N}E_{\boldsymbol{p}^{\prime}}+m_{N}^{2}-m_{X}^{2}\right)+m_{X}^{2}(E_{\boldsymbol{p}}-2m_{N})\right.$
$\displaystyle\left.+E_{\boldsymbol{p}}m_{N}(E_{\boldsymbol{p}}-m_{N})\right]\,,$
$\displaystyle\mathcal{B}_{\mathrm{Re}(b_{6})}(E_{\boldsymbol{p}},E_{\boldsymbol{p}^{\prime}})$
$\displaystyle=4m_{N}m_{X}^{2}\left(E_{\boldsymbol{p}}-E_{\boldsymbol{p}^{\prime}}\right)^{2}\left[E_{\boldsymbol{p}^{\prime}}\left(2E_{\boldsymbol{p}}m_{N}-m_{N}^{2}+m_{X}^{2}\right)\right.$
$\displaystyle\left.+E_{\boldsymbol{p}}(m_{N}-m_{X})(m_{N}+m_{X})+2m_{N}m_{X}^{2}\right]\,,$
$\displaystyle\mathcal{C}_{\mathrm{Re}(b_{6})}(E_{\boldsymbol{p}},E_{\boldsymbol{p}^{\prime}})$
$\displaystyle=16m_{N}^{2}m_{X}^{2}\left(E_{\boldsymbol{p}}-E_{\boldsymbol{p}^{\prime}}\right)^{2}\left(-m_{N}E_{\boldsymbol{p}^{\prime}}+E_{\boldsymbol{p}}m_{N}+2m_{X}^{2}\right)\,.$
(A.5)

  * •

For $\eta=\mathrm{Im}(b_{6})$, the model dependent functions in Eq. (A.2) can
explicitly be written as follows

$\displaystyle\mathcal{A}_{\mathrm{Im}(b_{6})}(E_{\boldsymbol{p}},E_{\boldsymbol{p}^{\prime}})$
$\displaystyle=8m_{N}^{2}\left(E_{\boldsymbol{p}}-E_{\boldsymbol{p}^{\prime}}\right)\left\\{E_{\boldsymbol{p}^{\prime}}\left[2E_{\boldsymbol{p}}^{2}m_{N}^{2}-2E_{\boldsymbol{p}}m_{N}^{2}E_{\boldsymbol{p}^{\prime}}+m_{N}m_{X}^{2}(2E_{\boldsymbol{p}}+m_{N})\right.\right.$
$\displaystyle\left.\left.-m_{X}^{4}\right]-E_{\boldsymbol{p}}m_{N}^{2}m_{X}^{2}+m_{X}^{4}(E_{\boldsymbol{p}}-2m_{N})\right\\}\,,$
$\displaystyle\mathcal{B}_{\mathrm{Im}(b_{6})}(E_{\boldsymbol{p}},E_{\boldsymbol{p}^{\prime}})$
$\displaystyle=2m_{N}\left(E_{\boldsymbol{p}}-E_{\boldsymbol{p}^{\prime}}\right)^{2}\left\\{E_{\boldsymbol{p}}^{2}m_{N}^{2}(E_{\boldsymbol{p}}-2m_{N})\right.$
$\displaystyle\left.+E_{\boldsymbol{p}^{\prime}}\left[-m_{N}^{2}E_{\boldsymbol{p}^{\prime}}\left(E_{\boldsymbol{p}^{\prime}}+E_{\boldsymbol{p}}+2m_{N}\right)+E_{\boldsymbol{p}}m_{N}^{2}(E_{\boldsymbol{p}}+4m_{N})\right.\right.$
$\displaystyle\left.\left.+2m_{N}m_{X}^{2}(2E_{\boldsymbol{p}}+m_{N})+2m_{X}^{4}\right]-2E_{\boldsymbol{p}}m_{N}^{2}m_{X}^{2}-2m_{X}^{4}(E_{\boldsymbol{p}}-2m_{N})\right\\}\,,$
$\displaystyle\mathcal{C}_{\mathrm{Im}(b_{6})}(E_{\boldsymbol{p}},E_{\boldsymbol{p}^{\prime}})$
$\displaystyle=-8m_{N}^{2}\left(E_{\boldsymbol{p}}-E_{\boldsymbol{p}^{\prime}}\right)^{2}\left[E_{\boldsymbol{p}}^{2}m_{N}^{2}+m_{N}^{2}E_{\boldsymbol{p}^{\prime}}\left(E_{\boldsymbol{p}^{\prime}}-2E_{\boldsymbol{p}}\right)-4m_{X}^{4}\right]\,.$
(A.6)

  * •

For $\eta=\mathrm{Re}(b_{7})$, we find

$\displaystyle\mathcal{A}_{\mathrm{Re}(b_{7})}(E_{\boldsymbol{p}},E_{\boldsymbol{p}^{\prime}})$
$\displaystyle=8m_{N}m_{X}^{2}\left(-m_{N}E_{\boldsymbol{p}^{\prime}}+E_{\boldsymbol{p}}m_{N}+2m_{X}^{2}\right)\left[E_{\boldsymbol{p}^{\prime}}\left(m_{N}E_{\boldsymbol{p}^{\prime}}+m_{N}^{2}-m_{X}^{2}\right)\right.$
$\displaystyle\left.+m_{X}^{2}(E_{\boldsymbol{p}}-2m_{N})+E_{\boldsymbol{p}}m_{N}(E_{\boldsymbol{p}}-m_{N})\right]\,,$
$\displaystyle\mathcal{B}_{\mathrm{Re}(b_{7})}(E_{\boldsymbol{p}},E_{\boldsymbol{p}^{\prime}})$
$\displaystyle=4m_{X}^{2}\left(E_{\boldsymbol{p}}-E_{\boldsymbol{p}^{\prime}}\right)\left(-m_{N}E_{\boldsymbol{p}^{\prime}}+E_{\boldsymbol{p}}m_{N}+2m_{X}^{2}\right)$
$\displaystyle\times\left[E_{\boldsymbol{p}^{\prime}}\left(2E_{\boldsymbol{p}}m_{N}-m_{N}^{2}+m_{X}^{2}\right)+E_{\boldsymbol{p}}(m_{N}-m_{X})(m_{N}+m_{X})\right.$
$\displaystyle\left.+2m_{N}m_{X}^{2}\right]\,,$
$\displaystyle\mathcal{C}_{\mathrm{Re}(b_{7})}(E_{\boldsymbol{p}},E_{\boldsymbol{p}^{\prime}})$
$\displaystyle=16m_{N}m_{X}^{2}\left(E_{\boldsymbol{p}}-E_{\boldsymbol{p}^{\prime}}\right)\left(-m_{N}E_{\boldsymbol{p}^{\prime}}+E_{\boldsymbol{p}}m_{N}+2m_{X}^{2}\right)^{2}\,.$
(A.7)

  * •

Finally, for $\eta=\mathrm{Im}(b_{7})$, we obtain

$\displaystyle\mathcal{A}_{\mathrm{Im}(b_{7})}(E_{\boldsymbol{p}},E_{\boldsymbol{p}^{\prime}})$
$\displaystyle=8m_{N}^{2}m_{X}^{2}\left(E_{\boldsymbol{p}}-E_{\boldsymbol{p}^{\prime}}\right)\left[E_{\boldsymbol{p}}^{2}m_{N}+E_{\boldsymbol{p}^{\prime}}\left(m_{N}E_{\boldsymbol{p}^{\prime}}+m_{N}^{2}+m_{X}^{2}\right)\right.$
$\displaystyle\left.-E_{\boldsymbol{p}}\left(m_{N}^{2}+m_{X}^{2}\right)\right]\,,$
$\displaystyle\mathcal{B}_{\mathrm{Im}(b_{7})}(E_{\boldsymbol{p}},E_{\boldsymbol{p}^{\prime}})$
$\displaystyle=4m_{N}^{2}m_{X}^{2}\left(E_{\boldsymbol{p}}-E_{\boldsymbol{p}^{\prime}}\right)^{2}\left[(2E_{\boldsymbol{p}}-m_{N})E_{\boldsymbol{p}^{\prime}}+E_{\boldsymbol{p}}m_{N}-2m_{X}^{2}\right]\,,$
$\displaystyle\mathcal{C}_{\mathrm{Im}(b_{7})}(E_{\boldsymbol{p}},E_{\boldsymbol{p}^{\prime}})$
$\displaystyle=16m_{N}^{2}m_{X}^{2}\left(E_{\boldsymbol{p}}-E_{\boldsymbol{p}^{\prime}}\right)^{2}\left(-m_{N}E_{\boldsymbol{p}^{\prime}}+E_{\boldsymbol{p}}m_{N}-m_{X}^{2}\right)\,.$
(A.8)

##  Appendix B Cross sections for relativistic dark matter-electron scattering

In the laboratory frame, the differential cross section for DM-electron
scattering can be written as follows

$\displaystyle\frac{{\rm
d}\sigma_{eX}(E_{\boldsymbol{p}},E_{\boldsymbol{k}^{\prime}})}{{\rm
d}E_{\boldsymbol{k}^{\prime}}}=\frac{1}{32\pi
m_{e}(E_{\boldsymbol{p}}^{2}-m_{X}^{2})}\overline{|\mathcal{M}_{eX}(E_{\boldsymbol{p}},E_{\boldsymbol{k}^{\prime}})|^{2}}\,,$
(B.1)

where $E_{\boldsymbol{k}^{\prime}}$ is the final state electron energy. As for
the case of DM-nucleon scattering, we evaluate Eq. (B.1) by the combined use
of FeynRules [30], CalcHEP [31] and analytical calculations for validation.
Here, we express the squared modulus of the scattering amplitude in Eq. (B.1)
as

$\displaystyle\overline{|\mathcal{M}_{eX}(E_{\boldsymbol{p}},E_{\boldsymbol{k}^{\prime}})|^{2}}=\frac{\eta^{2}h_{3}^{2}\,\mathcal{D}_{\eta}(E_{\boldsymbol{p}},E_{\boldsymbol{k}^{\prime}})}{3m_{X}^{4}\left[2m_{e}(E_{\boldsymbol{k}^{\prime}}-m_{e})+m_{A^{\prime}}^{2}\right]^{2}}\,.$
(B.2)

Below, we specify the model dependent function
$\mathcal{D}_{\eta}(E_{\boldsymbol{p}},E_{\boldsymbol{k}^{\prime}})$ for
different choices of coupling constant $\eta$:

$\displaystyle\mathcal{D}_{b_{5}}(E_{\boldsymbol{p}},E_{\boldsymbol{k}^{\prime}})$
$\displaystyle=8m_{e}\left[m_{e}^{2}(E_{\boldsymbol{k}^{\prime}}-m_{e})^{2}+2m_{e}m_{X}^{2}(E_{\boldsymbol{k}^{\prime}}-m_{e})+3m_{X}^{4}\right]$
$\displaystyle\times\left[2E_{\boldsymbol{p}}m_{e}(-E_{\boldsymbol{k}^{\prime}}+E_{\boldsymbol{p}}+m_{e})+m_{X}^{2}(m_{e}-E_{\boldsymbol{k}^{\prime}})\right]\,,$
$\displaystyle\mathcal{D}_{\mathrm{Re}(b_{6})}(E_{\boldsymbol{p}},E_{\boldsymbol{k}^{\prime}})$
$\displaystyle=8m_{e}^{2}m_{X}^{2}(E_{\boldsymbol{k}^{\prime}}-m_{e})\left\\{m_{e}\left[E_{\boldsymbol{k}^{\prime}}^{2}-E_{\boldsymbol{k}^{\prime}}(2E_{\boldsymbol{p}}+3m_{e})+2\left(E_{\boldsymbol{p}}^{2}+E_{\boldsymbol{p}}m_{e}+m_{e}^{2}\right)\right]\right.$
$\displaystyle\left.+m_{X}^{2}(E_{\boldsymbol{k}^{\prime}}-3m_{e})\right\\}\,,$
$\displaystyle\mathcal{D}_{\mathrm{Im}(b_{6})}(E_{\boldsymbol{p}},E_{\boldsymbol{k}^{\prime}})$
$\displaystyle=-8m_{e}^{2}(E_{\boldsymbol{k}^{\prime}}-m_{e})\left\\{-m_{e}m_{X}^{2}\left[-E_{\boldsymbol{k}^{\prime}}(2E_{\boldsymbol{p}}+m_{e})+2E_{\boldsymbol{p}}^{2}+2E_{\boldsymbol{p}}m_{e}+m_{e}^{2}\right]\right.$
$\displaystyle\left.+2E_{\boldsymbol{p}}m_{e}^{2}(m_{e}-E_{\boldsymbol{k}^{\prime}})(-E_{\boldsymbol{k}^{\prime}}+E_{\boldsymbol{p}}+m_{e})-m_{X}^{4}(E_{\boldsymbol{k}^{\prime}}-3m_{e})\right\\}\,,$
$\displaystyle\mathcal{D}_{\mathrm{Re}(b_{7})}(E_{\boldsymbol{p}},E_{\boldsymbol{k}^{\prime}})$
$\displaystyle=-8m_{e}m_{X}^{2}\left(-E_{\boldsymbol{k}^{\prime}}m_{e}+m_{e}^{2}-2m_{X}^{2}\right)\left\\{m_{e}\left[E_{\boldsymbol{k}^{\prime}}^{2}-E_{\boldsymbol{k}^{\prime}}(2E_{\boldsymbol{p}}+3m_{e})\right.\right.$
$\displaystyle\left.\left.+2\left(E_{\boldsymbol{p}}^{2}+E_{\boldsymbol{p}}m_{e}+m_{e}^{2}\right)\right]+m_{X}^{2}(E_{\boldsymbol{k}^{\prime}}-3m_{e})\right\\}\,,$
$\displaystyle\mathcal{D}_{\mathrm{Im}(b_{7})}(E_{\boldsymbol{p}},E_{\boldsymbol{k}^{\prime}})$
$\displaystyle=8m_{e}^{2}m_{X}^{2}(E_{\boldsymbol{k}^{\prime}}-m_{e})\left\\{m_{e}\left[E_{\boldsymbol{k}^{\prime}}^{2}-E_{\boldsymbol{k}^{\prime}}(2E_{\boldsymbol{p}}+3m_{e})+2\left(E_{\boldsymbol{p}}^{2}+E_{\boldsymbol{p}}m_{e}+m_{e}^{2}\right)\right]\right.$
$\displaystyle\left.+m_{X}^{2}(m_{e}-E_{\boldsymbol{k}^{\prime}})\right\\}\,.$
(B.3)

##  Appendix C Scattering cross sections for SIMP DM

For the DM-electron scattering cross section for the non-abelian SIMP model,
we find the expression

$\displaystyle\frac{{\rm d}\sigma}{{\rm d}E_{e}}$
$\displaystyle=\frac{1}{192\pi
m_{X}^{4}}\frac{\sin^{2}(2\theta^{\prime}_{X})e^{2}\varepsilon^{2}g_{X}^{2}}{E_{\boldsymbol{p}}^{2}-m_{X}^{2}}\Big{\\{}\mathscr{A}(E_{e})\left[E_{\boldsymbol{p}}^{2}-(E_{e}-m_{e})E_{\boldsymbol{p}}\right]-\mathscr{B}(E_{e})\Big{\\}}\mathscr{C}(E_{e})\,,$
(C.1)

where we collected terms depending on different powers of
$E_{\boldsymbol{p}}$, and introduced the two coefficients

$\displaystyle\mathscr{A}(E_{e})$
$\displaystyle=2(E_{e}-m_{e})^{2}m_{e}^{3}+10(E_{e}-m_{e})m_{e}^{2}m_{X}^{2}+24m_{e}m_{X}^{4}\,,$
(C.2)

and

$\displaystyle\mathscr{B}(E_{e})$
$\displaystyle=(E_{e}-m_{e})m_{X}^{2}\left[(E_{e}-m_{e})m_{e}^{3}+(3E_{e}-m_{e})m_{e}m_{X}^{2}+12m_{X}^{4}\right]\,.$
(C.3)

The overall factor $\mathscr{C}(E_{e})$ arises from the propagators due to
$\tilde{Z}^{\prime}$ and $\tilde{X}_{3}$ exchange, and it is given by

$\displaystyle\mathscr{C}(E_{e})=\left(\frac{1}{2E_{e}m_{e}-2m_{e}^{2}+m_{\tilde{Z}^{\prime}}^{2}}-\frac{1}{2E_{e}m_{e}-2m_{e}^{2}+m_{\tilde{X}_{3}}^{2}}\right)^{2}\,.$
(C.4)

In the non-relativistic limit, Eq. (C.1) reduces to Eq. (66) from [24] if
integrated from 0 to $2\mu^{2}v^{2}/m_{e}$, where $\mu$ is the DM-electron
reduced mass, while $v$ is the DM-electron relative velocity.

For the DM-nucleon scattering cross section for the non-abelian SIMP model, we
find the expression

$\displaystyle\frac{{\rm d}\sigma}{{\rm d}E_{\boldsymbol{p}^{\prime}}}$
$\displaystyle=\frac{1}{192\pi
m_{X}^{4}}\frac{\sin^{2}(2\theta^{\prime}_{X})e^{2}\varepsilon^{2}g_{X}^{2}}{E_{\boldsymbol{p}}^{2}-m_{X}^{2}}\left[F_{1}^{2}\mathcal{A}_{1}(E_{\boldsymbol{p}},E_{\boldsymbol{p}^{\prime}})+F_{2}^{2}\mathcal{A}_{2}(E_{\boldsymbol{p}},E_{\boldsymbol{p}^{\prime}})+2F_{1}F_{2}\mathcal{A}_{12}(E_{\boldsymbol{p}},E_{\boldsymbol{p}^{\prime}})\right]$
$\displaystyle\times\frac{\mathscr{D}(E_{\boldsymbol{p}},E_{\boldsymbol{p}^{\prime}})}{4}\,,$
(C.5)

where

$\displaystyle\mathcal{A}_{1}(E_{\boldsymbol{p}},E_{\boldsymbol{p}^{\prime}})$
$\displaystyle=4\Big{\\{}2E_{\boldsymbol{p}}m_{N}^{3}E_{\boldsymbol{p}^{\prime}}^{3}+2E_{\boldsymbol{p}^{\prime}}\Big{[}E_{\boldsymbol{p}}^{3}m_{N}^{3}+E_{\boldsymbol{p}}m_{N}^{2}m_{X}^{2}(5E_{\boldsymbol{p}}+m_{N})$
$\displaystyle+m_{N}m_{X}^{4}(15E_{\boldsymbol{p}}+m_{N})+6m_{X}^{6}\Big{]}$
$\displaystyle-
m_{N}E_{\boldsymbol{p}^{\prime}}^{2}\Big{[}4E_{\boldsymbol{p}}^{2}m_{N}^{2}+m_{N}m_{X}^{2}(10E_{\boldsymbol{p}}+m_{N})+3m_{X}^{4}\Big{]}$
$\displaystyle-
E_{\boldsymbol{p}}m_{X}^{2}\Big{[}E_{\boldsymbol{p}}m_{N}^{3}+m_{N}m_{X}^{2}(3E_{\boldsymbol{p}}+2m_{N})+12m_{X}^{4}\Big{]}\Big{\\}}\,,$
(C.6)

$\displaystyle\mathcal{A}_{2}(E_{\boldsymbol{p}},E_{\boldsymbol{p}^{\prime}})$
$\displaystyle=\left(E_{\boldsymbol{p}}-E_{\boldsymbol{p}^{\prime}}\right)\Big{\\{}E_{\boldsymbol{p}}^{3}m_{N}^{2}(E_{\boldsymbol{p}}-2m_{N})+2E_{\boldsymbol{p}}^{2}m_{N}m_{X}^{2}(2E_{\boldsymbol{p}}-5m_{N})$
$\displaystyle+E_{\boldsymbol{p}^{\prime}}\Big{[}6E_{\boldsymbol{p}}^{2}m_{N}^{3}+E_{\boldsymbol{p}^{\prime}}\Big{(}m_{N}E_{\boldsymbol{p}^{\prime}}(m_{N}E_{\boldsymbol{p}^{\prime}}+2m_{N}^{2}-4m_{X}^{2})-2E_{\boldsymbol{p}}m_{N}^{2}(E_{\boldsymbol{p}}+3m_{N})$
$\displaystyle-2m_{N}m_{X}^{2}(4E_{\boldsymbol{p}}+5m_{N})+10m_{X}^{4}\Big{)}+4m_{X}^{4}(7E_{\boldsymbol{p}}+9m_{N})$
$\displaystyle+4E_{\boldsymbol{p}}m_{N}m_{X}^{2}(2E_{\boldsymbol{p}}+5m_{N})\Big{]}+2E_{\boldsymbol{p}}m_{X}^{4}(5E_{\boldsymbol{p}}-18m_{N})-48m_{X}^{6}\Big{\\}}\,,$
(C.7)

$\displaystyle\mathcal{A}_{12}(E_{\boldsymbol{p}},E_{\boldsymbol{p}^{\prime}})$
$\displaystyle=-2\left(E_{\boldsymbol{p}}-E_{\boldsymbol{p}^{\prime}}\right)\left(-m_{N}E_{\boldsymbol{p}^{\prime}}+E_{\boldsymbol{p}}m_{N}+2m_{X}^{2}\right)\Big{\\{}E_{\boldsymbol{p}}^{2}m_{N}^{2}+m_{N}E_{\boldsymbol{p}^{\prime}}\big{[}m_{N}E_{\boldsymbol{p}^{\prime}}$
$\displaystyle-2\left(E_{\boldsymbol{p}}m_{N}+m_{X}^{2}\right)\Big{]}+2E_{\boldsymbol{p}}m_{N}m_{X}^{2}+12m_{X}^{4}\Big{\\}}\,,$
(C.8)

and

$\displaystyle\mathscr{D}(E_{\boldsymbol{p}},E_{\boldsymbol{p}^{\prime}})=\left(\frac{1}{m_{\tilde{Z}^{\prime}}^{2}+2m_{N}(E_{\boldsymbol{p}}-E_{\boldsymbol{p}^{\prime}})}-\frac{1}{m_{\tilde{X}_{3}}+2m_{N}(E_{\boldsymbol{p}}-E_{\boldsymbol{p}^{\prime}})}\right)^{2}\,.$
(C.9)

## References

  * [1] J. Cooley et al., _Report of the Topical Group on Particle Dark Matter for Snowmass 2021_ , 2209.07426.
  * [2] M. Battaglieri et al., _US Cosmic Visions: New Ideas in Dark Matter 2017; Community Report_ , _FERMILAB-CONF-17-282-AE-PPD-T_ (2017) [1707.04591].
  * [3] A. Mitridate, T. Trickle, Z. Zhang and K. M. Zurek, _Snowmass white paper: Light dark matter direct detection at the interface with condensed matter physics_ , _Phys. Dark Univ._ 40 (2023) 101221 [2203.07492].
  * [4] R. Essig, J. Mardon and T. Volansky, _Direct Detection of Sub-GeV Dark Matter_ , _Phys. Rev._ D85 (2012) 076007 [1108.5383].
  * [5] Planck collaboration, N. Aghanim et al., _Planck 2018 results. VI. Cosmological parameters_ , _Astron. Astrophys._ 641 (2020) A6 [1807.06209].
  * [6] B. W. Lee and S. Weinberg, _Cosmological Lower Bound on Heavy Neutrino Masses_ , _Phys. Rev. Lett._ 39 (1977) 165.
  * [7] LSND collaboration, L. B. Auerbach et al., _Measurement of electron - neutrino - electron elastic scattering_ , _Phys. Rev. D_ 63 (2001) 112001 [hep-ex/0101039].
  * [8] MiniBooNE DM collaboration, A. A. Aguilar-Arevalo et al., _Dark Matter Search in Nucleon, Pion, and Electron Channels from a Proton Beam Dump with MiniBooNE_ , _Phys. Rev. D_ 98 (2018) 112004 [1807.06137].
  * [9] LDMX collaboration, T. Åkesson et al., _Light Dark Matter eXperiment (LDMX)_ , 1808.05219.
  * [10] A. Berlin, N. Blinov, G. Krnjaic, P. Schuster and N. Toro, _Dark matter, millicharges, axion and scalar particles, gauge bosons, and other new physics with LDMX_ , _Physical Review D_ 99 (2019) .
  * [11] P. deNiverville and C. Frugiuele, _Hunting sub-gev dark matter with the $\mathrm{NO}\nu\mathrm{A}$ near detector_, _Phys. Rev. D_ 99 (2019) 051701.
  * [12] P. deNiverville, C.-Y. Chen, M. Pospelov and A. Ritz, _Light dark matter in neutrino beams: Production modeling and scattering signatures at MiniBooNE, t2k, and SHiP_ , _Physical Review D_ 95 (2017) .
  * [13] P. deNiverville, C.-Y. Chen, M. Pospelov and A. Ritz, _Light dark matter in neutrino beams: production modelling and scattering signatures at MiniBooNE, T2K and SHiP_ , _Phys. Rev. D_ 95 (2017) 035006 [1609.01770].
  * [14] B. Batell, R. Essig and Z. Surujon, _Strong constraints on sub-gev dark sectors from slac beam dump e137_ , _Phys. Rev. Lett._ 113 (2014) 171802.
  * [15] DELPHES 3 collaboration, J. de Favereau, C. Delaere, P. Demin, A. Giammanco, V. Lemaitre, A. Mertens et al., _DELPHES 3, A modular framework for fast simulation of a generic collider experiment_ , _JHEP_ 02 (2014) 057 [1307.6346].
  * [16] P. deNiverville, M. Pospelov and A. Ritz, _Observing a light dark matter beam with neutrino experiments_ , _Phys. Rev. D_ 84 (2011) 075020 [1107.4580].
  * [17] X. Chu, J. Hisano, A. Ibarra, J.-L. Kuo and J. Pradler, _Multipole vector dark matter below the GeV-scale_ , 2303.13643.
  * [18] P. Gondolo, I. Jeong, S. Kang, S. Scopel and G. Tomar, _Phenomenology of nuclear scattering for a WIMP of arbitrary spin_ , _Phys. Rev. D_ 104 (2021) 063018 [2102.09778].
  * [19] P. Gondolo, S. Kang, S. Scopel and G. Tomar, _Effective theory of nuclear scattering for a WIMP of arbitrary spin_ , _Phys. Rev. D_ 104 (2021) 063017 [2008.05120].
  * [20] R. Catena, K. Fridell and M. B. Krauss, _Non-relativistic Eff‘ective Interactions of Spin 1 Dark Matter_ , _JHEP_ 08 (2019) 030 [1907.02910].
  * [21] R. Catena, K. Fridell and V. Zema, _Direct detection of fermionic and vector dark matter with polarised targets_ , _JCAP_ 1811 (2018) 018 [1810.01515].
  * [22] R. Catena, D. Cole, T. Emken, M. Matas, N. Spaldin, W. Tarantino et al., _Dark matter-electron interactions in materials beyond the dark photon model_ , _JCAP_ 03 (2023) 052 [2210.07305].
  * [23] J. B. Dent, L. M. Krauss, J. L. Newstead and S. Sabharwal, _General analysis of direct dark matter detection: From microphysics to observational signatures_ , _Phys. Rev._ D92 (2015) 063515 [1505.03117].
  * [24] S.-M. Choi, H. M. Lee, Y. Mambrini and M. Pierre, _Vector SIMP dark matter with approximate custodial symmetry_ , _JHEP_ 07 (2019) 049 [1904.04109].
  * [25] T. Åkesson, A. Berlin, N. Blinov, O. Colegrove, G. Collura, V. Dutta et al., _Light dark matter experiment (ldmx)_ , 2018.  10.48550/ARXIV.1808.05219.
  * [26] R. Catena, J. Conrad and M. B. Krauss, _Compatibility of a dark matter discovery at XENONnT or LZ with the WIMP thermal production mechanism_ , _Physical Review D_ 97 (2018) .
  * [27] G. Arcadi, M. Dutra, P. Ghosh, M. Lindner, Y. Mambrini, M. Pierre et al., _The waning of the WIMP? A review of models, searches, and constraints_ , _Eur. Phys. J._ C78 (2018) 203 [1703.07364].
  * [28] A. Berlin, N. Blinov, G. Krnjaic, P. Schuster and N. Toro, _Dark Matter, Millicharges, Axion and Scalar Particles, Gauge Bosons, and Other New Physics with LDMX_ , 1807.01730.
  * [29] A. Berlin, P. deNiverville, A. Ritz, P. Schuster and N. Toro, _Sub-gev dark matter production at fixed-target experiments_ , _Phys. Rev. D_ 102 (2020) 095011.
  * [30] A. Alloul, N. D. Christensen, C. Degrande, C. Duhr and B. Fuks, _FeynRules 2.0 - A complete toolbox for tree-level phenomenology_ , _Comput. Phys. Commun._ 185 (2014) 2250 [1310.1921].
  * [31] A. Belyaev, N. D. Christensen and A. Pukhov, _CalcHEP 3.4 for collider physics within and beyond the Standard Model_ , _Comput. Phys. Commun._ 184 (2013) 1729 [1207.6082].
  * [32] G. Belanger, F. Boudjema, A. Pukhov and A. Semenov, _MicrOMEGAs 2.0: A Program to calculate the relic density of dark matter in a generic model_ , _Comput. Phys. Commun._ 176 (2007) 367 [hep-ph/0607059].
  * [33] P. Gondolo and G. Gelmini, _Cosmic abundances of stable particles: Improved analysis_ , _Nucl. Phys._ B360 (1991) 145.
  * [34] E. Izaguirre, G. Krnjaic and B. Shuve, _Discovering inelastic thermal relic dark matter at colliders_ , _Physical Review D_ 93 (2016) .
  * [35] Z. et al. (Particle Data Group), _Review of Particle Physics_ , _Progress of Theoretical and Experimental Physics_ 2020 (2020)  [https://academic.oup.com/ptep/article-pdf/2020/8/083C01/34673722/ptaa104.pdf].
  * [36] P. Ilten, Y. Soreq, M. Williams and W. Xue, _Serendipity in dark photon searches_ , _Journal of High Energy Physics_ 2018 (2018) .
  * [37] R. Catena, D. Cole, T. Emken, M. Matas, N. Spaldin, W. Tarantino et al., _Dark matter-electron interactions in materials beyond the dark photon model_ , _Journal of Cosmology and Astroparticle Physics_ 2023 (2023) 052.
  * [38] T. R. Slatyer, _Indirect dark matter signatures in the cosmic dark ages. i. generalizing the bound on $s$-wave dark matter annihilation from planck results_, _Phys. Rev. D_ 93 (2016) 023527.
  * [39] M. Dutra, M. Lindner, S. Profumo, F. S. Queiroz, W. Rodejohann and C. Siqueira, _MeV dark matter complementarity and the dark photon portal_ , _Journal of Cosmology and Astroparticle Physics_ 2018 (2018) 037.
  * [40] H. Liu, T. R. Slatyer and J. Zavala, _Contributions to cosmic reionization from dark matter annihilation and decay_ , _Physical Review D_ 94 (2016) .
  * [41] H. Liu, W. Qin, G. W. Ridgway and T. R. Slatyer, _Lyman- $\alpha$ constraints on cosmic heating from dark matter annihilation and decay_, _Phys. Rev. D_ 104 (2021) 043514.
  * [42] D. Banerjee et al., _Dark matter search in missing energy events with NA64_ , _Phys. Rev. Lett._ 123 (2019) 121801 [1906.00176].
  * [43] P. Schuster, N. Toro and K. Zhou, _Probing invisible vector meson decays with the NA64 and LDMX experiments_ , _Phys. Rev. D_ 105 (2022) 035036 [2112.02104].
  * [44] BaBar collaboration, J. P. Lees et al., _Search for Invisible Decays of a Dark Photon Produced in ${e}^{+}{e}^{-}$ Collisions at BaBar_, _Phys. Rev. Lett._ 119 (2017) 131804 [1702.03327].
  * [45] E. Kou et al., _The Belle II Physics Book_ , _Progress of Theoretical and Experimental Physics_ 2019 (2019)  [https://academic.oup.com/ptep/article-pdf/2019/12/123C01/32693980/ptz106.pdf].
  * [46] F. Kahlhoefer, K. Schmidt-Hoberg, T. Schwetz and S. Vogl, _Implications of unitarity and gauge invariance for simplified dark matter models_ , _JHEP_ 02 (2016) 016 [1510.02110].
  * [47] C. Chang, P. Scott, T. E. Gonzalo, F. Kahlhoefer and M. White, _Global fits of simplified models for dark matter with GAMBIT II. Vector dark matter with an $s$-channel vector mediator_, 2303.08351.