New production channels for light dark matter in hadronic showers

Hadronic showers transfer a relevant amount of their energy to electromagnetic subshowers. We show that the generation of"secondary"dark photons in these sub-showers is significant and typically dominates the production at low dark photon masses. The resulting dark photons are however substantially less energetic than the ones originating from mesons decay. We illustrate this point both semi-analytically and through Monte Carlo simulations. Existing limits on vector-mediator scenarios for light dark matter are updated with the inclusion of the new production processes.

One of the most compelling empirical arguments to search for extensions of the Standard Model (SM) of elementary particles is the need to explain the nature of dark matter (DM). In years past, theoretical and experimental efforts mainly catalysed around the hypothesis that DM corresponds to a Weakly Interacting Massive Particle (WIMP) with electroweak scale mass (for a recent review, see e.g. Ref. [1]). Such a hypothesis is certainly well grounded, given that in the early Universe WIMPs would be produced via thermal processes, and their subsequent annihilation with typical weak interaction rates would leave, almost independently of other details, a relic density of the correct size to match the observed cosmological amount of DM. However, null results of an extensive and long lasting search program that combined direct, indirect, and collider probes are presently triggering a waning of the WIMP paradigm [2]. While WIMP searches should certainly continue until all experimentally accessible corners of the parameter space are thoroughly probed, it is now timely and important to put no lesser vigor in exploring also other pathways.
One alternative scenario, well motivated in first place by the evidence that DM is reluctant to interact with ordinary matter, conjectures the existence of a new class of relatively light elementary particles not charged under the SM SU (3) C × SU (2) L × U (1) Y gauge group. After all, even in the SM some particles are uncharged under one or more gauge group factors, so that an extension to include a new sector blind to all SM interactions is not particularly exotic. In addition, given that in the SM there is no shortage of states with mass below, say, 1 GeV/c 2 , it is also rather natural to hypothesise that the same could be true for dark sector particles, the lightest of which would be stable, thus providing a light dark matter (LDM) candidate. On the other hand, the dark sector could well come equipped with its own set of interactions (to which SM particles should clearly be blind) and if this set also contains the simplest type of gauge force, corresponding to a U (1) gauge factor, then mixing between the dark spin-1 boson (often referred to as "dark photon" and denoted as V in this work) and the photon would naturally occur [3]. This would provide a portal though which the SM and the dark sector could communicate.
Recent years have witnessed a steadily growing interest towards LDM and its possible detection through the vector portal, and many studies have appeared deepening our understanding of the theoretical models and of their phenomenology, see for example Refs. [4][5][6][7][8][9][10]. Interestingly, besides promoting new experimental programs aiming to search both for the V and for LDM particles [11][12][13], the LDM paradigm also stimulated the reanalysis and reinterpretation of old data originally collected to search for other types of particles [5,[14][15][16]. Accelerator-based thick-target experiments at moderate beam energy (∼ 10÷100 GeV) are the ideal tool to probe the new hypothesis, since they have a very large discovery potential in a wide area of parameters space. Within this context, the main experimental techniques that have been considered so far are (1) missing energy/momentum/mass searches with electron and/or positron beams [16][17][18][19], (2) electron and proton thicktarget experiments searching for light new particles via their scattering in a downstream detector [4,20], and (3) decay of long-lived dark sector fields into SM particles [21][22][23][24][25][26][27][28][29].
Proton beam-dump experiments show an enhanced sensitivity to the dark sector. Thanks to the large beam energy and accumulated charge, typically higher than those in electrons and positrons counterparts, a large LDM signal yield is expected, usually at the price of a larger background [15,30,31]. The experimental intensity frontier is currently extremely active, and many new experiments will start to take data during the course of the next decade [11,32,33], making the accurate estimation of their potential reaches an important issue. This is particularly true for dark sector searches carried out at proton beam-dump experiments designed for neutrino physics [20,34], such as MiniBooNE [35], SBND [36], ICARUS [37] or DUNE [38], where the irreducible neutrino background calls for an even more careful evaluation of the expected LDM signal.
In the aforementioned vector portal scenario in which a light dark photon interacts with the SM sector via feeble gauge interactions, the main LDM production mechanism involved in a proton beam-dump experiment is the two-photons decay of light mesons (π 0 and η), where dark sector particles are produced thanks to the γ − V mixing. This production mechanism has been widely studied in the last decades. However, a proton-induced hadronic shower is always accompanied by an electromagnetic counterpart, which carries a significant fraction of the primary beam energy. This allows for a rich variety of electron-and positron-induced LDM production processes, incrementing the flux of LDM particles from the thick target, and thus the experimental sensitivities. An early attempt to consider this effect was presented in [39], considering only the V visible decay. In this work, for the first time we estimate the LDM production rate from the electromagnetic components of proton beam-dump experiments. We show that, in some cases, this is the dominant LDM production mechanism for a non-negligible region of the dark sector parameter space. Furthermore, we demonstrate that thanks to these new production processes proton beam-dump experiments can also probe non-minimal dark sector scenarios that were, so far, considered to be an unique prerogative of lepton-beam efforts, such as protophobic models [11] in which the dark photon coupling to quarks is strongly suppressed. The recently proposed protophobic fifth-force interpretation [40,41] of the observed anomalies in internal e ± pair creation in 8 Be and 4 He nuclear transitions [42,43] is an example of a particularly intriguing new physics scenario that the new production mechanisms allow to test also in proton beam thick-target experiments.
The paper is organised as follows. In Sec. II we describe the phenomenology of LDM production by secondary electrons and positrons in a ∼ 100 GeV proton beam-dump experiment, discussing both the main properties of proton-induced electromagnetic showers and the dominant LDM processes induced by e + and e − at this energy scale. In Sec. III we present the details of the numerical procedure that we developed to compute the enhanced sensitivity of proton beam-dump experiments, which results from taking into account the new production processes. Finally, in Sec. IV, after briefly reviewing the main features of a representative set of proton beamdump experiments, we present the corresponding exclusion limits and sensitivity curves, updated by including the new LDM production channels.

II. LDM PRODUCTION BY SECONDARY e ± IN PROTON BEAM-DUMP EXPERIMENTS
The production of dark-sector particles in a proton beam-dump experiment is a multi-step process involving the secondary particles produced in the thick target by the impinging hadron. Due to the variety of secondary particles being part of the developing hadronic shower, a large number of production mechanisms is possible. In this work, we include for the first time electronand positron-induced processes in the computation of the LDM yield of proton beam-dump experiments. In order to do so, we decouple the problem into two separate parts: the development of the EM shower which is controlled by SM physics, and the new physics processes generating the dark photon (which we assume to be strongly sub-dominant compared to the former). More in detail, we will first revisit the typical structure of the EM component of a proton-induced hadronic shower, for a primary beam energy in the 10÷100 GeV range. We will next discuss the main processes responsible for LDM production by electrons and positrons in this regime. Finally, we will focus on the production and detection of LDM in a typical proton-beam, thick-target experiment.
A. Production of e ± in proton-induced hadronic showers When a high-energy proton impinges on a thick target, a cascade of secondary hadrons with progressively degrading energy is produced, mostly containing protons, neutrons, and pions. Due to the isospin symmetry of hadron-induced reactions, approximately 1/3 of the lat-ter are π 0 . These immediately decay to high-energy γγ pairs, which in turn initiate an EM shower accompanying the hadronic one. A similar argument applies for η and η mesons, although their contribution to the EM shower is reduced, both because of the smaller production cross section, and because of the lower branching fraction for the γγ decay. On average, the fraction of the primary proton energy transferred to the EM component is of the order 50% for a 100 GeV impinging proton [44].
While a complete treatment based on numerical simulations will be presented in the next sections, a relatively good approximation of the energy distributions can be obtained from a semi-analytical approach. Starting from the typical differential number density of secondary neutral mesons n M0 (E) from a pN collision at the beam energy, with N a nucleus of the target material, the differential yield of mesons in the hadronic shower, per POT, can be estimated approximately by just considering the first interaction of the proton: where σ pN is the inelastic proton-nucleon cross section, A is the atomic mass of the target material, ρ the density, N A = 6.022×10 23 , L is the length of the active part of the target, and the second equality follows from the definition of the nuclear interaction length λ T . If the target is thick enough, L λ T , and in the approximation of only considering the first generation of secondary particles in the hadronic shower, we can set L λ T , which results in the simplified expression . Clearly, this approximation is expected to be more accurate for energies close to the beam energy, while for lower energies the actual number of neutral mesons would be underestimated. This effect is clearly visible in Fig. 1, where we show the differential π 0 yield from a 120 GeV proton beam impinging on a thick graphite target, comparing the approximate result from Eq. 1 (red curve) with that obtained from a full simulation of the hadronic shower made with Geant4 [45] (blue curve). We used the QGSPJETII software [46] to compute n M0 (E) -a similar calculation with the EPOS-LHC software [47] yielded the same conclusion.
While a thorough description of the EM shower development requires a complete Monte Carlo calculation, an approximate evaluation of the electrons and positrons track length can still be obtained with an analytical approach. We introduce the dimensionless shower depth parameter in unit of radiation length t ≡ d/X 0 and the shower age as function of the energy E and t [48,49]: where E γ is the energy of the photon inducing the shower. The differential distribution n e of electron/positron grows exponentially with s, corresponding to the power Comparison of the differential π 0 yield per proton on target from a 120 GeV proton impinging on a thick graphite target. Red curve: result obtained from Eq. 1 based on EPOS-LHC [47]. Blue curve: results of a full Geant4 -based simulation.
law scaling Ultimately, the lowest energy electrons/positrons (resp. photons) in the shower start interacting with the medium mostly via ionisation (resp. Compton scattering) and the shower stops developing. The critical energy c for which this happens is roughly defined as the energy for which the electron bremsstrahlung and ionisation rates are equal (in fact the energy at which the ionisation loss per X 0 is equal to the electron/positron energy). Denoting with Z the atomic number of the medium, the critical energy can be approximated by [50]: A direct consequence of the two energy loss mechanisms described above is that one typically expects that in a thick target the differential density distribution should be dominated by electron/positrons around the critical energy.
In order to describe more quantitatively the electromagnetic shower, we will broadly follow the approach of Rossi and Griesen [48] as reported in [49]. We refer to Appendix A for details. The first step is to obtain the differential energy spectra of both electrons and positrons, which in this formalism are treated on equal footing, neglecting initially the influence of ionisation or Compton where the factor 1/2 comes from the fact that the analytical approach does not distinguish between electrons and positrons and dNγ dE is the differential yield of photons from the mesons decay. Note that T ± (E) has dimension of GeV −1 .
We validate this approach in Fig. 2. In particular we show for reference the full result obtained from a Geant4 simulation [45], as is described in the next sections. We present the positrons track-length times energy squared distribution as function of the energy of the positrons. The semi-analytical approach carries an important uncertainty in that it does not account for the full dynamics of the hadronic shower, and it assumes instead that the initial proton interacts only once. Accordingly, the number of nuclei targets is set in Eq.(1) by what is assumed to be the "active" part of the target. We can either set L to the nuclear interaction length, or we can make the more conservative choice of setting L to the nuclear collision length, thus ensuring that the incoming proton would not loose energy before generating the shower. The results obtained for these two choices delimit the blue region in Fig. 2, which can be taken as a proxy for the typical uncertainty associated to the semianalytical procedure. In any case, we find a very good agreement with the full numerical approach. We further observe that the semi-analytical approach becomes more conservative with increasing proton beam energy (in SHiP for instance) as secondary mesons carry enough energy to generate sizeable sub-showers of their own.
One important comment is that this approach does not incorporate the angular distribution of the produced electrons/positrons. As can be readily inferred from the relative low energy of the peak of the spectrum in Fig. 2, the electrons/positrons angular distribution has a nonnegligible width. Depending on the geometry of the experiment (detector size and detector-dump distance), this effect can be critical, since it affects the angular distribution of the LDM particles produced, and thus the signal yield. In this work, we accounted for it by evaluating the double-differential track length T ± (E, Ω). As an example, Fig. 3 shows the angular distribution of positrons produced in the DUNE target by the 120 GeV Fermilab proton beam.

B. LDM production channels
The procedure described above is completely general and can be applied to any light new particle coupling to the electrons/positrons or to the light quarks (for instance axion-like particles and milli-charged particles). For concreteness, in this work we focused on the case of a LDM scenario where sub-GeV DM particles interact with the SM via a dark photon mediator V µ (with field strength F µν and dark gauge coupling g D ). The corresponding Lagrangian contains the following terms: where the parameter ε weights the kinetic mixing, B µν the hypercharge field strength, and J µ D is the dark gauge current, which depends on the details of the dark sector. After electroweak symmetry breaking, and after performing a standard redefinition of the photon field γ → γ−εV to diagonalise the kinetic term, the dark photon also acquires a ε-suppressed interaction with the SM electromagnetic current: Note that the dark photon mass m V can originate either from the Stueckelberg mechanism or from the VEV of a dark Higgs boson. The latter typically constitutes an important part of the phenomenology if it has the same mass as the dark matter candidate [23,26,51]; on the contrary, it can basically decouple if it is heavier than the dark photon. Here we will consider explicitly the second scenario. Finally, specifying the precise nature of the dark matter candidate χ is not critical for the scope of this work. In order to compare our result with the recent limits from the MiniBooNE collaboration, we considered a complex scalar dark matter candidate, although our conclusions also apply for other standard choices (Majorana dark matter, pseudo-Dirac dark matter with a small mass splitting, etc...) since their production and detection mechanisms are similar. For the case of a complex scalar the dark current is given by: As long as m V > 2m χ , the interaction in Eq. (8) leads to rapid dark photon decay into dark matter particles: this is the so-called invisible decay scenario on which we focus. Note that often in the literature an extra factor of 1/2 is included in the normalisation of the dark gauge current in Eq. (10). Thus, when relevant to carry out proper comparisons, we have rescaled the existing limits on the dark gauge coupling in Eq. (8) to account for the choice of normalisation. 1 For low mass dark sectors, the main production mechanisms for dark photon from the hadronic development of the shower are from the decay of light unflavored mesons. Depending on the mass of the dark photon, the dominant meson decay process are π 0 → γV , η, η → γV or ρ, ω → V → χχ * (in case the dark photon decays into dark sector particles). The typical branching ratio is given by In particular, note that there is no α em insertion so that this process is only mildly suppressed. On the other hand, hadronic showers develop a large electromagnetic component from the radiative decays of light neutral mesons π 0 , η. All relevant processes here depend on the density of the relevant targets (either nuclei for bremsstrahlung or atomic electrons for FIG. 4. Dominant processes for dark photon production during the electromagnetic development of a shower positron/photon processes). While the dominant production mechanism for vector mediators in electron beam dumps is mostly via electron bremsstrahlung, it was recently realised that production mechanisms based on secondary positrons can dominate in the low mass ranges [54]. Since in hadronic showers the yield of secondary electrons and positrons is almost the same, positron-related processes dominate the LDM production rate.
Denoting E ± the energy of the incoming positron/electron in the lab frame, the main processes responsible for dark photon production by secondary e + /e − are the following: • Bremsstrahlung of electrons and positrons off nuclei, e ± N → e ± N V with typical cross section: where ξ ∼ Z 2 is the effective flux of photon from the accelerated nuclei in the incoming electron/positron frame [55]. We observe that σ brem increases quadratically for small dark photon mass, but is, however, severely suppressed by α 3 em . Also, the emitted dark photons are typically very energetic, since they carry most of the energy of the initial e + /e − , with the median value for E V given by: This mechanism dominates the dark photon production in electron beam-dump experiments, due to the fact that it is enhanced for very energetic primary electrons (see the comparison with the resonant production mode in [54,56]). In the proton-shower induced environment, both the electrons and the positrons are secondary particles, and therefore they contribute equally to the bremsstrahlung production rate. We review in more detail this production mechanism and our numerical approach for this process in Appendix B.
• Direct resonant positrons annihilation on target atomic electrons, e + e − → V [57]. The cross section is given by: While this process can only occurs around the resonant energy (depending on the width of the dark photon, which is here relatively large due to the dark decay V → χχ * ), it is still important because it is only suppressed by α em [57]. Furthermore, given the restricted kinematics, the energies of the incoming positron and of the outgoing dark photon are related by This translates into a lower limit on the accessible dark photon masses: where E th is the experimental detection threshold.
• Associated production from positrons in the shower, e + e − → γV . In the limit where E + m e m 2 V , the cross section becomes which is typically α 2 em suppressed but is enhanced by a 1/m e factor. Note that, due to the presence of an additional photon in the final state, in this case the energy of the emitted dark photon can differ from m 2 V 2me . In fact, as shown in Appendix C, around half of the dark photons from associated production retain most of the energy of the incoming positron For sizeable dark gauge coupling, V decays to a χ * χ pair with near 100% branching ratio, thus allowing to easily derive the LDM production yield. In principle, χ * χ pairs can be also produced via exchange of an offshell V , and this process can be relevant especially when considering a large dark gauge coupling α D ∼ 0.1. Accounting for off-shell χ * χ production requires including in the resonant positrons annihilation the finite V width, and considering the full four-particle s-channel reaction e + e − → V * → χ * χ. More precisely, in the limit where the center-of-mass (CM) energy √ s is much larger than m V (particularly relevant for small, MeV-scale dark photon), the off-shell contribution can be estimated as: In particular, compared to the associated production (17), the log-enhanced α em log(2E + /m e ) term is replaced by the dark gauge coupling term α D . Therefore, in case α D 0.1, σ assoc is negligible with respect to σ off−shell . On the other hand, the experimental energy threshold tends to suppress both these processes with respect to bremsstrahlung. Finally, the electromagnetic shower further contain a significant number of photons, making the Compton-like scattering process γe − → e − V also a potentially relevant production channel. In the limit where E + m e m 2 V , the cross section becomes σ Compton σ assoc 2 which is also α 2 em suppressed. Note that this cross section falls much faster than that for associated production at larger dark photon mass. We present a thorough description of the impact of this channel, comparing it to the associated and bremsstrahlung production channels, in Appendix C. In conclusion, for all the experiments that we will consider in this paper, the two dominant new LDM production mechanisms are the resonant positron annihilation and the electron/positron bremsstrahlung.
In order to illustrate the respective importance of mesons decay process with respect to shower-induced ones, we present in Fig. 5 the corresponding dark photon production rates for the 8 GeV proton beam servicing the MiniBooNE experiment. We used the full Geant4 simulation described in the next section to obtain both the distribution of light mesons and the track length of secondary positrons. Interestingly, the secondary production strongly dominates in the lower mass regimes. This is both due to the fact that the meson production saturates in this regime and that the showers provide an abundant number of positrons and electrons with enough energy to produce such light dark photons. Both hadronic and shower-based processes have the same production rate for a dark photon mass around m cross ∼ 16 MeV. This "crossing" mass depends more generally on the energy available in the initial proton beam as well as on the material of the target. For instance, for the 120 The blue line corresponds to the rate for standard hadronic production processes. We have applied basic cuts on the Geant4 objects: their angle θ with respect to the beam axis is selected such that sin θ < 0.2, and their kinetic energy should be larger than 10 MeV.
GeV beam from Fermilab's main injector, which will be used by the DUNE experiment, m cross ∼ 20 MeV, while for the proposed SHiP experiment with access to the SPS 400 GeV beam and a high-Z material target, m cross ∼ 30 MeV.

C. Experimental LDM production and detection
The typical setup of a proton beam-dump experiment is shown in Fig. 6. The primary proton beam impinges on a thick target, where LDM particles are produced. These propagate straight towards a detector with cross size S placed at distance D downstream that reveals them. A sizeable amount of shielding material is placed between the dump and the detector to range out all other particles produced by the primary beam, except neutrinos.
In all the experimental setups considered in this paper, the distance between the target and the detector is much larger than the length of the target, so that the entire shower can be approximated as starting from the initial vertex. Therefore, the number of LDM particles emitted through a process characterised by a cross section σ can be computed as: where the X 0 is the radiation length of the material, ρ its mass density, A its atomic mass and N A = 6.022 × 10 23 . 2 Depending on the production process being considered, T − (E) and/or T + (E) should be used. Similarly, the differential yield dN dEχ can be obtained by replacing σ → dσ/dE χ .
As discussed before, this approach does not incorporate the angular distribution of the produced electrons/positrons. A rough estimate of the detector geometric acceptance is ε T ∼ S/(θ χ D) 2 , with θ χ being the average LDM emission angle.
If the χ couples diagonally with the V , the two main processes responsible for the interaction with the detector are the elastic scattering off electrons and the quasielastic scattering off nucleons. In the electron case, since m e m V , the electron carries most of the impinging χ energy and gives rise to an electromagnetic shower in the detector. In the nucleon case, instead, due to the nucleon larger mass, the recoil energy is typically lower, making the signal corresponding to this process more difficult to identify. For this reason, in this work we focus on the χ − e scattering process only. The differential cross section for χe → χe scattering with respect to the electron recoil energy E f in the laboratory frame is [15]: where E is the incoming χ energy and f and s stand for fermion and scalar χ respectively; The total 2 Note that the cross section has to be expressed in cm −2 .
signal yield can then be obtained analytically, convolving the differential cross section with the incoming LDM distribution and the cut efficiency for electron recoil detection. Note that while in this paper we consider the detection of LDM via its scattering in the detector, the main idea of secondary dark photon production is relevant also for other types of dark sector searches. We finally observe that, while in this work we focused on the case of LDM detection through the elastic scattering on atomic electrons, our idea also applies to LDM models predicting similar interaction mechanisms in the detector. For example, in inelastic dark matter scenarios (iDM) [21], if the splitting between the two dark χ 1 e + e − states is small with respect to the beam energy scale, the leptons-induced LDM yield in the beam dump would not change significantly. At the same time, provided m χ2 > m χ1 + 2m e , the expected signature in the detector would be either the direct decay χ 2 → χ 1 e + e − within the detector when the χ 2 state is sufficiently long-lived, or the non-diagonal scattering χ 1 N → χ 2 N , with N an atomic nucleus, followed by the decay χ 2 → χ 1 e + e − . In both cases, the result is a significant energy deposition in the detector. In particular, we note that in the limit where the heavy state χ 2 has a decay length much larger than the distance to the detector, the lower boost factor of secondary production events will enhance the detection prospects. We will investigate the effect of showerinduced iDM production in a future work (see e.g. [22][23][24][25][26][27][28][29] for recent works discussing the iDM physics case).

III. NUMERICAL EVALUATION
To re-evaluate the exclusion limits implied by existing proton beam-dump results when the lepton-induced sec-ondary production processes are properly included, and to estimate the sensitivity of planned experiments, the expected number of signal events within the detector has to be computed as a function of the model parameters, and compared with the background yield. We performed the calculation of the signal yield numerically, decoupling the evaluation of the LDM production in the beam-dump from the subsequent propagation and detector interaction as described below.
All the necessary numerical ingredients, including in particular the track length distributions used to describe the electrons and positrons from the sub-showers are available on the Zenodo online repository [58].

A. LDM production
The evaluation of the LDM production in the dump was further factorized into two independent steps: i) the calculation of the electrons and positrons track-length in the target, and ii) the computation of the LDM differential yield from e + interactions.
For each of the detector setups that we have considered in this work, and that are described in the next section, we have computed T ± (E, Ω), i.e. the electrons/positrons differential track-length distribution as a function of the particle energy and angle, by means of a Geant4 simulation. We have used the standard G4EmStandardPhysics physics list to describe EM interactions, and the FTFP BERT HP physics list to parameterize hadronic reactions. We have developed a custom class, inheriting from G4SteppingAction, that records, for each electron and positron step in the target, the corresponding particle energy and direction. The output of the simulation is the distribution T ± (E, Ω) for discrete bins of the two observables. For the energy, we have used a bin width ∆E corresponding to ∼ 0.1% of the primary proton beam energy. Since in the simulation, with default physics lists settings, the typical energy loss for each positron step inside the dump volume is already much smaller than ∆E, we did not include any explicit step limiter. Finally, to speed-up the calculation, we introduced for all particles an energy threshold equivalent to the detection threshold, discarding from the simulation all particles falling below this value. In order to make a fair comparison between the electron-and positroninduced production mechanisms with the "traditional" processes usually considered for proton beam-dump experiments involving neutral mesons decays, in the simulations we have also sampled the differential distribution n M 0 (E, Ω) for M 0 = π 0 , η.
The LDM yield in the target was then computed using the MADDUMP software [52] and a modified version of the Monte Carlo generator BdNMC [20] depending on the production process. The former is a plugin for the MadGraph5 aMC@NLO program [59,60] that allows to compute the differential yield of LDM particles in the target from the knowledge of T ± (E, Ω). In particular, we used MADDUMP to generate a list of outgoing dark matter momenta from all the leptonic production channels, including the s-channel e + e − → V ( * ) → χχ * , the associated production and the bremsstrahlung processes. For the latter, we adopted the nuclear form-factor parameterisation described in Ref. [55]. On the other hand the hadronic production processes were handled by BdNMC. For the production via light meson decays, we used the light neutral meson distributions including secondary mesons as given by Geant4 (instead of the build-in empirical distributions) and we have simulated their decay to dark matter via the vector portal. For completeness, we have further included the proton bremsstrahlung process and dark photon production via resonant vector meson mixing as it is implemented in BdNMC [20] (in particular, the timelike form factor used in the production rate is derived from [61] and hence incorporates the effect of ρ/ω meson production). We observe that, in the current version, both MADDUMP and BdNMC assume that all LDM particles are produced at the beginning of the target, neglecting the development of the EM shower in the corresponding volume. However, as already mentioned, this approximation is well justified by the much larger distance between the target and the detector.
The advantage of this dual approach, rather than handling together the description of the EM shower development and the production of LDM in a single simulation, is the fact that, for each considered experiment, T ± (E, Ω) and n M 0 (E, Ω) have to be computed only once, thus saving a significant amount of computation time. Only the evaluation of the LDM yield has to be repeated for different values of m V and m χ .
Finally, to account for the different materials in the target geometry, the procedure we adopted was to compute separately for each of them the e + /e − differential track length and the LDM yield using the procedure described before, summing the obtained results. To speed-up the calculation, only the materials with a non-negligible T ± relative weight ( 1%) were further considered.

