Dark fluxes from electromagnetic cascades

We study dark sector production in electromagnetic (EM) cascades. This problem requires accurate simulations of Standard Model (SM) and dark sector processes, both of which impact angular and energy distributions of emitted particles that ultimately determine flux predictions in a downstream detector. We describe the minimal set of QED processes which must be included to faithfully reproduce a SM cascade, and identify a universal algorithm to generate a dark sector flux given a Monte-Carlo simulation of a SM shower. We provide a new tool, $\texttt{PETITE}$, which simulates EM cascades with associated dark vector production, and compare it against existing literature and"off the shelf"tools. The signal predictions at downstream detectors can strongly depend on the nontrivial interplay (and modelling) of SM and dark sector processes, in particular multiple Coulomb scattering and positron annihilation. We comment on potential impacts of these effects for realistic experimental setups.

General ultraviolet (UV) completions of the Standard Model (SM) [1], the existence of dark matter [2,3], and neutrino masses [4,5] motivate the existence of feebly interacting particles [6,7].If the mass scale of these particles is below the electroweak scale, they must be singlets under the SM gauge group [8][9][10], dramatically restricting the possible set of interactions with the dark sector relevant at experimentally-accessible energies.This observation has lead to several benchmark "portal" interactions to be identified as the primary targets for experimental investigation [6,11,12].
Intensity frontier experiments with thin or thick targets [13][14][15][16][17][18][19][20][21][22][23][24][25][26][27] often provide the most sensitive probes of dark sectors, especially for particle masses below a few GeV [6,11].It is therefore important to quantify and maximize the impact of near term intensity frontier searches.Crucially, this requires a detailed understanding of both the production and detection of physics beyond the SM (BSM).Visible detector signatures involve an energy deposition in a detector via a decay into visible particles or via scattering.The former is easily evaluated, while the latter requires some atomic and/or nuclear physics input.For the purposes of dark sector searches a theoretical precision of O(10%) is acceptable, and this is achievable in many instances.
The determination of production of BSM particles in a thick target/beam stop is more challenging.High energy particles generate hadronic and electromagnetic (EM) cascades within the target leading to an exponentially growing multiplicity of final states with energies ranging over multiple decades [28,29].The production of dark states from these cascades has a complicated dependence upon the nature of the couplings to the dark sector, and on the distribution of energies and particle multiplicities in the hadronic and electromagnetic sectors.Particularly for low-mass dark states, there are many possible sources of emission such as meson decay and bremsstrahlung.Historically this complexity has been circumvented in one of two ways.First, one can consider BSM particle production only from the attenuated primary beam, since the attenuation can be estimated analytically for lepton beams [30].Alternatively, one can focus on production in the first interaction length, such that attenuation and effects of subsequent hadronic and EM cascades can be neglected [15,[31][32][33][34][35].Both approaches are conservative, systematically under-predicting the dark sector production rates.This guarantees that excluded regions of parameter space based on total rate analyses are robust against uncertainties, since additional production channels can only strengthen the bounds obtained from null observations.
For electrophilic gauge extensions of the SM, e.g., L e − L µ or L e − L τ , the dominant production mode for dark states is from the EM shower.Similarly, in certain axion models secondary electrons dominate axion production even in a proton beam dump [36].Conversely, if a dark photon couples democratically to charged leptons and hadrons, it has large production rates from meson decays [37,38] or proton bremsstrahlung [39,40].While we do not specialize to a particular UV model, we will focus on the former case where the EM shower is the dominant source of dark sector particles.Thus, in this work we study EM shower evolution and its contributions to dark sector particle production in thick target environments.While this problem has been discussed previously [41][42][43][44], a variety of different approximations have been used in the literature, some of which are mutually incompatible; there is currently no consensus on which predictions are reliable.With this in mind, our goals in this work are the following: cascade from dark sector production, enabling a self-contained treatment of each.The problem of dark sector production in an EM shower therefore represents a perturbatively calculable problem that can be computed via a series of controllable approximations.In fact, a similar approach could also be used for dark sector production in hadronic showers, but we do not explore this possibility here.
To distill the EM dynamics down to their minimal components, we have implemented a custom lightweight package, PETITE (Package for Electromagnetic Transitions In Thick-target Environments), specialized for the production of light dark sector particles (in addition to electrons, photons and positrons) in a beamdump or other fixed target experiment with ∼1−100 GeV beams. 1 We make use of well-documented and controlled approximations whenever possible, and model matter as a set of screened Coulomb sources, and a homogeneous gas of electrons at rest.We find that a good description of SM dynamics can be obtained with just three hard processes and a reasonable model of both continuous energy losses and multiple Coulomb scattering.We will show that some of these central processes are mismodelled in existing SM simulation tools.We also identify disagreements between existing results in the literature for BSM production in EM cascades, stemming from practical or conceptual issues.For example, at the practical level, stable numerical simulation of light particle bremsstrahlung requires careful MC sampling of phase space, which is difficult to achieve with many existing MC tools due to approximate phase space singularities.At the conceptual level, we find that certain methods for reprocessing a SM simulation of an EM cascade to produce dark sector states based on rescaling of cross sections are theoretically inconsistent.This method predicts incorrect kinematics of BSM particles, and can lead to order-of-magnitude errors in the signal rate for realistic experimental setups.Our simulation tool, PETITE, does not attempt to replicate or supplant existing tools such as GEANT-4 [47][48][49], however it can serve as a lightweight, fast MC and/or as a useful diagnostic tool.Important questions can be asked, and answered, using PETITE such as • How do SM uncertainties on cascade propagation influence BSM fluxes?
• How should one efficiently simulate and/or sample BSM events from a SM shower?
• Which methodologies are reliable?Do they agree?And if not, why?
Although far from exhaustive, we try to answer an interesting subset of these questions in Sec. 4. We focus on five concrete phenomenological issues which are relevant for dark sector production in EM showers: 1) the impact of multiple Coulomb scattering (MCS) modelling on geometric acceptances, 2) the treatment of resonant dark vector production and its interplay with continuous energy loss and MCS, 3) kinematic correlations in simulations of SM bremsstrahlung, 4) a comparison of dark sector bremsstrahlung in PETITE to an approximate method that instead re-purposes simulated SM bremsstrahlung kinematics, and 5) calculation of total fluxes at downstream detectors for a set of realistic experimental configurations, in a full PETITE simulation.
The goal of this paper is to provide a brief, but pedagogical, introduction to the physics of EM cascades, and to explain how BSM microphysics should be consistently implemented in any future simulation.It is organized as follows.In Sec. 2 we describe EM cascades in the SM and their implementation within PETITE.In Sec. 3 we discuss relevant production mechanisms for dark sector states, specializing to massive vector bosons.Our focus is on EM production, however we also highlight for the reader regions of parameter space in which hadronic production is known to dominate the flux in UV-complete models like the dark photon.In Sec. 4 we investigate the phenomenological questions listed above.In each case we focus on experimentally motivated distributions, but defer detailed analyses of specific experiments to future work.Finally in Sec. 5 we summarize our findings and comment on future directions.
We discuss several technical details in the appendices.In Appendix A we collect the known expressions for QED processes needed for simulating EM cascades and discuss various approximations that were used to obtain them.In Appendix B we describe the atomic and nuclear form factors that enter dark and SM bremsstrahlung, and pair production.The PETITE models for continuous energy loss and MCS are described in Appendix C. Fast and reliable simulation of SM processes necessitated a multitude of variable transformations and integrator tuning which we discuss in Appendix D. In Appendix E we validate the implementation of these processes by comparing PETITE to a semi-analytic shower model and to GEANT-4.In Appendix F we describe the numerical integration of dark sector bremsstrahlung and compare our implementation to other available programs.Finally, we provide explicit expressions for dark Compton scattering, and e + e − → γV in Appendix G. Our implementation of positron annihilation into dark vectors includes soft-and collinear photon resummation which is discussed in Appendix H.

Electromagnetic cascades
We first discuss the processes necessary to describe EM cascades involving electrons, positrons, and photons, before explaining how they are modelled in PETITE.These particles can interact with atomic electrons or nuclei in the material.In nuclear interactions, we treat the atoms as static and include a form factor which accounts for both atomic screening and the loss of nuclear coherence as the momentum exchanged with the target varies.These form factors are described in Appendix B. In what follows we denote nuclei by their atomic charge Z.Here we focus on the dominant components of EM showers, namely e ± and γ, but muons can be included straightforwardly.

Relevant processes
There are five processes which dominate the dynamics of an electromagnetic cascade.These can be classified as "discrete" and quasi-continuous.For discrete processes we generate kinematics of outgoing particles in individual events -expressions for the differential cross sections can be found in Appendix A. The three reactions of this type are 1.Bremsstrahlung on atomic nuclei: e ± Z → e ± Zγ.This process is dominated by small momentum transfers, Q, that satisfy QR ≪ 1 where R is the nuclear radius.We treat the nucleus as a static fixed Coulomb potential, use the Bethe-Heitler matrix element [50], and include form factors to account for screening and finite nuclear size effects at small and large momentum transfers respectively [51], see Appendix B. A minimum energy cut (whose default value is ω cut = 1 MeV), is applied to avoid infrared (IR) singularities from soft photons, see Appendix A for further details.While this cut determines the number of discrete "steps" taken in a cascade, BSM observables are not sensitive to this as long as the dark sector particle mass or detector threshold is larger than ω cut .
2. Bethe-Heitler pair production: γZ → e + e − Z.The matrix element for pair production is related to that of bremsstrahlung via crossing symmetry.The kinematics are similar, and we make use of the same small angle approximations, Thomas-Fermi screening form factor, and Bethe-Heitler matrix element.The process is IR finite and no minimum energy cut is required.
3. Scattering on atomic electrons : (γ, e ± ) e − bound → (γ, e ± )e − and e + e − bound → γγ.We use tree-level expressions for Bhabha, Møller, and Compton scattering as well as for e + e − annihilation, see Eqs. (A.1) to (A.4).Atomic electrons are modelled as free and at rest.These reactions tend to dominate at lower energies and are responsible for an asymmetry between electron and positron populations.Atomic electrons that are ionized above threshold in the simulation lead to an excess of e − over e + .For these processes, we enforce a minimum energy E min of the outgoing e ± .
In addition to these reactions, there are two processes that are much more common when e ± , γ transit material.They are so frequent that they must be modelled as quasi-continuous in practice; we will refer to these collisions as "soft" to contrast them with hard reactions above which are discrete in nature.The soft processes we consider are 1.Coulomb scattering on nuclei: e ± Z → e ± Z.This process is dominated by frequent small angle scattering which leads to a Gaussian angular distribution with power law tails that can be modelled using multiple scattering theory [52,53].We discuss the impact of different MCS models on dark sector fluxes in Sec.4.1 (the specific models implemented in PETITE are described in Appendix C.2) and highlight the importance of MCS in positron annihilation in Sec.4.2.

