Spin-1 thermal targets for dark matter searches at beam dump and fixed target experiments

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.


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.

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 → 2 processes and 2 → 2 forbidden annihilations, whereas in the simplified models we consider here the relic abundance is set by pair annihilation into visible particles.

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 µ , playing the role of DM, and a mediator particle described by the real vector field A ′µ .The following Lagrangian specifies the interactions between the complex and real vector fields X µ and A ′µ , respectively, and the Dirac spinors f associated with the electrically charged SM fermions, 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 ′µ .In order to make the analogy with the dark photon model explicit, one could identify h 3 with h 3 = eε, where, ε is the so-called kinetic mixing parameter, which enters the Lagrangian of the dark photon model via the term −(ε/2)F µν F ′µν , F µν and F ′µν 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: or Im[b 7 ].As far as the DM and mediator particle mass are concerned, we denote them by m X and m A ′ , 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 ′ ≥ 2m X .Under this mass hierarchy, and since g D ≫ eε, 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 ′ ≤ 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].

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) × U (1) Z ′ symmetry group under which none of the SM particles is charged.The corresponding gauge couplings are g X and g Z ′ .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) × U (1) Z ′ .The gauge bosons associated with the new symmetry group are denoted by X i,µ , i = 1, 2, 3 and Z ′ µ , or, equivalently, by ,µ and Z ′ µ , respectively.The SU (2) × U (1) Z ′ 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 µ and two linear combinations of Z ′ µ and X 3,µ , denoted here by Z′ µ and X3,µ , respectively.Their masses are [24], where with Here q S is the U Z ′ (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], 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 Z′ µ , X3µ and the charged SM fermions arise from a kinetic mixing term added to the Lagrangian of this SU X (2) × U (1) Z ′ gauge model.Specifically, they are given by [24] For example, for sin(2θ ′ X ) = −0.1 and α D ≡ g 2 X /(4π) = 0.5, Eq. (2.6) is always satisfied for perturbative values of g Z ′ .Consequently, when the mass hierarchy in Eq. (2.2) is realised, X µ is a stable DM candidate, while X3 and Z′ mediate the interactions between X µ and the SM fermions via Eq.(2.5).Finally, let us note that for cos(θ ′ X ) = 1, Eq. (2.4) reduces to the first line in Eq. (2.1) with b 5 = −2Im[b 6 ] = g X , if one integrates by part the second line in Eq. (2.4) and uses the equation of motion for X µ (i.e.Proca equation), which implies ∂ µ X µ = 0.

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.