B. Detector interaction and normalisation
We have used BdNMC to simulate the propagation and interaction of light dark matter with the detector. More precisely, we propagated the LDM particles to the detector and estimated their intersection with the detectors using the internal BdNMC routines. The scattering probability as a function of the dark matter nature (complex scalar or Dirac-fermion) was estimated using Eq. (21) (note that the complex scalar case was already present in the original BdNMC code). In order to simulate accurately the detector response, we added at the generator-level the selection cuts from the experiments. To speed-up the calculation, basic energy cuts were included directly in the cross section evaluation, while the more advanced ones (such as that on E e θ 2 e ) were applied after the scattering events had been simulated.
Finally, starting from the knowledge of the sensitivity of a given experiment in terms of signal yield, the corresponding reach curve was sampled as follows. To reduce the number of free parameters, we adopted the standard choice m V = 3m χ and α D = 0.1. Observing that in the scenario considered in this work all LDM particles are produced promptly in the beam dump, we can expect that for a given set of reduced model parameters the foreseen signal yield in the detector will scale as: where N 0 S (m χ ) is the signal yield corresponding to the kinetic mixing parameter ε 0 . We can thus obtain the limit for ε by inverting the previous relation.

IV. APPLICATIONS AND EXAMPLES
In this section, we present the revised exclusion limits and we discuss the estimates of the sensitivities that we have obtained for a representative selection of existing and planned proton beam-dump experiments, after the new positrons annihilation production mechanism is included in the evaluation of the LDM yield. For each case we briefly discuss the relevant experimental details, and the assumptions made in carrying out the analysis(see also Tab. I). Our results are summarised in Sec. IV E.