2.
Ionization: e ± e − bound → e ± e − free .Atomic ionization dominates energy loss for electrons at energies below the critical energy, O(10 − 100 MeV) depending on Z [54].Furthermore, at higher energies where bremsstrahlung is more important (and which we simulate explicitly) ionization can still influence the shower dynamics.We treat ionization energy losses with a constant dE/dx = 2 MeV cm2 /g, i.e., independent of energy (this number is only weakly dependent on the type of atom [54]).
In addition to these two processes there are contributions to energy loss coming from soft photons emitted in bremsstrahlung where the outgoing photon energy is below ω cut .However, for the default value of the cut that separates hard from soft, ω cut = 1 MeV, the contribution from below-threshold bremsstrahlung will be subdominant to the energy loss from ionization, see Appendix C for more details.
The small number of processes, along with the applicability of perturbation theory makes EM cascades well suited to a systematic investigation of secondary production.We make further approximations for some of the reactions above that we now detail before moving on to discuss PETITE itself.
Pair production is subject to collinear divergences, while bremsstrahlung has overlapping soft and collinear divergences.However a further singularity exists in both processes which actually dominates the phenomenology and this is the low-Q 2 singularity of Coulomb scattering.This makes the phase space for bremsstrahlung and pair production very peaked at small Q 2 .We make use of small-angle approximations that treat the collinear singularities analytically.At the same time this approach simplifies the treatment of the Coulomb singularity.The relevant expressions for differential cross-sections can be found in [55] and also in Appendix A.

Implementation within PETITE
We now focus on the explicit implementation of the above processes in PETITE.The design of PETITE involves a Python class called Shower.A Shower is initiated with a material (e.g., graphite with Z = 6), a minimum energy cut (E min ≥ ω cut ) below which particles are not tracked dynamically, and a dictionary of auxiliary information (e.g. the differential cross sections as a function of energy).
A Shower is a collection of particles produced in a cascade, each of which is represented by the Particle class.This class contains all relevant information for a given particle, including its four-momentum, position, and a dictionary of additional information.This dictionary contains a unique ID number, the Particle Data Group identifying number (PDG-ID), its parent's PDG-ID and the process and particle from which it was generated, its generation number (the number of hard interactions before the particle was created), and weight (for weighted events).Other info such as the mass and stability of the particle are also stored.
An EM cascade is generated by calling the Shower method generate_shower which takes as an input an initial Particle that seeds the shower, e.g., an electron or a photon with four momentum (E, p x , p y , p z ) entering at position (0, 0, 0).By instantiating a Shower with arguments describing the initial particle, material, E min , and path of pre-computed VEGAS integrators, PETITE sets up the necessary infrastructure to calculate the showers.The shower itself is generated by calling generate_shower, which loops through all Particle in the Shower.At each iteration, Particles are propagated using the method propagate_particle which determines the step size x step and whether a discrete/hard interaction takes place in that step; the probability of such an interaction is P = 1 − e −xstep/λ(E) with λ(E) the mean free path calculated using the tabulated cross-sections stored in Shower. 2 The four-momentum and position of the Particle are updated at each step, taking into account multiple scattering and energy loss (see above and Appendix C.1). Steps are repeated until a reaction takes place (as determined by an accept-reject algorithm using the probability P ), or until the energy of the particle falls below E min .
Having conditioned on a reaction occurring, the reaction type is randomly selected according to the relative cross-sections ("branching ratios") available to the current Particle.The functions used in the shower generation that encode the physical processes (e.g., bremsstrahlung, Compton scattering, etc), can be found in all_processes.py and radiative_return.py,while MCS implementations are found in moliere.py.The most computationally intensive step is generating the final state kinematics for the selected process.In PETITE this is achieved using importance sampling from a pre-trained VEGAS integrator at a nearby (slightly larger) energy chosen such that the phase space of the pre-trained integrator is a strict superset of the phase space of interest, see Appendix D. The nearby energy is used as an estimator, but the matrix element used for importance sampling is evaluated at the energy of interest.This enables us to correctly (and relatively quickly) sample the full multiparticle phase space in each reaction.The generate_shower loop continues until the energy of all active particles has dropped below E min .Note that while PETITE comes with pre-computed VEGAS integrators, the user may want to generate others, for example a coarser grid of generators or new integrators for different dark vector masses.To this end, PETITE comes with a set of utilities to facilitate such a task.
Parent particles are kept in the event record, along with their daughters which are created at every vertex.For example, after a bremsstrahlung event e + Z → e + γZ the event record will contain the initial parent positron e + , the daughter positron e + , and the daughter photon γ.All the information of the parent positron is kept for bookkeeping, but it is not propagated further.Daughter particles propagate until interacting.Thus, each Shower contains the full shower history, and the user can then plot distributions, and/or event displays using the provided utilities.
A code-focused guide for PETITE, along with some example applications, is available on the associated GitHub repository .A visualization of the output from PETITE is shown in Fig. 1.A comparison of the performance of PETITE and GEANT-4 is provided in Appendix E.

