Light dark Higgs boson in minimal sub-GeV dark matter scenarios

Minimal scenarios with light (sub-GeV) dark matter whose relic density is obtained from thermal freeze-out must include new light mediators. In particular, a very well-motivated case is that of a new"dark"massive vector gauge boson mediator. The mass term for such mediator is most naturally obtained by a"dark Higgs mechanism"which leads to the presence of an often long-lived dark Higgs boson whose mass scale is the same as that of the mediator. We study the phenomenology and experimental constraints on two minimal, self-consistent dark sectors that include such a light dark Higgs boson. In one the dark matter is a pseudo-Dirac fermion, in the other a complex scalar. We find that the constraints from BBN and CMB are considerably relaxed in the framework of such minimal dark sectors. We present detection prospects for the dark Higgs boson in existing and projected proton beam-dump experiments. We show that future searches at experiments like Xenon1T or LDMX can probe all the relevant parameter space, complementing the various upcoming indirect constraints from astrophysical observations.


Introduction
Among the many puzzles facing the Standard Model (SM) of particle physics, the issue of dark matter (DM) is certainly one of the most pressing. While the prime candidate of the last decades has been the Weakly Interactive Massive Particle (WIMP, see, e.g. [1,2] for the latest reviews), direct, indirect and collider searches have so far failed to give an uncontroversial signal of such particles. Among the alternative ideas for dark matter that have emerged over the years, sub-GeV dark matter is gaining momentum, thanks both to a rich upcoming experimental program and to the fact that, similarly to the WIMP, it relies on the robust, UV-insensitive, thermal freeze-out mechanism to achieve the correct relic density (see [3] and [4] for reviews). These dark matter scenarios typically involve a dark matter candidate interacting with SM particles through a light mediator. In this article we shall focus on a specific class of models where the mediator is a new gauge boson, V , corresponding to a spontaneously broken new abelian gauge group U (1) D , because of their viability in providing a light thermal dark matter as well as because of many experimental searches devoted to such models. We will refer to this new gauge boson as the dark photon in the following.
Since the new U (1) D gauge group can mix with the Standard Model U (1) Y gauge group, the dark photon acts as a proper mediator between the dark and visible sectors. Such dark gauge groups have been particularly often used in a dark matter context due to their interesting properties and experimental prospects for detection (some very recent examples are, e.g. [5][6][7][8][9]). For instance, they can give rise to simple Self-Interacting Dark Matter models (SIDM, see [10] for the latest review) which could lead to better agreement between numerical simulations and the astrophysical observations.
The simplest way to generate the dark photon mass perturbatively is through a "dark Higgs mechanism". This assumes the presence of an additional dark Higgs boson which gives the dark photon its mass through a Vacuum Expectation Value (VEV) v S . Thus, a complete, self-consistent "dark sector" contains a dark matter candidate, the dark photon and the dark Higgs boson. Crucially, the dark Higgs boson mass is also proportional to v S , so that it should naturally be at the same mass scale as the dark photon and not significantly heavier.
Paradoxically, most of the literature on the field either focused on the dark Higgs boson, with or without the dark photon, or assumed that it decouples from the rest of the spectrum and concentrated on the dark matter and the dark photon only (one of the recent exceptions is [11] with a focus on the relic density constraint). In contrast, we present in this paper two minimal, self-consistent and perturbative models for the dark sector and systematically study the large part of the parameter space where the dark Higgs boson is light. In this case the dark photon, dark matter and dark Higgs boson must all be considered simultaneously.
As we will see below, the most important characteristic of a light dark Higgs boson is the fact that its lifetime is typically of order of one second or longer. Indeed, when the dark Higgs boson is lighter than the dark photon and of twice the dark matter mass, its decay is particularly suppressed as it can only proceed through a loop-induced coupling to light Standard Model particles. Such long-lived dark Higgs boson have been studied independently for several years and have been shown to possibly leave a signal in long baseline neutrino experiments and more generally in so-called "beamdump" experiments (see, e.g. [12][13][14][15]). Light dark Higgs boson originating for instance from the decay of a light meson can travel through the shielding of these beam-dump experiments and subsequently decay in the downstream detector. 1 We will re-evaluate this particular search strategy for detecting dark Higgs bosons and show that they are currently not sensitive enough to reach the thermal value target in our two minimal models.
The second main result is that the relic density calculation is thoroughly modified by the presence of new dark matter annihilation channels involving a dark Higgs boson. In particular, the long lifetime of the dark Higgs boson implies that the thermal freeze-out mechanism proceeds as in a two-component dark matter scenario. However, its presence also opens up new additional s-wave annihilation channels for dark matter at the time of recombination and leads therefore to severe bounds from CMB observations [19,20].
Finally, a long-lived dark Higgs boson is constrained by Big Bang Nucleosynthesis (BBN) related data [21,22], especially given that its metastable density obtained from thermal freeze-out can be larger than that of the dark matter. Nonetheless, we will show that light dark Higgs bosons in our two minimal dark sector models have metastable density substantially smaller than the Higgs portal case and can alleviate significantly the bounds presented in [21].
The paper is organized as follows. We first present in Sec. 2 two models of the dark sector framework, as well as existing constraints on the dark photon from various experiments. We then focus in Sec. 3 on the dark matter candidate and the effect of the presence of the dark Higgs boson on its relic density and on the constraints from CMB. Section 4 discusses detection prospect for the dark Higgs boson in beam-dump experiments as well as constraints related to BBN. Finally, in Sec. 5 we summarize our results and conclude. The appendix contains additional details on the calculation of dark Higgs boson production cross section from light mesons decay.