A. MiniBooNE
MiniBooNE is a proton beam-dump experiment at Fermilab, originally designed to measure short-baseline neutrino oscillations [65]. The MiniBooNE detector is a 6 m radius spherical tank, filled with 818 tons of mineral oil [35]. It is installed approximately 540 m downstream of a beryllium neutrino production target, where the 8 GeV proton beam from the Fermilab Booster impinges on.
Recently, a dedicated LDM measurement was performed by the MiniBooNE-DM collaboration using data corresponding to 1.86 · 10 20 protons on target [62]. Since neutrino interactions in the detector represent an irreducible background for the LDM measurement, the experiment was performed by steering the primary proton beam in an "off-target" configuration, to avoid neutrino production in the target, and to impinge directly on the steel beam-dump installed 50 m downstream. This resulted in a neutrino background reduction of a factor 30. The experiment considered both the nucleon and the electron scattering channel to detect LDM, with the latter providing the most stringent limits. After employing a sophisticated set of selection cuts to discriminate between the LDM signal and the residual backgrounds, zero events were observed in the signal region. This allowed the collaboration to set a 90% CL limit on the LDM parameters space, corresponding to 2.3 expected signal events.
To compute the LDM flux in MiniBoone, we described the beam dump in Geant4 as a 4 m long steel block. Since this correspond to approximately 24 hadronic interaction lengths, we ignored any further downstream material. Also, we did not include any material upstream the thick target. In this work, we only considered the χ − e − scattering process. We reproduced the MiniBooNE-DM analysis following the same strategy adopted in Ref. [20]. We parametrized the MiniBooNE-DM response with the following selection cuts, E e > 75 MeV and cos(θ e ) > 0.99, where E e and θ e are, respectively, the scattered electron energy, and the angle measured with respect to the primary beam direction. The validity of this parametrisation can be assessed from Fig. 7, where we compare the sensitivity for the "traditional" LDM production as reported by the MiniBooNE-DM collaboration (dashed orange line) with that obtained applying the aforementioned selection cuts (solid rust line) observing a very good agreement.

