11institutetext: Department of Physics and Astronomy, University of Victoria,
Victoria, BC V8P 5C2, Canada22institutetext: Department of Physics and
Astronomy, York University, Toronto, Ontario, M3J 1P3, Canada33institutetext:
Particle Theory Department, Fermi National Accelerator Laboratory, Batavia, IL
60510, USA44institutetext: Theoretical Physics Department, CERN, Esplande des
Particules, 1211 Geneva 23, Switzerland55institutetext: Department of Physics
and Astronomy, Mitchell Institute for Fundamental Physics and Astronomy, Texas
A&M University, College Station, TX 77843, USA66institutetext: Walter Burke
Institute for Theoretical Physics, California Institute of Technology,
Pasadena, CA 91125, USA

# Dark fluxes from electromagnetic cascades

Nikita Blinov 3 Patrick J. Fox 4,5 Kevin J. Kelly 3 Pedro A.N. Machado 6 Ryan
Plestid<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

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, 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.

††preprint: CALT-TH/2024-001, FERMILAB-PUB-24-0003-T, CERN-TH-2024-006, MI-
HET-825

##  1 Introduction

General ultraviolet (UV) completions of the Standard Model (SM) Craig:2022cef
, the existence of dark matter Bertone:2004pz ; Battaglieri:2017aum , and
neutrino masses Kajita:2016cak ; McDonald:2016ixn motivate the existence of
feebly interacting particles Beacham:2019nyx ; Antel:2023hkf . If the mass
scale of these particles is below the electroweak scale, they must be singlets
under the SM gauge group L3:2001xsz ; CMS:2019hsm ; ATLAS:2023sbu ,
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 Beacham:2019nyx ; Gori:2022vri ;
Harris:2022vnx .