Minimal dark sector models and bounds on the dark photon
We present in this section two minimal, self-consistent dark sector models for a sub-GeV dark matter. As was discussed in the Introduction, such dark sectors typically include three types of fields: • an extra gauge boson (called "dark photon" in the following) V corresponding to "dark" gauge group U (1) D with a gauge constant g V ; • a complex scalar S with charge q S , called henceforth "dark Higgs boson". It spontaneously breaks the dark gauge group through a VEV, v S ; • a dark matter particle χ with charge q χ . We will consider both a complex scalar and a Majorana fermion dark matter candidate. As usual, we will assume that a discrete Z 2 symmetry protects the dark matter from decaying.

Lagrangian, masses and lifetimes
The gauge and matter content that we are considering implies that the dark sector can be coupled to the SM either through kinetic mixing between the two abelian gauge groups or by mixing between the SM Higgs H and the dark Higgs boson S. While both portals are a priori open, in this article we will focus on the vector portal. We will furthermore argue below that this is the most natural choice given the sub-GeV mass domain we are interested in. The kinetic mixing can in principle arise from loops of heavy fields charged under both gauge groups. We will assume in the following that they are safely decoupled at the energy scale that we consider. Given that the dark matter candidates must be charged under the new gauge group U (1) D , care must be taken when choosing them such that the dark gauge group remains anomaly-free. In particular, this excludes one single Majorana dark matter candidate, albeit a non-minimal scenario with a second heavier Majorana field canceling the anomaly is still possible. Consequently, we will consider in this paper two minimal, self-consistent, models for the dark matter candidates: • model pDF : the pseudo-Dirac fermion case, where a Dirac fermion χ = (χ L , χ † R ) dark matter acquires additional Majorana masses from its Yukawa interactions with the dark Higgs boson; • model CS : the complex scalar dark matter case, where we also denote the dark matter field by χ.
The simplest charge assignment in the pDF case is a U (1) D charge +2 for the dark Higgs boson S and ±1 for the two dark matter fermions χ L and χ R . In the CS case, we assign a charge +1 to the dark Higgs boson S and +1 to the complex scalar dark matter χ.
The effective Lagrangian for the dark photon vector and the dark Higgs boson fields in these two minimal dark sector models is then given by while the DM field is introduced either as a scalar or a fermion through the Lagrangian where V pDF and V CS describe the mixing of the DM particle with the dark Higgs boson S. We parametrize them as If µ 2 S > 0, and in the relevant limit where λ SH λ S , λ H , we can solve the tadpole equations for the VEVs of the SM Higgs v H and of the dark Higgs boson v S , leading to where µ 2 H and λ H are respectively the SM Higgs mass term and self quartic coupling. At zeroth order in v S /v H , the dark Higgs boson mass M S and dark photon mass M V are where we have introduced the dark Higgs boson U (1) D charge q S . In particular, the dark Higgs boson is lighter than the dark photon when This case will be of particular interest since the dark Higgs boson is then long-lived, as we will see in the next section. Notice that for typical SM-like values λ ∼ 0.1 and g V ∼ 0.5, the dark Higgs boson is indeed lighter than the dark photon. Furthermore, when the dark gauge coupling is chosen near its perturbativity bound with α D ≡ g 2 V /4π of order 0.5, then having a dark Higgs boson heavier than the dark photon leads to λ S > 1.25q 2 S and therefore possible non-perturbative behavior in the dark sector. For large values of α D , assuming the dark Higgs boson to be heavy enough to completely decouple from the rest of the dark sector is hence impossible in a minimal perturbative setup.
The kinetic mixing parameter should be small enough to avoid various experimental bounds discussed in the following sections. In a Grand Unified Theory context, the required small values for ε could be obtained from loops of heavy particles charged under both the SM hypercharge U (1) Y and the new U (1) D gauge group [23], with values between 10 −2 and 10 −5 depending on whether the mixing is generated at one or two-loops. 2 Finally, in the pDF case, the dark Higgs boson VEV leads to Majorana mass terms for the left-handed and right-handed components of χ. After diagonalizing the mass matrix, the lightest eigenstate χ 1 becomes our dark matter candidate. Notice that in principle y SL = y SR so that gauge coupling of schematic form χ 1 χ 1 V and χ 2 χ 2 V are a priori generated (albeit typically suppressed compared to χ 1 χ 2 V term).

Dark Higgs boson lifetime
When the tree-level decay of dark Higgs boson to dark matter is kinematically forbidden and its mixing with SM Higgs boson is negligible, the only decay mode available is through a triangular diagram of the form given in Fig. 1a. Furthermore, when M S < 2m µ the dominant decay mode is S → e + e − with the dark Higgs boson width given by [25] where the loop function is expressed as The above expressions also apply to the decay to muons by just replacing m e with m µ . In particular, in the limit m e M S , M V , we have The corresponding lifetime τ S is presented in Fig. 1b as a function of M S /M V , from which one can recover the exact value for any set of parameters using the previous scaling relations.
As an order of magnitude estimate, we then have where f are the kinematically accessible SM fermions, α em is the electromagnetic finestructure constant. In particular for M S below the dimuon mass threshold we find (2.12) In principle, the mixing between the Standard Model Higgs and the dark Higgs boson through the mixing quartic coupling λ SH could lead to additional decay channels. However, since the Higgs boson VEV contributes at tree level to the dark Higgs boson mass by λ SH v 2 we need for a dark Higgs boson mass between 10 and 100 MeV. If the dark Higgs boson could only decay to e + e − through its mixing with the SM Higgs, its lifetime τ S,Hmix would then be parametrically given by (2.14) This implies that, unless one is prepared to significantly tune the theory to ensure a light dark sector while keeping a large λ SH , the decay through SM Higgs mixing should be significantly smaller than the loop-induced one.
In the rest of this article, we will therefore neglect the Higgs-portal related effects (which includes the quartic λ SH , but also for simplicity, the dark matter/Higgs quartic λ χH in the CS case).

