{'fraction_non_alphanumeric': 0.03011755284157538, 'fraction_numerical': 0.014100620656122186, 'mean_word_length': 4.328253581225236, 'pattern_counts': {'":': 0, '<': 15, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Signatures of dark Higgs boson in light fermionic dark matter scenarios

Thermal dark matter scenarios based on light (sub-GeV) fermions typically require the presence of an extra dark sector containing both a massive dark photon along with a dark Higgs boson. The latter generates both the dark photon mass and an additional mass term for the dark sector fermions. This simple setup has both rich phenomenology and bright detection prospects at high-intensity accelerator experiments. We point out that in addition to the well studied pseudo-Dirac regime, this model can achieve the correct relic density in three different scenarios, and examine in details their properties and experimental prospects. We emphasize in particular the effect of the dark Higgs boson on both detection prospects and cosmological bounds.


INTRODUCTION
Weakly interacting massive particle (WIMP) as a dark matter (DM) candidate with a sub-GeV mass scale has attracted a great deal of attention (for recent reviews see e.g. [1][2][3]).
Much recent effort has gone to explore the nature and properties of DM of sub-GeV mass interacting through a light mediator of the same scale, hence reproducing WIMP scenarios at the MeV scale. Such scenarios involve the presence of a MeV scale vector boson charged under an extra U (1) D , usually referred to as dark photon, that acts as the light mediator for dark matter interactions. The kinetic mixing between the dark photon and the Standard Model (SM) photon is then invoked to introduce interactions between the dark sector and SM particles (see [4,5] for reviews). We focus on a model containing a complete, self consistent "dark sector" built from a Dirac fermion DM candidate, a dark photon and a dark Higgs boson. Interestingly, bounds from cosmic microwave background (CMB) observations can be avoided by introducing Yukawa couplings between the dark Higgs boson and the DM field which lead to an additional Majorana mass term once the dark Higgs boson acquires its vacuum expectation value.
In this article, we build on our previous work [6] and show that this model exhibits a very rich DM phenomenology. Indeed, we have identified four distinct scenarios which achieve the correct relic density while avoiding the bounds from CMB: (i) an inelastic DM (iDM ) regime, (ii) a Majorana DM (mDM ) regime, (iii) a secluded regime and (iv) a forbidden regime. The iDM regime is defined as the region of the parameter space with pseudo-Dirac DM composed of two Majorana states with a small mass splitting. This scenario has been well studied in the past, and the dark Higgs boson does not play any significant role there.
The mDM regime has larger mass splitting between the two Majorana mass eigenstates, of the order of the DM mass, and markedly different detection prospects in accelerator-based experiments. The secluded regime refers to the region where the DM mass is close to the dark Higgs boson mass such that the annihilation channel into dark Higgs bosons becomes critical in reaching the thermal target. Note that this definition of "secluded" region is different from the one used in previous literature [7]. And finally, the forbidden regime is the region of parameter space where the DM mass is close to that of the dark photon but less than that of the dark Higgs boson. Thus, the annihilation channel of DM into dark photon pair or a dark Higgs boson and a dark photon is forbidden kinematically, but nonetheless contribute to relic density due to thermal effects.
While direct detection experiments are not typically sensitive to our model due to various suppression mechanisms [5] (since scattering processes are either inelastic or velocitysuppressed), detection prospects in accelerator-based experiments are bright and have been enthusiastically studied in recent years for a variety of contexts and experimental projects (some very recent examples are for instance [8][9][10][11][12][13][14][15][16]). However, bounds for sub-GeV fermionic DM were mainly obtained in the pseudo-Dirac regime, either by searching for the decay of the long-lived heavy component of the pseudo-Dirac DM, or by focusing instead on the signatures from DM scattering, as advertised in, e.g. [17][18][19][20][21] (a important exception is of course search focusing on the dark photon itself, which remain relatively agnostic about the details of the dark sector, see [4,5] for a review).
Our main point is that the accelerator-based bounds differ significantly in all of the four regimes identified above, especially when the parameter space is restricted to the parameter space reaching the thermal target. In order to illustrate this fact, we focus on two old experiments: LSND [22] that uses a proton beam and a relatively short beamline and the electron beam dump experiment E137 [23], and recast their null results into constraints for each of the four regimes. On top of providing the current experimental status of each regime, this shows that the pseudo-Dirac regime is not the only region of the parameter space which can be significantly constrained by accelerator experiments. Furthermore, we show that when long-lived, the dark Higgs boson can significantly alter the accelerator phenomenology of our model, by either providing sizable new constraints (mainly in the forbidden regime), or on the contrary by reducing dramatically their reach. In particular, while searches for dark Higgs boson had already been advertised long ago in, e.g. [17,[24][25][26] and considered in our previous work [6], we complement these analyses by including the dark Higgs boson decay into a dark photon and an e + e − pair.
The paper is organized as follows. We begin by describing the light fermionic DM model studied here in Sec. 2. We also describe the different regions of parameter space mentioned above in more detail, and briefly review the cosmological constraints on the dark Higgs boson. In Sec. 3, we present the constraints from beam dump experiments like LSND and E137 for each of the scenarios mentioned above. In Sec. 4 we summarize our results and conclude.

A self-consistent minimal framework
Our goal is to build a model of thermal fermionic DM with mass in the sub-GeV range while relying essentially on simple SM-like building blocks. Arguably the simplest and least constrained way of doing this is to rely on a vector mediator mechanism, where the portal between the dark sector and the SM is provided by the kinetic mixing between a new, spontaneously-broken, abelian gauge group U (1) D (under which the SM is neutral) 1 and the SM U (1) Y . In this approach, two main constraints need to be factored in: gauge anomaly cancellation for the dark gauge group, and indirect detection bounds. The former implies that at least two fermionic fields with opposite U (1) D charge must be added to the theory, and the latter that their mass terms must be at least partially of Majorana type to avoid an s-wave annihilation through an off-shell dark photon. Interestingly, both constraints can be straightforwardly satisfied by including a Yukawa coupling between the dark Higgs boson and these new fermionic states.
More precisely, we define a Dirac fermion χ = (χ L ,χ R ) DM with charge 1 under the dark U (1), which will acquire additional Majorana masses from its Yukawa interactions with the dark Higgs boson of U (1) D charge q S = +2. The effective Lagrangian for the dark photon vector (V ) and the dark Higgs boson (S) fields in this minimal dark sector model is then given by where F µν is the corresponding field tensor for U (1) D , and we have introduced a kinetic mixing term parametrized by ε. This term is an invariant of both gauge groups which, even if not present at tree-level in the theory, can be generated by loop corrections of heavy vectorlike fermion charged under both U (1)s [28]. Furthermore, we will fix the value of the gaugepreserving quartic coupling mixing the dark Higgs boson with the SM Higgs to zero, since, as was shown in [6], its effects are negligible for values compatible with a natural splitting 1 Note that building models where part of the SM fields are charged under the new dark gauge group is possible but non-trivial, see e.g. [27].
between dark sector and the SM (see also [29] where this coupling is generated directly by supersymmetry, as is shown to have small values). The DM field is then introduced through the Lagrangian After the dark Higgs boson acquires a Vacuum Expectation Value (VEV), the mass of the dark Higgs boson and of the dark photon become correlated and are given by In term of Weyl fermions, the interaction and mass terms for the DM fields are where the DM mass matrix is The lightest eigenstate has typically a negative mass and the splitting between the two mass eigenstates χ 1 and χ 2 is As usual in this case, there are typically two possible conventions, one can either keep the negative mass but ensure that the rotation matrix in the DM sector Z X is real, or else ensure that both eigenstates have positive masses, at the expense of introducing an imaginary rotation matrix. Following [30], most authors considering the pseudo-Dirac limit have used the latter choice as this leads to a very simple form for the gauge interaction between mass eigenstates χ 1 and χ 2 . We will make throughout this paper the opposite choice to have potentially a negative mass for χ 1 but real rotation matrices as this is easier to implement numerically while not in the pseudo-Dirac limit. In the following, we will also denote by χ 1 and χ 2 the Majorana fermions corresponding to the previously defined Weyl fermion states for notational simplicity. We also use the notation M χ (e.g in plots) to refer to the absolute value of the DM mass M χ = |M χ 1 | and when M χ 1 > 0, we define Defining the effective gauge interaction between χ 1 and χ 2 by g D, 12 , we have which reduces to the dark gauge coupling g V when the Dirac mass dominates, as expected.
On the other hand, the DM candidate χ 1 has an effective gauge coupling given by (2.10) which vanishes in the limit where y L = y R and is very small for y L , y R 1 or y L ∼ y R .
A key element of this model is the fact that there are only two independent mass scales: the fermion Dirac mass and the dark Higgs boson VEV whose interplay will determine most of the phenomenology of the model. As we will see later, the DM mass is typically restricted by CMB bounds and relic density evaluation to be smaller than the dark photon one: M χ < M V . Freedom in the spectrum of the model is therefore mainly determined by the position of the dark Higgs boson mass and the four regimes we will identify below largely depends on it.
Analytically, we find that the dark Higgs boson is lighter than the dark photon when A similarly simple expression can be found to assess whether or not the dark Higgs boson is lighter than the splitting between the two fermion states χ 1 and χ 2 . When M χ 1 < 0, one has M S < ∆ χ when While the Higgs quartic and the dark gauge coupling are a priori free parameters of our model, we can place an upper bound by requiring them to remain perturbative at least up to the TeV scale. 2 Considering the gauge coupling alone leads to the usual result α D ≡ g 2 V /(4π) 0.5. On the other hand, and especially in the pseudo-Dirac regime where y L , y R 1, the Higgs quartic beta function is given by 12) and is positive. In the naive assumption of constant gauge coupling we get λ 1TeV . In practice, perturbativity of the dark Higgs quartic hence typically limits α D 0.1 − 0.15 depending on the size of the negative contributions from the Yukawa coupling y L , y R and of the initial value of λ S . In the rest of this paper, we make sure that our parameter space satisfies perturbativity by evaluating the spectrum obtained by using SPheno [31,32] code created by SARAH [33][34][35]. In particular, we run the initial couplings up to the TeV scale to check that they remain perturbative, then evaluate all masses at the sub-GeV scale we are interested at tree-level.

Simplified models and dark Higgs boson
Our numerical results presented in the rest of this Section are based on a scan of the parameter space of our model (see Appendix A for the input parameters and the chosen range) to identify all relevant DM regions. We used MultiNest [36] to direct the scanning process toward relic density values compatible with the result from the Planck Collaboration [37] Ωh 2 = 0.1188 ± 0.0010. The code BayesFITS is used to interface all the public codes used, including a slightly modified version of MicrOMEGAs v.4.3.5 [38], SPheno, and a heavily modified version of the code BdNMC from [20].
In this simple model, the correct relic density of DM is typically obtained through three main diagrams shown in Figure 1. Crucially, each of them suppresses the annihilation of DM at the time of CMB recombination, as required by the stringent bounds on the annihilation process [39]. In particular, the s-channel annihilation corresponding to the diagram of of p-wave type and thus velocity-suppressed. Note that the annihilation into SV can also lead to the correct relic density in a similar fashion as the V V channel, albeit in a smaller region of the parameter space (we include this region into the larger "forbidden regime" one in the following).
Based on these processes, we have identified four typical scenarios where one can obtain the correct DM relic density while avoiding CMB bounds: the inelastic DM regime (iDM), the Majorana DM regime (mDM), the secluded regime and the forbidden regime. These four scenarios are represented in the M S /M χ plane in Figure 2. As was underlined in [6], when the dark Higgs boson is long-lived (typically M S < 2M χ and M S < M V ), the relic density must be obtained by considering a two-component DM-like scenario to obtain simultaneously the correct relic density and the dark Higgs boson metastable density. Let us review the properties of each scenario in turn: • Inelastic DM regime (iDM): This scenario has been the most studied one in the literature. It roughly corresponds to the mass spectrum M χ < M S , M V with M χ 1 < 0 in our convention and ∆ χ < 0.5M χ as can be seen in Figure 2. In this regime the role of the dark Higgs boson is typically neglected and an approximate number symmetry is introduced so that y R , y L 1. Under this assumption, DM is of the pseudo-Dirac type and the mechanism to obtain the correct relic density is a coannihilation through an off-shell dark photon. This scenario has the advantage of achieving the correct relic density for a wide range of parameters and leading to strong signatures in accelerator (2.13) with couplings of typical order g 11 ∼ 0, g 12 ∼ g V and y DM 1. Note that a priori the dark Higgs boson mass is a free parameter in this scenario, as long as M χ < M S . We show in Figure 3 the transition between the two regimes : M χ < M S and M χ > M S (secluded regime, described below). An interesting observation is that one should strictly speaking also impose M χ 2 < M S since for M χ < M S < M χ 2 the dark Higgs boson already modifies significantly the relic density by depleting the abundance of the heavy state χ 2 and hence reducing the efficiency of the coannihilation mechanism.
For larger values of ∆ χ (above 0.1), this scenario is in fact already quite constrained by accelerator-based experiments if the χ 1 is assumed to make for the whole DM relic density, as we will see in the next sections. • Secluded regime: When the dark Higgs boson mass becomes of the same order or lighter than the DM particle M S ∼ M χ , the relic density is mainly fixed by t-channel 3 As described in the Appendix A, we then have to fix y L and y R so as to be consistent with this mass ratio.
annihilation into a pair of dark Higgs bosons (Figure 1(c)), corresponding to the lefthand side of Figure 3. Since the annihilation proceeds through an unsuppressed tchannel process, low values y R , y L 1 are typically preferred to obtain the correct relic density, which in turn leads to small splitting between χ 1 and χ 2 . Notice that our "secluded' regime differs slightly from the one of [5] since we consider annihilation into the scalar dark Higgs boson and not into the dark photon. Accelerator signatures of this regime are therefore similar to the iDM regime, but the relic density is independent of the kinetic mixing, leading to a larger range of accessible parameter space. A typical mass ratio for this scenario is For longer lifetimes (corresponding to ε 10 −6 ), bounds from CMB deformation (see, e.g. [43]) due to dark Higgs boson decays during the recombination era becomes relevant and can probe a dark Higgs boson metastable density up to Ωh 2 S ∼ 10 −12 . Furthermore, when the DM mass is around 5 MeV and below, CMB bounds on the effective number of neutrinos arise [44], with some model dependency on the precise value of the limit. In the same mass range and for strong gauge coupling α D 0.5, constraints from self-interaction of the DM may also become relevant [5].
A core-collapsed supernova like SN1987A gives rise to a hot neutron star or a protoneutron star environment. The creation of new particles that can interact with SM particles and their subsequent escape can lead to cooling of the interior of the proto-neutron star.
The core of the proto-neutron star is around T ∼ 30 MeV and a significant loss of energy in the form of dark sector particles in the MeV range can result in the supernova cooling at a faster rate than expected. In previous studies, the dark photon production inside the supernova core through bremsstrahlung process has been used to constrain the kinetic mixing parameter ε. A recent update of this constraint includes a DM candidate in addition to the dark photon [57]. In the context of the model studied here, the dark sector includes a further addition of the dark Higgs boson.
For SN1987A constraints on DM in our model, the inelastic DM case when ∆ = 0 is the relevant case considered in [57]. Note that the constraints for ∆ > 0 are not directly applicable here, because of the model restrictions used in [57]. In particular, the presence of elastic scattering χ 1 p → χ 1 p is neglected, unlike the case for our model where it becomes significant when the splitting ∆ χ increases (see Eq. (2.10)).
In principle the effect of this constraint on dark Higgs boson should be mitigated by two reasons: the production of dark Higgs boson being suppressed compared to dark photon and DM, and the ambiguity of the distribution of DM and dark photon inside the supernova particularly for the large values of ε and other relevant couplings considered here, for which the dark Higgs boson can scatter off dark photon and/or DM at a significantly large rate.
Obtaining a proper lower limit on the parameter space will require a calculation involving full simulation of the dark sector population inside the supernova.

DARK HIGGS BOSON AT ACCELERATOR EXPERIMENTS
We present in this section the bounds on the four scenarios presented above that arise from accelerator-based experiments. The fact that the dark sector in our model not only contains a dark photon and a DM particle, but also a dark Higgs boson and a heavy dark sector state χ 2 will significantly enhance the prospect of its experimental verification.
Most generic searches typically focus on a dark photon and either assume that it decays invisibly and search for missing-energy signatures [58,59] -the strongest bounds using this strategy are currently from BaBar [60] and NA64 [61]. Alternatively, one can also try to observe its decay as bumps in a dilepton invariant spectrum such as in NA-48/2 [62], BaBar [63] and LHCb [64]. If the dark photon decays visibly but with a longer lifetime then a range of long-baseline experiments (see [4] for a review) further probe the range ε 10 −4 − 10 −5 .
One can also search for scattering of DM particles produced through on-shell dark photon decay in various beam dump experiment. This has been an active field of research for a decade [17-21, 65, 66], and we have included the case of scattering in the LSND experiment as an example of the current sensitivity. 5 In this work, while including the previous bounds, we focus instead on stronger but more model-dependent bounds arising from the decay of either the heavy dark sector state χ 2 or of the dark Higgs boson. In particular, we complement our previous work [6,67] in several directions by including more production channels for the dark Higgs boson and by considering also its decay channel through an off-shell dark photon, which becomes dominant when M S > M V . Furthermore, we have also included the bounds from the decay of the longlived heavy dark sector state χ 2 → χ 1 e + e − , studied in [29] in a supersymmetric context and later by [65] in a setup similar to our iDM regime. Being obtained with a completely different code, our results can be considered as an independent cross-check for the iDM regime.

Dark Higgs boson production and decay
In accelerator-based experiment, the dark Higgs boson is mainly produced through either the chain decay of some heavy dark sector states (for example χ 2 , as shown in Figure 4(a)), when this is kinematically allowed, or through dark Higgstrahlung after an excited dark photon is produced as shown in Figure 4(b). In general, the former has higher rates but at the expense of a stronger model dependence. Both processes can happen either in proton beam dump experiments, for instance from meson decay, or in electron beam-dump experiment through bremsstrahlung production of dark photon (dark bremsstrahlung).
For the dark bremsstrahlung production (in our case relevant for the electron beam dump experiment E137) we have used the public code Madgraph [68]. Some details about numerical setup and the target form factors used in the calculation are given in Appendix B.
For the meson decay cases, we use a modified version of the code BdNMC [20] where we have directly included the expression for the various differential branching ratios relevant to both the dark Higgs boson and heavy dark sector states production. We refer to [6] for the details of the dark Higgstrahlung process in the case of proton beam dumps. In this work, we will focus on the remaining expressions which are needed for the dark sector state production.
In order to treat the chain decay case in proton beam dump, we have recalculated the production rates for DM in meson decays through an on-shell or off-shell dark photon including the proper rotation matrices for the dark sector states χ 1 and χ 2 . This allows one to smoothly interpolate between the various regimes described above (in particular, we recover the standard results from ref. [65] for the iDM regime). As in [6], we are interested in the differential decay rate dBR π 0 →γχ i χ j /dsdθ V , where θ V denotes the angle between χ i and the mediating dark photon V * in the rest frame of the latter and s denotes the four-momentum squared of V * . A key parameter is the coupling g ij between χ i χ j and V , which is already described analytically in Eq. (2.9) and Eq. (2.10). In term of the rotation matrix Z X relating the DM gauge eigenstates to the mass eigenstates, it is given by 1) and the vertex includes a γ µ γ 5 factor in our conventions (real mixing matrices but potentially negative mass for χ 1 ). We obtain where S is a symmetry factor equal to 1/2 if i = j and 1 otherwise, Γ V is the width of the dark photon and the triangular function λ is defined as As usual, this straightforwardly applies to the case of the η meson by replacing m π 0 by m η and BR π 0 →γγ by BR η→γγ = 0.394. The subsequent decay of the heavy dark sector state χ 2 = χ 1 S is then considered to occur instantaneously and with a branching ratio of one if it is kinematically accessible. On the other hand if M V + 2m e < M S < 2M χ , the decay channel S → V e + e − becomes kinematically accessible. Interestingly, for relatively small splitting, the decay width takes a simple form which can be parametrically expressed as As it was already noticed in [26] in a related context, the dark Higgs boson lifetime is now significantly shorter and scales only as ε 2 , similarly to the heavy dark sector decay considered in [65]. This allows one to probe extremely small values of ε since a sizable fraction of the dark Higgs bosons produced in accelerator-based experiments will decay in the detector. 6 6 It is interesting to compare this result with the expression for the lifetime of the long-lived heavy dark

Accelerator experiments constraints
In the recent years, there has been a surge in interest in examining the potential reach of upcoming neutrino-related experiments in probing light DM sectors. In this section, we rather take the approach of considering the currently applicable bounds from the existing experiments LSND [22] and E137 [23] to explore the current limits on the four scenarios described above. In particular, we want to point out in which case the presence of the dark Higgs boson affects strongly the constraints and phenomenology of the model. We want to study in detail the signatures from the secluded, forbidden and Majorana DM regime which have not been studied in details so far. We refer to the appendix for more details on the Monte-Carlo simulation used to get the expected number of events.
a. Inelastic DM regime (iDM). This regime has been already thoroughly analyzed in the context of both existing experiments and for several prospective ones. Furthermore, as shown in [6], in this case the dark Higgs boson is extremely long-lived, so that the bounds from dark Higgs boson-related processes are not relevant. On the other hand, the decay of long-lived heavy dark sector states give strong signatures in many existing and upcoming fixed-target and accelerator-based experiments. In order to cross-check our numerical tools, we have reproduced in Figure 5 the current constraints which arise from the LSND and E137 experiments and found very good agreement with the existing literature.
b. Majorana regime (mDM). When the mass splitting between χ 1 and χ 2 increases, the heavy state can decay before reaching the detector. This typically happens for values of the kinetic mixing parameter below the missing energy search limit when ∆ χ M χ for LSND and even lower values in E137 as can be seen in Figure 6. Consequently, the limits from, e.g., LSND present now an upper bound for these types of process. We illustrate this in Figure 6 for two values of the splitting parameter ∆ χ = 0.8 and ∆ χ = 2. Interestingly, potentially upcoming experiments such as JSNS 2 [69] or upgraded SeaQuest [14] will have a shorter beam-line and hence can reduce this blind spot. Additionally, the typical relic density prediction is now modified with respect to the standard iDM case and depends sector state χ 2 which can be easily derived from the results of [65], (3.6) A second interesting aspect of this regime occurs when M S > ∆ χ . In this case, the heavy dark sector state χ 2 can decay instantaneously by emitting a dark Higgs boson. Given that the dark Higgs boson lifetime is several orders of magnitude larger, the expected reach is drastically reduced. We illustrate this effect in Figure 7.
c. Secluded regime. The accelerator bounds in this regime are very similar to the iDM case. The main difference is the fact that for a given dark Higgs boson mass, a lower bound on the kinetic mixing parameter is given by the BBN bounds discussed in the previous section. We show this in Figure 8 for two typical mass ratios. Note that this bound is weakened in the region of the secluded regime with M χ M S M V as the metastable density of dark Higgs boson is depleted by its annihilation into DM. Overall this scenario is 7 For the Figure 6, we have fixed y L and y R such that the effective Yukawa coupling y DM is given by consistent with the bounds described in Appendix A.   It is therefore long-lived and typically decays into electrons for the parameter range of interest in this paper. In the absence of any other dark sector states, the typical bounds on the dark photon decay exhibit a mass dependent gap between 10 −5 − 10 −3 (see e.g [5] for an up-to-date summary). Interestingly, the presence of additional dark sector states can strongly help in closing this blind spot. Indeed, when the dark Higgs boson is heavier than the dark photon, its lifetime is short enough to lead to a sizable number of events in both  It is important to note that these bounds also applies to generic Higgsed dark photon scenario for M S > M V even without the presence of DM.

SUMMARY AND CONCLUSIONS
In this paper, we have considered a simple self consistent model of light (sub-GeV) fermionic DM, whose relic density is fixed by freeze-out. This model is often used when studying pseudo-Dirac DM since large part of its parameter space leads to this paradigm.
Our main point is that the same model, when all particles of the spectra are properly ac-  We present in this appendix the numerical setup used to obtain the number of events in the experiments E137 and LSND, whose main characteristics are described in Table II.

Production at electron and proton beam dump
In proton beam dump experiment (and more particularly LSND), the dominant production mechanism for a beam energy in the tens of GeV is the decay of light mesons. In this paper, we build upon the numerical tools developed in [6] for this setup and based on a modified version of the code BdNMC [20]. Focusing on LSND, we start by simulating the kinematic distribution of the light meson π 0 based on a weighted Burman-Smith distribution in order to account for the various target material (water, then high-Z metal) used over the experiment lifetime. The (differential) branching ratios for π 0 decay into the relevant final states χ i χ j γ and SV γ are then calculated analytically and used to sample the kinematics of the final state and the total number of produced events. More details can be found in [6], including the expression for the differential decay rate for π 0 → SV γ, while the rate for π 0 → χ i χ j γ is shown in Eq 3.2. When the heavy dark sector state's (χ 2 ) fast decay into Sχ 1 is kinematically available, we decay it by assuming the process to occur instantaneously.
Depending on the final search channel we are interested in (dark matter scattering, χ 2 decay or S decay) the relevant particle kinematics are then stored and passed to the detector simulation part of the code.
We proceed differently for electron beam dump in that we instead rely on the public code Madgraph [68] to produce the events. In our calculation, we use the following form factor for the target atom with an elastic and inelastic term given by [24,[70][71][72] G el 2 = a 2 t 1 + a 2 t where a = 111Z −1/3 /m e with m e being the electron mass, d = 0.164 GeV 2 A −2/3 with A being the number of nucleons in the target, a = 773Z −2/3 /m e with Z being the number of protons in the target, m p is the proton mass and µ p = 2.79. The atom-dark photon vertex then becomes (P i + P f ) µ √ G 2 , where P i and P f are initial and final momenta of the atom respectively, while G 2 = G el 2 +G in 2 . We implement the above atomic and nuclear form factors using the in-built procedure of Madgraph [73].
In order to account for the energy damping of the electron beam in E137 as it goes through the aluminum target, we run the previous code for three different energies of 18 GeV, 10.5 GeV and 3.5 GeV. The resulting cross-section is then weighted by the usual energy distribution from Ref. [72]: where we have directly integrated over the thick target length (approximated to infinity), E 0 = 20 GeV and E e = 18, 10.5, 3.5 GeV in our case. 8 Finally we reconstruct an unweighted sample of final states by adding events from our three energy bins proportionally to the bin's weight. This unweighted sample of final states' kinematic is then passed to the detector simulation part of the code, along with the total number of produced events.

Decay and detector simulations
This part of our simulation, also loosely based on BdNMC [20], takes as input the kinematics information of final states and the total number of produced events (regardless of the electron or proton beam dump origin of the events). The scattering detection simulation is imported from the original BdNMC code and is described in [20]. In this section, we will focus on the detector simulation for a decaying χ 2 or S particles.
We start by analytically calculating the lifetime of the particle and its decay probability within a given decay volume consisting of the detector itself for LSND and a fraction of the open space before the detector for E137 (we will comment on this particular aspect at the end of this section). We then sample the decay product kinematics using the hardcoded expression for the differential decay rate for the processes χ 2 → χ 1 e + e − , S → e + e − or S → V e + e − . At this point we determine whether any more decay occurs before the detector (e.g dark photon decay for the S → V e + e − case) then pass the decay products to the detector simulation.
First, the simulation of the LSND detector response to the event is based on the search in [74], so that we require the electrons in final state to either be reconstructed as a single track (namely with an angular spread of the e + e − pair below 12 • ) or that only one of the electron is reconstructed (using an electron detection efficiency of 19%). Overall, the reconstructed object is then required to have an energy between 18 and 50 MeV and be forward-oriented with cos θ b > 0.9, where θ b is the angle with the beam-line in the laboratory frame, in order to fake the signature of an elastic events as searched for in [74]. Following [65], we then use a 55-events limit to derive the bounds from this experiment.
Second, for the E137 experiment, we require that at least one of the electrons crosses the detector with an angle to the beam line smaller than 30 mrad and an energy E > 1 GeV. No such events were observed by E137 (see, e.g Figure 9 of Ref. [75]) so that we place an upper limit of 3-events. An important comment regarding E137 is that the decay volume in front of the detector was an open space at atmospheric pressure. The typical radiation length for electrons in the air is of 304m [76], but perhaps more importantly the corresponding Moliere radius is 73m. This implies that after a distance of order meter, the energy of the electromagnetic shower is spread out in a transverse direction significantly larger than the spread of an event as seen by the E137 (see Figure 8 of Ref. [75]). While a proper modeling of the electromagnetic shower will be probably needed to go beyond the order of magnitude evaluation, we will take the simpler approach of considering only those events for which the