Relic density
Let us denote the DM particle associated with X µ 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], 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 socalled 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 ε ∼ 10 −4 (10 −6 ) and m X ∼ 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 → 2 processes and 2 → 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 X3 mediator mass ratio determines whether the DM relic abundance is set by 3 → 2 processes or 2 → 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.

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 X3 is in kinetic equilibrium with the SM thermal bath during the DM freeze-out.Interestingly, the in-equilibrium decay of X3 into SM particles, combined with the effective scattering of DM with the X3 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θ ′ X ), m X and m X3 , the value of ε required for X3 and, consequently, for the DM particles to be in kinetic equilibrium at the freeze-out temperature T f can be estimated from where H(T f ) is the Hubble rate at T f , and the equilibrium densities for X and X3 are given by 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 ≡ m X /T .In Eq. (2.8), Γ X3 is the total decay rate of X3 , for which the expression is given in [24].

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 ′ particles are produced through the kinetic mixing term in Eq. (2.1), either via dark bremsstrahlung or via meson decays.These A ′ s are then expected to decay into a DM particle/anti-particle pair for values of ε 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 X3 and the Z′ mediator can in principle be produced in fixed target and beam dump experiments.In our calculations we will focus on sin(2θ ′ X ) = −0.1, which implies that only Z′ particles can significantly be produced via the interaction in Eq. (2.5), and that, once produced, these Z′ particles will dominantly decay invisibly via gauge bosons self-interactions, as long as ε is small compared to g X .Consequently, for sin(2θ ′ 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.

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 ′ mediators (in the case of simplified models) or Z′ 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, in the case of simplified models.For the case of SIMP DM, we use Eq.(2.10) but with and A ′ → Z′ .Here, F 1 and F 2 are nuclear form factors, and we list them in appendix A.

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.

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 → 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 ′ > 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 ≈ m µ , where DM annihilation into muons become kinematically accessible.The resonance features appear because of resonances in the cross sections for DM annihilation into hadrons.
In the case of vector SIMP DM, we rely on previus results from [24].For each given α D = g 2 X /(4π) (g X is one of the gauge couplings introduced in Sec.(2.2)) and sin(2θ   40) from [24].
extract the value of the m X3 /m X ratio that gives the correct DM relic density from Fig. 3 in [24].While for α D = 0.1 the correct relic abundance can be obtained for DM masses above approximately 10 MeV, for α 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 X3 /m X ratio, focusing on the case in which 3 → 2 are the dominant production mechanism.As one can see from this figure, for m X3 /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 X3 /m X > 2, the 3 → 2 process is kinematically not allowed.

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 Re[b 7 ] models.We take the 90% C.L. exclusion limits from the work of [37].

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 Im[b 7 ] ̸ = 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-α observations on our pwave 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 ′ /m X ) and models.

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 [7,16] and MiniBooNE [8], a proton beam is incident on a target and the following chain of interactions occur producing DM: pp → Xπ 0 ; π 0 → γA ′ ; A ′ → DM DM , where in the inclusive process, pp → Xπ 0 , X denotes an unspecified/unmeasured set of particles, whereas π 0 can also be η 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 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.

Monophoton searches
We also include searches for single photon events at e + e − colliders, where DM is produced through the process e + e − → γA ′ , A ′ → DM DM .

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 ε 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 ε vs m X for phase 3 of Belle-II reported in [45], to obtain the expected exclusion limits on our spin-1 DM model.

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, where i → i = Xe − → Xe − , J = 0, and

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.

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 = ε 2 α D (m X /m A ′ ) 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 ̸ = 0 for both m A ′ /m X = 3 and m A ′ /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 ′ /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 ̸ = 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 p /m X ) 2 ≫ 1 relative to the case of the complex scalar DM model, where E p is the relativistic energy of the incoming DM particle in the rest frame of the downstream  cancellations between contributions from the m Z′ and m X′ 3 mediators.In Fig. 6, we assume α D = 0.5, m A ′ /m X = 3, and m X3 /m X ≃ 2. Fig. 7 reports the results of an analogous analysis where we assume a different combination of parameters, namely α D = 0.1, m A ′ /m X = 3, and m X3 /m X ≃ 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 α D in the two figures.

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 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 → 2 processes and forbidden annihilations.It also exhibits a lower bound on the DM particle mass arising from the relic density constraint.The larger α 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. and (C.9)

and (h 3 ,
Im[b 7 ]) ̸ = 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 α D = g 2 D /(4π), where g D is one of b 5 , Re[b 6 ], Im[b 6 ], Re[b 7 ]

Figure 1 :
Figure1: 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 ′ = 3m X , and α D = g 2 D /(4π) = 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.

Figure 2 :
Figure 2: DM relic abundance, Ωh 2 , as a function of m X and m X3 /m X for SIMP DM thermally produced by the freeze-out of 3 → 2 processes.Here we take α D = 0.5.The horizontal dashed line corresponds to the observed value of Ωh 2 .Coloured lines have been obtained by evaluating Eq. (40) from[24].

Figure 4 : 4 . 1 m
Figure 4: Same as Fig. 3, now for the Re[b 6 ], Im[b 6 Re[b 7 ], and Im[b 7 ] models with α D = 0.5 and m A ′ = 3m χ .The models with Re[b 7 ] ̸ = 0 and Im[b 6 ] ̸ = 0 are strongly constrained by beam dump experiments, while the model with Im[b 7 ] ̸ = 0 is ruled out by CMB observations.The 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.