Constraints on the dark photon
If the dark matter is heavier than half of the dark photon mass, the dark photon decays mainly into a pair of leptons. This minimal scenario is mostly constrained by searches for bumps in the di-lepton invariant spectrum at NA-48/2 [26] and BaBar [27], setting bounds for ε 10 −3 . Slightly less competitive bounds also arise from rare meson  decays. One can also obtain bounds from electron beam-dump experiments like E137, E141 or E774, for small kinetic couplings leading to a long-lived dark photon decaying to visible sector. These searches give a lower bound on the kinetic mixing for a dark photon with mass in the tens of MeV range (see, e.g. [4] for a summary of the current bounds). The most relevant case for the parameter space considered here is when the dark photon decay channel to dark matter is kinematically open. The dark photon thus decays invisibly with the strongest bounds being set by searches at BaBar [28] and NA64 [29]. More precisely, the BaBar analyses searches for narrow peaks in the distribution of missing mass arising from e + e − → γV events. Their limit excludes the region ε > 10 −3 for the dark photon mass range we consider, which in particular rules out the dark photon explanation for the (g − 2) µ excess. Secondly, the NA64 Collaboration recently released bounds on the decay V → invisible. Their limits significantly exceed the one set by BaBar for M V 100 MeV, reaching ε < 10 −4 below 10 MeV. An explicit visualization of these bounds will be shown below in Fig. 4 in Sec. 3.2. Note that the projected bounds from the LDMX proposal (see, e.g. [3]) will cover almost all of the parameter space consistent with the relic density thermal value target as shown in Fig. 4.
In the following and for all the numerical results, we used the code MultiNest [30] to direct the scanning procedure, based on the dark matter relic density. All data points presented in this paper are therefore compatible with the result from the Planck Collaboration [31] Ωh 2 = 0.1188 ± 0.0010 at 95% CL. The interfaces with the various public codes used here is done with the help of the private code BayesFITS. We use a slightly modified version of MicrOMEGAs v.4.3.5 [32] (and of its two-component dark matter module). We evaluate the spectrum from the non-SUSY SPheno [33,34] code generated by SARAH (see Refs. [35][36][37]). We use renormalization group evolution of the hidden sector parameters to ensure their perturbativity up to the electroweak scale, and evaluate all masses at tree-level due to the light scale considered. Finally, the estimation of the number of events in beam-dump experiments is obtained from a substantially modified version of BdNMC from [18] (more particularly, we have used the original code to extract the distributions of initial mesons and expanded its routines to the production and detection processes relevant for the dark Higgs boson).
In the following, we will restrict our analysis to the case where the dark Higgs boson is below the dimuon threshold, so that it can only decay to an e + e − pair. The dark photon is also considered to be lighter than around 500 MeV, so that the leptonic decay channels still dominate its decay width compared to hadronic ones (see [38]). We summarize the independent parameters and their scanned ranges and priors in Table 1. Note that we do not vary the SM Higgs parameters. In particular, we take advantage of the relation (2.8) to trade v S for M V as an input parameter, so that we vary g V , ε, M S , M V , m χ , y SL and y SR in the pDF model and g V , ε, M S , M V , m χ , λ Sχ and λ χ in the CS model.

Light DM phenomenology
In this section we discuss the phenomenology of the light DM candidate in our two minimal dark sector models. We focus particularly on the relic density constraints and on the bounds from CMB power spectrum for s-wave annihilation processes occurring during the recombination era. We begin with a discussion of relic density for the pseudo-Dirac fermion (pDF case) and complex scalar (CS case) DM candidates. In the following, the dark matter mass is denoted by M χ , which hence refers to the mass of lightest mass eigenstate χ 1 in the pDF case.