Dark sector fluxes from electromagnetic cascades
We now focus on the flux of dark sector particles that can emerge from an EM cascade in a thick target.Some of the earliest discussions of light particle production in beam dumps occurred in the context of searches for fractional charges [56], heavy leptons [57] and neutrino-like particles [58,59].In these experimental works, analytic approximations of the secondary spectra of photons and electrons, due to Tsai and Whitis [30], were used to estimate signal yields in thick targets.Later Tsai's analytic results were applied to axion bremsstrahlung from an attenuated electron beam [36,60].This method has been adopted to reinterpret old experimental results in the context of other dark sector particles like the dark photon [61,62].A different signal prediction technique was used in analyzing the electron beam dump E137 data [63]: simulated SM shower spectra were convolved with cross-sections for multiple axion production mechanisms, including bremsstrahlung, annihilation, and Primakoff-type reactions.To our knowledge this was one of the first works to systematically consider a full set of processes possible in a EM shower for a dark sector model.In Ref. [64,65] the secondary production of muons is considered since the focus of those works is on muonphillic interactions.More recently in Refs.[42,43], it was realized that for generic dark sector couplings, resonant e + e − → X production can be extremely efficient for m X ≲ 50 MeV, possibly supplying the dominant flux of new physics particles at experiments with primary or secondary positrons.This qualitative conclusion was supported by an independent study in Ref. [44].Secondary fluxes of photons have also been demonstrated to be of interest in the context of axion like particles [36].
Recently, several approaches for beam dump dark sector production have been developed.In Ref. [44] a scheme was proposed to model dark vector fluxes by reweighing the distribution of SM photons produced in a GEANT-4 simulation; crucially this procedure does not account for possible differences in angular distributions of visible and dark states.Two GEANT-4 native solutions to the production of dark sector states in medium  were presented in Refs.[66,67].In these works, SM cascade and dark sector production are simulated simultaneously by using an unphysically large dark sector coupling.As emphasized in Ref. [43], this approach is inefficient in two distinct ways.First, the dark sector coupling cannot be made too large; otherwise the properties of the SM cascades are altered, ultimately leading to inaccurate predictions of dark sector fluxes.
The second disadvantage is that one cannot reuse the expensive part of the simulation (the SM cascade) as one scans over BSM particle masses and portal types.The authors of Ref. [43] have implemented a solution to the above problems by separately simulating the SM cascade in GEANT-4 and BSM production in MadDump [41] or BdNMC [31].By linking multiple independent codes, it becomes difficult to assess systematic uncertainties in SM or dark sector production modelling.Our goal is to develop an algorithm to predict dark sector production from all processes that can occur in a given model, which can address some of the issues highlighted above.As a concrete example we will study the case of a dark sector vector V with mass m V which couples to predominantly to the electron vector current.There are a multitude of UV-complete models which contain this, and other couplings, in their low energy limit.In general, at the energies we are concerned with, the IR couplings could be to baryons, leptons, and mesons: In this work we will primarily focus on the case where the coupling to electrons is the dominant coupling, c e ≡ εe and c p,n,π ≈ 0. This can occur for electrophilic gauge extensions of the SM, e.g., L e −L µ or L e −L τ , and EM shower-induced production dominates.Similarly, in certain axion models secondary electrons dominate axion production even in a proton beam dump [36], although through a different Lorentz structure.Alternatively, a dark photon couples democratically to charged leptons and hadrons (|c p | = |c e | = c π = εe), giving large production rates from meson decays [37,38] or proton bremsstrahlung [39,40].We will not explicitly discuss hadronic production modes and refer the interested reader to the relevant literature, e.g., [15,[31][32][33][34][35].

Beyond the Standard Model fluxes from electromagnetic cascades
EM production in a thick target stemming from Eq. (3.1) is dominated by three main mechanisms: 1. Dark bremsstrahlung: e ± Z → e ± ZV .In the static limit (i.e., the nucleus is taken as infinitely heavy), this process is kinematically allowed for E e ≥ m V .Due to the non-zero dark vector mass, the matrix element for dark bremsstrahlung contains collinear, but not soft divergences.The collinear divergences only appear in the limit where both the electron and vector boson can be treated as approximately massless.This results in a matrix element that prefers E V ≫ m V .In turn this means dark bremsstrahlung dominantly produces hard dark vectors [61,[68][69][70][71].This highly peaked differential cross-section presents challenges for fast and accurate MC sampling.In Appendix F we discuss our approach, using the adaptive VEGAS algorithm [72,73] to sample the differential cross-sections presented in [70,71].

Dark Compton scattering: γe
This process is operational provided E γ ≥ m 2 V /(2m e ).The corresponding differential cross section can be found in Appendix G.

Positron annihilation: e
bound → V (nγ).This channel can be the largest contribution to the dark vector flux when it is open.The kinematic threshold is set by s = s th with s th = m 2 V .For s > s th the emission of initial state radiation of soft and collinear photons brings the energy onto the resonant peak.This is the radiative return effect seen in e + e − collisions.At threshold the tree-level expression for e + e − → V γ diverges like σ ∝ 1/(s − s th ).When integrated over an incoming positron energy distribution, this produces a logarithmic singularity.One-loop radiative corrections are need to ensure the cross section is finite.This issue can be dealt with at tree level by using an IR cutoff on the outgoing photon energy [43].Here we follow the lepton-collider literature [74][75][76][77] and make use of resummed QED distribution functions for the electrons and positrons in the collision.This removes the singularity in the vicinity of the reaction threshold, and resums large double-logarithms (see Sec. 4.2 and Appendix H for more details).

Generating dark vectors from a Shower using PETITE
We now describe a model-independent algorithm for the production of new weakly-coupled physics in a SM cascade.The small coupling between SM and the dark sector ensures that production of dark states is a small perturbation on the evolution of the shower.Thus, the dark production and the SM shower development can be simulated separately at leading non-trivial order in the small coupling.While we focus on the specific implementation within PETITE, the general set-up outlined here can be reproduced in other codes.Importantly, because the dark sector production is entirely separate from the generation of the SM cascade, any standard shower simulation tool can be used to generate the initial SM shower.This can then be "dressed" with BSM production.
Whether generated by PETITE, or in another simulation, let us assume that we have a list of SM particle momenta, ⃗ p, and production positions, ⃗ x, associated with a SM shower.Our goal is then to produce a possible BSM interaction position and probability for each particle in this list.For each dark sector process p available to the SM particle, the probability of interaction dw (p) occurring between z and z + dz is proportional to the probability for the BSM process to occur and to the probability of no SM process happening up to that point: where z = 0 corresponds to the SM particle's starting position ⃗ x.The cross section σ BSM and mean free path λ SM MFP depend on z through the particle's energy E(z); for photons no continuous energy loss occurs and so the energy is fixed.We assume a homogeneous target with number density n (p) , which is process-dependent.Note that the mean free path is determined by SM processes in the small dark vector coupling limit, i.e., λ MFP ≈ λ SM MFP .Sampling z from this distribution therefore gives an interaction energy E(z) at which the differential cross-section for p is sampled to generate dark vector kinematics.This event is assigned the weight w (p) = dz dw (p) /dz.The collection of dark sector particle momenta produced in this way represents a set of possible shower histories which produced a single dark sector particle each.
For most processes the distribution in Eq. (3.2) is relatively flat, but, as discussed above, the cross section for positron annihilation, e + e − → V (nγ), is highly energy dependent with a peak at s = m 2 V .Due to the multiple Coulomb scattering, there can be significant deflection of the parent particle between its initial energy E 0 = E(z = 0) and the energy at which the dark vector is produced.It is therefore crucial for resonant e + e − → V (nγ) production to properly sample from Eq. (3.2), Sec.4.2 highlights the impact of this sampling procedure.
The PETITE algorithm for producing dark states from SM cascade is summarized below.
PETITE algorithm shower = existing Standard Model shower (e.g. as generated by PETITE) for particle in shower do if particle = e + then processlist = {Dark bremsstrahlung, Dark annihilation} else if particle = e − then processlist = {Dark bremsstrahlung} else if particle = γ then processlist = {Dark Compton} end if for p in processlist do if energy of particle is above energy threshold for process p then • determine the interaction position/energy by sampling from Eq. (3.2).
• generate dark sector kinematics of the outgoing BSM particle by drawing from dσ • record the weight w (p) position, momentum, and particle species.end if end for end for We will occasionally refer to this procedure as resampling.It generates a list of weighted dark sector emissions that could occur in a SM shower and it can be repeated for several distinct SM showers to obtain average predictions for dark sector production.Note that the procedure outlined above produces an equal number of weighted events for each available dark process.Processes can instead be sampled based on their relative importance within PETITE.
This algorithm is implemented in PETITE in the DarkShower class.A DarkShower is initiated by passing it a Shower object which contains a pre-simulated SM shower.This can be performed either by specifying an external file location or passing a Shower object created by PETITE.The BSM events are produced by running the method generate_dark_shower, which loops through all particles in the shower producing weighted events.

Algorithm for long cascades
In the discussion above we have assumed that less than one dark sector particle is produced per electromagnetic shower.The multiplicity of particles in a shower grows rapidly with incident energy, which can result in O(few) or greater dark sector particles being produced within a single shower.In this case one must take care to account for the correlations among the descendants of a given parent particle.For instance, suppose our parent particle is an electron with energy E 0 .Let us suppose that the probability of dark sector emission can be estimated by P (DS|e − ) ∼ O(ε 2 ).If the typical number of SM daughters in a shower is ⟨N SM ⟩ and then one must consider multiple dark emissions within a single progenitor's shower.To do this one would draw the number of daughters from a Poisson distribution with the appropriate mean, then sample each daughter sequentially, altering the subsequent dynamics of the SM shower along each leg.
If ⟨N D ⟩ is so large that the ultra-weak coupling limit is not valid, then the algorithm discussed above can be viewed as a subroutine.In this more general case one first draws the number of dark sector particles produced, executes PETITE algorithm, then modifies the shower dynamics along the appropriate leg.The procedure is repeated until the O(few) dark sector events are produced after which the SM simulation is set back to its initial configuration and the entire procedure repeated.

Comparison to previous work
PETITE differs from previous implementations discussed in the introduction of this section in four key ways.First, having a minimal SM reaction network, using PETITE we have been able to pinpoint the impact of SM simulation details on downstream fluxes.In practice, this means that we can easily turn off different SM processes, modify their modelling, and study their effects on dark sector yields.We provide some examples of these studies in Sec. 4. Second, in contrast to Refs.[66,67] we can re-purpose the same SM shower for many different BSM parameter choices allowing for efficient parameter scans.Third, in contrast to Ref. [44] we sample BSM events from their true underlying differential cross section.This avoids issues that stem from very different differential distributions of SM and BSM bremsstrahlung.This problem is explored in more detail in Sec.4.4.Finally, in contrast to Ref. [43] we work at the level of MC event record, rather than energy-angle distributions.This is advantageous for tracking uncertainties and allows greater flexibility for interfacing with different event generators.Implementation differences aside, we find reasonable, although not perfect, agreement of dark vector yields between PETITE and Ref. [43] in Sec.4.5.

Key processes and their impact on phenomenology
We now discuss applications of PETITE to experimentally-relevant questions, and compare its predictions to existing methods for dark particle production in a thick target.The typical experimental setup we are interested in consists of a beam impinging a thick target and a downstream detector at a small solid angle from the target, such as the MiniBooNE, DUNE and BDX experiments.While specific detector signatures rely on the details of a given experimental setup, we focus on the predictions at the level of the dark sector particle flux.Several experiments of interest have a proton beam, which can initiate a EM cascade through π 0 production and decay, and through bremsstrahlung.Since PETITE does not simulate these hadronic interactions, we obtain the EM showers either from GEANT-4, which simulates the full hadronic and EM showers, or from PYTHIA [78], which provides the primary π 0 whose decay photons are showered with PETITE.Overall, we find only marginal differences between GEANT-4 and PYTHIA when using them to predict relevant observables.
In what follows we consider both SM and BSM processes whose detailed implementation can potentially impact experimental results.First, we examine the modelling of multiple Coulomb scattering (MCS) and its effect on the angular spread of SM cascades and of the resulting dark vectors in Sec.4.1.MCS turns out to have an interesting interplay with dark sector production from resonant positron annihilation, which we discuss in Sec.4.2.Next we describe the challenges of simulating SM bremsstrahlung and compare PETITE to other software in Sec.4.3; we find that the sampled phase space distributions can differ significantly between different codes, but these differences can be washed out in thick targets where, e.g., MCS plays an important role in determining the shower angular spread.In Sec.4.4 we highlight the qualitative differences between SM and dark sector bremsstrahlung which make it difficult to map SM photons into dark sector particle emissions, as was done previously.Finally in Sec.4.5 we use PETITE to estimate dark sector signal yields for a realistic set of experimental parameters.See also Appendix E for validations of PETITE standard EM showers against GEANT-4.

Multiple scattering and downstream fluxes
There are a variety of Standard Model effects that can conceivably influence downstream fluxes.We will start by studying how the modelling of MCS can change the angular distribution during a shower's development, thereby affecting downstream acceptances.Modelling differences can arise due to the specific implementation in a code (e.g., choices of particle step size) or variations in the microscopic input (e.g., the probability distribution of scattering angles).Here we will focus on the latter effect.
The goal of MCS theory is to compute the distribution of scattering angles f (θ, t) after a charged particle has traversed a thickness t of material.Analytic calculations of this quantity rely on approximating the atomic potential by the Thomas-Fermi model [52,79].Comparisons of MCS models implemented in GEANT-4 find ∼10% agreement with empirical measurements [80,81].The general form of this distribution is Gaussian with power law tails, which result from the central limit of a large number of small angle deflections and rare large angle scatters, respectively.In many practical applications, only the Gaussian approximation is used with different choices of its width [53].We will illustrate the potential impact of MCS modelling by calculating the angular spread of dark vectors resulting from EM showers simulated with two different MCS models.
We start with the method of Molière [79], which was improved by Bethe [52].We will summarize Bethe's results to explain the different approximations used in computing MCS.The number of electrons in an angular interval dθ after traversing a material thickness t is f (θ, t)θdθ.Molière theory allows an asymptotic expansion of f (θ, t) in a parameter, B, which encodes the depth traversed in a given material.Specifically, B is the solution of the transcendental equation where χ 2 c ≡ 4πN tα 2 Z(Z + 1)z 2 /(pβ) 2 is a reference scattering angle that characterizes the probability of Rutherford scattering, with N the atomic number density, Z the atomic number, and z, p and β the charge, momentum and velocity of the propagating particle; the screening angle χ 2 a = (1.13+ 3.76α 2 )λ 2 /a 2 encodes the suppression of Rutherford scattering, where α = zZα/β, and λ = 1/p is the de Broglie wavelength of the electron and a = 0.885Z −1/3 /(αm e ) is the Fermi radius of the atom.The screening angle can be determined by calculating the single scattering by a Thomas-Fermi potential.More details are given in Ref. [52], and in Appendix C. The ratio χ 2 c /χ 2 a is approximately the mean number of scatters in a thickness t [52,53], so unless t is so small that the Molière theory would no longer be applicable.The expansion developed by Bethe can be summarized as where f i are given as integrals containing Bessel functions.In particular, f 0 is a Gaussian and the higher f i lead to a broadening of the tails of this Gaussian.The point here is that this approximation, while not trivial to compute, can be systematically improved, as Bethe argues that B typically assumes values between 5 and 20, and the approximation can be valid even for large angles.
The second model we consider is a Gaussian-only approximation to MCS due to Lynch and Dahl [53].The width of the Gaussian is also characterized by χ 2 c and χ 2 a as detailed in Appendix C, with numerical  Typical photon energies from π 0 decays from a 120 GeV proton beam range from 1-30 GeV.(Right) The same distribution below the cut p T /p Z ≤ 2.5/574.Results are obtained from 3000 independent cascades simulated with PETITE.The annihilation probability is roughly four orders of magnitude higher than the probability of dark bremsstrahlung and the curves have been re-scaled so as to fit on the same axes.factors fit to a Monte Carlo simulation of many discrete Coulomb scatters.This prescription is more accurate than the slightly simpler Gaussian distribution given in the PDG [54].
In the left panel of Fig. 2 we show the impact of the choice of MCS model on the angular distribution of dark vectors of mass 5 MeV produced in EM showers initiated by a 10 GeV photon impinging a carbon target.Note that the typical photon energy from the decay of π 0 produced when a 120 GeV proton beam hits a target is 1 − 30 GeV.The green curve labelled "Bethe-Molière" is obtained by retaining the first two terms in the expansion discussed above, f 0 and f 1 .The magenta curves correspond to the Gaussian Lynch-Dahl MCS model.
To better appreciate the impact of MCS modelling, we show a DUNE-like acceptance as a vertical blue line (2.5 m/574 m = 4.3 mrad), and present a "zoomed in" plot of the accepted events in the right panel.The more complete Bethe-Molière theory predicts a systematic shift to larger values of p T , with effects on the order of ∼10%, and can lead to some geometrical acceptance impact on downstream detectors.For example, the geometric acceptance for dark vectors with mass m V = 5 MeV produced via resonant annihilation changes from ϵ geo ≈ 0.49 when using the Lynch-Dahl implementation to ϵ geo ≈ 0.43 when using the more accurate Bethe-Moliere theory.

Resonant positron annihilation
Positron annihilation to dark sector states via e + e − → V (nγ) has been shown to be the most efficient production mechanism available in EM cascades for light dark sector masses m V ≲ 100 MeV.This is true in a general EM cascade [42], and even in proton beam dump experiments [43,44].Recent results from the NA64 collaboration explicitly make use of resonant production [14,82].
Despite being the dominant production mechanism for certain dark vector masses, theoretical treatments of e + e − → V (nγ) in the literature do not account for several important effects.For experiments with a detector far away from the target the most relevant effect is the interplay of annihilation and the modelling of MCS (the other effect, modelling of radiative return, is discussed in Appendix H).Because the annihilation  cross-section is strongly peaked at E res = m 2 V /2m e , typical positrons with E e + > E res in the cascade must propagate through the material to lose energy through ionization and soft bremsstrahlung until they reach this resonance.During the propagation the e + undergoes MCS which modifies its direction; the direction of the positron at the interaction point then determines the direction of the outgoing V , since the initial electron is nearly stationary and any initial state radiation is assumed to have negligible p T .In previous work it was instead assumed that the initial positron direction determines the V direction [43].Below we will show that the effect of MCS on the positron as it propagates between its birth and resonant annihilation into V can have significant impact on experimental acceptances and must be treated carefully.
In Fig. 3, we present the differential flux of resonantly produced dark vectors in a thick graphite target struck with 120 GeV protons for m V = 3 MeV (left) and 6.2 MeV (right).We simulate π 0 production from pp collisions with PYTHIA and inject them into PETITE, which decays π 0 → γγ and generates the EM shower.The fluxes are given as a function of p T /p z , effectively the angular spread of the beam; vertical dashed lines represent the approximate angular acceptance of ArgoNeuT and a DUNE-like near detector, at 0.2 and 4.3 mrad.To illustrate the importance of propagating the positron down to the interaction point, we present two fluxes.In green, we assume that the positron direction when it produces the dark vector through e + e − → V (nγ) is given by its direction when it was created upstream in the shower, pi , leading to a dark vector direction pV = pi .The orange curve demonstrates the impact of taking the V direction as the direction of the positron when it annihilates into V (i.e., it is sampled after propagation and energy loss), namely, pV = psamp .This propagation accounts for MCS which induces a non-negligible angular spread for the positron.We can clearly see that compared to the previous case V tends to be created with significantly larger opening angles.In particular, properly accounting for MCS can reduce the number of dark vectors within the geometrical acceptance of DUNE or ArgoNeuT by roughly an order of magnitude.Taking DUNE as a concrete example, the geometric acceptance differs by a factor of roughly 4 for m V = 6.2 MeV, and by a factor of roughly 15 for m V = 3 MeV.
Besides properly modelling the angular spread, PETITE also includes a careful treatment of singularities from soft and collinear photon emission by making use of QED parton distribution functions that resum QED corrections with leading-log accuracy.We do not attempt a "matching and merging" with fixed order results, leaving this for future work.The resulting cross section prediction differs from naive tree-level estimates by as much as ∼20%, see Appendix H for further details.Last, we note that we do not include atomic binding in our current simulations.

Standard Model bremsstrahlung
Both bremsstrahlung and pair production are difficult processes to simulate due to their non-trivial phase space singularity structures.The bremsstrahlung cross section is singular for collinear and soft photon emission, and in the low momentum transfer limit (see Appendices A and D for more details).Collinear singularities are regulated by the finite electron mass, while Coulomb singularities are regulated either by kinematic cuts or by atomic screening.The IR divergence due to soft photon emission is effectively cut-off by the minimum energy threshold of the shower.The difficulties of sampling from such a highly singular distribution, and the need for fast efficient code has lead the developers of popular tools for EM cascades such as GEANT-4, EGS-5 [83], and FLUKA [84,85], to approximate kinematic samplings by default. 3We focus on bremsstrahlung although default implementations of pair production, which is related to bremsstrahlung by a crossing symmetry, have the same issues.In particular, in GEANT-4, a full five-dimensional sampling of phase space for e + e − pair production is available using the Bethe-Heitler matrix element via the method G4BetheHeitler5DModel [86].
EGS-5 and GEANT-4 sample approximated distributions of photon energy and angle, and neglect momentum transfer to the nucleus, which leads to a spurious correlation between electron and photon transverse momentum (see Sec. 2.7.1 of [83] and Sec.10.2 of [87]).FLUKA uses the Berger-Seltzer parameterization which is two dimensional [84,88].This reduction from a four-dimensional to two-dimensional phase space makes sampling much faster.The electron phase space, however, is then not properly correlated with the outgoing photon.To illustrate this point, we compare the correlation between electron and photon transverse momenta for the first splitting of a 10 GeV electron traversing graphite predicted by PETITE, GEANT-4 and MadGraph in Fig. 4. In GEANT-4, the final p T of the electron is set equal and opposite to the p T of the photon, essentially constraining the nuclear recoil to be vanishing (middle panel).In reality, these variables are not strongly correlated as can be seen in the left panel.MadGraph (right panel) lacks coverage of the low p T phase space, which is not a problem for predicting observables in high energy colliders such as the LHC, but in principle could affect dark sector production in EM showers.We note that predicting events at such low p T region is not the raison d'être of MadGraph, and in order to produce this panel, we have turned off all cuts in available in the run card.Overall, bremsstrahlung is an unusually difficult process to simulate, due to its three-fold singularity structure.
Although these differences are quite large for individual bremsstrahlung events, the effect on a bulk shower is much less pronounced due to angular spreading from MCS.We have checked that the downstream effects of mismodelled bremsstrahlung p T distributions are not dramatic in thick targets but may be relevant for thin target experiments such as LDMX.

Dark vs. standard bremsstrahlung
PETITE uses the exact tree-level dark bremsstrahlung differential cross-section to produce kinematics and to normalize rates.This allows us to highlight the differences of this process to its SM counterpart.The matrix element for SM photon bremsstrahlung, e ± Z → e ± Zγ, is strongly peaked towards small momentum transfers, leading to soft and collinear divergences.The outgoing photons tend to be produced with small energy fraction and with small p T /p z (the characteristic opening angles are θ γ ∼ m e /E e , c.f. Eq. (A.5)).For dark vector bremsstrahlung, on the other hand, nonzero m V regulates these divergences, leading to qualitatively different kinematics for m V > m e : the dark vectors tend to carry a significant fraction of the incoming e ± energy and their emission angle is correlated with this fraction.While technical details can be found in Appendix A, we illustrate this point in a concrete example.Let us take an experimental setup similar to the BDX experiment at Jefferson Lab [18].We simulate an EM shower from a 10.6 GeV electron beam impinging on a thick target.In the left panel of Fig. 5, we show the kinematics (p T /p z versus energy) of the photons produced via bremsstrahlung, as this is the dominant production mechanism.The preference for soft and collinear photon emission can be clearly seen.To compare to dark vector kinematics, we simulate the production of m V = 100 MeV dark photons in PETITE.For this mass, production due to hadronic reactions (e.g., π 0 → γV ), positron annihilation and dark Compton scattering are all subdominant, leaving dark bremsstrahlung as the dominant channel.In the middle panel of Fig. 5, we present the kinematics for the dark vectors produced via dark bremsstrahlung.We can see two main differences with respect to photon kinematics: the V spectrum is much harder, and there is a strong correlation between E V and its opening angle.For reference, we plot a horizontal line corresponding to the BDX detector geometrical acceptance, a 50 cm × 50 cm target placed 10 m downstream from the back of the target [18].
Due to these kinematic differences, one cannot estimate the dark bremsstrahlung flux by simply reweighting the photon bremsstrahlung flux in EM showers by the ratio of total cross sections.To appreciate this, we present the dark bremsstrahlung flux obtained using a re-weighting procedure from Ref. [44] in the right panel of Fig. 5.This scheme uses the distribution of photons produced in a GEANT-4 simulation and applies a re-weighting factor bin-by-bin in energy and angle.The number of photons, N γ (θ i , E j ), produced through bremsstrahlung in the (i, j)-th bin centred around angle θ i and energy E j is converted to the number of dark vectors, N V (θ i , E j ), in the same bin using the (angle-independent) ansatz In the right panel of Fig. 5 we can see that the re-weighting procedure leads to mismodeling of the dark bremsstrahlung fluxes, with a much stronger preference for harder and more collinear dark vectors.As we will see later, this can have significant impact on dark vector flux predictions in specific experimental setups.The re-weighting procedure can overestimate the dark vector flux by 2 or more orders of magnitude at MiniBooNE or DUNE.For example, for the BDX setup we consider in Fig. 5 the re-weighting procedure overestimates the dark vector geometric acceptance (i.e. the fraction of the binned histogram below the dashed black line) by 35% (> 99% acceptance vs. 74% acceptance).At the level of the overall rate we find that Ref. [43] over-predicts the yield by roughly 5000 for m V = 100 MeV (we believe this comes from the normalization of f (x) being off for large m V ).

Impact on experimental yields
In this subsection we compare existing results from the literature with the output from PETITE, including all effects discussed above.We focus on the simplest figure of merit, which we take to be the total number of dark sector particles passing through a detector in experimental setups similar to MiniBooNE and DUNE.In this subsection (and this subsection only) we consider a coupling of the dark vector V µ to the hadronic electromagnetic current J µ with the same coupling as electrons, i.e. |c p | = |c e | = c π = εe in Eq. (3.1).This arises, for example, when considering kinetic mixing portals with dark photons.We compare dark vector production in the Booster Neutrino Beam (BNB) for dark vector with parent particles' opening angles up to sin θ beam < 0.2 from PETITE and Ref. [43] in the left panel of Fig. 6.The other two panels show PETITE and Ref. [44] predictions for the dark vector flux at MiniBooNE (middle panel) and DUNE (right panel).The output of PETITE for dark bremsstrahlung agrees well with Ref. [43], as can be see in the left panel.Although our treatment of MCS in positron annihilation is different from Ref. [43] its impact appears to be negligible for the large angular cut considered in that panel; for more realistic angular selection these differences are more apparent as we see in the middle and right panels.The difference between the meson decay production curves (blue lines) in the centre and right panel are due to heavier mesons (η, η ′ ) and proton bremsstrahlung, which we have not accounted for (we only considered π 0 decays).We have not been able to reproduce the direct hadronic production yield curve of Ref. [43] as shown in the left panel (PETITE is roughly 2× as large), but find that the photons which dominate the EM cascade are reasonably well modelled.We therefore find qualitative, but not perfect quantitative, agreement with Ref. [43].While secondary meson  The flux emanating from the Booster Neutrino Beam with parent particles satisfying a cut of sin θ beam < 0.2 on their angle relative to the beam axis.This configuration does not correspond to a specific detector and was chosen for direct comparison with Ref. [43].Flux of dark vectors passing through the MiniBooNE (Centre) and DUNE (Right) detectors, which have been chosen for direct comparison of PETITE with Ref. [44].
production affects the total rate of dark vector production for large opening angles, it is a subdominant effect when considering detectors far downstream from the target.The middle and right panels reveal discrepancies with Ref. [44].As we discussed in Sec.4.4, the reweighting procedure for dark bremsstrahlung production (orange lines) used in Ref. [44] overestimates the dark vector production in the forward direction.By simulating the dark vector bremsstrahlung production using the appropriate matrix element, we find a much lower geometrical acceptance.An immediate consequence of this is that dark bremsstrahlung is typically subdominant compared to other production mechanisms for masses between 1 MeV and 1 GeV in proton fixed target experiments.Note that we still expect dark bremsstrahlung to be relevant in electron fixed target experiments for heavier masses where positron annihilation (green) is kinematically forbidden, particularly for hadrophobic models such as L e − L µ where dark bremsstrahlung can compete with primary production from meson decays.
Another major finding of our study, is that positron annihilation rates (green lines) are modified substantially for detectors with small angular acceptances (see Sec. 4.2).We reaffirm the findings of Refs.[41][42][43][44] that dark annihilation is an important production channel, and find reasonable quantitative agreement with their estimates for heavier dark vector masses.At low masses, we find discrepancies as large as a few orders of magnitude, which stem from the proper sampling of p T from MCS.We are able to reproduce the results of [44] both using their methodology and a modified version of PETITE that uses the initial positron momentum to determine the vector direction (c.f.Sec.4.2).Overall, due to small angular acceptances, we find that a proper simulation of dark bremsstrahlung and a proper sampling of dark vector p T in positron annihilation can change detector yields in realistic experiment setups by up to several orders of magnitude, depending on the mass of the dark vector.This directly affects estimated experimental sensitivities to a wide range of models, including dark photons, L e − L µ and other flavour combinations, and B − L.

Discussion
The prediction for dark sector fluxes emanating from a thick target is a highly non-trivial task due to the number and complexity of processes involved in realistic models.In this work we have tackled the problem of secondary production in an electromagnetic (EM) cascade.We have developed PETITE (Package for Electromagnetic Transitions In Thick-target Environments), a Python program which includes a lightweight EM shower generator, as well as a novel and universal algorithm for dressing (or resampling) generic SM Monte Carlo (MC) simulations with rare dark sector events.While we focused on estimating the production of dark vector particles in EM showers, our algorithm can be used for hadronic cascades, other BSM models or even for pure SM applications.In the latter case the resampling technique enables the simulation of rare but important events, providing a simple alternative to the MC biasing method [89].As resampling can be applied to any MC data (or implemented in any program), it can be easily interfaced with experimental simulation pipelines.
We have performed comparisons between PETITE and other methods for predicting SM and BSM fluxes, extensively validating our approach and revealing several improvements enabled by PETITE.For example, PETITE fully samples SM bremsstrahlung phase space, in contrast to existing tools such as GEANT-4, EGS-5, and FLUKA which predict an incorrect correlation of the outgoing photon and electron momenta.The significance of this effect depends on the target thickness, as it can be washed out in thick targets.We also addressed outstanding disagreements in predictions of detector fluxes from BSM particle bremsstrahlung, finding that our first-principles-based simulation can correct potential signal yields by orders of magnitude.Finally, we have identified a nontrivial interplay between the multiple Coulomb scatterings (MCS) and dark vector production via positron annihilation.Due to the angular spread induced by MCS and the small geometrical acceptance of realistic experimental setups, properly accounting for MCS in positron annihilation can change downstream yields by more than an order of magnitude.
Taking it all together, different qualitative conclusions can be gleaned from our work depending upon the primary incident particle on target.Focusing on the specific scenario of a dark photon for concreteness, we conclude the following.
• Protons on target: At lower dark photon masses, positron annihilation dominates the dark vector flux, but MCS broadening makes the flux lower than previously estimated.For higher masses, production is mostly from meson decays.Dark bremsstrahlung rarely dominates the flux.
• Electrons on target: Dark bremsstrahlung is important, and thus should be properly simulated accounting for the full phase space.A substantial positron flux emerges after a few radiation lengths, and positron annihilation is relevant.As expected, hadronic production is negligible.
We highlight that our findings are relevant for other dark sector models, such as B − L or L e − L µ .In all cases, we have substantially improved theoretical control over the flux predictions of dark sector particles originating in EM showers.Other improvements can be implemented within PETITE, such as the inclusion of the Landau-Pomeranchuk-Migdal effect relevant at high energies, more realistic models of ionization, and MCS, other BSM scenarios, and standard hadronic shower processes to extend our treatment to all showers.Besides, PETITE can also be useful to simulate rare standard model processes (e.g., µ + µ − production) via the resampling technique developed here.Looking forward, it will be important to perform fully realistic simulations of dark sector particle production in upcoming experiments, such as the SBN detectors, DUNE, DarkQuest, and BDX.Furthermore, although PETITE offers a fast and lightweight solution to the generation of EM showers, it is meant to complement rather than replace more comprehensive generators such as GEANT-4 or FLUKA.Therefore, it will be useful to interface the resampling algorithm presented herein with these tools, together with possible improvements on certain microphysics, e.g., SM bremsstrahlung.are supported by Fermi Research Alliance, LLC under Contract No. DE-AC02-07CH11359 with the U.S. Department of Energy, Office of Science, Office of High Energy Physics.KJK is supported in part by DOE grant DE-SC0010813.NB was supported in part by NSERC, Canada.We used Package-X [90] and FeynCalc [91] to derive and check several analytic results in this work.Our numerical work was enabled by VEGAS [72,73], numpy [92], scipy [93] and matplotlib [94].

A Review of QED processes
Processes involving two-to-two kinematics in QED are well known.Here, for completeness, we supply the relevant differential cross sections for Compton scattering (e ± γ → e ± γ), electron-positron annihilation (e − e + → γγ), Bhabha scattering (e + e − → e + e − ), and Møller scattering (e − e − → e − e − ): We have taken the electrons, photons, and positrons to be unpolarized in all scattering processes.Note that s, t, and u are defined according to the momentum ordering written for each process, p 1 p 2 → p 3 p 4 , where Bremsstrahlung and pair-production in the field of a nucleus were first computed by Bethe and Heitler [50]; textbook discussions with modern notation can be found in Itzykson and Zuber [95], and (in the greatest detail) in Landau and Lifshitz [55].In the small-angle regime (which dominates) the differential cross section for bremsstrahlung process of e ± Z → e ± Zγ is given by [55] Here the photon energy, in the lab frame, is ω and the incoming (outgoing) lepton energy, also in the lab frame, is ϵ (ϵ ′ = ϵ − ω).We cut off the IR range of the photon energy, but do not alter the upper bound from the maximum allowed kinematically, so ω cut < ω < ϵ − m e .The variables δ (δ ′ ) are proxies for θ (θ ′ ) which are the angles between the outgoing photon and the incoming (outgoing) lepton.They are defined as The angles δ and e /m e ]; in the limit m e /E (′) e → 0 this tends to the positive real axis [0, ∞).The azimuthal angle between the incoming and outgoing lepton is denoted ϕ and lies in the range 0 ≤ ϕ < 2π.The momentum transfer to the nucleus, Q 2 = −q µ q µ = q 2 , can be expressed as The total phase space d Π = dωdδdδ ′ dϕ.We also explicitly include an atomic form factor F (Q 2 ) 2 [51,96] which screens interactions at small Q 2 -it has the behaviour that for small Q 2 , F (Q 2 ) 2 ∝ Q 4 (eliminating the Q −4 singularity present in Eq. (A.5)) and for large Q 2 , it scales like Z 2 (the unscreened case).The effect of screening is most apparent for large-Z targets and large incident ϵ.See Appendix B for explicit definitions of these form factors and the appropriate simplifications we take for SM processes.Pair-production is related to bremsstrahlung by crossing symmetry and therefore it may be described using a similar set of approximations.It is the dominant process for sufficiently high-energy photons propagating through a thick target [54].The differential cross-section in the small-angle approximation is given by [55] Here, the incoming photon has energy ω, and the outgoing electron (positron) has energy ϵ − (ϵ + = ω − ϵ − ).Similar to Eq. (A.7), q is the momentum transferred to the nucleus and in this approximation, the expression for Q 2 is The integration phase space is d Π = dϵ + dδ + dδ − dϕ.The energy ϵ + spans the range [m e , ω − m e ].The other three variables have the same range as for electron bremsstrahlung, and δ ± are related to the outgoing positron/electron angles θ ± through a similar relation as Eq.(A.6).The atomic form factor F (Q 2 ) 2 is the same as the above, as well.For both pair production and bremsstrahlung a nuclear form factor can be similarly included, but its effect is minimal because of the dominance of small Q 2 where F nucl. (Q 2 ) ≈ 1.This is not necessarily true for dark photon production, which we treat separately in Appendix F. For high-Z nuclei, Coulomb corrections modify the cross section [51,97,98].We have not included this effect within PETITE, which means that the PETITE mean free paths are somewhat shorter than in nature.The advantage of not including Coulomb corrections is that the theory is entirely self-consistent (i.e., the total and differential cross sections derive from the same formulae).

B Form factors
Depending on the momentum exchanged when a particle interacts with material, the scattering may be coherent off the entire atom, nucleus, a part thereof, or off an atomic electron.The form factor F (Q 2 ), accounts for both atomic screening at low Q 2 and the loss of nuclear coherence at large Q 2 .It contains both a coherent (elastic) contribution [51,99], which scales as the square of the nuclear charge, and an incoherent (inelastic) contribution which scales linearly with Z Where the two contributions are [61] In each case the form factor is a product of terms accounting for atomic and hadronic effects.For G el these factors encode atomic screening of nuclear charge and the finite size of nucleus, while for G inel they describe atomic excitation and finite proton size.The atomic parameters are given by The proton's magnetic moment is µ p = 2.79, and m p is the mass of the proton.Notice that the second bracketed term in Eq. (B.3) is not squared, see Footnote 4 of Ref. [43].Atomic screening is important for SM bremsstrahlung and pair production.When considering pair production of electrons, one almost never encounters Q 2 large enough that the nuclear form factor, or incoherent scattering contributes.Therefore, in PETITE for SM interactions, G inel is set to zero, and only the atomic form factor in included in G el (Q 2 ) (i.e.d → ∞).Appropriate modifications to include the loss of nuclear coherence and incoherent scattering could be easily added if so desired (this would be important, for instance, if one wanted to consider γZ → µ + µ − Z).

C.1 Energy loss
In PETITE we have not attempted a stochastic treatment of energy loss. 4We instead treat energy loss with a deterministic function ⟨dE/dx⟩ which includes contributions from soft bremsstrahlung (photon emissions below the cut E min -see Sec.2.2) and ionization.Because we take E min ∼ MeV, well below the critical energy for most materials O(10 − 100 MeV) [54], ⟨dE/dx⟩ is approximately an energy-independent constant since hard bremsstrahlung (photon emissions above E min ) is explicitly simulated as a discrete process.
The propagation of a particle between its birth and hard interaction takes place in a series of steps, with each step much smaller than the mean free path.After each step the energy is adjusted using the average energy loss; if the energy falls below E min , the particle is terminated and will not be propagated anymore.

C.2 Multiple Coulomb scattering
MCS is crucial for the proper modelling of shower dynamics, especially its transverse distribution.PETITE has two different implementations of MCS available, and the ability to tune the typical scattering angle from MCS.In what follows all formulae are written for particles with unit charge.In both the Lynch-Dahl and Bethe-Molière formalisms [52,53], rescaling the outputted MCS angle can be interpreted as re-scaling the characteristic scattering angle χ c .
Bethe-Molière The slower implementation in PETITE is based on the theory of multiple scattering first developed by Molière [79], and then simplified and improved by Bethe [52].In practice, MCS in this formalism is implemented using an asymptotic expansion in a parameter B, defined by the transcendental equation where with the characteristic scattering, χ c , and screening, χ a , angles being defined below Eq. (4.1).Here t is the thickness of material traversed in g/cm 2 , β the particle velocity, Z the atomic number and A the atomic weight in g/mol.Ω is the mean number of scatters, so typically Ω ≫ 1 unless t < 10 −3 X 0 [53].
The Molière distribution is given by The probability distribution function for scattering through an angle θ is then where ϑ = θ/ Bχ 2 c .Note that f i for i > 0 correspond to non-Gaussian corrections that come from rare, relatively large-angle deflections.Because of this non-Gaussianity, numerical sampling of Bethe-Molière distribution is relatively slow.
Lynch-Dahl (Gaussian) In many applications a Gaussian approximation to MCS can be used, enabling very fast sampling.The default settings in PETITE use the Lynch-Dahl parametrization where the angles are sampled from a normal distribution with variance [53], where Ω is given in Eq. (C.2) and we take F = 0.98 in PETITE.

D Sampling and VEGAS
For the sake of integration and event sampling of SM bremsstrahlung and pair production, we perform a change of coordinates relative to d Π (defined below Eq. (A.7) and Eq.(A.9)).For bremsstrahlung, the new coordinates are defined in Table 1.Since the integrand in Eq. (A.5) is singular for small q 2 and scales like 1/ω, it peaks when δ ≈ δ ′ and when δ, δ ′ , ω → 0. Also, by converting to dimensionless integration variables, we are more readily able to recycle trained VEGAS integration maps for new incident electron/positron energies.
Table 1: Integration variables x i used for training VEGAS on SM bremsstrahlung emission for electrons and positrons.The right column provides the mapping between these and the kinematics variables used in Eq. (A.5).
VEGAS Integration Variable -SM Bremsstrahlung Mapping to Eq. (A.5) We also perform a parameter redefinition for pair production to improve the efficiency and accuracy of event sampling, with the mapping provided in Table 2.This has the same advantages as in the bremsstrahlung case due to the q −4 nature of Eq. (A.8).
Finally, the inclusion of atomic screening tames the low-q 2 Coulomb singularity.This renders pairproduction and bremsstrahlung numerically stable.The inclusion of atomic screening is also physical, and modifies the cross section and average angles of emission for both processes.
Table 2: Integration variables x i used for training VEGAS on SM pair production.The right column provides the mapping between these and the kinematics variables used in Eq. (A.5).

VEGAS Integration Variable -SM Pair Production
Mapping to Eq. (A.8)

E Validation of SM showers in PETITE
In this section we compare SM showers generated with PETITE to predictions from analytic theory [30] and from GEANT-4.

E.1 Comparison with analytic shower theory
The analytic shower theory developed by Tsai and Whitis in Ref. [30] included only bremsstrahlung and pair production with approximate differential rates in order to obtain tractable integro-differential equations for particle intensities as a function of depth and energy.These equations were solved iteratively for sequential generations of particles starting from an initial electron beam.
In Fig. 7 we turn off all processes except for bremsstrahlung and pair production in PETITE and compare track lengths to analytic theory generation by generation, for an initial 10 GeV electron on a graphite target.The track length for particle type i = {e − , e + , γ} is where t is the depth inside the target in radiation lengths, I i (t, E) is the differential particle intensity, i.e. dEI i (t, E) is the total current passing through a surface at depth t.The track length is a useful quantity in practice because the total yield of dark sector particles is proportional to dET i (E)σ(E).We can therefore think of dET i (E) as the cumulative distance travelled by all particles of type i with energy in the interval [E, E + dE] in the entire cascade (i.e., it takes into account not only the actual path length of individual particles, but also their multiplicity in the shower).The track length is extracted numerically from PETITE following the prescription from Ref. [100]: where I i (t j , E) is the "measured" intensity at depth t j (the depths at which the intensities are calculated are spaced by ∆t).We used 10 4 equally-spaced t j 's but the numerical results are stable to variations of this number by a factor of ∼10.
The main features of Fig. 7 are that the analytic theory provides an excellent description of 0th generation electrons (the beam electrons that have undergone any number of bremsstrahlung events) and 1st generation   [30].All processes except for bremsstrahlung and pair production have been turned off in PETITE to closely match the assumptions of Ref. [30].The PETITE and analytic theory curves are labelled by the generation number, which counts the number of parent photons a particle has had: generation 0 corresponds to the primary beam that lost energy due to bremsstrahlung, generation 1 includes photons produced from the primary beam and e ± they generate.The analytic theory describes well the straggling of the primary beam and first generation photons, but it systematically over-predicts e ± pair production at lower energies.The bars on the PETITE results are statistical uncertainties for 10 3 showers.photons (the photons produced in those same bremsstrahlung events).However, it significantly over-predicts first generation electrons and positrons (e ± that are pair produced from 1st generation photons).We have tracked this issue to the "complete screening" approximation employed in Ref. [30].This approximation replaces the energy-dependence of the cross-section σ(E), by its asymptotic value at large energies, σ asym , which is absorbed into the definition of the radiation length.Tsai and Whitis scale all their length scales to radiation lengths, whereas in PETITE each step is determined using the full energy-dependent mean free path λ(E).Since σ(E) ≤ σ asym we have λ(E) ≥ λ asym .This results in a suppression of the normalized track length distribution T (E)/X 0 that is energy dependent.The over-prediction of secondary e ± pairs compared to simulations is also evident in previous work, see, e.g., Fig. 4 of Ref. [100].
It is also interesting to consider the impact of all other processes that are present in PETITE, but not in Ref. [30]: these include annihilations, Compton, Møller and Bhabha scattering, ionization energy loss and all photon emissions (without generation restrictions).In Fig. 8 we compare PETITE with all of these processes ("full PETITE") enabled to the analytic theory.We see that the additional physics does not qualitatively impact the conclusions of previous paragraph; the main difference is that because PETITE computes all γ emissions, the photon track length is larger than that predicted by Ref. [30] (this is already evident in the right panel of Fig. 7).

E.2 Comparison with GEANT-4
In many experiments, dark sector signal yield is maximized when the parent SM particle has a small transverse momentum and a large energy.Thus, modelling this part of the SM cascade is crucial.In Fig. 9 we show a comparison of energy-weighted shower distributions (in transverse and longitudinal distances r T and z) between PETITE and GEANT-4, taking an incoming beam of 10 GeV photons in a graphite target.The statistical significance of differences between the two distributions is estimated by computing a signed χ 2 for 0.1 average particle energies (E) in each r T and z bin: where σ 2 E is the variance of E. The particle number in each bin is sensitive to the minimum tracking energy E min ; smaller tracking energies imply more low-energy particles and greater sensitivity to low-energy physics.Therefore in Fig. 9 we show two different choices of E min .The trend is that earlier parts of the showers (smaller r T and z which are usually the most relevant for dark sector production) are consistent between the two programs since generally χ 2 ≪ 1 in each bin.Statistically-significant differences emerge at large transverse and longitudinal distances when we lower E min .We have identified the likely culprit of these differences to be modelling of multiple Coulomb Scattering; GEANT-4 implements separate MCS models for low-and high-energy particles, and it always includes rare high-deflection events, whereas PETITE uses a single model for MCS.We have checked, for example, that turning off MCS completely in both simulations (an unphysical limit) or using the Bethe model in PETITE (which includes the rare high-angle scatterings but is noticeably slower than a simple Gaussian model) improves agreement.We emphasize that for most experimental set-ups that feature a distant, on-axis detector these differences are immaterial, since the signal events will be generated by the earlier, energetic and forward-going part of the cascade.The PETITE and GEANT-4 distributions in this region of phase space are in very good agreement.

F Dark vector bremsstrahlung
The mass of the dark sector particles qualitatively changes the energy and angular distribution of its emission compared to its SM counterpart.These important features are determined by the soft and collinear singularity structure of the amplitude.This means that accurate simulation of this interaction can be numerically challenging despite the simplicity of the underlying process.In this Appendix we validate our implementation of dark bremsstrahlung and compare it to existing tools.
Neglecting mass-suppressed diagrams where the dark vector couples to the nucleons, the tree-level amplitude for e(p) + N (P ) → e(p ′ ) + N where q = P ′ − P and Q 2 = −q 2 .Spin-summing and averaging the squared the matrix element, one obtains an integral expression for the differential cross section [70,71] Here I ϕ is the result of integrating the matrix element squared over the azimuthal angle of the exchanged momentum ⃗ q and E 0 is the initial electron's energy in the nucleus rest frame.The full expressions for I ϕ , Q min,max are presented in Ref. [70,71].
The dark sector and beam particle particle masses, m V and m e , regulate the soft and collinear divergences, which can be seen explicitly in the kinematic regime with q 2 ≈ −Q 2 min [61], where the minimum photon virtuality is This approximate form is valid when E ≫ m V .This exercise also shows that for m V > m e , the amplitude is enhanced for x ≈ 1 as stressed in Ref. [61].MeV.We show MadGraph, MadGraph with settings that improve numerical stability, CalcHEP (which achieves good stability with default settings), the treatment in PETITE, and the Weizsacker-Williams approximation commonly employed in the dark photon literature [61].
The q 2 ≈ −Q 2 min regime mentioned above is relevant because the photon propagator ∼ 1/q 2 contributes another apparent singularity at small Q 2 .This divergence is regulated in two ways: by the kinematic cutoff Q 2 min , and by the atomic and nuclear form-factor 2) and Eq.(B.4)).The latter effect is due to the nuclear charge being screened by the atomic electrons in this low momentum transfer limit.The kinematic cutoff is dominant when m V ≳ 3 MeVZ 1/6 E/GeV, which is the case for most of the parameter space we are interested in (especially if experimental detection places a ∼ GeV energy threshold on signal events, leading to a similar lower bound on the initial electron energy E).At large energies and small masses, however, the form factor suppression is relevant.
To summarize, the bremsstrahlung amplitude has apparent singularities that are regulated by particle masses or atomic screening.This leads to a highly peaked differential cross-section that can be difficult to sample from.This is particularly challenging if the peaks occur in nontrivial surfaces in phase space (i.e., if they do not align with the integration variables).We find that the adaptive algorithm of VEGAS [72,73] is able to correctly integrate and sample from the partially integrated differential cross-section from Refs.[69][70][71].We validate our implementation by comparing the total cross-section, energy and angular distributions of the emitted dark sector states to other calculations using CalcHEP [102], MadGraph [103] and the Weizsäcker-Williams approximation [69].
In Fig. 10 we compare the total cross section as a function of beam energy for a fixed m V = 0.1 GeV and a carbon target evaluated using different methods.We see that our implementation is in agreement with existing MC tools, CalcHEP and MadGraph.We note that the default MadGraph settings (as of version 3.4.0)are inappropriate for this particular process with no cuts on the emitted particles.This issue has been noted in, e.g., Refs.[43,104].In the former, the atomic form factor was set to 0 for very low photon virtualities to avoid this problem.The issue seems to be in the channel decomposition automatically performed by MadGraph [105], which by default looks at the (singular) propagator 1/q 2 instead of the (regular) combination F (Q 2 )/q 2 .We found that this results in a biased event sample, leading to inconsistent results in terms of cross-section and distributions.This problem is completely addressed by changing the single diagram enhancement strategy to sde_strat 1 and t channel integration strategy to t_strat 2 (see Ref. [105] for a detailed description of these options).Both CalcHEP and our implementation simply use VEGAS on the entire differential cross-section which avoids these complications.In Fig. 11 we show that angular and energy distributions of the emitted vector particles are in perfect agreement between PETITE and MadGraph (with sde_strat 1 and t_strat 2) within MC-statistical uncertainties.