B. NOνA
NOνA is a neutrino experiment at Fermilab studying the oscillation of muon neutrinos to electron neutrinos [66]. The experiment measures neutrinos produced in the NuMI target facility by the 120 GeV proton beam from the FNAL Main Injector [67]. The NOνA near detector (NOνA-ND) is located 990 m downstream from the target, at 14.6 mrad angle from the primary beam direction. Such off-axis configuration was chosen to optimise the neutrino energy distribution for the oscillation measurement. The detector is a large volume of plastic (PVC) extrusions filled with liquid scintillator (active volume), followed by a muon detector made of alternating steel planes and scintillator planes. The active volume is a high-granularity sampling calorimeter, characterised by enhanced PID and tracking capabilities. The corresponding mass is approximately 193 · 10 3 kg, for a total volume of 3.9 × 3.9 × 12.67 m 3 [68].
A first estimate of the NOνA-ND sensitivity to LDM was discussed in [34] where, however, only the χ − e − scattering channel was considered. This result was based on a preliminary report of the ν − e elastic scattering analysis performed by the collaboration [69], for a total exposure of 2.97 × 10 20 POT. Both the elastic neutrino scattering signal (120 expected events) and the corresponding backgrounds (40 expected events) were treated as an irreducible background for the LDM search, for a 90% CL exclusion limit of 16.4 LDM events.
In this work, we computed the NOνA-ND sensitivity to LDM by simulating electrons-and positrons-induced production processes in the NUMI target. We implemented the official Geant4 description of the target geometry and materials, as was used to measure fundamental neutrino properties [66], and that was provided to us by the NOνA  (54) TABLE I. Beam, target, and detector main characteristics for the experiments considered in this work, along with the total number of protons on target (PoT) and the typical lower energy cut. The distance (to the centre of the experiment) D and the typical detector dimensions (length L and cross-area S) are also indicated. Note that the SHiP design is not final. We list the number of events corresponding to a 90% confidence level exclusion limit. In the DUNE-PRISM case, we considered both an on-axis and an off-axis configuration (see text for details). The references in the first column refer either to a published analysis in the case of existing constraints, or to projected bounds in the case of planned experiments.
collaboration. We considered the NOνA-ND active volume described before, with an average electron number density n e 3 · 10 23 cm −3 . Finally, we parameterized the detector response to the scattered electron with the following selection cuts: E e > 500 MeV, E e ·θ 2 e < 5 MeV, where θ e is measured with respect to the impinging particle direction.