Relic density
The relic density of DM in the standard freeze-out scenario is obtained by solving the following Boltzmann equation where n χ is the density of the DM species and σv is the thermally averaged annihilation rate of DM. The thermally averaged annihilation rate is given by where K 1 and K 2 are modified Bessel's functions. A useful parametrization of the annihilation rate is in terms of s-wave and p-wave annihilations like σv ≡ σ 0 x −n , with x = m χ /T . Here n = 0 for s-wave and n = 1 for p-wave annihilation. In this parametrization, x at freeze-out is given by [39] x f = ln 0.038(n + 1) where, following the notation of [39], we note that g represents the DM degrees of freedom, while g * and g * ,s represent the number of relativistic degrees of freedom at freeze-out. With the above expression for x f one can write the approximate expression for relic density as In the two minimal dark sector scenarios we consider, the dark matter particle can be either a pseudo-Dirac fermion (pDF case) or a complex scalar (CS case). In both cases, including the dark Higgs boson field leads to several new annihilation channels in a similar manner to the usual supersymmetric WIMP. The usual behavior considered by the previous literature corresponds to the case when the dark Higgs boson is significantly heavier than the dark matter candidate so that annihilation into dark Higgs boson is suppressed even with thermal effects included. The dominant process is a s-channel annihilation to SM particles through an off-shell dark photon with the annihilation cross-section, for instance in the CS case given by [40,41] (3.5) Using the above expression for σ 0 in Eq. (3.3) and Eq. (3.4) we arrive at the following estimate for the relic density (3.6) On the other hand, for M χ ∼ M S , dark matter annihilation into final states involving dark Higgs boson become relevant. They proceed either through a t-channel exchange of a dark matter particle, or through a dark Higgs boson s-channel. These new mechanisms alone could explain the current relic density for a dark sector coupling between dark Higgs boson and dark matter in the range we consider, and therefore have to be taken into account. A key complication of this setup is that the dark Higgs boson is a metastable particle with lifetime above 0.01 s in all of our parameter space. Consequently, thermal freeze-out proceeds akin to a two-component dark matter scenario. This is especially relevant when the mass of the dark matter and of the dark Higgs boson are of the same order, so that both χχ → SS and SS → χχ processes are occurring at the time of dark matter freeze-out. Furthermore, in the case of the pDF model, the fact that we consider the Yukawa couplings to the two Weyl components to be different in general (i.e, y SL = y SR in contrast with [42]) implies that the annihilation channels χ 1 χ 1 → e + e − and χ 1 χ 2 → SS are also available.
In Fig. 2 we represent the relevant annihilation channels that contribute to the relic density in the CS case and the pDF case. We see from the figure that in the  10 50 Figure 2: Points satisfying the dark matter relic density constraints in the M χ − M S plane, sorted according to the dominant annihilation channels at freeze-out in the pDF case (a) and the CS case (b). In (b), the region with M χ > M S , which is excluded by CMB bounds, has been indicated.
CS case, when M χ M S the SS channel dominates, while as M S approaches 2M χ the Sγ channel becomes dominant. When M S > 2M χ there are no S final states with only the e + e − channel remaining available to achieve the correct relic density, while the M S < M χ region is entirely excluded by CMB. For the pDF case the picture is more complicated due to the presence of coannihilation channels. 3 However, one still sees the e + e − channel being predominant in the M S > 2M χ region, and the SS channel dominating when M χ M S . In the regions where M χ < M S < 2M χ and M S < M χ , the Sγ channel is mostly responsible for achieving the correct relic density. In the very low mass region (M χ < 10 MeV), our choice of parameter range (and in particular m χ > 10 MeV) implies that most of our points have a very large splitting between both dark matter components χ 1 and χ 2 . 4 In particular, the channel χ 1 χ 2 → γV is sufficiently thermally-suppressed to lead to the correct relic density for some our points.
Thus the presence of dark Higgs bosons changes significantly the relic density evaluation in our two models. However, it also leads to two additional difficulties. First, the presence of a large metastable density of dark Higgs boson after thermal freeze-out may lead to strong tensions with BBN. We will explore this aspect in Sec. 4.2. Second, as we will see in the next section, the presence of the new annihilation channels, while significantly reducing the constraints arising from the relic density, may on the other hand be in strong tension with indirect bounds from CMB power spectrum.  Figure 3: Constraints from the CMB power spectrum on s-wave annihilation processes from [19,20] as a function of the dark matter mass M χ , with each point satisfying the dark matter relic density constraint for the pDF model (a) and for the CS case (b). The points have been sorted according to the dominant s-wave dark matter annihilation channels at recombination time. The red dashed line corresponds to projected indirect detection bounds from dark matter annihilating into γγ from [43]. It therefore only applies to the γγ channel, with bounds on annihilation into dark Higgs boson and its subsequent decay to e + e − pairs being weaker.

Direct and indirect detection bounds
The CMB power spectrum has been measured with high precision and as such can impose stringent constraints on the nature of DM. In particular DM that injects energy in the form of electromagnetically interacting SM particles in the inter-galactic medium (IGM) can significantly alter the recombination history of the universe by ionizing and heating the IGM gas. Such injections from DM annihilation can be parametrized as p ann = f σv /M χ , where f denotes the efficiency with which the energy injected by DM annihilations is transferred to the IGM. Usually the constraints from s-wave DM annihilations which do not depend on velocity of DM can be very stringent and virtually rule out most models with m χ < 10 GeV [19]. Since electrons and photons are the most efficient at ionizing the IGM, the annihilation channels that are most severely constrained produce e − s and γ-rays in their final states. For the pDF model, when λ Sχ L = λ Sχ R annihilation into an e + e − pair as χ 1 χ 1 → V * → e + e − becomes accessible. It is however safely suppressed by mixing matrices elements and the off-shell nature of the V in all our parameter space. When M χ < M V , M S , the only other relevant annihilation channel is χ 1 χ 1 → γγ, which is suppressed as ε 4 . If M χ > M S , then the annihilation χ 1 χ 1 → Sγ dominates. This channel is only suppressed as ε 2 and may lead to excluding some points at the boundary of our parameter space as shown in Fig. 3.
The situation is very different in the CS model, as t-channel annihilation into dark Higgs boson χχ → SS is completely unsuppressed when M S < M χ . Hence CMB bounds essentially rule out this portion of the parameter space as we present in Fig. 3.  Figure 4: Constraints on the dark photon mass M V and on the kinetic mixing parameter ε from NA64 and BaBaR missing energy searches, with the projected bounds from one year exposure of the annual modulation signals at Xenon1T according to [46] and from the full dataset of NA64 (corresponding to 10 11 electrons on target, see [3]). Points satisfying the dark matter relic density and relevant BBN and CMB constraints are shown for the pDF model (a) and for the CS case (b). The dashed red line represents the projected LDMX bound [3].
The only other channel which is not velocity-suppressed is then χχ → γγ, which is, as for the pDF case, suppressed as ε 4 . Notice that in both cases, if M χ > M V , other annihilation channels involving the dark photon open up which could lead to more severe bounds. However, they typically also significantly reduce the relic density, limiting the possibility to reach the thermal value target for the range of ε we consider in this paper.
The bounds from CMB quoted in Fig. 3 depend in principle on the annihilation products (in particular they have been calculated for the e + e − , γγ and e + e − via S decays). However, they do not differ significantly in the dark matter mass region we are interested in, so that we chose to use the bounds from χχ → γγ. 5 In the very low dark matter mass region of our plot (M χ 10 MeV), BBN-related bounds from energy injection from dark matter annihilation at freeze-out could become relevant [44]. However, they are model dependent and may in particular be modified due to the presence of a potentially long-lived dark Higgs.
Furthermore, indirect dark matter searches which focus on the MeV mass range for dark matter as considered in [43] also have the potential to significantly constrain our parameter space. This is also shown in Fig. 3 where we have represented the projected bounds from [43] for a one-year effective exposure from Galactic center observations of an experiment with the same specifications as e-ASTROGAM [45].
Finally, in the case of complex scalar dark matter CS, direct detection experiments searching for annual modulation signals through electron recoil are also relevant. Such searches were performed at XENON10 and XENON100 [46,47], leading to the following bounds [48] σ SI where the last inequality is the derived XENON10/XENON100 bound L(M χ ) which depends on the precise dark matter mass (see [46]). We present in Fig. 4 the corresponding bound as function of the dark matter mass for all points of our scans satisfying the relic density constraint. The projected bound from XENON1T can probe almost all of the parameter space where we found the correct relic density. Note that while we focused on xenon-based experiments, upcoming semi-conductor ones will also probe all of the parameter space of this model [49].