Intensity frontier experiments with thin or thick targets Berlin:2018bsc ;
NA64:2023wbi ; deNiverville:2011it ; Batell:2014mga ; SHiP:2015vad ;
BDX:2016akw ; deGouvea:2018cfv ; DeRomeri:2019kic ; Machado:2019oxb ;
Harnik:2019zee ; ArgoNeuT:2019ckq ; DUNE:2020fgq ; Breitbach:2021gvv ;
Apyan:2022tsd ; Coloma:2023adi often provide the most sensitive probes of dark
sectors, especially for particle masses below a few GeV Beacham:2019nyx ;
Gori:2022vri . 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
$\mathcal{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 Bhabha:1937xxx ; Gaisser:1990vg . 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 Tsai:1966js . 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
deNiverville:2011it ; deNiverville:2016rqh ; Magill:2018tbb ; Magill:2018jla ;
Berryman:2019dme ; Plestid:2020kdm . 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_{\mu}$ or
$L_{e}-L_{\tau}$, 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 Tsai:1986tx . Conversely, if a dark
photon couples democratically to charged leptons and hadrons, it has large
production rates from meson decays Batell:2009di ; Reece:2009un or proton
bremsstrahlung Blumlein:2013cua ; Foroughi-Abari:2021zbm . 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 Buonocore:2018xjk ; Nardi:2018cxi ; Celentano:2020vtu ;
Capozzi:2021nmp , 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:

  1. 1.

To provide a (SM) generator-independent solution for producing a dark sector
flux from a SM cascade.

  2. 2.

To describe and implement a principled, flexible, and systematically
improvable framework for producing a Monte Carlo (MC) event record of dark
particles emanating from EM showers.

  3. 3.

To compare existing results on BSM particles produced in EM showers, and to
identify which methods are reliable and to diagnose the cause of any
disagreements.

EM showers are largely driven by a small set of reactions all of which are
perturbatively calculable and systematically improvable if needed. We adopt
the phenomenologically-motivated assumption that the impact of dark sector
emissions on the dynamics of EM showers is negligible. This allows us to
“factorize” the SM 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 beam-dump or other fixed target experiment with ${\sim}1{-}100$ GeV
beams.111Strictly speaking the approximations we use are valid for electron,
positron, and photon energies greater than $\sim 1~{}{\rm MeV}$, and less than
$100~{}{\rm GeV}$ where the Landau-Pomeranchuk-Migdal effect substantially
modifies bremsstrahlung cross sections Landau:1953um ; Migdal:1956tc . 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 GEANT4:2002zbu ;
Allison:2006ve ; Allison:2016lfl , 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^{-}\rightarrow\gamma V$ in Appendix G. Our
implementation of positron annihilation into dark vectors includes soft- and
collinear photon resummation which is discussed in Appendix H.

##  2 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^{\pm}$ and $\gamma$, but muons can be included straightforwardly.

###  2.1 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. 1.

Bremsstrahlung on atomic nuclei: $e^{\pm}Z\rightarrow e^{\pm}Z\gamma$. This
process is dominated by small momentum transfers, $Q$, that satisfy $QR\ll 1$
where $R$ is the nuclear radius. We treat the nucleus as a static fixed
Coulomb potential, use the Bethe-Heitler matrix element Bethe:1934za , and
include form factors to account for screening and finite nuclear size effects
at small and large momentum transfers respectively Tsai:1973py , see Appendix
B. A minimum energy cut (whose default value is $\omega_{\rm
cut}=1\,\mathrm{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 $\omega_{\rm cut}$.

  2. 2.

Bethe-Heitler pair production: $\gamma Z\rightarrow 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. 3.

Scattering on atomic electrons : $(\gamma,e^{\pm})~{}e^{-}_{\rm
bound}\rightarrow(\gamma,e^{\pm})e^{-}$ and $e^{+}e^{-}_{\rm
bound}\rightarrow\gamma\gamma$. We use tree-level expressions for Bhabha,
Møller, and Compton scattering as well as for $e^{+}e^{-}$ annihilation, see
Eqs. A.7 to A.10. 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_{\rm min}$ of the
outgoing $e^{\pm}$.

In addition to these reactions, there are two processes that are much more
common when $e^{\pm},\gamma$ 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. 1.

Coulomb scattering on nuclei: $e^{\pm}Z\rightarrow e^{\pm}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 Bethe:1953va ; Lynch:1990sq . 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 Sec. C.2) and highlight the importance
of MCS in positron annihilation in Sec. 4.2.

  2. 2.

Ionization: $e^{\pm}e^{-}_{\rm bound}\rightarrow e^{\pm}e^{-}_{\rm free}$.
Atomic ionization dominates energy loss for electrons at energies below the
critical energy, $\mathcal{O}(10-100\;\mathrm{MeV})$ depending on $Z$
ParticleDataGroup:2022pth . 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 $\differential E/\differential x=2~{}{\rm MeV}~{}{\rm
cm}^{2}/{\rm g}$, i.e., independent of energy (this number is only weakly
dependent on the type of atom ParticleDataGroup:2022pth ).

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 $\omega_{\mathrm{cut}}$. However, for the default value of the
cut that separates hard from soft, $\omega_{\mathrm{cut}}=1\,\mathrm{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 Berestetskii:1982qgu and also in Appendix A.

###  2.2 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_{\rm min}\geq\omega_{\rm 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_{\rm 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_{\rm step}$ and
whether a discrete/hard interaction takes place in that step; the probability
of such an interaction is $P=1-{\rm e}^{-x_{\rm step}/\lambda(E)}$ with
$\lambda(E)$ the mean free path calculated using the tabulated cross-sections
stored in Shower.222The step size $x_{\rm step}\approx\lambda(E)/10$ is chosen
small enough such that $\lambda(E)$ does not change significantly over the
step, and large enough such that the probability of interaction is not too
small. The four-momentum and position of the Particle are updated at each
step, taking into account multiple scattering and energy loss (see above and
Sec. 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_{\rm 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_{\rm
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\rightarrow e^{+}\gamma Z$ the event record will contain the initial
parent positron $e^{+}$, the daughter positron $e^{+}$, and the daughter
photon $\gamma$. 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.

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~{}{\rm 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.

##  3 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 Bellamy:1967qf , heavy leptons Barna:1968su and neutrino-like
particles SLAC:E56 ; Rothenberg:1972pt . In these experimental works, analytic
approximations of the secondary spectra of photons and electrons, due to Tsai
and Whitis Tsai:1966js , were used to estimate signal yields in thick targets.
Later Tsai’s analytic results were applied to axion bremsstrahlung from an
attenuated electron beam Donnelly:1978ty ; Tsai:1986tx . This method has been
adopted to reinterpret old experimental results in the context of other dark
sector particles like the dark photon Bjorken:2009mm ; Andreas:2012mt . A
different signal prediction technique was used in analyzing the electron beam
dump E137 data Bjorken:1988as : 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.
Marsicano:2018vin ; Rella:2022len the secondary production of muons is
considered since the focus of those works is on muonphillic interactions. More
recently in Refs. Nardi:2018cxi ; Celentano:2020vtu , it was realized that for
generic dark sector couplings, resonant $e^{+}e^{-}\rightarrow X$ production
can be extremely efficient for $m_{X}\lesssim 50~{}{\rm 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. Capozzi:2021nmp . Secondary fluxes of photons
have also been demonstrated to be of interest in the context of axion like
particles Tsai:1986tx .

Recently, several approaches for beam dump dark sector production have been
developed. In Ref. Capozzi:2021nmp 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. Bondi:2021nfp ; Eichlersmith:2022bit . In these works, SM cascade and
dark sector production are simulated simultaneously by using an unphysically
large dark sector coupling. As emphasized in Ref. Celentano:2020vtu , 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. Celentano:2020vtu have implemented a solution to
the above problems by separately simulating the SM cascade in GEANT-4 and BSM
production in MadDump Buonocore:2018xjk or BdNMC deNiverville:2016rqh . 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:

$\mathcal{L}\supset
V_{\mu}\left(c_{e}\bar{e}\gamma^{\mu}e+c_{p}\bar{p}\gamma^{\mu}p+c_{n}\bar{n}\gamma^{\mu}n\right)+c_{\pi}\frac{e}{32\pi^{2}f_{\pi}}\pi^{0}\epsilon^{\mu\nu\alpha\beta}F_{\mu\nu}\partial_{\alpha}V_{\beta}+\ldots$
(3.1)

In this work we will primarily focus on the case where the coupling to
electrons is the dominant coupling, $c_{e}\equiv\varepsilon e$ and
$c_{p,n,\pi}\approx 0$. This can occur for electrophilic gauge extensions of
the SM, e.g., $L_{e}-L_{\mu}$ or $L_{e}-L_{\tau}$, and EM shower-induced
production dominates. Similarly, in certain axion models secondary electrons
dominate axion production even in a proton beam dump Tsai:1986tx , although
through a different Lorentz structure. Alternatively, a dark photon couples
democratically to charged leptons and hadrons
($|c_{p}|=|c_{e}|=c_{\pi}=\varepsilon e$), giving large production rates from
meson decays Batell:2009di ; Reece:2009un or proton bremsstrahlung
Blumlein:2013cua ; Foroughi-Abari:2021zbm . We will not explicitly discuss
hadronic production modes and refer the interested reader to the relevant
literature, e.g., deNiverville:2011it ; deNiverville:2016rqh ; Magill:2018tbb
; Magill:2018jla ; Berryman:2019dme ; Plestid:2020kdm .

###  3.1 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. 1.

Dark bremsstrahlung: $e^{\pm}Z\rightarrow e^{\pm}ZV$. In the static limit
(i.e., the nucleus is taken as infinitely heavy), this process is
kinematically allowed for $E_{e}\geq 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}\gg m_{V}$. In turn this
means dark bremsstrahlung dominantly produces hard dark vectors Bjorken:2009mm
; Liu:2016mqv ; Liu:2017htz ; Tsai:1989vw ; Gninenko:2017yus . 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 Lepage:1977sw ; Lepage:2020tgj to sample the differential cross-
sections presented in Tsai:1989vw ; Gninenko:2017yus .

  2. 2.

Dark Compton scattering: $\gamma e^{-}_{\rm bound}\rightarrow Ve^{-}_{\rm
free}$. This process is operational provided $E_{\gamma}\geq
m_{V}^{2}/(2m_{e})$. The corresponding differential cross section can be found
in Appendix G.

  3. 3.

Positron annihilation: $e^{+}e^{-}_{\rm bound}\rightarrow V(n\gamma)$. This
channel can be the largest contribution to the dark vector flux when it is
open. The kinematic threshold is set by $s=s_{\rm th}$ with $s_{\rm
th}=m_{V}^{2}$. For $s>s_{\rm 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^{-}\rightarrow V\gamma$ diverges like
$\sigma\propto 1/(s-s_{\rm 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 Celentano:2020vtu . Here we follow the lepton-collider
literature Kuraev:1985hb ; Nicrosini:1986sm ; Karliner:2015tga ; Greco:2016izi
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).

###  3.2 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, $\vec{p}$, and production positions,
$\vec{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 $\differential w^{(p)}$ occurring between $z$ and $z+\differential
z$ is proportional to the probability for the BSM process to occur and to the
probability of no SM process happening up to that point:

$\frac{\differential w^{(p)}}{\differential z}=n^{(p)}\sigma^{(p)}_{\rm
BSM}(z)~{}\exp\left(-\int_{0}^{z}\frac{\differential z^{\prime}}{\lambda_{\rm
MFP}^{\rm SM}(z)}\right),$ (3.2)

where $z=0$ corresponds to the SM particle’s starting position $\vec{x}$. The
cross section $\sigma^{(p)}_{\rm BSM}$ and mean free path $\lambda_{\rm
MFP}^{\rm SM}$ 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., $\lambda_{\rm MFP}\approx\lambda_{\rm MFP}^{\rm
SM}$. 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)}=\int\differential z\;\differential w^{(p)}/\differential z$. 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^{-}\to
V(n\gamma)$, is highly energy dependent with a peak at $s=m_{V}^{2}$. 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^{-}\rightarrow V(n\gamma)$ 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=\\{\texttt{Dark\ bremsstrahlung, Dark annihilation}\\}$

else if $particle=e^{-}$ then

$processlist=\\{\texttt{Dark\ bremsstrahlung}\\}$

else if $particle=\gamma$ then

$processlist=\\{\texttt{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
$\differential\sigma_{\rm BSM}^{(p)}$.

• compute weight $w^{(p)}$ by integrating Eq. 3.2.

• 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.

###  3.3 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
$\mathcal{O}({\rm 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({\rm
DS}|e^{-})\sim\mathcal{O}(\varepsilon^{2})$. If the typical number of SM
daughters in a shower is $\langle N_{\rm SM}\rangle$ and

$\varepsilon^{2}\times\langle N_{\rm SM}\rangle\sim\mathcal{O}(1),$ (3.3)

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 $\langle N_{D}\rangle$ 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 $\mathcal{O}({\rm
few})$ dark sector events are produced after which the SM simulation is set
back to its initial configuration and the entire procedure repeated.

###  3.4 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. Bondi:2021nfp ; Eichlersmith:2022bit we can re-
purpose the same SM shower for many different BSM parameter choices allowing
for efficient parameter scans. Third, in contrast to Ref. Capozzi:2021nmp 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. Celentano:2020vtu 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. Celentano:2020vtu in Sec. 4.5.

##  4 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 $\pi^{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 Bierlich:2022pfr , which
provides the primary $\pi^{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.

###  4.1 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(\theta,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 Moliere:1947zza ; Bethe:1953va .
Comparisons of MCS models implemented in GEANT-4 find ${\sim}10\%$ agreement
with empirical measurements Gottschalk:1993xxx ; Makarova:2017xxx . 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
Lynch:1990sq . 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 Moliere:1947zza , which was improved by
Bethe Bethe:1953va . We will summarize Bethe’s results to explain the
different approximations used in computing MCS. The number of electrons in an
angular interval $\differential\theta$ after traversing a material thickness
$t$ is $f(\theta,t)\theta\differential\theta$. Molière theory allows an
asymptotic expansion of $f(\theta,t)$ in a parameter, $B$, which encodes the
depth traversed in a given material. Specifically, $B$ is the solution of the
transcendental equation

$B-\ln B=\ln\frac{\chi_{c}^{2}}{1.167\chi_{a}^{2}},$ (4.4)

where $\chi_{c}^{2}\equiv 4\pi Nt\alpha^{2}Z(Z+1)z^{2}/(p\beta)^{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 $\beta$ the charge, momentum and velocity of the propagating
particle; the screening angle
$\chi_{a}^{2}=(1.13+3.76\tilde{\alpha}^{2}){\mkern 0.75mu\mathchar
22\relax\mkern-9.75mu\lambda}^{2}/a^{2}$ encodes the suppression of Rutherford
scattering, where $\tilde{\alpha}=zZ\alpha/\beta$, and ${\mkern
0.75mu\mathchar 22\relax\mkern-9.75mu\lambda}=1/p$ is the de Broglie
wavelength of the electron and $a=0.885Z^{-1/3}/(\alpha 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.
Bethe:1953va , and in Appendix C. The ratio $\chi_{c}^{2}/\chi_{a}^{2}$ is
approximately the mean number of scatters in a thickness $t$ Bethe:1953va ;
Lynch:1990sq , so $\chi_{c}^{2}/\chi_{a}^{2}\gg 1$ (so
$B\sim\ln\chi_{c}^{2}/\chi_{a}^{2}$) unless $t$ is so small that the Molière
theory would no longer be applicable. The expansion developed by Bethe can be
summarized as

$f(\theta,t)=\frac{1}{\chi_{c}^{2}B}\left[f_{0}(\theta)+B^{-1}f_{1}(\theta)+B^{-2}f_{2}(\theta)+\ldots\right],$
(4.5)

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 Lynch:1990sq . The width of the Gaussian is also characterized
by $\chi_{c}^{2}$ and $\chi_{a}^{2}$ as detailed in Appendix C, with numerical
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 ParticleDataGroup:2022pth .

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 $\pi^{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\;\mathrm{m}/574\;\mathrm{m}=4.3\;\mathrm{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 ${\sim}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~{}{\rm MeV}$ produced via
resonant annihilation changes from $\epsilon_{\rm geo}\approx 0.49$ when using
the Lynch-Dahl implementation to $\epsilon_{\rm geo}\approx 0.43$ when using
the more accurate Bethe-Moliere theory.

Figure 2: Impact of MCS on BSM fluxes for a $10~{}{\rm GeV}$ incident photon,
on graphite, producing dark vectors with mass $m_{V}=5~{}{\rm 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\;\mathrm{m}/574\;\mathrm{m}=4.3~{}\mathrm{mrad})$. Typical photon
energies from $\pi^{0}$ decays from a 120 GeV proton beam range from 1-30 GeV.
(Right) The same distribution below the cut $p_{T}/p_{Z}\leq 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.

###  4.2 Resonant positron annihilation

Positron annihilation to dark sector states via $e^{+}e^{-}\rightarrow
V(n\gamma)$ has been shown to be the most efficient production mechanism
available in EM cascades for light dark sector masses $m_{V}\lesssim
100~{}{\rm MeV}$. This is true in a general EM cascade Nardi:2018cxi , and
even in proton beam dump experiments Celentano:2020vtu ; Capozzi:2021nmp .
Recent results from the NA64 collaboration explicitly make use of resonant
production Andreev:2021fzd ; NA64:2023wbi .

Despite being the dominant production mechanism for certain dark vector
masses, theoretical treatments of $e^{+}e^{-}\rightarrow V(n\gamma)$ 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_{\rm res}=m_{V}^{2}/2m_{e}$, typical
positrons with $E_{e^{+}}>E_{\rm 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
Celentano:2020vtu . 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 $\pi^{0}$ production from $pp$
collisions with PYTHIA and inject them into PETITE, which decays
$\pi^{0}\rightarrow\gamma\gamma$ 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^{-}\to V(n\gamma)$ is given by its
direction when it was created upstream in the shower, $\hat{p}_{i}$, leading
to a dark vector direction $\hat{p}_{V}=\hat{p}_{i}$. 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, $\hat{p}_{V}=\hat{p}_{\rm samp}$. 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~{}{\rm MeV}$, and by a factor of roughly
$15$ for $m_{V}=3~{}{\rm MeV}$.

Figure 3:  Comparison of $p_{T}/p_{z}$ (defined relative to the beam axis) for
two different sampling procedures for positron annihilation for
$m_{V}=3~{}{\rm MeV}$ (left) and $m_{V}=6.2~{}{\rm MeV}$ (right). The green
curve take the dark vector direction, $\hat{p}_{V}$, as the initial positron’s
direction, $\hat{p}_{V}=\hat{p}_{i}$, 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,
$\hat{p}_{V}=\hat{p}_{\mathrm{samp}}$. The left and right dashed vertical
lines are representative of angular acceptances at ArgoNeuT (0.2 mrad) and
DUNE (4.3 mrad), respectively.

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 ${\sim}20\%$,
see Appendix H for further details. Last, we note that we do not include
atomic binding in our current simulations.

###  4.3 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 Hirayama:2005zm , and FLUKA Ferrari:2005zk ; Bohlen:2014buj
, to approximate kinematic samplings by default.333It should be noted that
all three programs provide a more precise treatment of the total but not
differential cross section than PETITE (e.g., by including Coulomb
corrections). PETITE uses the tree-level Bethe-Heitler expressions for both
the differential and total cross section. We 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 Bernard:2018hwf .

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
Hirayama:2005zm and Sec. 10.2 of GEANT_physicsmanual ). FLUKA uses the Berger-
Seltzer parameterization which is two dimensional SELTZER198595 ;
Ferrari:2005zk . 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.

Figure 4: Transverse momentum distribution for daughter electrons and photons
from bremsstrahlung produced by the first splitting of a $E_{e}=10~{}{\rm
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
Berestetskii:1982qgu . The artificially tight correlation between $p_{T}$ of
the daughter $\gamma$ and daughter $e$ in GEANT-4 is readily seen, as is the
lack of phase space coverage in MadGraph.

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.

###  4.4 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^{\pm}Z\to e^{\pm}Z\gamma$, 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 $\theta_{\gamma}\sim m_{e}/E_{e}$, c.f. Eq. A.11). 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^{\pm}$ 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 BDX:2016akw . 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.,
$\pi^{0}\to\gamma 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 $\times$ 50 cm target placed 10 m downstream from the back of the target
BDX:2016akw .

Due to these kinematic differences, one cannot estimate the dark
bremsstrahlung flux by simply re-weighting 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. Capozzi:2021nmp 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_{\gamma}(\theta_{i},E_{j})$, produced through bremsstrahlung in the
$(i,j)$-th bin centred around angle $\theta_{i}$ and energy $E_{j}$ is
converted to the number of dark vectors, $N_{V}(\theta_{i},E_{j})$, in the
same bin using the (angle-independent) ansatz

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.6 (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.

$\begin{split}N_{V}(\theta_{i},E_{j})&=\varepsilon^{2}N_{\gamma}(\theta_{i},E_{j})\times
f\quantity(\tfrac{m_{V}}{E_{e}^{\rm est}(E_{j})})~{},\\\ f(x)&={\rm
min}\quantity(1,1154\exp[-24.32x^{0.3174}])~{},\\\ E_{e}^{\rm
est}(E_{\gamma})&=1.0773E_{\gamma}+13.716~{}{\rm MeV}~{}.\end{split}$ (4.6)

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.
Celentano:2020vtu over-predicts the yield by roughly 5000 for
$m_{V}=100~{}{\rm MeV}$ (we believe this comes from the normalization of
$f(x)$ being off for large $m_{V}$).

###  4.5 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_{\mu}$ to the hadronic electromagnetic
current $J_{\mu}$ with the same coupling as electrons, i.e.
$|c_{p}|=|c_{e}|=c_{\pi}=\varepsilon 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\theta_{\rm beam}<0.2$ from PETITE
and Ref. Celentano:2020vtu in the left panel of Fig. 6. The other two panels
show PETITE and Ref. Capozzi:2021nmp predictions for the dark vector flux at
MiniBooNE (middle panel) and DUNE (right panel).

Figure 6:  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\theta_{\rm 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. Celentano:2020vtu . Flux of dark
vectors passing through the MiniBooNE (Centre) and DUNE (Right) detectors,
which have been chosen for direct comparison of PETITE with Ref.
Capozzi:2021nmp .

The output of PETITE for dark bremsstrahlung agrees well with Ref.
Celentano:2020vtu , as can be see in the left panel. Although our treatment of
MCS in positron annihilation is different from Ref. Celentano:2020vtu 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 ($\eta$, $\eta^{\prime}$) and proton bremsstrahlung,
which we have not accounted for (we only considered $\pi^{0}$ decays). We have
not been able to reproduce the direct hadronic production yield curve of Ref.
Celentano:2020vtu as shown in the left panel (PETITE is roughly $2\times$ 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. Celentano:2020vtu . While secondary meson 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. Capozzi:2021nmp .
As we discussed in Sec. 4.4, the re-weighting procedure for dark
bremsstrahlung production (orange lines) used in Ref. Capozzi:2021nmp
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_{\mu}$ 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. Buonocore:2018xjk ;
Nardi:2018cxi ; Celentano:2020vtu ; Capozzi:2021nmp 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
Capozzi:2021nmp 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_{\mu}$ and other
flavour combinations, and $B-L$.

##  5 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 Ghinescu:2021sjm . 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_{\mu}$. 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., $\mu^{+}\mu^{-}$ 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.

###### Acknowledgements.

We thank Peter LePage for help with the Python implementation of VEGAS. We
thank Adrian Thompson for helpful correspondences regarding the implementation
of dark vector production in Capozzi:2021nmp . We thank Bhaskar Dutta, Wes
Ketchum, Shirley Li, and Stephan Meighen-Berger for useful discussions, and
Brian Batell for feedback on the manuscript. RP is supported by the Neutrino
Theory Network under Award Number DE-AC02-07CHI11359, the U.S. Department of
Energy, Office of Science, Office of High Energy Physics under Award Number
DE-SC0011632, and by the Walter Burke Institute for Theoretical Physics. PJF
and PANM 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 Patel:2015tea and
FeynCalc Shtabovenko:2020gxv to derive and check several analytic results in
this work. Our numerical work was enabled by VEGAS Lepage:1977sw ;
Lepage:2020tgj , numpy harris2020array , scipy 2020SciPy-NMeth and matplotlib
Hunter:2007 .

##  Appendix 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^{\pm}\gamma\rightarrow e^{\pm}\gamma$), electron-positron
annihilation ($e^{-}e^{+}\to\gamma\gamma$), Bhabha scattering
($e^{+}e^{-}\rightarrow e^{+}e^{-}$), and Møller scattering
($e^{-}e^{-}\rightarrow e^{-}e^{-}$):

$\displaystyle\derivative{t}\sigma_{\rm Comp.}$
$\displaystyle=\frac{8\pi\alpha^{2}}{(s-m_{e}^{2})^{2}}\quantity[\quantity(\frac{m_{e}^{2}}{s-m_{e}^{2}}+\frac{m_{e}^{2}}{u-m_{e}^{2}})^{2}+\quantity(\frac{m_{e}^{2}}{s-m_{e}^{2}}+\frac{m_{e}^{2}}{u-m_{e}^{2}})-\frac{1}{4}\quantity(\frac{s-m_{e}^{2}}{u-m_{e}^{2}}+\frac{u-m_{e}^{2}}{s-m_{e}^{2}})]~{},$
(A.7) $\displaystyle\derivative{t}\sigma_{\rm Ann.}$
$\displaystyle=\frac{8\pi\alpha^{2}}{s(s-4m_{e}^{2})}\quantity[\quantity(\frac{m_{e}^{2}}{t-m_{e}^{2}}+\frac{m_{e}^{2}}{u-m_{e}^{2}})^{2}+\quantity(\frac{m_{e}^{2}}{t-m_{e}^{2}}+\frac{m_{e}^{2}}{u-m_{e}^{2}})-\frac{1}{4}\quantity(\frac{t-m_{e}^{2}}{u-m_{e}^{2}}+\frac{u-m_{e}^{2}}{t-m_{e}^{2}})]~{},$
(A.8) $\displaystyle\derivative{t}\sigma_{\rm Bhabha}$
$\displaystyle=\frac{4\pi\alpha^{2}}{u(u-4m_{e}^{2})}\quantity[\frac{1}{t^{2}}\quantity[\tfrac{1}{2}(s^{2}+u^{2})+4m_{e}^{2}(t-m_{e}^{2})]+\frac{2}{tu}\quantity(\tfrac{1}{2}s-m_{e}^{2})(\tfrac{1}{2}s-3m_{e}^{2})+t\leftrightarrow
u]~{},$ (A.9) $\displaystyle\derivative{t}\sigma_{\text{M\o ller}}$
$\displaystyle=\frac{4\pi\alpha^{2}}{s(s-4m_{e}^{2})}\quantity[\frac{1}{t^{2}}\quantity[\tfrac{1}{2}(s^{2}+u^{2})+4m_{e}^{2}(t-m_{e}^{2})]+\frac{2}{tu}\quantity(\tfrac{1}{2}s-m_{e}^{2})(\tfrac{1}{2}s-3m_{e}^{2})+t\leftrightarrow
u]~{}.$ (A.10)

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}\to p_{3}p_{4}$, where
$s=(p_{1}+p_{2})^{2}=(p_{3}+p_{4})^{2}$,
$t=(p_{3}-p_{1})^{2}=(p_{4}-p_{2})^{2}$, and
$u=(p_{3}-p_{2})^{2}=(p_{4}-p_{1})^{2}$.

Bremsstrahlung and pair-production in the field of a nucleus were first
computed by Bethe and Heitler Bethe:1934za ; textbook discussions with modern
notation can be found in Itzykson and Zuber Itzykson:1980rh , and (in the
greatest detail) in Landau and Lifshitz Berestetskii:1982qgu . In the small-
angle regime (which dominates) the differential cross section for
bremsstrahlung process of $e^{\pm}Z\rightarrow e^{\pm}Z\gamma$ is given by
Berestetskii:1982qgu

$\begin{split}\frac{\differential\tilde{\sigma}}{\differential\tilde{\Pi}}=\frac{8}{\pi}\left\lvert
F(Q^{2})\right\rvert^{2}\alpha^{3}\frac{\epsilon^{\prime}m_{e}^{2}}{\omega\epsilon
Q^{4}}\delta\delta^{\prime}\times\bigg{\\{}\frac{\delta^{2}}{(1+\delta^{2})^{2}}+\frac{\delta^{\prime
2}}{(1+\delta^{\prime
2})^{2}}&+\frac{\omega^{2}}{2\epsilon\epsilon^{\prime}}\frac{\delta^{2}+\delta^{\prime
2}}{(1+\delta^{2})(1+\delta^{\prime 2})}\\\ &\hskip
43.36464pt-\quantity(\frac{\epsilon^{\prime}}{\epsilon}+\frac{\epsilon}{\epsilon^{\prime}})\frac{\delta\delta^{\prime}\cos\phi}{(1+\delta^{2})(1+\delta^{\prime
2})}\bigg{\\}}~{}.\end{split}$ (A.11)

Here the photon energy, in the lab frame, is $\omega$ and the incoming
(outgoing) lepton energy, also in the lab frame, is $\epsilon$
($\epsilon^{\prime}=\epsilon-\omega$). We cut off the IR range of the photon
energy, but do not alter the upper bound from the maximum allowed
kinematically, so $\omega_{\rm cut}<\omega<\epsilon-m_{e}$. The variables
$\delta$ ($\delta^{\prime}$) are proxies for $\theta$ ($\theta^{\prime}$)
which are the angles between the outgoing photon and the incoming (outgoing)
lepton. They are defined as

$\delta^{(\prime)}=\frac{\epsilon^{(\prime)}\theta^{(\prime)}}{m_{e}}.$ (A.12)

The angles $\delta$ and $\delta^{\prime}$ span $[0,\pi\times
E_{e}^{(\prime)}/m_{e}]$; in the limit $m_{e}/E_{e}^{(\prime)}\rightarrow 0$
this tends to the positive real axis $[0,\infty)$. The azimuthal angle between
the incoming and outgoing lepton is denoted $\phi$ and lies in the range
$0\leq\phi<2\pi$. The momentum transfer to the nucleus,
$Q^{2}=-q_{\mu}q^{\mu}=\mathbf{q}^{2}$, can be expressed as

$Q^{2}=m_{e}^{2}\left[\left(\delta^{2}+\delta^{\prime
2}-2\delta\delta^{\prime}\cos\phi\right)+m_{e}^{2}\left(\frac{1+\delta^{2}}{2\epsilon}+\frac{1+\delta^{\prime
2}}{2\epsilon^{\prime}}\right)^{2}\right].$ (A.13)

The total phase space
$\differential\tilde{\Pi}=\differential\omega\differential\delta\differential\delta^{\prime}\differential\phi$.
We also explicitly include an atomic form factor $\left\lvert
F(Q^{2})\right\rvert^{2}$ Schiff:1951zza ; Tsai:1973py which screens
interactions at small $Q^{2}$ – it has the behaviour that for small $Q^{2}$,
$\left\lvert F(Q^{2})\right\rvert^{2}\propto Q^{4}$ (eliminating the $Q^{-4}$
singularity present in Eq. A.11) 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 $\epsilon$. 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 ParticleDataGroup:2022pth . The differential cross-section in the
small-angle approximation is given by Berestetskii:1982qgu

$\begin{split}\frac{\differential\tilde{\sigma}}{\differential\tilde{\Pi}}=\frac{8}{\pi}\left\lvert
F(Q^{2})\right\rvert^{2}\alpha
r_{e}^{2}\frac{m_{e}^{4}\epsilon_{+}\epsilon_{-}}{\omega^{3}Q^{4}}\delta_{+}\delta_{-}\times\bigg{\\{}-\frac{\delta_{+}^{2}}{(1+\delta_{+}^{2})^{2}}-\frac{\delta_{-}^{2}}{(1+\delta_{-}^{2})^{2}}&+\frac{\omega^{2}}{2\epsilon_{+}\epsilon_{-}}\frac{\delta_{+}^{2}+\delta_{-}^{2}}{(1+\delta_{+}^{2})(1+\delta_{-}^{2})}\\\
&+\quantity(\frac{\epsilon_{+}}{\epsilon_{-}}+\frac{\epsilon_{-}}{\epsilon_{+}})\frac{\delta_{+}\delta_{-}\cos\phi}{(1+\delta_{+}^{2})(1+\delta_{-}^{2})}\bigg{\\}}.\end{split}$
(A.14)

Here, the incoming photon has energy $\omega$, and the outgoing electron
(positron) has energy $\epsilon_{-}$ ($\epsilon_{+}=\omega-\epsilon_{-}$).
Similar to Eq. A.13, $q$ is the momentum transferred to the nucleus and in
this approximation, the expression for $Q^{2}$ is

$Q^{2}=m_{e}^{2}\left[\left(\delta_{+}^{2}+\delta_{-}^{2}+2\delta_{+}\delta_{-}\cos\phi\right)+m_{e}^{2}\left(\frac{1+\delta_{+}^{2}}{2\epsilon_{+}}+\frac{1+\delta_{-}^{2}}{2\epsilon_{-}}\right)\right].$
(A.15)

The integration phase space is
$\differential\tilde{\Pi}=\differential\epsilon_{+}\differential\delta_{+}\differential\delta_{-}\differential\phi$.
The energy $\epsilon_{+}$ spans the range $[m_{e},\omega-m_{e}]$. The other
three variables have the same range as for electron bremsstrahlung, and
$\delta_{\pm}$ are related to the outgoing positron/electron angles
$\theta_{\pm}$ through a similar relation as Eq. A.12. The atomic form factor
$\left\lvert F(Q^{2})\right\rvert^{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_{\rm nucl.}(Q^{2})\approx 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 Bethe:1954zz
; Davies:1954zz ; Tsai:1973py . 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).

##  Appendix 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 Tsai:1973py ;
PhysRevD.8.3109 , which scales as the square of the nuclear charge, and an
incoherent (inelastic) contribution which scales linearly with $Z$

$|F(Q^{2})|^{2}=Z^{2}G_{\rm el}(Q^{2})+ZG_{\rm inel}(Q^{2})~{}.$ (B.16)

Where the two contributions are Bjorken:2009mm

$\displaystyle G_{\rm el}(Q^{2})$
$\displaystyle=\quantity(\frac{a^{2}Q^{2}}{1+a^{2}Q^{2}})^{2}\quantity(\frac{1}{1+Q^{2}/d})^{2}~{},$
(B.17) $\displaystyle G_{\rm inel}(Q^{2})$
$\displaystyle=\quantity(\frac{a^{\prime 2}Q^{2}}{1+a^{\prime
2}Q^{2}})^{2}\quantity(\frac{1+\frac{Q^{2}}{4m_{p}^{2}}(\mu_{p}^{2}-1)}{(1+Q^{2}/\Lambda_{H}^{2})^{4}})~{}.$
(B.18)

In each case the form factor is a product of terms accounting for atomic and
hadronic effects. For $G_{\rm el}$ these factors encode atomic screening of
nuclear charge and the finite size of nucleus, while for $G_{\rm inel}$ they
describe atomic excitation and finite proton size. The atomic parameters are
given by

$a=\frac{111}{Z^{1/3}m_{e}}\mbox{\quad
and\quad}a^{\prime}=\frac{773}{Z^{2/3}m_{e}}~{},$ (B.19)

while the hadronic parameters are

$\Lambda_{H}=0.84~{}{\rm GeV}\mbox{\quad
and\quad}d=\left(0.4A^{-1/3}\,\mathrm{GeV}\right)^{2}~{}.$ (B.20)

The proton’s magnetic moment is $\mu_{p}=2.79$, and $m_{p}$ is the mass of the
proton. Notice that the second bracketed term in Eq. B.18 is not squared, see
Footnote 4 of Ref. Celentano:2020vtu .

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_{\rm inel}$ is set
to zero, and only the atomic form factor in included in $G_{\rm el}(Q^{2})$
(i.e. $d\to\infty$). 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 $\gamma
Z\rightarrow\mu^{+}\mu^{-}Z$).

##  Appendix C Continuous processes

###  C.1 Energy loss

In PETITE we have not attempted a stochastic treatment of energy loss.444In
principle such a description would involve proper sampling from the Landau-
Vavilov distribution. This is computationally expensive and is therefore not
implemented in standard codes such as GEANT-4, which instead employ toy models
of ionization energy loss. We instead treat energy loss with a deterministic
function $\langle\differential E/\differential x\rangle$ which includes
contributions from soft bremsstrahlung (photon emissions below the cut $E_{\rm
min}$ – see Sec. 2.2) and ionization. Because we take $E_{\rm
min}\sim\mathrm{MeV}$, well below the critical energy for most materials
$\mathcal{O}(10-100\;\mathrm{MeV})$ ParticleDataGroup:2022pth ,
$\langle\differential E/\differential x\rangle$ is approximately an energy-
independent constant since hard bremsstrahlung (photon emissions above $E_{\rm
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_{\rm 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 Bethe:1953va ; Lynch:1990sq ,
rescaling the outputted MCS angle can be interpreted as re-scaling the
characteristic scattering angle $\chi_{c}$.

#### Bethe-Molière

The slower implementation in PETITE is based on the theory of multiple
scattering first developed by Molière Moliere:1947zza , and then simplified
and improved by Bethe Bethe:1953va . In practice, MCS in this formalism is
implemented using an asymptotic expansion in a parameter $B$, defined by the
transcendental equation

$B+\ln B=\ln\frac{\Omega}{1.167}~{},$ (C.21)

where

$\Omega=\frac{\chi_{c}^{2}}{\chi_{a}^{2}}=\frac{7800~{}t~{}(Z+1)Z^{1/3}}{\beta^{2}A\quantity[1+3.34\quantity(\frac{Z\alpha}{\beta})^{2}]}$
(C.22)

with the characteristic scattering, $\chi_{c}$, and screening, $\chi_{a}$,
angles being defined below Eq. 4.4. Here $t$ is the thickness of material
traversed in g/cm2, $\beta$ the particle velocity, $Z$ the atomic number and
$A$ the atomic weight in g/mol. $\Omega$ is the mean number of scatters, so
typically $\Omega\gg 1$ unless $t<10^{-3}X_{0}$ Lynch:1990sq .

The Molière distribution is given by

$\displaystyle f(x)$ $\displaystyle=f_{0}(x)+B^{-1}f_{1}(x)+\ldots$ (C.23)
$\displaystyle f_{0}(x)$ $\displaystyle={\rm e}^{-x}~{},$ (C.24)
$\displaystyle f_{1}(x)$ $\displaystyle={\rm
e}^{-x}(1-x)\quantity[\log(x)-{\rm Ei}(x)+(1-2{\rm e}^{-x})]~{}.$ (C.25)

The probability distribution function for scattering through an angle $\theta$
is then

$P(\theta)~{}\theta\differential\theta=f\quantity(\vartheta^{2})~{}\vartheta\differential\vartheta,$
(C.26)

where $\vartheta=\theta/\sqrt{B\chi_{c}^{2}}$. 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 Lynch:1990sq ,

$\displaystyle\sigma^{2}$
$\displaystyle=\frac{\chi_{c}^{2}}{1+F}\quantity[\frac{1+v}{v}\log(1+v)-1]~{},$
(C.27) $\displaystyle v$ $\displaystyle=\frac{\Omega}{2(1-F)}~{}.$ (C.28)

where $\Omega$ is given in Eq. C.22 and we take $F=0.98$ in PETITE.

##  Appendix 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
$\differential\tilde{\Pi}$ (defined below Eq. A.13 and Eq. A.15). For
bremsstrahlung, the new coordinates are defined in Table 1. Since the
integrand in Eq. A.11 is singular for small $q^{2}$ and scales like
$1/\omega$, it peaks when $\delta\approx\delta^{\prime}$ and when
$\delta,\delta^{\prime},\omega\rightarrow 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.11. VEGAS Integration Variable – SM Bremsstrahlung |  Mapping to Eq. A.11
---|---
$x_{1}\in[0,1]$ | $\omega=E_{\gamma}^{\rm min.}+x_{1}\left(\epsilon-m_{e}-E_{\gamma}^{\rm min.}\right)$
$x_{2}\in[0,2]$ | $\delta=\frac{\epsilon}{2m_{e}}\left(x_{2}+x_{3}\right)$
$x_{3}\in[-2,2]$ | $\delta^{\prime}=\frac{\epsilon}{2m_{e}}\left(x_{2}-x_{3}\right)$
$x_{4}\in[0,1]$ | $\phi=2\pi\left(x_{4}-\frac{1}{2}\right)$

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.14.

Finally, the inclusion of atomic screening tames the low-$q^{2}$ Coulomb
singularity. This renders pair-production 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.11. VEGAS Integration Variable – SM Pair Production |  Mapping to Eq. A.14
---|---
$x_{1}\in[0,1]$ | $\epsilon_{+}=m_{e}+x_{1}\left(\omega-2m_{e}\right)$
$x_{2}\in[0,2]$ | $\delta_{+}=\frac{\omega}{2m_{e}}\left(x_{2}+x_{3}\right)$
$x_{3}\in[-2,2]$ | $\delta_{-}=\frac{\omega}{2m_{e}}\left(x_{2}-x_{3}\right)$
$x_{4}\in[0,1]$ | $\phi=2\pi x_{4}$

##  Appendix E Validation of SM showers in PETITE

In this section we compare SM showers generated with PETITE to predictions
from analytic theory Tsai:1966js and from GEANT-4.

###  E.1 Comparison with analytic shower theory

The analytic shower theory developed by Tsai and Whitis in Ref. Tsai:1966js
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\;\mathrm{GeV}$ electron on a graphite
target. The track length for particle type $i=\\{e^{-},e^{+},\gamma\\}$ is

$T_{i}(E)=X_{0}\int_{0}^{T}\differential tI_{i}(t,E),$ (E.29)

where $t$ is the depth inside the target in radiation lengths, $I_{i}(t,E)$ is
the differential particle intensity, i.e. $\int 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 $\int\differential ET_{i}(E)\sigma(E)$. We can therefore think
of $\differential ET_{i}(E)$ as the cumulative distance travelled by all
particles of type $i$ with energy in the interval $[E,E+\differential E]$ 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. Marsicano:2018krp :

$T_{i}(E)=X_{0}\Delta t\sum_{j}I_{i}(t_{j},E),$ (E.30)

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 $\Delta 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 ${\sim}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 photons (the
photons produced in those same bremsstrahlung events). However, it
significantly over-predicts first generation electrons and positrons
($e^{\pm}$ that are pair produced from 1st generation photons). We have
tracked this issue to the “complete screening” approximation employed in Ref.
Tsai:1966js . This approximation replaces the energy-dependence of the cross-
section $\sigma(E)$, by its asymptotic value at large energies, $\sigma_{\rm
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
$\lambda(E)$. Since $\sigma(E)\leq\sigma_{\rm asym}$ we have
$\lambda(E)\geq\lambda_{\rm 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^{\pm}$ pairs compared to simulations is
also evident in previous work, see, e.g., Fig. 4 of Ref. Marsicano:2018krp .

It is also interesting to consider the impact of all other processes that are
present in PETITE, but not in Ref. Tsai:1966js : 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
$\gamma$ emissions, the photon track length is larger than that predicted by
Ref. Tsai:1966js (this is already evident in the right panel of Fig. 7).

Figure 7: Track length distributions for $e^{-}$, $e^{+}$ and $\gamma$ for an
initial $10\;\mathrm{GeV}$ $e^{-}$ beam on a thick graphite target. Each panel
compares predictions of PETITE to the analytic calculations of Tsai and Whitis
Tsai:1966js . All processes except for bremsstrahlung and pair production have
been turned off in PETITE to closely match the assumptions of Ref. Tsai:1966js
. 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^{\pm}$
they generate. The analytic theory describes well the straggling of the
primary beam and first generation photons, but it systematically over-predicts
$e^{\pm}$ pair production at lower energies. The bars on the PETITE results
are statistical uncertainties for $10^{3}$ showers. Figure 8: Same as Fig. 7,
but with all processes enabled in PETITE.

Figure 9: Signed $\chi^{2}$ test statistic for PETITE and GEANT-4 predictions
of $dE/dr_{T}dz$, the energy-weighted single-shower particle distribution in
transverse ($r_{T}$) and longitudinal distance ($z$) from the beam axis, for a
10 GeV photon impinging on a graphite target. Each column corresponds to
$e^{-}$, $e^{+}$ and $\gamma$ respectively. The two rows use different values
of $E_{\mathrm{min}}$, the minimum particle energy tracked in the shower. The
$\chi^{2}$ values of $\mathcal{O}(1)$ indicate similarity of distributions
(i.e., within statistical uncertainties). The vertical and horizontal dashed
lines indicate estimated $95\%$ shower energy containment in the longitudinal
and transverse directions; they correspond to $\sim 14.4X_{0}$ and
$2r_{\mathrm{Moliere}}$, respectively Fabjan:2003aq .

###  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 $\chi^{2}$
for average particle energies ($\overline{E}$) in each $r_{T}$ and $z$ bin:

$\chi^{2}(r_{T},z)=\operatorname{sgn}\left(\overline{E}(\texttt{PETITE})-\overline{E}(\texttt{GEANT-4})\right)\frac{(\overline{E}(\texttt{PETITE})-\overline{E}(\texttt{GEANT-4}))^{2}}{\sigma_{E}^{2}(\texttt{PETITE})+\sigma^{2}_{E}(\texttt{GEANT-4})}$
(E.31)

where $\sigma_{E}^{2}$ is the variance of $\overline{E}$. The particle number
in each bin is sensitive to the minimum tracking energy $E_{\mathrm{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_{\mathrm{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
$\chi^{2}\ll 1$ in each bin. Statistically-significant differences emerge at
large transverse and longitudinal distances when we lower $E_{\mathrm{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.

##  Appendix 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)\to
e(p^{\prime})+N(P^{\prime})+V(k)$ is Liu:2017htz

$\mathcal{M}=e^{2}\varepsilon\frac{F(Q^{2})}{-Q^{2}}\varepsilon_{\mu}(k)^{*}\bar{u}(p^{\prime})\left[(\not{P}+\not{P}^{\prime})\frac{{\rm
i}(\not{p}-\not{k}+m_{e})}{m_{V}^{2}-2p\cdot
k}\gamma^{\mu}+\gamma^{\mu}\frac{{\rm
i}(\not{p}^{\prime}+\not{k}+m_{e})}{m_{V}^{2}+2p^{\prime}\cdot
k}(\not{P}+\not{P}^{\prime})\right]u(p)~{}.$ (F.32)

where $q=P^{\prime}-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 Tsai:1989vw ; Gninenko:2017yus

$\frac{\differential\sigma}{\differential
E_{V}\differential\cos\theta}=\left(\frac{\varepsilon^{2}\alpha^{3}|\vec{k}|}{8M_{N}^{2}|\vec{p}||\vec{k}-\vec{p}|}\right)\int_{Q^{2}_{\rm
min}}^{Q^{2}_{\rm max}}\frac{\differential
Q^{2}}{Q^{4}}|F(Q^{2})|^{2}I_{\phi}(t,E_{V},\cos\theta;M_{N},m_{V},E_{e})~{}.$
(F.33)

Here $I_{\phi}$ is the result of integrating the matrix element squared over
the azimuthal angle of the exchanged momentum $\vec{q}$ and $E_{0}$ is the
initial electron’s energy in the nucleus rest frame. The full expressions for
$I_{\phi},Q_{\rm min,max}$ are presented in Ref. Tsai:1989vw ;
Gninenko:2017yus .

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}\approx-Q^{2}_{\mathrm{min}}$ Bjorken:2009mm ,
where the minimum photon virtuality is

$Q_{\mathrm{min}}\approx\frac{m_{V}^{2}}{2E}=0.5\;\mathrm{keV}\;\;\bigg{(}\frac{m_{V}}{10\;\mathrm{MeV}}\bigg{)}^{2}\bigg{(}\frac{100\;\mathrm{GeV}}{E}\bigg{)}~{}.$
(F.34)

This approximate form is valid when $E\gg m_{V}$. This exercise also shows
that for $m_{V}>m_{e}$, the amplitude is enhanced for $x\approx 1$ as stressed
in Ref. Bjorken:2009mm .

The $q^{2}\approx-Q^{2}_{\mathrm{min}}$ regime mentioned above is relevant
because the photon propagator $\sim 1/q^{2}$ contributes another apparent
singularity at small $Q^{2}$. This divergence is regulated in two ways: by the
kinematic cut-off $Q^{2}_{\mathrm{min}}$, and by the atomic and nuclear form-
factor $F(Q^{2})\propto a^{2}Q^{2}$ for $a^{2}Q^{2}\lesssim 1$ with
$1/a\approx 4.6Z^{1/3}\times 10^{-6}$ GeV (see Eq. B.17 and Eq. B.19). 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}\gtrsim 3\;\mathrm{MeV}Z^{1/6}\sqrt{E/\mathrm{GeV}}$,
which is the case for most of the parameter space we are interested in
(especially if experimental detection places a ${\sim}$ 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 Lepage:1977sw ; Lepage:2020tgj is able to
correctly integrate and sample from the partially integrated differential
cross-section from Refs. Tsai:1989vw ; Liu:2017htz ; Gninenko:2017yus . 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 Belyaev:2012qa , MadGraph Alwall:2014hca and the Weizsäcker-
Williams approximation Liu:2017htz .

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. Celentano:2020vtu ;
Forbes:2022bvo . 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 Mattelaer:2021xdr ,
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. Mattelaer:2021xdr 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.

Figure 10: 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 Bjorken:2009mm .

Figure 11: Energy (left) and angle (right) distributions of the emitted vector
particle in $e^{-}+N\to e^{-}+N+V$ with an beam energy of $10$ GeV and
$m_{V}=0.1$ GeV.

##  Appendix 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\gamma\rightarrow eV$.
The differential cross section for dark Compton scattering is closely related
to that for Compton scattering, Eq. A.7:

$\derivative{t}\sigma_{\rm Dark\
Comp.}=\varepsilon^{2}\derivative{t}\sigma_{\rm
Comp.}-\frac{4\varepsilon^{2}\pi\alpha^{2}m_{V}^{2}}{(s-m_{e}^{2})^{4}(u-m_{e}^{2})^{2}}\left[stu-
m_{e}^{2}m_{V}^{2}(s+u)+m_{e}^{2}\left(2su+m_{e}^{2}\left(t-2m_{e}^{2}\right)\right)\right]~{},$
(G.35)

where the Mandelstam variables in $\sigma_{\rm Comp.}$ are redefined and
satisfy the relation $s+t+u=2m_{e}^{2}+m_{V}^{2}$.

In the $m_{e}\to 0$ limit, the differential cross sections for dark Compton
scattering and annihilation $e^{+}e^{-}\to\gamma V$

$\displaystyle\derivative{t}\sigma_{\rm Dark\;Comp.}$
$\displaystyle=\frac{2\varepsilon^{2}\pi\alpha^{2}}{s^{2}}\left[\frac{2m_{V}^{2}(s+u)-2m_{V}^{4}-s^{2}-u^{2}}{su}\right]~{},$
(G.36) $\displaystyle\derivative{t}\sigma_{\rm Dark\;Ann.}$
$\displaystyle=\frac{2\varepsilon^{2}\pi\alpha^{2}}{s^{2}}\left[\frac{tu(2m_{V}^{4}-2m_{V}^{2}(t+u)+t^{2}+u^{2})}{(t-m_{e}^{2})^{2}(u-m_{e}^{2})^{2}}\right]~{},$
(G.37)

where the Mandelstam variables are $s=(p_{\gamma}+p_{e^{-}})^{2}$,
$u=(p_{e^{-}}-k_{V})^{2}$ and $t=(p_{\gamma}-k_{V})^{2}$ for Compton
scattering. For the process of $e^{+}e^{-}\rightarrow\gamma V$ the Mandelstams
are $s=(p_{e^{+}}+p_{e^{-}})^{2}$, $u=(k_{V}-p_{e^{-}})^{2}$,
$t=(k_{\gamma}-p_{e^{-}})^{2}$. 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.

##  Appendix 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. Bjorken:1988as in
the context of axions and in Ref. Marsicano:2018krp for dark photons. At tree
level, there are two processes $e^{+}e^{-}\to V$ and $e^{+}e^{-}\to\gamma 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. Marsicano:2018krp avoid
this issue by putting an (arbitrary) minimum energy cut on the SM photon.
Additionally, the process $e^{+}e^{-}\to\gamma V$ contains a collinear
singularity (regulated by $m_{e}$) when the emitted photon is aligned with
$e^{\pm}$; 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_{V}^{2}$, where $x_{\pm}$ are the longitudinal momentum
fractions of $e^{\pm}$. The cross-section for this radiative return is given
by

$\sigma_{rr}(s)=\left[\frac{4\pi^{2}\alpha\varepsilon^{2}}{s}\beta_{f}\left(\frac{3-\beta_{f}^{2}}{2}\right)\right]\int_{m_{V}^{2}/s}^{1}\frac{dx}{x}f_{e}(x,s)f_{e}\left(\frac{m_{V}^{2}}{xs},s\right),$
(H.38)

where $\beta_{f}^{2}=1-4m_{e}^{2}/m_{V}^{2}$ 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. Nicrosini:1986sm .
We will use the Kuraev-Fadin Kuraev:1985hb expression555As illustrated in Ref.
Greco:2016izi , this PDF works best as $x\to 1$; we use it because of its
simplicity.

$f_{e}(x,s)=\frac{\beta}{16}\left((8+3\beta)(1-x)^{\beta/2-1}-4(1+x)\right)+\dots,$
(H.39)

where $\beta=(2\alpha/\pi)(\ln s/m^{2}_{e}-1)$ and the ellipsis stands for
higher order terms in $\beta$. Here we see the impact of the collinear
singularities in factors of $\ln s/m_{e}^{2}$, and the soft divergence as
$(1-x)^{1-\beta/2}$. For the typical GeV-scale energies $\log s/m_{e}^{2}$ is
$\mathcal{O}(10)$. The resummation therefore has important effects both at
large incoming particle energies (where these logarithms grow), and also near
threshold $s\approx m_{V}^{2}$. In particular, evaluating Eq. H.38 as $s\to
m_{V}^{2}$, one finds that $\sigma_{rr}\propto(s-m_{V}^{2})^{-1+\beta}$,
compared to $(s-m_{V}^{2})^{-1}$ in the tree-level approximation for
$e^{+}e^{-}\to\gamma 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.38) to generate
momentum fractions $x$ and $m_{V}^{2}/(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_{V}^{2}/(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\%$.

Figure 12: Ratio of the radiative return cross-section $\sigma_{rr}$ to the
(leading order) tree-level annihilation cross-section $\sigma_{LO}$ for
several dark vector masses, as a function of $\sqrt{s}$. The cross section is
suppressed by an $O(1)$ 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
$\sqrt{s}$, because radiative return allows the parton to access the
resonance.

## References

  * (1) N. Craig et al., “Snowmass Theory Frontier Report”, 2211.05772.
  * (2) G. Bertone, D. Hooper and J. Silk, “Particle dark matter: Evidence, candidates and constraints”, _Phys. Rept._ 405 (2005) 279 [hep-ph/0404175].
  * (3) M. Battaglieri et al., “US Cosmic Visions: New Ideas in Dark Matter 2017: Community Report”, in _U.S. Cosmic Visions: New Ideas in Dark Matter_ , 7, 2017, 1707.04591.
  * (4) T. Kajita, “Nobel Lecture: Discovery of atmospheric neutrino oscillations”, _Rev. Mod. Phys._ 88 (2016) 030501.
  * (5) A. B. McDonald, “Nobel Lecture: The Sudbury Neutrino Observatory: Observation of flavor change for solar neutrinos”, _Rev. Mod. Phys._ 88 (2016) 030502.
  * (6) J. Beacham et al., “Physics Beyond Colliders at CERN: Beyond the Standard Model Working Group Report”, _J. Phys. G_ 47 (2020) 010501 [1901.09966].
  * (7) C. Antel et al., “Feebly Interacting Particles: FIPs 2022 workshop report”, in _Workshop on Feebly-Interacting Particles_ , 5, 2023, 2305.01715.
  * (8) L3 collaboration, “Search for heavy neutral and charged leptons in $e^{+}e^{-}$ annihilation at LEP”, _Phys. Lett. B_ 517 (2001) 75 [hep-ex/0107015].
  * (9) CMS collaboration, “Search for vector-like leptons in multilepton final states in proton-proton collisions at $\sqrt{s}$ = 13 TeV”, _Phys. Rev. D_ 100 (2019) 052003 [1905.10853].
  * (10) ATLAS collaboration, “Search for third-generation vector-like leptons in $pp$ collisions at $\sqrt{s}=13\,\text{TeV}$ with the ATLAS detector”, _JHEP_ 07 (2023) 118 [2303.05441].
  * (11) S. Gori et al., “Dark Sector Physics at High-Intensity Experiments”, 2209.04671.
  * (12) P. Harris, P. Schuster and J. Zupan, “Snowmass White Paper: New flavors and rich structures in dark sectors”, in _Snowmass 2021_ , 7, 2022, 2207.08990.
  * (13) A. Berlin, N. Blinov, G. Krnjaic, P. Schuster and N. Toro, “Dark Matter, Millicharges, Axion and Scalar Particles, Gauge Bosons, and Other New Physics with LDMX”, _Phys. Rev. D_ 99 (2019) 075001 [1807.01730].
  * (14) NA64 collaboration, “Search for Light Dark Matter with NA64 at CERN”, 2307.02404.
  * (15) P. deNiverville, M. Pospelov and A. Ritz, “Observing a light dark matter beam with neutrino experiments”, _Phys. Rev. D_ 84 (2011) 075020 [1107.4580].
  * (16) B. Batell, R. Essig and Z. Surujon, “Strong Constraints on Sub-GeV Dark Sectors from SLAC Beam Dump E137”, _Phys. Rev. Lett._ 113 (2014) 171802 [1406.2698].
  * (17) SHiP collaboration, “A facility to Search for Hidden Particles (SHiP) at the CERN SPS”, 1504.04956.
  * (18) BDX collaboration, “Dark Matter Search in a Beam-Dump eXperiment (BDX) at Jefferson Lab”, 1607.01390.
  * (19) A. de Gouvêa, P. J. Fox, R. Harnik, K. J. Kelly and Y. Zhang, “Dark Tridents at Off-Axis Liquid Argon Neutrino Detectors”, _JHEP_ 01 (2019) 001 [1809.06388].
  * (20) V. De Romeri, K. J. Kelly and P. A. N. Machado, “DUNE-PRISM Sensitivity to Light Dark Matter”, _Phys. Rev. D_ 100 (2019) 095010 [1903.10505].
  * (21) P. A. Machado, O. Palamara and D. W. Schmitz, “The Short-Baseline Neutrino Program at Fermilab”, _Ann. Rev. Nucl. Part. Sci._ 69 (2019) 363 [1903.04608].
  * (22) R. Harnik, Z. Liu and O. Palamara, “Millicharged Particles in Liquid Argon Neutrino Experiments”, _JHEP_ 07 (2019) 170 [1902.03246].
  * (23) ArgoNeuT collaboration, “Improved Limits on Millicharged Particles Using the ArgoNeuT Experiment at Fermilab”, _Phys. Rev. Lett._ 124 (2020) 131801 [1911.07996].
  * (24) DUNE collaboration, “Prospects for beyond the Standard Model physics searches at the Deep Underground Neutrino Experiment”, _Eur. Phys. J. C_ 81 (2021) 322 [2008.12769].
  * (25) M. Breitbach, L. Buonocore, C. Frugiuele, J. Kopp and L. Mittnacht, “Searching for physics beyond the Standard Model in an off-axis DUNE near detector”, _JHEP_ 01 (2022) 048 [2102.03383].
  * (26) A. Apyan et al., “DarkQuest: A dark sector upgrade to SpinQuest at the 120 GeV Fermilab Main Injector”, in _Snowmass 2021_ , 3, 2022, 2203.08322.
  * (27) P. Coloma, J. López-Pavón, L. Molina-Bueno and S. Urrea, “New Physics searches using ProtoDUNE and the CERN SPS accelerator”, 2304.06765.
  * (28) H. J. Bhabha and W. Heitler, “The passage of fast electrons and the theory of cosmic showers”, _Proceedings of the Royal Society of London. Series A - Mathematical and Physical Sciences_ 159 (1937) 432.
  * (29) T. K. Gaisser, _Cosmic rays and particle physics_. 1990.
  * (30) Y.-S. Tsai and V. Whitis, “Thick target bremsstrahlung and target consideration for secondary particle production by electrons”, _Phys. Rev._ 149 (1966) 1248.
  * (31) P. deNiverville, C.-Y. Chen, M. Pospelov and A. Ritz, “Light dark matter in neutrino beams: production modelling and scattering signatures at MiniBooNE, T2K and SHiP”, _Phys. Rev. D_ 95 (2017) 035006 [1609.01770].
  * (32) G. Magill, R. Plestid, M. Pospelov and Y.-D. Tsai, “Millicharged particles in neutrino experiments”, _Phys. Rev. Lett._ 122 (2019) 071801 [1806.03310].
  * (33) G. Magill, R. Plestid, M. Pospelov and Y.-D. Tsai, “Dipole Portal to Heavy Neutral Leptons”, _Phys. Rev. D_ 98 (2018) 115015 [1803.03262].
  * (34) J. M. Berryman, A. de Gouvea, P. J. Fox, B. J. Kayser, K. J. Kelly and J. L. Raaf, “Searches for Decays of New Particles in the DUNE Multi-Purpose Near Detector”, _JHEP_ 02 (2020) 174 [1912.07622].
  * (35) R. Plestid, V. Takhistov, Y.-D. Tsai, T. Bringmann, A. Kusenko and M. Pospelov, “New Constraints on Millicharged Particles from Cosmic-ray Production”, _Phys. Rev. D_ 102 (2020) 115032 [2002.11732].
  * (36) Y.-S. Tsai, “Axion bremsstrahlung by an electron beam”, _Phys. Rev. D_ 34 (1986) 1326.
  * (37) B. Batell, M. Pospelov and A. Ritz, “Exploring Portals to a Hidden Sector Through Fixed Targets”, _Phys. Rev. D_ 80 (2009) 095024 [0906.5614].
  * (38) M. Reece and L.-T. Wang, “Searching for the light dark gauge boson in GeV-scale experiments”, _JHEP_ 07 (2009) 051 [0904.1743].
  * (39) J. Blümlein and J. Brunner, “New Exclusion Limits on Dark Gauge Forces from Proton Bremsstrahlung in Beam-Dump Data”, _Phys. Lett. B_ 731 (2014) 320 [1311.3870].
  * (40) S. Foroughi-Abari and A. Ritz, “Dark sector production via proton bremsstrahlung”, _Phys. Rev. D_ 105 (2022) 095045 [2108.05900].
  * (41) L. Buonocore, C. Frugiuele, F. Maltoni, O. Mattelaer and F. Tramontano, “Event generation for beam dump experiments”, _JHEP_ 05 (2019) 028 [1812.06771].
  * (42) E. Nardi, C. D. R. Carvajal, A. Ghoshal, D. Meloni and M. Raggi, “Resonant production of dark photons in positron beam dump experiments”, _Phys. Rev. D_ 97 (2018) 095004 [1802.04756].
  * (43) A. Celentano, L. Darmé, L. Marsicano and E. Nardi, “New production channels for light dark matter in hadronic showers”, _Phys. Rev. D_ 102 (2020) 075026 [2006.09419].
  * (44) F. Capozzi, B. Dutta, G. Gurung, W. Jang, I. M. Shoemaker, A. Thompson et al., “Extending the reach of leptophilic boson searches at DUNE and MiniBooNE with bremsstrahlung and resonant production”, _Phys. Rev. D_ 104 (2021) 115010 [2108.03262].
  * (45) L. D. Landau and I. Pomeranchuk, “Limits of applicability of the theory of bremsstrahlung electrons and pair production at high-energies”, _Dokl. Akad. Nauk Ser. Fiz._ 92 (1953) 535.
  * (46) A. B. Migdal, “Bremsstrahlung and pair production in condensed media at high-energies”, _Phys. Rev._ 103 (1956) 1811.
  * (47) GEANT4 collaboration, “GEANT4–a simulation toolkit”, _Nucl. Instrum. Meth. A_ 506 (2003) 250.
  * (48) J. Allison et al., “Geant4 developments and applications”, _IEEE Trans. Nucl. Sci._ 53 (2006) 270.
  * (49) J. Allison et al., “Recent developments in Geant4”, _Nucl. Instrum. Meth. A_ 835 (2016) 186.
  * (50) H. Bethe and W. Heitler, “On the Stopping of fast particles and on the creation of positive electrons”, _Proc. Roy. Soc. Lond. A_ 146 (1934) 83.
  * (51) Y.-S. Tsai, “Pair Production and Bremsstrahlung of Charged Leptons”, _Rev. Mod. Phys._ 46 (1974) 815.
  * (52) H. A. Bethe, “Moliere’s theory of multiple scattering”, _Phys. Rev._ 89 (1953) 1256.
  * (53) G. R. Lynch and O. I. Dahl, “Approximations to multiple Coulomb scattering”, _Nucl. Instrum. Meth. B_ 58 (1991) 6.
  * (54) Particle Data Group collaboration, “Review of Particle Physics”, _PTEP_ 2022 (2022) 083C01.
  * (55) V. B. Berestetskii, E. M. Lifshitz and L. P. Pitaevskii, _Quantum Electrodynamics_ , vol. 4 of _Course of Theoretical Physics_. Pergamon Press, Oxford, 1982.
  * (56) E. H. Bellamy, R. Hofstadter, W. L. Lakin, M. L. Perl and W. T. Toner, “A Search For Fractionally Charged Particles”, _Phys. Rev._ 166 (1968) 1391.
  * (57) A. Barna, J. Cox, F. Martin, M. L. Perl, T. H. Tan, W. T. Toner et al., “A Search For New Particles Produced By High-Energy Photons”, _Phys. Rev._ 173 (1968) 1391.
  * (58) M. Schwartz et al., “Electrodynamics at very low and very high energies”, .
  * (59) A. F. Rothenberg, “A Search for unknown sources of neutrino like particles”, other thesis, Stanford University, 5, 1972.
  * (60) T. W. Donnelly, S. J. Freedman, R. S. Lytel, R. D. Peccei and M. Schwartz, “Do Axions Exist?”, _Phys. Rev. D_ 18 (1978) 1607.
  * (61) J. D. Bjorken, R. Essig, P. Schuster and N. Toro, “New Fixed-Target Experiments to Search for Dark Gauge Forces”, _Phys. Rev. D_ 80 (2009) 075018 [0906.0580].
  * (62) S. Andreas, C. Niebuhr and A. Ringwald, “New Limits on Hidden Photons from Past Electron Beam Dumps”, _Phys. Rev. D_ 86 (2012) 095019 [1209.6083].
  * (63) J. D. Bjorken, S. Ecklund, W. R. Nelson, A. Abashian, C. Church, B. Lu et al., “Search for Neutral Metastable Penetrating Particles Produced in the SLAC Beam Dump”, _Phys. Rev. D_ 38 (1988) 3375.
  * (64) L. Marsicano, M. Battaglieri, A. Celentano, R. De Vita and Y.-M. Zhong, “Probing Leptophilic Dark Sectors at Electron Beam-Dump Facilities”, _Phys. Rev. D_ 98 (2018) 115022 [1812.03829].
  * (65) C. Rella, B. Döbrich and T.-T. Yu, “Searching for muonphilic dark sectors with proton beams”, _Phys. Rev. D_ 106 (2022) 035023 [2205.09870].
  * (66) M. Bondi, A. Celentano, R. R. Dusaev, D. V. Kirpichnikov, M. M. Kirsanov, N. V. Krasnikov et al., “Fully Geant4 compatible package for the simulation of Dark Matter in fixed target experiments”, _Comput. Phys. Commun._ 269 (2021) 108129 [2101.12192].
  * (67) T. Eichlersmith, J. Mans, O. Moreno, J. Muse, M. Revering and N. Toro, “Simulation of dark bremsstrahlung in Geant4”, _Comput. Phys. Commun._ 287 (2023) 108690 [2211.03873].
  * (68) Y.-S. Liu, D. McKeen and G. A. Miller, “Validity of the Weizsäcker-Williams approximation and the analysis of beam dump experiments: Production of a new scalar boson”, _Phys. Rev. D_ 95 (2017) 036010 [1609.06781].
  * (69) Y.-S. Liu and G. A. Miller, “Validity of the Weizsäcker-Williams approximation and the analysis of beam dump experiments: Production of an axion, a dark photon, or a new axial-vector boson”, _Phys. Rev. D_ 96 (2017) 016004 [1705.01633].
  * (70) Y.-S. Tsai, “Production Of Neutral Bosons By An Electron Beam”, _Phys. Rev. D_ 40 (1989) 760.
  * (71) S. N. Gninenko, D. V. Kirpichnikov, M. M. Kirsanov and N. V. Krasnikov, “The exact tree-level calculation of the dark photon production in high-energy electron scattering at the CERN SPS”, _Phys. Lett. B_ 782 (2018) 406 [1712.05706].
  * (72) G. P. Lepage, “A New Algorithm for Adaptive Multidimensional Integration”, _J. Comput. Phys._ 27 (1978) 192.
  * (73) G. P. Lepage, “Adaptive multidimensional integration: VEGAS enhanced”, _J. Comput. Phys._ 439 (2021) 110386 [2009.05112].
  * (74) E. Kuraev and V. S. Fadin, “On Radiative Corrections to e+ e- Single Photon Annihilation at High-Energy”, _Sov. J. Nucl. Phys._ 41 (1985) 466.
  * (75) O. Nicrosini and L. Trentadue, “Soft Photons and Second Order Radiative Corrections to e+ e- —$>$ Z0”, _Phys. Lett. B_ 196 (1987) 551.
  * (76) M. Karliner, M. Low, J. L. Rosner and L.-T. Wang, “Radiative return capabilities of a high-energy, high-luminosity $e^{+}e^{-}$ collider”, _Phys. Rev. D_ 92 (2015) 035010 [1503.07209].
  * (77) M. Greco, T. Han and Z. Liu, “ISR effects for resonant Higgs production at future lepton colliders”, _Phys. Lett. B_ 763 (2016) 409 [1607.03210].
  * (78) C. Bierlich et al., “A comprehensive guide to the physics and usage of PYTHIA 8.3”, _SciPost Phys. Codeb._ 2022 (2022) 8 [2203.11601].
  * (79) G. Moliere, “Theorie der Streuung schneller geladener Teilchen I. Einzelstreuung am abgeschirmten Coulomb-Feld”, _Z. Naturforsch. A_ 2 (1947) 133.
  * (80) B. Gottschalk, A. M. Koehler, R. J. Schneider, J. M. Sisterson and M. S. Wagner, “Multiple Coulomb scattering of 160 MeV protons”, _Nuclear Instruments and Methods in Physics Research B_ 74 (1993) 467.
  * (81) A. Makarova, B. Gottschalk and W. Sauerwein, “Comparison of geant4 multiple coulomb scattering models with theory for radiotherapy protons”, _Physics in Medicine & Biology_ 62 (2017) 5959.
  * (82) Y. M. Andreev et al., “Improved exclusion limit for light dark matter from e+e- annihilation in NA64”, _Phys. Rev. D_ 104 (2021) L091701 [2108.04195].
  * (83) H. Hirayama, Y. Namito, A. F. Bielajew, S. J. Wilderman and W. R. Nelson, “The EGS5 code system”, .
  * (84) A. Ferrari, P. R. Sala, A. Fasso and J. Ranft, “FLUKA: A multi-particle transport code (Program version 2005)”, .
  * (85) T. T. Böhlen, F. Cerutti, M. P. W. Chin, A. Fassò, A. Ferrari, P. G. Ortega et al., “The FLUKA Code: Developments and Challenges for High Energy and Medical Applications”, _Nucl. Data Sheets_ 120 (2014) 211.
  * (86) D. Bernard, “A 5D, polarised, Bethe-Heitler event generator for $\gamma\to e^{+}e^{-}$ conversion”, _Nucl. Instrum. Meth. A_ 899 (2018) 85 [1802.08253].
  * (87) GEANT4 collaboration, “Physics reference manual”, tech. rep., 2022\.
  * (88) S. M. Seltzer and M. J. Berger, “Bremsstrahlung spectra from electron interactions with screened atomic nuclei and orbital electrons”, _Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms_ 12 (1985) 95.
  * (89) S. Ghinescu, B. Döbrich, E. Minucci and T. Spadaro, “A biased MC for muon production for beam-dump experiments”, _Eur. Phys. J. C_ 81 (2021) 767 [2106.01932].
  * (90) H. H. Patel, “Package-X: A Mathematica package for the analytic calculation of one-loop integrals”, _Comput. Phys. Commun._ 197 (2015) 276 [1503.01469].
  * (91) V. Shtabovenko, R. Mertig and F. Orellana, “FeynCalc 9.3: New features and improvements”, _Comput. Phys. Commun._ 256 (2020) 107478 [2001.04407].
  * (92) C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau et al., “Array programming with NumPy”, _Nature_ 585 (2020) 357.
  * (93) P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau et al., “SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python”, _Nature Methods_ 17 (2020) 261.
  * (94) J. D. Hunter, “Matplotlib: A 2d graphics environment”, _Computing in Science & Engineering_ 9 (2007) 90.
  * (95) C. Itzykson and J. B. Zuber, _Quantum Field Theory_ , International Series In Pure and Applied Physics. McGraw-Hill, New York, 1980.
  * (96) L. I. Schiff, “Energy-Angle Distribution of Thin Target Bremsstrahlung”, _Phys. Rev._ 83 (1951) 252.
  * (97) H. A. Bethe and L. C. Maximon, “Theory of Bremsstrahlung and Pair Production. 1. Differential Cross Section”, _Phys. Rev._ 93 (1954) 768.
  * (98) H. Davies, H. A. Bethe and L. C. Maximon, “Theory of Bremsstrahlung and Pair Production. 2. Integral Cross Section for Pair Production”, _Phys. Rev._ 93 (1954) 788.
  * (99) K. J. Kim and Y.-S. Tsai, “Improved weizsäcker-williams method and its application to lepton and $w$-boson pair production”, _Phys. Rev. D_ 8 (1973) 3109.