C. SHiP
SHiP is a proposed beam-dump experiment at CERN SPS to search for weakly interacting long lived particles [63]. The SHiP detector, currently being designed, foresees two complementary apparatus, to investigate the hidden sector exploiting both the visible decay signature of hidden particles and the recoil signal from the scattering on atomic electrons and nuclei. In particular, the SHiP Scattering and Neutrino Detector (SND) is a hybrid apparatus consisting of alternating layers of an absorber, nuclear emulsion films and fast electronic trackers, characterized by a very low detection threshold and enhanced PID capability. The detector is located approximately 40 m from the production target where the 400 GeV proton beam impinges on.
A first estimate of the SHiP experiment sensitivity to LDM was discussed in [20] considering both the χ − e − and the χ − N scattering processes. More recently, the SHiP collaboration presented an updated limit for the χ − e − channel, based on a robust evaluation of the irreducible neutrino background and on a realistic parameterization of the foreseen detector response, for a total exposure of 2 · 10 20 POT [53].
In this work, we evaluated the SHiP sensitivity to LDM as follows. We computed the LDM flux due to positrons annihiliation in the beam dump with Geant4, implementing the current target geometry and material composition that were provided to us by the collaboration. We parameterized the SND active volume as a 90 × 75 × 320 cm 3 volume, located 38 m from the beam dump, with a fiducial mass of 10 ton. The following selection cuts were applied to the scattered electron kinematics, 1 GeV < E e < 20 GeV, 10 mrad < θ e < 20 mrad, with θ e measured with respect to the impinging LDM particle direction. Within this signal region, we assumed an irreducible neutrino background of 800 events [53]. This corresponds to a 90% CL exclusion limit of 38 events.