Light dark Higgs boson
We now turn to the second light state of our dark sector: the dark Higgs boson. As we have shown in Sec. 2.2, this particle is long-lived in most of our parameter space. We explore in this section two consequences of this long lifetime: the detection prospects at proton beam-dump experiments, and the constraints from BBN-related observables.

Beam dump experiments
Fixed target experiments are well suited for the detection of light dark sector particles. They typically involve a high-intensity, but relatively low-energy proton or electron beam impacting the target, producing a shower of secondary particles, which are later disposed off in a large shielding. Long-lived or stable dark matter particles are produced at a low rate in the target, but since they interact very weakly with the shielding, they travel to a downstream detector which can subsequently detect them. In particular, when the dark photon decays into dark matter particles, it effectively produces a "dark matter beam" and the possible scattering of dark matter in the detector can then be estimated [12,[16][17][18]. In particular, a case comparable to our fermion dark matter scenario pDF has been studied in [42]. In this section, we will focus on three proton beam-dump experiments: LSND [50], miniBooNE [51] and the proposed SBND experiment at Fermilab [52]. The details of the experimental setups are presented in Table 2. These three experiments rely on proton beams with relatively low energy so that we expect dark sector production through bremsstrahlung and direct production to be sub-dominant compared to the meson decay mechanism [18].
Notice that electron beam-dump experiments, like E137 [53], can also lead to dark sector beams through dark photon production by bremsstrahlung. However, the bounds on the kinetic mixing parameter ε from the latter were found in [54] (in a context roughly similar to ours -albeit in a supersymmetric model) to be always significantly weaker than the current missing energy bound ε < 10 − 3. The case studied in [42], which we will considered in more details at the end of this section, is a notable exception.    Table 3: Summary of the relevant characteristics of mesons productions in the experiments considered. Note that the lower energy at LSND prevents the production heavier mesons.

Dark Higgs boson production through meson decay
Proton beam-dump experiments could be practically seen as light meson factories, with around one neutral pion created for each proton on the target. We furthermore include the production of heavier η, ρ and ω mesons. The relevant number of mesons produced in each experiment is given in Table 3 based on [41,55]. We simulate their kinematic distribution by using a weighted Burman-Smith distribution to account for the different target material used by the LSND experiment over its lifetime (water, then high-Z metal) and an averaged π + and π − Sanford-Wang distribution for MiniBooNE and SBND.
The produced meson has a tiny chance of decaying into dark sector particles. In this decay, dark Higgs boson can be produced from an excited dark photon through a "dark" Higgstrahlung mechanism. The processes for the scalar meson decay are π 0 , η → γV * , V * → SV, and for the vector meson case The corresponding Feynman diagrams are shown in Fig. 5.
Focusing on the first process, we can write the branching ratio for a neutral pion as where dΠ π 0 →γV * and dΠ V * →V S represent the usual two-body decay phase space, |M| 2 is the squared, averaged amplitude, s is the squared momentum of the excited dark photon and is integrated between (M V + M S ) 2 and m 2 π 0 . The relevant quantity for our Monte-Carlo simulation is the differential decay rate dBR π 0 →γSV dsdθ , where θ is the angle between the dark Higgs boson and the excited dark photon in the rest frame of the latter. We find (see Appendix A for details) where q S is the dark Higgs boson U (1) D charge, Γ V is the width of the dark photon (which can be neglected in practice) and λ is given by The case of the η meson is completely similar, with the replacement m π 0 → m η and BR π 0 →γγ → BR η→γγ = 0.394. We have also checked agreement with the integrated standard results of [12]. The second process, corresponding to vector meson decays, is a simpler two-body decay. The branching ratio is given by and similarly for ω mesons. While the processes described above are typically suppressed compared to the onshell production of dark matter particles from dark photon decay, the dark Higgs boson on the other hand is easier to detect as one can search directly for its decay products. Note that due to the absence of gauge vertices between two dark Higgs bosons and the dark photon, the only scattering process available is through dark Higgs boson mixing with the SM Higgs boson and is therefore negligible here.

Dark Higgs boson decay and detection
As discussed before, when the decay of a dark Higgs boson into two dark photons or dark matter particles is kinematically forbidden, it is long-lived and can only decay to an e + e − pair through the loop-diagram process shown in Fig. 1a. This is in principle a very distinctive signature compared to dark matter scattering. In practice however, most of the existing experimental bounds are derived from neutrino-electron scattering signal, which consist of only one charged track. In detail, for each of the considered experiments, we have: • LSND: We choose to use the search [56] for electron neutrino ν e via the inclusive charged-current reaction ν e + C → e − + X. 6 Following [15], we will consider that the outgoing e + e − pair is interpreted as a single electron event satisfying the energy cut, 60 MeV < E e + + E e − < 200 MeV and use the electron detection efficiency of around 10%. Given the uncertainties presented in [56] (see especially Fig 29 and the Tables IV and V), and the fact that the energy distribution of our process would not have been uniform, we will consider that 25 events should have been observed and draw our contours accordingly. As was already pointed out in [15] for dark photon searches, a re-analysis of the LSND data focused on pair of e + e − events and increasing the energy threshold would significantly improve the limit from this experiment.
• MiniBooNE: We concentrate on the "off-target" dataset used in [58] for dark matter searches, and therefore require the electron and positron tracks to satisfy cos α > 0.99 where α is the angle to the beam axis and have energy in 50 MeV < E e ± < 600 MeV. The efficiency for detecting leptons is taken to be 35% from [59]. Following [42], we will require that both leptons are sufficiently separated so that miniBooNE could resolve both tracks (with a angular resolution of 2 • ). Since no such search has been yet released, we can only give projections.
• SNBD: We will conservatively apply the same lepton detection efficiency and cut cos α > 0.99 as in the MiniBooNE analysis for this experiment, as this is enough to significantly extend the reach of MiniBooNE.
Once the dark Higgs bosons have been produced, they will travel through the shielding before decaying into the detector. The probability of a decay event happening within the detector is simply given by where γ is the Lorentz factor, L d is the distance between the target and the entry point of the dark Higgs in the detector and L cr is the length of the intersection of the dark Higgs trajectory with the detector. In the limit where γvτ S L d , L cr the probability reduces to The number of dark Higgs bosons detected scales as ε 6 α 2 D so that even a tiny modification of the kinetic mixing will lead to drastic changes in the detection signature. We show in Fig. 6 the number of events expected in all three experiments considered here as a function of M S . We have chosen ε = 0.001 and α D = α em but the expected number of events for any other values of these parameters can be recovered from the previously mentioned scaling relations. In particular, notice that SBND will improve on the miniBooNE bound by one order of magnitude, provided a suitable search strategy is implemented.
Compared with the standard bounds from dark matter searches in this experiment, as in, e.g. [18,42,58], our expected number of events is even more sensitive to the kinetic mixing parameter ε. In both of our models, we found that the thermal value target is still out of reach of beam dump experiment as shown in Fig. 7, where we have shown the projected number of events at SBND. The cases of LSND and miniBooNE are similar, with no points compatible with the relic density constraint leading to more than a few expected events.
Hence the situation for dark Higgs boson search at proton beam-dump experiments is relatively similar to the one for the dark matter scattering searches in the same detectors, with the thermal value target out of reach of current experiments [18]. One interesting exception in the pDF case was pointed out in [42]. When dark matter is produced from dark photon decay, the heaviest mass eigenstate χ 2 can only decay to χ 1 through an off-shell dark photon, in the process χ 2 → χ 1 e + e − , which leads to a long lifetime of order constraints from BBN-related observables, one can search for the e + e − pair produced by the decay. The reach is then significantly stronger as we show in Fig. 8. However, their bounds depends significantly on ∆ and is rapidly not competitive for lower values.
In [42] the thermal value target was almost systematically excluded for dark matter masses in our range of interest, however the fact that the dark Higgs boson opens several new annihilation channels modifies strongly this prediction, as we show in Fig. 8. Furthermore, the presence of a light dark Higgs boson modifies even more significantly the phenomenology when M χ 2 − M χ 1 > M S . Indeed, the structure of our Lagrangian (namely the possibility of different Yukawa couplings between the righthanded and left-handed part of the original dark matter Dirac field) allows for the unsuppressed decay χ 2 → χ 1 S. Hence, in this particular regime, the previous search channel is no longer open, but should be replaced by a search for dark Higgs boson decay as described above. This production mechanism should however be orders of magnitude larger than the Higgstrahlung, as it proceeds completely on-shell, leading to much stronger bounds than the one from Fig. 7. Thus, it would be interesting to re-run the search presented in [42] while including the effect of a light dark Higgs boson. We save this analysis for future work.

BBN constraints
Bounds on dark Higgs bosons from BBN can be surprisingly strong, limiting lifetime to be as small as 0.1 s for sub-GeV dark Higgs when mixing with the SM Higgs boson is considered, as shown in [21]. As we will show in this section, these constraints   Figure 8: An example of bounds on the parameter y ≡ ε 2 α D (M χ /M V ) 4 as a function of the DM mass M χ for ∆ ∼ 0.1, α D ∼ 0.1 from [42]. We show the points from our scans satisfying all our constraints as well as ∆ < 0.1, α D < 0.1 and M V ∼ 3M χ 1 (a) and M V ∼ 10M χ 1 (b). Since the bound is weaker for smaller ∆ and α D , the represented lines are the strongest possible bounds from the analysis of [42] for both sets of points.
will be mitigated in our case due to two factors. First, due to its small mass, the dark Higgs boson decays only leptonically during BBN, and second, the annihilation mechanisms for our U (1) D -charged Higgs boson are significantly more effective than the one in [21], so that the metastable density of dark Higgs boson after freeze-out is orders of magnitude smaller.
The decay products of long lived particles like the dark Higgs boson during the evolution of the Universe can distort the agreement between the standard BBN predictions and experimental observations of primordial abundances of light nuclei, in particular 3 He and D. However, the annihilation of dark Higgs bosons during freeze-out provides a mechanism for depletion that can in turn ameliorate this potential disagreement. The energy injections from the decay of such long lived particles can be at early or late time. Here, early time refers to the early stages of BBN when t 10 s, wherein decays from a long lived particle could affect the neutron to proton ratio, n/p, or the effective number of neutrino species, N eff . Late time refers to the later stages of BBN when t 100 s which affects the final primordial abundances of light nuclei. We shall discuss constraints from both early as well as late time energy injection from S decays.
First we consider constraints from energy injection at early time. In particular, hadronic decays of dark Higgs boson, like for example mesons, occurring in the early universe could significantly alter the n/p ratio. Similarly, the direct production of neutrons and protons through quarks and gluons when the dark Higgs boson is sufficiently heavy can also give rise to stringent constraints on the lifetime of the dark Higgs boson [22,[60][61][62][63]. However, in our case we restrict the dark Higgs boson mass M S to be less than the dimuon threshold. This also means that there are no hadronic modes available for the dark Higgs boson decay and the only possible decay mode involves electrons. As a result we avoid stringent constraints from hadronic injections and in- The definition of effective number of neutrino species assumes that the three neutrino species instantaneously decouple giving a definite neutrino-photon temperature ratio T ν /T γ = (4/11) 1/3 , so that 6) where N ν = 3 is the number of neutrino species. Since the energy injected by the S decays will lead to a reheating of the electron-photon bath with respect to the neutrinos, this decreases T ν /T γ . And as can be seen from Eq. (4.6), this leads to a lowering of N eff . We use the result obtained in ref. [21] where an approximate analytical approach was adopted to calculate N eff assuming a neutrino decoupling temperature of 1.4 MeV. The 2σ lower bound from PLANCK [64] requires N eff > 2.71. We show in Fig. 9 the exclusion limit on dark Higgs boson lifetime, τ S , as a function of Ω S h 2 , the relic density the dark Higgs boson would have had today if it was stable. We see that most of the parameter space survives as a result of efficient annihilation channels of dark Higgs boson, particularly when M S > M χ and the dark Higgs boson can annihilate into dark matter thereby decreasing its abundance substantially. However, when M S < M χ this annihilation channel is not so efficient and the metastable abundance can be quite large and some of the parameter space especially above τ ∼ 100 s is ruled out. Next we consider the effect of late time energy injections at t 100 s. Such late energy injections can potentially destroy light nuclei through dissociation thereby al-tering their abundances. When the long lived particle primarily decays to hadrons the resulting hadro-dissociation can be very effective in reducing the primordial abundances even at relatively early times t ∼ 100 s. But once again in the dark Higgs boson scenario considered here the only viable decay mode is e + e − which leads to constraints only from electromagnetic showers. The absence of hadronic showers means that we avoid severe constraints from measurement of primordial abundances. The constraints from electromagnetic showers arise through photo-dissociation of light nuclei which become significant at t 10 4 s. At t ∼ 10 4 − 10 6 s, the photo-dissociation of deuterium, while at t 10 6 the over production of D and 3 He through the photo-dissociation of 4 He lead to the most stringent constraints [22]. The choice of parameters in this case as mentioned in the next section, leads to a lifetime in the range of 1 − 10 5 s. In this range of dark Higgs boson lifetime the bounds from N eff are the most stringent up to ∼ 10 4 s, and above 10 4 s the bounds from D/H and 3 He/D are the most stringent as far as BBN is concerned. In the dark Higgs boson scenario, however, there can be additional annihilation channels which can reduce the metastable abundance as mentioned in Sec. 3.1. For example, the production of a dark matter pair from dark Higgs boson annihilation can be significant for M S M χ , thereby potentially avoiding constraints from BBN. In Fig. 9 we show the exclusion limits from D/H and 3 He/D abundances. We see that most of the parameter space above 10 4 s is ruled out by these bounds, however one could still have substantial annihilation into dark matter which may allow a few points in the parameter space especially in the pDF case.
Bounds on the lifetime translate almost directly into a lower bound for the kinetic mixing parameter from Eq. (2.12). When the dark Higgs boson metastable density is large as no effective annihilation into dark matter is possible, then we have the rough bound ε 10 −4 . When M S M χ the metastable density is suppressed by the annihilation process SS → χχ which dominates over the reverse process, and the bounds are significantly weakened. Nonetheless, notice that since the dominant dark matter annihilation channel in this regime, χχ → γS, depends quadratically on ε, most points still have τ S 10 4 s as can be seen in Fig. 9. However, this is not a strong bound and some more fine-tuned points can have longer lifetime, of order 10 6 s. For such high values of τ S , the mixing with the SM Higgs boson (which we neglected following the discussion of Sec. 2.2) should become competitive to mediate the dark Higgs boson decay.

Summary and conclusions
We have argued that in models with a massive, but light, dark vector mediator, the spectrum should naturally contain a light dark Higgs boson, whose presence can substantially modify the predictions of the two models considered in this paper. In the plane M χ −M S we can identify four regions, as shown in Fig. 10, each with very distinct phenomenologies: High S aboundance Low S aboundance Scalar DM pseudo-Dirac DM Figure 10: The four phenomenologically distinct regions described in the text shown in the dark matter mass M χ versus the dark Higgs boson mass M S plane with data points from the CS case (blue) and the pDF case (orange) which satisfy all the constraints considered in this analysis. The grey-shaded regions mark the two regimes where the phenomenology does not significantly differ from previous studies.
• The short-lived dark Higgs boson regime corresponding to relatively heavy dark Higgs boson. This is the case considered in most of the previous literature, most notably recently in [42] for the pDF model. Dark Higgs bosons tend to decay instantaneously into a dark matter pair, leaving little new imprint, both in beamdump experiments and in cosmological observables.
• The long-lived dark Higgs boson regime in which the dark Higgs boson is light enough so that it cannot decay into dark photon or dark matter particles. Its decay products can then be observed in beam-dump experiments, even though the corresponding bounds are often weaker than the missing energy searches by BaBar and NA64. Depending on whether or not one has M S M χ , this regime divides into two sub-regions: -The low abundance region, M χ < M S < 2M χ , where the process SS → χχ is effective. The metastable density of dark Higgs bosons after freeze-out is therefore suppressed, so that the bounds from BBN are weakened. On the other hand, the relic density is realized mainly from χχ → Sγ processes, implying that smaller values of the kinetic mixing parameter will in general lead to overabundant dark matter, hence leading to a lower bound on ε for most of our points.  Figure 11: Summary of the various relevant bounds considered in this analysis, with the points satisfying the relic density constraint in the M S /τ S plane for the pDF model (a) and the CS case (b) overlaid. We have restricted the points to the region of longlived dark Higgs boson where the phenomenology is distinctively different from the previous studies of these models.
lation process for the dark Higgs boson. The consequent high metastable density of dark Higgs bosons translates into relatively strong bounds from BBN-related observables. Furthermore, the dark matter annihilation channel χχ → SS is kinematically open. Being an s-wave process for the model CS, this region is ruled out by CMB constraints, as can be seen in Fig. 10.
In Table 4, we give benchmark points for these regions satisfying all the constraints considered in this article.
In this paper, we have focused on the long-lived dark Higgs boson regime, as the secluded and short lived dark Higgs boson scenarios had already been covered extensively. We found that while the dark Higgs boson can in principle be produced and detected in proton beam-dump experiments, the thermal value target is out of reach of the experiments considered here. This conclusion should however be mitigated by several comments. First, as has been already advocated by many previous papers, it would be very interesting to make a re-analysis of the LSND data, possibly raising the energy threshold for the detected electrons and looking eventually for e + e − pair directly as this will significantly increase the reach of this experiment. Second, our conclusion regarding the reach of beam-dump experiments only applies to low-energy beam experiments, where the dominant production mechanism is meson decay. For more energetic beam experiments, or for electron beam dumps, different production channels for the dark Higgs boson, like direct production, should be considered. Finally, we expect that dark Higgs boson decay in proton beam-dump experiments could set stronger bounds than those of missing energy searches in the case of the pDF model when the process χ 2 → χ 1 S is available. Table 4: Benchmark points for the models analyzed in this work. Mass-related quantities are given in MeV and MeV 2 , cross-section times velocity are in cm 3 /s. On the other hand, the cosmology of the two models considered is significantly modified, with additional annihilation channels leading to various constraints. The bounds from the CMB arising from the fact that some of the new dark matter annihilation channels were unsuppressed at recombination time have been presented, excluding in particular completely the region M χ > M S in the CS case. Furthermore, BBN-related observables which arise as a consequence of the long lifetime of the dark Higgs boson were found to be relevant, but weaker than could have been expected from previous works. We have summarized the main constraints on both the CS and pDF models in Fig. 11.
Regarding earth-based experiments, the most promising discovery channels for these models seem to be the missing-energy searches as they exclude already large portion of the parameter space. In the case of the pDF model, direct detection in beamdump experiment of the decay of the heaviest dark matter field as advocated in [42] is a promising strategy which can be further combined with the search for dark Higgs boson from the χ 2 → χ 1 S channel when it is kinematically accessible. It would be interesting to study other types of cosmological probes for an extremely long lived dark Higgs boson, as for example possible supernovae-related constraints or possible signatures from dark Higgs boson (or dark photon) production in DM annihilation in the sun as was already studied in [14,65] for heavier dark matter candidates.
In the long run, the experimental prospects for both our models are bright. Almost all of the parameter space which meets the thermal value target will be independently probed by the next generation of projected electron beam-dump experiments (for instance LDMX), by direct detection experiments such as XENON1T for the CS model, and by indirect detection experiments searching for current dark matter annihilation in the MeV mass range.
In order to simulate the decay chains in our Monte-Carlo simulation, we need the differential branching ratio d 2 BR M →γSV dsdθ . Writing the two-body phase space differential element dΠ 2 , we can use the recursion relations to decompose the three-body phase space into two-body ones combined with an extra integral over the excited dark photon squared momentum s, leading to where we used the short-hand notation λ ≡ λ 1, M 2 V /s, M 2 S /s .