G Dark Compton scattering and resonant annihilation
Dark Compton scattering involves an incident photon which scatters on an electron and produces an electron and dark vector, eγ → eV .The differential cross section for dark Compton scattering is closely related to that for Compton scattering, Eq. (A.For the dark annihilation differential cross section we retained m e in the denominator which regulates collinear divergences.Note that we do not use this expression in PETITE; rather we implemented the results of a resummation of all soft and collinear SM emissions as described in Appendix H. Since the phase space of these processes is one-dimensional, sampling of the final state kinematics is fast.This is not so for dark bremsstrahlung, which is discussed in detail in Appendix F.

H Positron annihilation to dark vectors
The importance of secondary positron annihilation for the production of dark sector states in thick targets was first pointed out in Ref. [63] in the context of axions and in Ref. [100] for dark photons.At tree level, there are two processes e + e − → V and e + e − → γV that can be called annihilation (termed "resonant" and "non-resonant").The distinction between them is not obvious since in the second process, the SM photon can be arbitrarily soft.The authors of Ref. [100] avoid this issue by putting an (arbitrary) minimum energy cut on the SM photon.Additionally, the process e + e − → γV contains a collinear singularity (regulated by m e ) when the emitted photon is aligned with e ± ; this leads to the presence of large logarithms in the rate, which make perturbative calculations worse.To address these issues, we consider an alternative approach to simulating these processes.The emission of an arbitrary number of soft or collinear photons can be accounted for by using electron and positron "parton distribution functions" (PDFs), f e (x, Q).The PDFs calculate the probability for an electron/positron to carry a fraction x of longitudinal momentum of the beam particle when probed at scale Q.The largest annihilation rate occurs when the soft and collinear photon emissions bring the centre-of-mass energy of the colliding "partons" on resonance with the mass of the dark sector state: x + x − s = m 2 V , where x ± are the longitudinal momentum fractions of e ± .The cross-section for this radiative return is given by where β 2 f = 1 − 4m 2 e /m 2 V and we used the narrow width approximation.The PDFs are calculated by solving the Gribov-Lipatov equation, which resums virtual, soft, and collinear emissions.A useful compilation of various approximations for f e (x) can be found in Ref. [75].We will use the Kuraev-Fadin [74]  where β = (2α/π)(ln s/m 2 e −1) and the ellipsis stands for higher order terms in β.Here we see the impact of the collinear singularities in factors of ln s/m 2 e , and the soft divergence as (1 − x) 1−β/2 .For the typical GeV-scale energies log s/m 2 e is O (10).The resummation therefore has important effects both at large incoming particle energies (where these logarithms grow), and also near threshold s ≈ m 2 V .In particular, evaluating Eq. (H.1) as s → m 2 V , one finds that σ rr ∝ (s − m 2 V ) −1+β , compared to (s − m 2 V ) −1 in the tree-level approximation for e + e − → γV .In other words, the soft photon singularity becomes integrable after resummation.This is important when the cross-section is integrated with a distribution of incoming positron energies, as effectively happens in a EM shower in a beam dump.The resummed result gives a finite answer, while the tree-level cross-section would give infinity (without a minimum energy cut).
In order to simulate the radiative return production of dark sector states, we sample the parton luminosity function (the integrand in Eq. (H.1)) to generate momentum fractions x and m 2 V /(xs) carried by leptons in the CM frame.The direction of the outgoing V is taken to be in the same direction of the incoming e + and its energy is approximately m 2 V /(2xm e ).While this adequately models the dominant part of the phase space (the one experiencing a collinear enhancement), one may be interested in events where a photon is emitted with a significant transverse momentum.Here the tree-level calculation would be appropriate since one can use it to explicitly generate the outgoing photon kinematics.These two parts of phase space can be merged by selecting cuts on the maximum p T = Q in the PDF calculation to correspond to the minimum p T in the tree-level calculation.For simplicity, we only focus on the collinearly-enhanced phase space in this work.In Fig. 12 we compare radiative return and tree-level cross section calculations.We find that radiative return modifies the cross section at the level of 10 − 20%.fraction in the vicinity of the resonance because emission of soft-collinear photons take the parton below the resonance threshold.The cross section is enhanced for the same reason at higher √ s, because radiative return allows the parton to access the resonance.

3 D 4 F
Dark sector fluxes from electromagnetic cascades 3.1 Beyond the Standard Model fluxes from electromagnetic cascades 3.2 Generating dark vectors from a Shower using PETITE 3.3 Algorithm for long cascades 3.4 Comparison to previous work 4 Key processes and their impact on phenomenology 4.1 Multiple scattering and downstream fluxes 4.2 Resonant positron annihilation 4.3 Standard Model bremsstrahlung 4.4 Dark vs. standard bremsstrahlung 4.5 Impact on experimental yields Sampling and VEGAS E Validation of SM showers in PETITE E.1 Comparison with analytic shower theory E.2 Comparison with GEANT-Dark vector bremsstrahlung G Dark Compton scattering and resonant annihilation H Positron annihilation to dark vectors 1 Introduction

Figure 1 :
Figure 1: The projection in x − z (top) and y − z (bottom), where z is along the beam axis, of a typical shower development as simulated using PETITE.In this case the shower is initiated by a 10 GeV photon incident on a thick graphite target.Photons are shown as green wavy lines, electrons as blue straight lines, and positrons as red straight lines.Particles have been tracked until their energy drops below 3 MeV.

Figure 2 :
Figure 2: Impact of MCS on BSM fluxes for a 10 GeV incident photon, on graphite, producing dark vectors with mass m V = 5 MeV: (Left) Distribution of p T /p z for dark bremsstrahlung from positrons and electrons, and dark annihilation of positrons on atomic electrons.Distributions are shown for two different MCS implementations.A vertical blue line shows a rough angular acceptance for the DUNE near detector (2.5 m/574 m = 4.3 mrad).Typical photon energies from π 0 decays from a 120 GeV proton beam range from 1-30 GeV.(Right) The same distribution below the cut p T /p Z ≤ 2.5/574.Results are obtained from 3000 independent cascades simulated with PETITE.The annihilation probability is roughly four orders of magnitude higher than the probability of dark bremsstrahlung and the curves have been re-scaled so as to fit on the same axes.

Figure 3 :
Figure3: Comparison of p T /p z (defined relative to the beam axis) for two different sampling procedures for positron annihilation for m V = 3 MeV (left) and m V = 6.2 MeV (right).The green curve take the dark vector direction, pV , as the initial positron's direction, pV = pi , when it was created in the shower.The orange curve implements the direction of V as the direction of the positron accounting for propagation and energy loss, pV = psamp .The left and right dashed vertical lines are representative of angular acceptances at ArgoNeuT (0.2 mrad) and DUNE (4.3 mrad), respectively.

1 Figure 4 :
Figure 4: Transverse momentum distribution for daughter electrons and photons from bremsstrahlung produced by the first splitting of a E e = 10 GeV incident electron.Distributions are computed with PETITE (left), GEANT-4 (middle) and MadGraph with custom settings (right), see Appendix F. The implementation in PETITE reproduces the analytic expressions in [55].The artificially tight correlation between p T of the daughter γ and daughter e in GEANT-4 is readily seen, as is the lack of phase space coverage in MadGraph.

1 Figure 5 :
Figure 5: Comparison of reweighting vs resampling as a method for producing dark vector fluxes.The SM bremsstrahlung distribution from a 10.6 GeV e − is shown (Left), the dark vector distributions produced by reweighting the SM distribution using Eq.(4.3) (Centre), and the distribution of dark vectors produced via PETITE's resampling procedure (Right).A dashed line is drawn at p T /p Z = 0.5/10 to illustrate the BDX detector's angular acceptance.Each panel's colourscale is normalized to its distribution's highest-flux bin.

Figure 6 :
Figure6: Total yields of dark photons from electromagnetic secondary and hadronic primary production for three different experimental configurations.(Left) The flux emanating from the Booster Neutrino Beam with parent particles satisfying a cut of sin θ beam < 0.2 on their angle relative to the beam axis.This configuration does not correspond to a specific detector and was chosen for direct comparison with Ref.[43].Flux of dark vectors passing through the MiniBooNE (Centre) and DUNE (Right) detectors, which have been chosen for direct comparison of PETITE with Ref.[44].

PETITE
Track Length vs Analytics (Bremsstrahlung and Pair Production Only)

Figure 7 :
Figure7: Track length distributions for e − , e + and γ for an initial 10 GeV e − beam on a thick graphite target.Each panel compares predictions of PETITE to the analytic calculations of Tsai and Whitis[30].All processes except for bremsstrahlung and pair production have been turned off in PETITE to closely match the assumptions of Ref.[30].The PETITE and analytic theory curves are labelled by the generation number, which counts the number of parent photons a particle has had: generation 0 corresponds to the primary beam that lost energy due to bremsstrahlung, generation 1 includes photons produced from the primary beam and e ± they generate.The analytic theory describes well the straggling of the primary beam and first generation photons, but it systematically over-predicts e ± pair production at lower energies.The bars on the PETITE results are statistical uncertainties for 10 3 showers.

Figure 10 :
Figure10: Total cross section for dark vector bremsstrahlung of electrons on a carbon nucleus for m V = 100 MeV.We show MadGraph, MadGraph with settings that improve numerical stability, CalcHEP (which achieves good stability with default settings), the treatment in PETITE, and the Weizsacker-Williams approximation commonly employed in the dark photon literature[61].

Figure 11 :
Figure 11: Energy (left) and angle (right) distributions of the emitted vector particle in e − +N → e − +N +V with an beam energy of 10 GeV and m V = 0.1 GeV.