D. DUNE
DUNE is a large-scale experiment under construction in the US conceived for neutrino and proton decay studies [70]. DUNE will consist of a near detector, that will record interactions near the source of the beam, and of a much larger far detector, located underground 1,300 km downstream of the source. DUNE will detect neutrinos produced by the primary 120 GeV proton beam of the Fermilab accelerator complex impinging on a graphite target.
In a recent work it was shown that, despite the abundant neutrino background, a dedicated analysis with the DUNE near detector data will be able to explore unknown territories in the LDM parameters space, exploiting the χ − e − scattering channel [64]. In this work, we adopted the same description for the DUNE near detector geometry used in Ref. [64], considering a 3x4x5 m 3 liquid argon detector located 574 m downstream from the target. We described the target as a thin, 220-cm long graphite cylinder [71]. We parameterized the detector response with the following cuts on the scattered electron kinematics: E e θ 2 e < 2m e , E e > 50 MeV, with θ e measured with respect to the impinging χ direction.
To derive the DUNE near detector exclusion limits for LDM, we considered a total accumulated charge of 1.1 · 10 21 POT/year, and a 7-years long measurement. We observe that, as discussed in Ref. [64], the DUNE near detector sensitivity to LDM can be significantly enhanced by performing multiple measurements at different off-axis locations, to exploit the different angular spectra of the LDM signal and the neutrino background (DUNE-PRISM detector concept). In this work, for simplicity we performed a first estimate of the DUNE sensitivity to LDM produced by secondary e + considering both a single on-axis and a single off-axis measurement (at the maximum transverse distance of 36 m), leaving a more comprehensive evaluation for the future. We estimated the irreducible neutrino background for the on-axis (offaxis) measurement to be 71 · 10 3 ( 1500) events, assuming an equal experiment run time in neutrino and anti-neutrino mode [64]. This corresponds to a 90% CL exclusion limit of 350 (54) signal events.

E. Results
In this section we present our results for the limits and for the projected sensitivities of the four experiments described above, assuming that LDM is a complex scalar particle, that is for the model discussed in Sec. II B. In order to consistently compare the dark matter production via meson decay and via resonant production in the electromagnetic shower, we have used for the former the π 0 and η meson yields from the Geant4 simulation described in the previous sections. Similarly the assumptions on the detectors geometry, signal response, and backgrounds have been applied to both type of production.
In the following, limits denoted as ε e lim are based only on e + /e − processes, that is they are derived considering only dark photon interactions with the leptons. They can therefore also be used to constrain protophobic dark matter scenarios, for which proton beam-dump experiments are usually believed to have no sensitivity. For the coupling g e of a dark photon interacting dominantly with the leptons and with suppressed couplings to hadrons, the limits on the couplings are given by the simple relation: In the following figures, this "lepton-only" limit ε e lim is represented as a solid green line. Note that being electron-based experiments, the limits from NA64 and BaBar also apply in this case. We first considered the reach of the MiniBooNE experiment. As can be seen in Fig. 7, we find excellent agreement between our simulation using light meson production (orange dashed line) and the original limit from the collaboration [62] (rust solid line). This confirms the robustness of our calculations. The dotted and dashed green lines correspond, respectively, to the limits from bremsstrahlung, and from positron-induced production, including both resonant e + e − → V → χχ and associated e + e − → γV → γχχ processes. They contribute significantly to the total number of expected events for m V ∼ 20 MeV, thus significantly enhancing the full Mini-BooNE exclusion limits compared with those from the NA64 collaboration. The mass range where the pure resonant process is active is clearly visible in the plot. In particular, the lower bound at M χ1 ∼ 3 MeV (m V ∼ 10 MeV) is due to the fact that, following Eq. (16), a dark photon resonantly produced at this low mass does not transfer enough energy to the LDM particle (and ultimately to the scattered electron) to pass the E th se-  [16] and NA64 [17] collaborations. The dashed orange line corresponds to the sensitivity as extracted from [62], the rust solid line is our estimate based on hadronic processes only, the solid green line is our estimate based on secondary production processes only, and the thick black line is the combination of the two.
lection cut. For dark photon masses below this threshold, the dominant production processes are thus the dark bremsstrahlung from electrons and positrons and the associated dark photon production. Note that the limit from secondary production is conservative in that we do not include dark photon production via the Compton-like process γe − → V e − . 3 For the lowest dark photon mass, as discussed in Appendix B and C, the cross-section for bremsstrahlung increases quadratically with the inverse of the dark photon mass, while associated production saturates.
The impact of the energy threshold on the limit is further visible in Fig. 8, where we plot the expected sensitivity of the SHiP experiment. Also in this case, the comparison between our calculation (rust dashed line) and the results of the collaboration (orange dashed line) for light mesons LDM production show a relatively good agreement (notice that we did not include possible detection efficiencies in our estimate). Even if the experiment will use the 400 GeV SPS proton beam, leading in principle to high-energy electromagnetic showers, due to the high detection threshold (∼ 1 GeV) electrons-and positrons-induced processes represent only a small fraction of the final events.

SHiP collaboration
, and p brem. production secondary production Total FIG. 8. Projected reach of the SHiP experiment. The grey region represents the exclusion bounds from the BaBaR [16] and NA64 [17] collaborations. The dashed orange line is the limit extracted from [72], the rust line our estimate based on hadronic processes only, the solid green line our estimate based on the secondary production processes only, and the thick black line is the combination of both.
We show in more detail in Fig. 9 the LDM energy distribution for the different production mechanisms, for the specific choice m V = 30 MeV. The energy distribution for the leading mesons decay channel peaks as the highest energies, as expected since it originates from mesons from the primary hadronic shower. The secondary production from electrons/positrons bremsstrahlung retains a significant fraction of the energy of the shower and peaks just above the GeV. As shown in Appendix B, this is due to both the fact that bremsstrahlung dark photons typically retain all the energy of the incoming e + /e − and that the bremsstrahlung process itself is effective at large centerof-mass energy. Finally, LDM production through resonant positrons annihilation is peaked at a lower energy below the GeV, around half the energy of the outgoing dark photon E V = m 2 V /(2m e ) ∼ 0.9 GeV. In the case of the NOνA experiment, the large energy threshold E th = 0.5 GeV also limits significantly the contribution of electromagnetic shower-induced processes, with a corresponding lower mass threshold around M χ1 ∼ 10 MeV (m V ∼ 30 MeV), as seen in Fig. 10. Note that the relatively large energy threshold as well as the large distance between the beam dump and the experiment tends to reduce the contribution from the showergenerated events, since they are typically both less collimated and less energetic than their hadronic-generated counterparts. We illustrate the effect of lowering the energy threshold for the NOνA and SHiP experiments in Fig. 11. In this case, we did not combine the hadronic and leptonic limits as for the other plots, to illustrate that the background level are likely to be significantly modified, so that the proposed reaches should also be rescaled accordingly. On the other hand, it is clear that the ra-   [16] and NA64 [17] collaborations. The rust line is our estimate based on hadronic processes only, the solid green line our estimate based on secondary production processes only, and the thick black line is the combination of both.
tios between both production modes is not significantly modified by this change. In particular, in the case of the NOνA experiment, the small geometric acceptance of the experiment suppresses naturally the shower-induced events.
Finally, we present in Fig. 12 the long term prospect based on the near detector of the DUNE experiment. This experiment will adopt a much lower energy threshold than NOνA and SHiP. Consequently, we observe that the leptonic-induced events play an important role in the final production rates, particularly at small dark matter  Low , = . , = SHiP, , and p brem. prod. SHiP, secondary prod. NO A, , and p brem. prod. NO A, secondary prod.   [16] and NA64 [17] collaborations. The rust line represents our estimate based on hadronic processes only, the solid green line our estimate based on secondary production only, and the thick black line is the combination of both. masses. A particularity of the proposed DUNE-PRISM near detector concept is that it can be physically moved off-axis up to 36 m to reduce the overall background. While we did not performed a complete analysis like the one carried out in Ref. [64], we present in Fig. 13 the possible reach of the DUNE near detector in case it will be moved at the maximal off-axis distance, considering the same run parameters as the nominal on-axis mode. In- terestingly, the wide emission cone of the leptons-induced dark matter candidate enhances their importance with respect to the standard mesons decay processes.

V. CONCLUSIONS AND OUTLOOKS
When a high-energy proton beam impinges on a thick target, a large fraction of the primary energy is transferred to the electromagnetic component of the developed particles shower, resulting into an abundant production of photons, electrons, and positrons. In this work, starting from this observation, we have discussed for the first time the role of electrons-and positrons-induced processes in proton beam-dump experiments in relation to LDM searches. We have shown that LDM production from shower induced electromagnetic processes, that was so far overlooked, must be accounted for to properly assess the sensitivity of forthcoming proton-beam dump experiments, and to derive limits on the LDM parameter space from the analysis of existing data.
A numerical procedure, based on the MADDUMP and BdNMC simulations codes was developed to generate LDM particles, and, starting from the e + /e − differential track length in the target computed with a Geant4-based simulation, to propagate them into a downstream detector. We considered a representative set of proton thick-target experiments (MiniBooNE, NOνA, SHiP, and DUNE), finding that for each of them the new production mechanism results into a non-negligible increment of the sensitivity to LDM. For some regions of the parameters space, the e + /e − -induced processes actually represent the dominant production mechanism for LDM, and can lead to signal rates on par with the standard results. Due to the typically softer spectrum of LDM particles generated from e + /e − secondaries with respect to those originat-ing from mesons decays, this effect is more important for experiments characterised by low detection threshold on the scattered electron.
Before concluding, it should be emphasised that, while we focused on the case of a dark photon mediator, our analysis can be easily extended to any other LDM model. Given that the increase in the LDM particle yield that we obtain only depends on the inclusion of new production channels, our results can be relevant also for LDM searches based on detection strategies different from the simple χe − → χe − scattering considered here, as for example measurements of energy deposition in the detector from visible decays of long-lived dark sector states. Finally, while we concentrated on proton beam-dump experiments, it would also be important to properly account for the new processes analysed in this work for projected LHC-based intensity frontier experiments, such as FASER(ν) [73,74], MATHUSLA [75], Codex-b [76], ANUBIS [33] or MilliQan [77]. The extremely high energy available at LHC interaction points may actually lead to an even stronger production of dark sector particles from processes induced by electromagnetic showers. We thus believe that it would be particularly important for these experiments to consider carefully also showerbased dark sector productions, and not only to estimate correctly their sensitivity reach, but also to optimise the choice of the detection energy thresholds for the physics run.
Acknowledgments AC and LD warmly thanks L. Buonocore for very helpful discussions on MADDUMP. AC thanks P. Snopok and A. Habig for their support in implementing the official NUMI target geometry in our Monte Carlo simulation, W. Bonivento and T. Ruf for their support in implementing the official SHIP target geometry in our Monte Carlo simulation, and K. Kelly for his help concerning the neutrino background expected in the DUNE-PRISM measurement. LD and EN are supported by the INFN Iniziativa Specifica Theoretical Astroparticle Physics (TAsP-LNF).

Appendix A: Analytical treatment of EM showers
In this Appendix we describe the technical details of the analytical shower modelling. Our treatment is based on the study of the development of high energy cosmic ray showers in the atmosphere presented in Ref. [49], which in turn is based on the Rossi and Griesen approach [48].
The idea is to solve first the equations coupling the differential density of electrons/positrons n e (E, t) and of photons n 0 γ (E, t) as function of the depth parameter t (expressed in unit of radiation length), that read: where dσ b dx and dσ p dx are respectively the differential cross section for bremsstrahlung photon production and for e ± pair production, and σ p = 1 0 dx dσ p dx is the integrated pair production cross section. The two differential cross sections are given by: The first two terms in Eq. (A1) represent respectively the fraction of e ± of energy E which loose energy by bremsstrahlung, and the fraction of higher energy e ± which end up with energy E following a bremsstrahlung. The last term accounts for e ± produced via photon conversion. The two terms in Eq. (A2) represent, respectively, the photons lost to pair-production and the photons produced via bremsstrahlung. The parameter x represents the energy ratio E e /E γ between the incident e ± and the outgoing photon for bremsstrahlung, while it represents the opposite ratio for e ± pair production. The effective parameter b Z can be expressed as function of the atomic number Z of the medium as As was worked out long ago by Rossi and Greisen [48], it is possible to obtain an analytical solution for the above set of coupled equations valid for the later stage of shower development, i.e. when t, E/E γ 1. For a shower induced by a photon of energy E γ the solution reads: where we have used the primed notation for the derivatives with respect to s. The auxiliary function G γ→e (s) is defined as: while the two functions λ 1,2 read: We have used the following cross-sections momenta: which can also be straightforwardly expressed as (lengthy) expressions involving polylogarithm functions [49]. Once this un-cut distribution is estimated, the approach of Rossi and Griesen is to add a "loss" term in Eq. (A1) by replacing where c is the critical energy defined in Eq. (4). Approximate solutions to the new system of equations can be searched for in the form: n e (E, s) = n 0 e (E, s) × p 1 (E/ c , s) .
In general n 0 e (E, s) on the right-hand-side of this equation should be multiplied by a cut-off function p that can in principle be obtained by replacing n 0 e × p in the system of differential equations. In our paper, we are using for simplicity the interpolation of p, estimated at the shower maximum, that is p → p 1 (x = E/ c , s = 1) as given in [49]. Note that a good analytical interpolation in x = E/ c is given by: Finally, since the original hadronic shower produces a large number of photons with different energy, the resulting track-length distribution for the full electromagnetic shower is obtained by integrating over the initial differential distribution of photons, as shown in Eq. (7).
Appendix B: Numerical approach to bremsstrahlung processes Bremsstrahlung production of dark photons is traditionally the dominant production mechanism considered in electron beam-dump experiments. We give in this Appendix a few details about our estimation of this process via MadGraph5 aMC@NLO , starting from a brief summary of the analytical approach based on the Weizsacker-Williams approximation [55,78,79]. We present the result for the case of an incoming electron, but note that it also applies for the case of an incoming positron.
We consider the process where N is a nucleus with atomic number Z. For simplicity, we focus on the case of a monochromatic impinging beam (the extension to the realistic case through a track length approach is straightforward). We follow the notations and summarising the discussion of [55]. We define as E 0 (E V ) the energy of the incoming electron (outgoing dark photon) in the lab frame, and we introduce the ratio x ≡ E V /E 0 . As was noted in [78], the photons mediating the process are only very mildly virtual so that their interaction with the electron are dominated by their transverse polarisation. It is then possible to decompose the cross section into a real photon-electron scattering, e(p)γ(q) → e(p )V (k) where the photon has the (small) virtual momentum q ≡ P i − P f , and a form factor for the emission of the photon from the nucleus. Let us define t ≡ −q 2 (not to be confused with the the depth parameter t introduced in the previous Appendix) and call θ V the angle of the outgoing dark photon with respect to the incoming electron in the lab frame. The full cross section can be written [78]: with β V ≡ 1 − m 2 V /E 2 0 . Importantly, the cross section for the 2 → 2 process is estimated at the minimum virtuality t = t min . The term αemF π describes the effective photon flux integrated from t = t min to the total center of mass (CM) energy t max = s. It can be obtained by integrating the nuclear and atomic form factors over the virtuality: with G 2 (t) = G el 2 + G in 2 defined by with µ p = 2.79 and the proton mass m p = 0.938 GeV. Interestingly, we see that the form factors disfavour very soft or very hard photon exchanges due to either the screening from the electrons in the atomic cloud when or from the finite nuclear size in the other limit As pointed out by [55], all values of t contribute equally to the integral -in particular, the integral it is not dominated by t ∼ t min . Indeed, while the virtual photon propagator squared, 1/t 2 , is maximum at t = t min , the phase-space numerator balances it in the integral. The minimum value of t is given by where at t ∼ t min Following [55], the cross section for the 2 → 2 process at t ∼ t min can be written up to terms in m 2 e as: Putting everything together and neglecting the θ V dependence in F, the cross section can be integrated once yielding (note that the original expression from [55] missed a factor of 1/2 [79]). It is clear that this differential cross section has an approximate singularity for x ∼ 1, regulated by the electron mass at , where the subscript c 1 labels a first cutoff point. As remarked in [55], the approximation also breaks down if the virtuality is too large, yielding a second cutoff (1− The total cross section finally reads: where An important feature that can be read out from this formula is that the cross section is actually only mildly dependent on the incoming electron energy, either via the logarithm term (which saturates when the me m V contribution dominates), or via the form-factor contribution, which also saturates at high energy due to the atomic electrons screening.
We have simulated this process in MadGraph5 aMC@NLO using an effective N N γ interaction with form factor G 2 . This implies that we did not use the Weizsacker-Williams approximation for the cross section, but we directly estimated the 2 → 4 process with dark matter final states. Furthermore, in order to regulate the numerical divergence which arises for large electron energies when the exchanged photon is very soft, we have modified the form factor G 2 (t). In particular, due to the screening effects occurring when a 2 t 1, we know that this part of the phase space is sub-dominant in the final production rate. We therefore implemented a regularisation cut by setting the form factor to 0 in the "screened" region: G r 2 (t) = G 2 (t) for a 2 t > 1/3 0 for a 2 t < 1/3 .
We have explicitly checked that the value of the final cross section is not modified by varying the cut between a 2 t < 1 and a 2 t < 0.05, and agrees with the analytical expression developed above. Furthermore, we have verified that the differential distribution in angles and energy are also not affected by this regularisation procedure.

Appendix C: Associated and Compton-like process
We give in this Appendix more details about the associated production and Compton-like scattering which complement the pure resonant production of light dark matter.
The differential cross section for both processes peaks forward at θ ∼ 0, with θ the V production angle in the CM frame (although the associated production process is also enhanced in the opposite direction, θ ∼ π). For small angles and in the limit √ s m V , m e , the following similar expressions hold: In particular, both differential cross sections saturate at very small angle, when sθ 2 < 4m 2 e . The total cross sections are also equivalent, with where the factor of 2 is compensated by the fact that the associated production also generate efficiently events with a very forward photon, with the same rate as in the forward dark photon region. Hence both processes lead to similar production rates of energetic dark photons, and since the cross section does not depend on m V , we expect these rates to saturate in the light dark photon limit. Finally, note that we have considered for both processes the atomic electrons to be free (i.e. described by a plane wave wavefunction) and in particular we neglected the target electron motion [57]. Furthermore, we observe that in an electromagnetic shower, the distribution of photons actually follows relatively closely the one of the positron/electron as long as the energy is above the critical energy. One has T γ ∼ (1.3 − 1.5) · (T e + + T e − ) in most of the shower development -see for example the discussion in Ref. [49]. All in all, we therefore expect the production of very forward dark photons in the electromagnetic sub-shower to be a factor of 2 larger for the Compton-like production than for the associated production, albeit with very similar kinematics.
We have simulated the associated production process in MadGraph5 aMC@NLO using the positron track length estimated via Geant4 . As can be seen from the differential cross section Eq. (C1), the process has an approximate collinear divergence regulated by the electron mass which leads to a logarithmic enhancement of the total cross section. We numerically-regulated this divergence in MadGraph5 aMC@NLO by adding a generator-level cut on θ as θ > 10 −5 rad. Since this value is safely below the saturation value for the differential cross section 2m e / √ s in the whole range of energies considered in this work, the effect of this cut on the magnitude of the cross section is negligible. Furthermore, the associated cross section also presents an infrared divergence from soft photon emission when √ s ∼ m V , which is not present in the above formula since we assumed √ s m V . This second divergence formally cancels against the infrared divergence of the virtual 1-loop correction to the resonant production process, and represents therefore an higher order effect. That is, formally the events with a soft photon represent a QED radiative correction to the resonantly-produced dark photon. Since we are already simulating the treelevel resonant process, we imposed √ s > m V /0.95 at the generator-level, independently of the emission angle θ, to ensure that only events with sufficiently hard photons are simulated.
Summarising, we observe that in all the experiments we have considered the associated production process is FIG. 14. Production cross section for the associated e + e − → γV and Compton-like process γe − → e − V as function of the energy of the incoming particle (either a photon Eγ, a positron or an electron with energy E e ± ) in the laboratory frame. We chose ε = 0.001, mV = 10 MeV.
often sub-dominant compared to the bremsstrahlung or to the resonant production mechanism. Indeed, the latter strongly dominates due to its 2 α em scaling when enough positrons with adequate energy m 2 V /(2m e ) are produced in the showers. On the other hand, as can be seen in Figure 14, the bremsstrahlung cross section saturates at high incoming energy, while both associated and Compton-like process decrease due to their 1/s dependence. This implies that, even for very small dark photon masses where there is a 1/m 2 V enhancement, bremsstrahlung production gets contributions from positrons in the full range of energies available in the shower. Furthermore, in the opposite limit of large dark photon masses, where the resonant dark photon energy E res V , see Eq. (15), is larger than the beam energy and resonant production cannot occur, both associated and Compton-like processes are also forbidden. In this case, the bremsstrahlung process has access to a larger CM energy since it corresponds to an interaction with the nucleus, and can be effective up to E ± ∼ m V . We have included in our numerical evaluation the associated production rate, while we leave for future refinements the estimation of the LDM signal arising from Compton-like dark photon production. As pointed out above, we expect this process to be sizeable only in the limited region where the dark photons are massive enough to suppress bremsstrahlung, but light enough so that resonant production is not available due to the experimental energy thresholds.