Lepton Jets from Radiating Dark Matter

The idea that dark matter forms part of a larger dark sector is very intriguing, given the high degree of complexity of the visible sector. In this paper, we discuss lepton jets as a promising signature of an extended dark sector. As a simple toy model, we consider an $\mathcal{O}(\text{GeV})$ DM fermion coupled to a new $U(1)'$ gauge boson (dark photon) with a mass of order GeV and kinetically mixed with the Standard Model photon. Dark matter production at the LHC in this model is typically accompanied by collinear radiation of dark photons whose decay products can form lepton jets. We analyze the dynamics of collinear dark photon emission both analytically and numerically. In particular, we derive the dark photon energy spectrum using recursive analytic expressions, using Monte Carlo simulations in Pythia, and using an inverse Mellin transform to obtain the spectrum from its moments. In the second part of the paper, we simulate the expected lepton jet signatures from radiating dark matter at the LHC, carefully taking into account the various dark photon decay modes and allowing for both prompt and displaced decays. Using these simulations, we recast two existing ATLAS lepton jet searches to significantly restrict the parameter space of extended dark sector models, and we compute the expected sensitivity of future LHC searches.


I. INTRODUCTION
For decades, nature has been keeping us completely in the dark about the nature of dark matter. In view of this, the scope of theoretical and phenomenological studies has become significantly broader than it used to be some years ago. An often discussed possibility nowadays is a dark sector featuring new gauge interactions. Phenomenologically, such new gauge interactions can for instance lead to Sommerfeld enhancement of the dark matter (DM) annihilation cross section [1][2][3], to DM self-interactions [4,5], which may affect small scale structures in the Universe such as dwarf galaxies, or to DM bound states [6][7][8]. In this context, it is worth noting that there are long-standing discrepancies between theory and small scale structure observations, which could be explained either by dark matter self interactions [4,5] or by baryonic feedback [9,10]. Also recent observations of the galaxy cluster Abell 3827 could point towards self-interacting DM [11] (see, however, [12]).
Here, we discuss the prospects of detecting new gauge interactions in the dark sector at the LHC. We work in a simplified model with a dark matter fermion χ and a U (1) gauge boson A (dark photon), kinetically mixed with the photon [3,13] (see e.g. [14] for a review). This model, which we dub radiating dark matter, serves here as a proxy for a large class of more complicated scenarios, to which our conclusions can be applied as well. If the dark matter and dark photon masses are relatively small compared to the typical partonic center of mass energy at the LHC, emission of dark photons from a dark matter particle receives a strong enhancement in the collinear direction. This is analogous to the emission of collinear photons in QED or the emergence of a II. MODEL

A. The Lagrangian
The results of this paper can be applied to any model in which a 10 GeV DM particle interacts with a new 10 GeV gauge boson, as long as the typical energy at which DM particles are produced at the LHC is much higher than their masses. To illustrate our main points, we will use a specific toy model of radiating DM, in which the dark sector consists of a fermionic DM particle χ and a massive U (1) gauge boson A . The coupling strength between χ and A is described by a dark fine structure constant α A ≡ g 2 A /(4π). Moreover, the dark photon A is kinetically mixed with the SM photon. The Lagrangian for the dark sector is thus FIG. 1: Feynman diagram for dark matter pair production at the LHC through an s-channel Z resonance, followed by "dark radiation", i.e. emission of several-mostly soft or collinear-dark photons.
where m A and m χ are the masses of A and χ, respectively. We denote by F µν and F µν the field strength tensors of the SM photon and the dark photon A , respectively. The dimensionless coupling constant describes the strength of the kinetic mixing between the A and the photon. Due to SU (2) invariance, the A should also mix kinetically with the Z, but effects of this mixing are suppressed by factors of m 2 A /M 2 Z and are therefore neglected in the following. Since the kinetic mixing between A and the photon is too small to lead to significant DM production at the LHC, we assume an additional DM-SM coupling. For definiteness, we take this coupling to be through a heavy s-channel vector resonance Z with couplings to all quark flavors and to dark matter. Since the dynamics of dark radiation does not depend on the primary DM production mechanism but only on the production cross section and to some extent on the DM energy spectrum, our results will apply to any model in which DM particles can be produced in significant numbers at the LHC, including for instance models with contact interactions, t-channel mediators, or Higgs portal interactions. The relevant terms in the Lagrangian of our toy model are where q f is the SM quark field of flavor f and g q , g χ are the Z couplings to quarks and DM particles, respectively. For simplicity, we have assumed the coupling to quarks to be flavor universal. Fig. 1 illustrates DM pair production at the LHC via s-channel Z exchange, followed by radiation of several dark photons. Due to the assumed lightness of χ and A , the emission is enhanced in the collinear direction. Due to this enhancement, with a moderate dark fine structure constant α A ∼ O(0.1), we typically expect a few dark photons to be radiated in each DM pair production process. We will discuss the dynamics of this "dark radiation" in sec. III using a formalism analogous to parton showers in QCD.

B. Benchmark points
We define two benchmark points in the parameter space of our toy model, which we will use to illustrate our main points in secs. III and IV. In sec. IV, we will also discuss in detail how departing from these benchmark points affects our results. The two benchmark points A and B are summarized in table I, together with several phenomenological observables derived from them.
In both cases, we assume a Z mass of 1 TeV. We choose g q and g χ such that the resonant Z production cross section is about 1 pb at the 8 TeV LHC for benchmark point A, and about a factor of 10 smaller for benchmark point B. In both cases, the Z has a branching ratio ∼ 85% for the decay Z →χχ. We choose both m χ and m A to be of order GeV or below. Together with a moderately large dark fine structure constant α A = 0.2, this leads to the radiation of a significant number of A bosons when DM is produced at the LHC. Note that the number of m Z g q g χ m χ m A α A cτ σ 7 (Z ) σ 8 (Z ) σ 13 (Z ) Γ Z BR(Z → χχ) 2 n A 8 2 n A 13 [TeV] [GeV] [GeV] [mm] [ I: Values of the model parameters m Z (heavy mediator mass), g q , g χ (heavy mediator couplings), m χ (DM mass), m A (dark photon mass), α A (dark fine structure constant) and cτ (dark photon decay length) at our two benchmark points A and B. We also show the resulting values for several derived quantities, in particular the kinetic mixing corresponding to the given cτ and m A , the resonance cross sections σ 7 (Z ), σ 8 (Z ) and σ 13 (Z ) for Z production pp → Z at the 7 TeV, 8 TeV and 13 TeV LHC, respectively, the total decay width Γ Z of the heavy mediator, its branching ratio to DM pairs, and the average numbers 2 n A 8 (2 n A 13 ) of radiated dark photons in each DM pair production event at the 8 TeV (13 TeV) LHC.
radiated A bosons is almost independent of the collider center of mass energy because the boost of the primary DM particles, which determines the amount of radiation, is mostly dictated by the s-channel Z resonance. We choose the dark photon decay length cτ such that at benchmark point A, most dark photons will decay far from their production vertex, while at benchmark point B, the displacement of the decay vertex is much smaller. Thus, benchmark point A is in the realm of lepton jet searches requiring displaced decay vertices, while benchmark point B is most easily detectable in searches for prompt lepton jets. Note that in table I we treat cτ as a fundamental parameter and as a derived quantity. This will only be relevant when we vary the model parameters around their benchmark values in sec. IV (figs. 6 and 7). To make the discussion there more transparent, we choose to keep cτ rather than fixed when varying the other model parameters. Of course, as long as all other model parameters are fixed, there is a one-to-one correspondence between cτ and .
Our parameter choices are consistent with existing collider constraints: In particular, dijet searches do not exclude Z bosons with the parameters chosen here [38,39]. For instance, the ATLAS dijet search in 20.3 fb −1 of 8 TeV LHC data [38] restricts the product σ(Z ) × BR(Z → qq) × A of the Z production cross section, the branching ratio to a quark-antiquark final state, and the experimental acceptance factor A ∼ 0.6 to be 1 pb at 95% confidence level (CL). The corresponding CMS search [39] does not explicitly quote a limit at m Z = 1 TeV, but noting that no spectral feature is observed at this mass, the limits given on heavier Z bosons can be extrapolated down, indicating that the CMS limit is stronger than the ATLAS limits by perhaps a factor of two. The acceptance factor is similar in ATLAS and CMS. Limits on dijet resonances are also provided by the CDF experiment at the Tevatron, yielding σ(Z )×BR(Z → qq)×A 0.11 pb [40]. Since CDF requires both jets to have a rapidity |y| < 1, we conservatively estimate the acceptance in CDF is not larger than in CMS. Including an acceptance factor of 0.6, we predict at our benchmark point A σ(Z ) × BR(Z → qq) × A = 0.08 pb at the 8 TeV LHC and σ(Z ) × BR(Z → qq) × A = 0.001 pb at the Tevatron, which is clearly consistent with the above constraints. The pp → Z →χχ production cross sections at benchmark point B are even smaller.
DM production through heavy resonances is also constrained by monojet searches [41,42]. ATLAS sets a lower limit m Z / √ g q g χ 2 TeV [41], whereas we have m Z / √ g q g χ ∼ 3.2 TeV at benchmark point A and m Z / √ g q g χ ∼ 10 TeV at benchmark point B. Both benchmark points thus satisfy the limit.
We also need to consider cosmological constraints on the dark sector of our toy model. The most important question here is whether χ can be a thermal relic and account for all of the DM in the Universe. If there is no particle-antiparticle asymmetry in the dark sector, the answer is no.
The reason is that at our benchmark points the annihilation cross section [43] is much larger than the thermal relic cross section of few×10 −26 cm 3 /sec. Reducing the annihilation cross section would either require a significantly smaller α A , precluding significant dark radiation at the LHC, or a larger A mass, in particular m A > m χ . In the latter case, only the annihilation channelχχ →qq would remain open. The cross section for this channel is [44,45] σv χχ→qq where N f is the number of kinematically accessible quark flavors. Since σv χχ→qq scales as g 2 q g 2 χ m 2 χ /m 4 Z , it is much smaller than the thermal relic value for the small χ masses that we are interested in here. It would thus lead to overclosure of the Universe. For symmetric DM models, we therefore conclude that χ must be a subdominant component of the DM in the Universe 1 . This, of course, does not change the LHC phenomenology that we are interested in here, it only means that our toy model cannot provide a complete description of the dark sector. The fact that for symmetric DM models a discovery of χ through a dark radiation signature is only possible if χ is a subdominant component of DM is in line with the general statement that, in models with multiple thermally produced DM components, it is likely that the subdominant DM components are discovered at the LHC first due to their larger interaction strengths. For asymmetric DM, i.e. DM with a primordial particle-antiparticle asymmetry generated at some high scale, χ could account for all of the DM in the Universe. In this case, the large annihilation cross section to A pairs is not a problem but a feature because it allows for efficient annihilation of the symmetric component of the primordial DM soup. However another opposite charge particle might be needed due to gauge invariance for such light A [47]. Let us emphasize that, in this paper, we will indiscriminately call χ the dark matter, even though the above discussion shows that it does not need to be the dominant DM component.
Next, let us comment on constraints from direct and indirect DM searches. Without a primordial χ-χ asymmetry, both types of constraints are easily satisfied due to the required small relic abundance of χ. For asymmetric DM, indirect searches are also insensitive because χχ annihilation is impossible in such models. The sensitivity of direct detection experiments is limited due to the small m χ at our benchmark points. For benchmark point B with m χ = 0.4 GeV, the DM mass is obviously below the direct detection threshold, while for benchmark point A with m χ = 4 GeV, one might worry about limits from low-threshold experiments. The best limit in this mass range comes from CDMSlite and requires the spin-independent χ-nucleon scattering cross section σ χN to be below 1.5 × 10 −40 cm 2 if χ accounts for all of the DM in the Universe [48]. Our prediction at benchmark point A is σ χN 6.5 × 10 −42 cm 2 , avoiding the bound by more than an order of magnitude.
A different set of bounds on our model comes from the fact that χ particles have self-interactions, mediated by the A boson. Such DM self-interactions are constrained by observations of colliding galaxy cluster (most notably the Bullet Cluster) [49] and by measurements of the ellipticity of DM halos of groups of galaxies [50]. The resulting limit on the DM self-scattering cross-section is σ χχ /m χ 1 cm 2 /g = 1.78 × 10 −24 cm 2 /GeV .  Obviously, this constraint applies only if χ forms the bulk of the DM in the Universe. As explained above, this is only possible in our scenario if there is a particle-antiparticle asymmetry in the dark sector. In the perturbative regime, where α A m χ /m A 1, a simple calculation shows that the predicted non-relativistic self-scattering cross section is [5] σ χχ /m χ 5 × 10 −31 For our benchmark point A (B) from table I, this leads to σ χχ /m χ 10 −30 cm 2 /GeV (10 −29 cm 2 /GeV), well within the observational limit.
Regarding the dark photon mass m A and its kinetic mixing with the SM photon, strong constraints exist from low energy experiments and from astrophysics. These constraints require 10 −10 for m A 10 MeV, but relax to 10 −3 for 10 MeV m A 10 GeV, the mass range that we are interested in here. Our benchmark points from table I are thus consistent with limits from dark photon searches. We will discuss these limits in more detail in sec. IV, see in particular fig. 8.

C. Dark Photon decay
After dark photons have been radiated from pair-produced χ particles at the LHC, they decay to SM particles via their kinetic mixing with the photon (see eq. (1)). For m A 2 GeV and m A < 2m χ , the dominant decay modes include A → e + e − , µ + µ − , as well as decays to various hadronic resonances. The partial decay width to leptons of a specific flavor is The hadronic decay width can be calculated using measurements of hadron production cross sections at e + e − colliders, in particular the ratio R(s) ≡ σ(e + e − → hadrons)/σ(e + e − → µ + µ − ), where s is the center-of-mass energy. At s = m 2 A , the ratio between the partial A decay width into a given hadronic final state and the partial decay with to the µ + µ − final state should be equal to R(s): We use hadronic cross-sections from refs. [51,52] to calculate the A decay branching ratios. We include 19 decay channels, shown in fig. 2. In the calculation, care must be taken to avoid double counting of hadronic degrees of freedom. For example, the decay A → π + π − π 0 π 0 receives contributions both from the direct production of pions and from the decay A → ωπ 0 , followed by ω → π + π − π 0 . We choose to preferentially include channels with pions and kaons. For heavier mesons, we only include decay channels in which there is no double counting. In principle, there is also a double counting problem between kaon and pion final states because of the fully hadronic kaon decay channels. Fortunately, these channels have very small branching ratios and can be neglected.
For m A 2 GeV, hadronic A decays are more conveniently described in QCD language as decays to quark pairs. The partial decay width to a quark pair of a given flavor q f is where Q q f is the electric charge of q f and N c = 3 is a color factor. At intermediate A masses ∼ 2 GeV, a very large number of hadronic resonances appears, so that neither the QCD description in terms of quark final states nor the low energy effective theory involving just a few hadronic states are applicable. We therefore do not consider this mass range in the following, but we will see in sec. IV that our results on the LHC sensitivity can be smoothly interpolated between the QCD regime and the low-energy regime.

III. DARK PARTON SHOWER
Before presenting our numerical results, let us provide some analytic discussion of dark photon radiation. Part of the following discussion is based on refs. [53,54]. We first note that factorization theorems ensure that the cross section can be factorized into short range (hard process) and long range (shower) contributions [55]. Thus, we can treat the dark photon shower as independent of the primaryχχ production process. The relevant processes in the development of an A shower are χ → χ + A and A →χχ. In the following, we include only the first process, which is shown also in fig. 3. Its probability is much larger than the probability for A →χχ due to the structure of the respective splitting kernels. In particular, the splitting kernel for A →χχ is not divergent in the soft limit, in analogy to photon splitting in QED. The radiation of gauge bosons from fermions, on the other hand, is formally divergent in the soft and collinear phase space regions. In our scenario, the dark matter and dark photon masses regularize the radiation probability in these regions by providing a natural infrared cut-off. From a phenomenological point of view, we are more interested in collinear emission than in soft emission because the decay products of very soft dark photons would not be detectable at the LHC. We work under the assumption of strongly ordered emission, i.e. the virtuality of the incoming particle in a splitting process is assumed to be much larger than the virtuality of the outgoing particles. The latter are therefore treated as on-shell. Within the strongly ordered emission assumption, we can consider each splitting χ → χ + A as an isolated process, since a secondary splitting would be a small perturbation with respect to the first one.
The differential probability for a collinear splitting can be written as where t is the virtuality of the incoming DM particle and x is the fraction of its energy that is transferred to the outgoing DM particle. The quantity P χ→χ (x, t) is the splitting kernel that encodes the model-dependent details of the splitting probability. We call attention to the fact that we treat α A as independent of t, i.e. we neglect renormalization group running. Since the running is logarithmic, it will not play a major role in the development of the dark photon shower. Moreover, the running would depend on the full particle content of the dark sector and is thus highly model dependent. Eq. (10) can be understood as the squared matrix element of the splitting process, multiplied by the propagator of the incoming DM particle and integrated over the phase space of the outgoing particles. In the limit t → 0 (highly collinear splitting), the propagator factor 1/t in eq. (10) diverges as expected. It is prevented from going all the way to zero thanks to the non-zero dark photon mass.

A. Kinematics and notation
We write the off-shell 4-momentum of the incoming DM particle in a given splitting process as (see ref. [54] for a slightly different approach) Notice that in general E = √ŝ /2. For instance, if the second χ particle, which we call the spectator, is on-shell and does not radiate, energy momentum conservation requires its three-momentum to be (0, 0, −p). This fixes the energy of the spectator to be E 2 s = p 2 + m 2 χ , and hence E ≥ E s + m A if at least one dark photon is radiated in the event. In the following, we will nevertheless assume E = E s = √ŝ /2, the so called collinear approximation. Naturally, we expect it to fail when m A + m χ ∼ E, or when the opening angle between the splitting products is large, so the splitting is not collinear any more.
The on-shell 4-momenta of the outgoing DM particle and dark photon after the splitting are and respectively. Since the energies of the outgoing particles must be larger than their masses and because energies and momenta have to be positive, the following constraints must be fulfilled: The first relation defines in particular the allowed range of the variable x in the splitting probability eq. (10). Unfortunately, this range depends on E and thus on all preceding splitting processes.
This would preclude us from treating the dark photon shower analytically. Therefore, we will in the following take the minimum and maximum values of x to be where E 0 is the initial energy of the DM particle χ before emission of the first dark photon. Since most emitted dark photons are soft, we expect E to be close to E 0 for most splittings. The relations (15) would even be exact for m χ = m A = 0, but since we have m χ , m A E 0 , we expect them to be very good approximations also in our case.
The virtuality t is defined as and the splitting kernel is given by [56] We see that P χ→χ (x, t) becomes large for x → 1, where the emitted dark photon is very soft.
Since we are dealing with massive dark photons, x is kinematically required to always be < 1 (see eqs. (14) and (15)). In the following, we will neglect the mass dependent term in P χ→χ and call the thus approximated splitting kernel P χ→χ (x). Like the approximation used for x min and x max above, also this approximation is necessary to make different splittings completely independent of each other, a precondition for our analytic calculations. In order to also regularize the collinear t → 0 divergence in eq. (10), care must be taken also when choosing the integration limits for t. The lower bound t min for t is obtained when k t → 0, that is while the upper bound is obtained when k t is maximal, within the constraints (14), that is, with Interestingly, when x takes its maximum or minimum value allowed by eq. (14), k t,max becomes 0, which leads t min = t max and the splitting probability goes to 0. As in eq. (15), we will in the following take E → E 0 also when evaluating t min and t max . This approximation is necessary to make consecutive splitting processes completely independent and thus treat the dark photon shower analytically.

B. Number of emitted dark photons
A DM particle is produced at the LHC with a certain maximal virtuality t max and then radiates in general several dark photons. In each splitting, its virtuality is reduced, until it finally reaches the infrared cutoff t min . The expectation value for the number of radiated dark photons is Note that this compact expression is valid only if we neglect the t-dependence of the splitting kernel and of the integration boundaries x min , x max , t min , t max defined in the previous section. The probability that a DM particle χ radiates exactly m dark photons is given by a Poisson distribution [53]: To see this, note that eq. (10) implies that the probability for no splitting to occur between t max and t min is given by the Sudakov factor with n A defined in eq. (21). The probability for exactly one splitting is then In the first line, the factors ∆(t min , t) and ∆(t, t max ) give the probabilities for no splittings to happen in the intervals [t min , t) and (t, t max ], respectively. The splitting kernel P χ→χ (x) describes the probability that a splitting happens at virtuality t. For two splittings, we have in analogy In the second line, we have extended the integration region of the t integral to the full range [t min , t max ] and included a factor 1/2! to compensate for the double counting induced that way. This is the place where the approximations discussed in sec. III A are crucial: extending the integration interval is only possible if we can assume that x min , x max , t min and t max are the same for all splittings and that P χ→χ is independent of t. Eq. (25) straightforwardly generalizes to larger numbers of splittings, yielding the Poisson probabilities p m in eq. (22).

C. Recursion relation for the χ and A energy spectra
Let us now study the energy distributions of χ and A in events with radiation of multiple A bosons. Our starting point is the observation that, in the collinear limit, each splitting process is independent. As explained in sec. III A, we will make the crucial approximation that the integration boundaries x min , x max , t min , t max are the same for each emission process and are given by eqs. (15), (18) and (19) (with E → E 0 ). We moreover neglect again the t-dependence of the splitting kernel P χ→χ . These approximations will allow us to compute the energy distributions recursively.
If we call the final energy of the on-shell DM particle after final state radiation E χ and the energy fraction retained by the DM particle after all splittings X ≡ E χ /E 0 , 2 we can write the final DM energy distribution as where f χ,m (X) is the energy distribution of χ in events with exactly m emitted dark photons. f χ,m (X) obeys the recursion relation (see also [57]) where is the energy fraction for dark radiation showers with only a single splitting and Θ( · ) denotes a window function that is 1 when the condition in parentheses is satisfied and 0 otherwise. Note that for m = 1, the definitions of x as the fraction of energy the DM particle retains in a single splitting and of X as the fraction of energy the DM particle retains after all splittings are identical. For When evaluating f χ (X) numerically, we truncate the recursive series at m = 10, which is justified by the fact that the Poisson probability p m is very small at m n A . The distribution of the energy E A of radiated A bosons can be obtained in a similar way. With the notation Z = E A /E 0 , it is given by where f A ,k (Z) is the energy distribution of the k-th emitted A boson. Note that f A ,k (Z) is not a function of the total number of radiated dark photons m because each splitting process is independent. The recursion relation for f A ,k (Z) is

D. Energy spectra from an inverse Mellin transform
A particularly elegant way of computing the energy spectrum f χ (X) of dark matter particles and the spectrum f A (Z) of dark photons is by first computing the moments X s and Z s of these distributions, and then applying an inverse Mellin transform to recover f A (Z) itself. The Mellin transform of a function f (X) is defined as [58] and the inverse transform is given by The Mellin transform and its inverse are closely related to the Fourier and Laplace transforms [58] and can therefore be evaluated efficiently using the FFT (Fast Fourier Transform) algorithm [59]. For events in which a single A is emitted, the moments of the DM energy distribution after the emission, weighted by the probability of having exactly one emission, are given by In the last line, we have defined As in the previous section, we again take x min , x max , t min and t max to have the same value for each splitting process, and we neglect the t-dependence of P χ→χ . If two A bosons are emitted, we obtain analogously In the second line, we have again used the approximations from sec. III A to extend the integration region of the t integral (see the discussion below eq. (25)). Eq. (36) can be easily generalized to the case of m emitted A bosons: Summing over m, we find the moments of the overall A energy distribution: An inverse Mellin transform of ϕ(s) then yields the dark matter energy spectrum f χ (X). As can be seen from eq. (33), this requires evaluating ϕ(s) at complex values of s. A numerically stable way of doing this is by using the FFT algorithm to compute X s (eq. (35)), then using eq. (38), and afterwards apply another FFT to evaluate the inverse Mellin transform. The procedure for obtaining the A spectrum is alike and yields for the contribution of showers with m dark photons to the moment s, weighted by the probability of emitting exactly m dark photons, with Summing over all emissions we obtain  TABLE II: Characteristics of the dark photon shower for various choices of the DM mass m χ and the dark photon mass m A . The rows labeled "A" and "B" correspond to the two benchmark points from table I. In all cases, we assume χ pair production at a center of mass energy √ŝ = 1 TeV and we take α A = 0.2. We show the predicted number 2 n A of dark photons per event (with the factor of 2 coming from the fact that we consider χ pair production), the average fraction X that a χ particle retains of its initial energy after showering, and the average fraction Z of the initial χ energy that each dark photon receives. We compare the predictions of our semi-analytic calculations, eqs. (21), (38) and (41) (with s = 1), to the results of a Monte Carlo simulation in Pythia. As expected, the results satisfy the relation X + n A Z = 1.

E. Comparison of analytic results to Monte Carlo simulations
We have also simulated the radiation of A bosons fully numerically in Pythia 8 [35][36][37], using the "hidden valley" model implemented therein. Note that, like our analytic calculations, also the dark photon shower implemented in Pythia includes only dark photon radiation, χ → χA and χ →χA , not dark photon splitting, A →χχ.
To compare the χ and A energy spectra obtained using the different approaches, we assume χ pair production at a center of mass energy of √ŝ = 1 TeV. In Pythia, this is achieved by simulating the process e + e − → Z →χχ, followed by the dark photon shower. The results for several characteristic parameters of the dark parton shower, namely n A , X and Z , are shown in table II. Comparing Pythia results to the results obtained using eqs. (21), eq. (38) (with s = 1) and eq. (41) (with s = 1), we find excellent agreement. We attribute the remaining differences to the approximations we had to make in the analytical calculation, in particular the assumption that the upper and lower integration limits in x and t do not depend on x and t themselves.
A more detailed comparison between our analytic results and Pythia is shown in figs. 4 and 5. In fig. 4, we have plotted the distribution of the number of emitted dark photons per event (i.e. perχχ pair) at our two benchmark points from table I. Fig. 5 shows the energy spectra of χ and A after the dark photon shower has developed.
As we can see in fig. 4, Pythia results confirm that dark photon emission is approximately a Poisson process. We attribute the small deviations to the approximations we made in deriving eq. (22) in sec. III B.
Regarding the energy spectra plotted in fig. 5, we find excellent agreement between both of our analytic approaches (red dashed and blue dot-dashed curves) and also with Pythia results (black dashed curves). Note that the recursive method is numerically less stable than the Mellin transform approach because numerical errors-especially those incurred close to the upper integration boundary x max at which the splitting kernel has a regularized singularity-can accumulate in the recursive formulas (27) and (31). We have also carried out simple leading order simulations of the processes e + e − → Z →χχ and e + e − → Z →χχA in CalcHEP [60] (gray solid curves in fig. 5) and find that, as expected, this simulation reproduces our analytic result for the specific process e + e − → Z →χχA very well, but is insufficient as a proxy for the full dark photon spectrum. Since the full dependence of the matrix element on m χ and m A is taken into account in CalcHEP, it also illustrates that the approximations discussed in sec. III A are justified at this level. In k=1 p m f A ,k (Z) from events with a particular fixed number of dark photons to the overall energy spectra f χ (X) and f A (Z) (dashed rainbow-colored curves). Note that events with several emitted A bosons can dominate over events with a single A because 2 n A > 1 at our benchmark points. This illustrates that the A spectrum would be very difficult to obtain from a full matrix element calculation and necessitates a parton shower algorithm.
In summary, we have developed in this section two analytic methods for computing the energy spectrum of χ and A , and we have verified them by comparing to Monte Carlo simulations.

IV. COLLIDER SEARCHES FOR LEPTON JETS FROM DARK RADIATION
To constrain radiating DM models using existing LHC lepton jet searches, we numerically simulated the expected signal and background event spectra at center of mass energies of √ s = 7 TeV and √ s = 8 TeV. As in sec. III E, signal events are generated using the implementation of a "hidden valley" model in Pythia 8 [35][36][37]. We compare our simulation results to data from a 7 TeV ATLAS search for prompt lepton jets [15] and from an 8 TeV ATLAS search for displaced lepton jets [16]. We have also simulated events at √ s = 13 TeV to predict the sensitivity of future searches for prompt and displaced lepton jets.

A. Prompt lepton jets
Our prompt lepton jet analysis follows the 7 TeV, 5 fb −1 ATLAS search in ref. [15], but including only muonic lepton jets. The reason for discarding lepton jets with electrons is that it is very difficult to implement the corresponding selection criteria without a full detector simulation. The trigger requirement for muonic lepton jets is the presence of either three muons with transverse momenta p T > 6 GeV or one muon with p T > 18 GeV. Muons are required to have pseudorapidity |η| < 2.5, and they are required to have a track in the inner detector associated with them. We In all cases, we assume χ pair production at a center of mass energy √ŝ = 1 TeV. We compare the results from the recursion formulas described in sec. III C, the Mellin transform method developed in sec. III D, the dark photon shower simulation in Pythia, and a simple leading order simulation of e + e − → Z →χχA in CalcHEP. For the Mellin transform method, we also show the result separated according to the number of A bosons emitted in eachχχ pair production event. implement the last requirement by demanding that the parent A of the muon decays within a radius of 122.5 mm (corresponding to the position of the last silicon pixel layer) from the beam axis. Finally, muon tracks are required to have a transverse impact parameter |d 0 | < 1 mm with respect to the primary vertex.
We cluster muons into lepton jets using the following algorithm: we collect all muons with a cone of radius ∆R = 0.1 around the highest-p T muon in the event to form the first lepton jet. We then iterate this procedure, starting, in each iteration, with the highest-p T muon not associated with a lepton jet yet. Ref. [15] then defines two classes of events that are considered in the analysis: Double muon jet events require at least two muonic lepton jets, each of which must contain at least two muons with p T > 11 GeV. If the event was selected only by the single muon trigger, the leading muon in the lepton jet is in addition required to have p T > 23 GeV. Moreover, the two muons closest in p T within each lepton jet are required to have an invariant mass m µµ < 2 GeV.
To reject non-isolated lepton jets, the variable is defined, where the sum in the numerator runs over the E T of all calorimeter deposits within ∆R = 0.3 of any muon in the lepton jet, but excluding deposits within ∆R < 0.05 around any muon. The denominator is the p T of the lepton jet. The isolation cut imposes ρ < 0.3. Single muon jet events require only one muonic lepton jet, but this lepton jet must contain at least four muons with p T > 19 GeV, 16 GeV, 14 GeV for the three leading muons and p T > 4 GeV for all other muons. The same invariant mass cut as for the double muon jet analysis is imposed, and the isolation variable (42) has to satisfy ρ < 0.15.
Since we find that double muon jet events are much more sensitive to our radiating DM model than single muon jet events, we do not include the latter in the following analysis.
The dominant background to the muonic lepton jet search is from misidentified QCD multijet events, which contribute 0.5 ± 0.3 events in the signal region [15]. For the 7 TeV analysis, we can use this number at face value, while for the extrapolation to √ s = 13 TeV, we assume the background cross section after cuts to be a factor of 3 larger. This factor of 3 corresponds roughly to the scaling of the cross section for QCD jet production [61]. We assume the relative error on the background to remain the same as at 7 TeV. Note that, to be able to use this simple scaling of the background cross section, we have to use the same cuts at 7 TeV and at 13 TeV.

B. Displaced lepton jets
Our search for displaced lepton jets, based on the ATLAS 8 TeV, 20.3 fb −1 analysis [16], is somewhat more involved than the search for prompt lepton jets described in the previous section. First, we need to be more careful in simulating detector effects, which we do at particle level. We pay particular attention to the displaced A decay vertices, which lead to unusual detector signatures for each particle. The coordinates of the A decay vertex and the momenta of the A decay products are smeared by Gaussians with widths proportional to 1/L xy , where L xy is the transverse distance of the A decay vertex from the beam pipe. This way, we in particular take into account the fact that the momenta of charged particles that travel through more material are affected more strongly by the smearing. All parameters of our detector simulation are calibrated against the lepton jet reconstruction efficiencies shown in fig. 6 of ref. [16]. Since this plot is for a sample of dark photons with flat p T and η distributions in the range p T ∈ [10, 100] GeV and η ∈ [−2.5, 2.5], we generated a similar event sample for the comparison. To further improve the agreement between our predicted efficiencies and the results from ref. [16], we apply a fudge factor in each L xy bin.
Following [16], muons are selected from the event sample if they are produced inside the sensitive ATLAS detector volume, i.e. up to a radius of 7 m. Only loose requirements on the transverse momentum and pseudo-rapidity of reconstructed muons are applied: their p T must be above the trigger threshold of 6 GeV and the pseudorapidity |η| must be < 2.5. Since the goal of the search in ref. [16] was to constrain long-lived dark photons with displaced decay vertices, the background from prompt muons is reduced by requiring stand-alone muons, i.e. muons that are detected in the muon spectrometer but cannot be matched to a track in the inner detector. This essentially means that the dark photon decay in which the muon was produced must have happened outside the last pixel layer of the inner detector (ID), located at a radius of 122.5 mm from the beam pipe. To suppress cosmic ray muons, the reconstructed muon trajectory must still be associated with the primary interaction vertex. In particular, it must fulfill requirements on the impact parameters |z 0 | < 270 mm (longitudinal) and |d 0 | < 200 mm (transverse).
Hadronic jets (or, more precisely, calorimeter jets) are based on energy deposits in the electromagnetic calorimeters (ECAL) and hadronic calorimeters (HCAL), clustered using the anti-k T algorithm [62] implemented in FASTJET [63] with a cone radius of R = 0.4. Calorimeter jets must fulfill p T > 20 GeV and |η| < 2.5.
Three different types of lepton jets (LJ) are defined: A type-0 or "muonic" LJ arises from dark photon showers in which all A bosons decay to muons. It consists of at least two displaced and collimated muons with no additional activity nearby. If however, such muons are accompanied by a single calorimeter jet, we call the sum of these objects a type-1 or "mixed" LJ. A type-2 LJ or "calorimeter" LJ is defined as an isolated calorimeter jet that fulfills additional requirements to distinguish it from QCD background. Note that, especially for type-1 and type-2 LJs, the term lepton jet used in [16] is somewhat misleading because these objects can also arise from hadronic A decays. The clustering algorithm for all three LJ types is described in the following.
We collect all muons and calorimeter jets inside a cone of radius ∆R = 0.5 around the highestp T muon in the event. If at least one other muon and no calorimeter jet is found inside this cone, those muons are combined to a muonic LJ. If at least one muon and exactly one calorimeter jet is found inside the cone, all objects are combined into a mixed LJ. If no second muon is found within ∆R = 0.5 of the highest-p T muon, the muon is discarded. The algorithm then continues with the highest-p T muon that has not been associated with a lepton jet or discarded yet. After all muons have been processed, the remaining calorimeter jets are reconstructed as calorimeter LJs if their electromagnetic (EM) fraction is lower than 0.1. This is for example true for dark photons decaying to charged pions, but also for a decays to electrons happening inside the hadronic calorimeter. Since the pseudo-rapidity region 1.0 < |η| < 1.4 tends to create a fake low EM fraction due to the transition between the barrel and endcap EM calorimeters (ECAL), this region is excluded. Additionally the width W of a calorimeter LJ, defined by needs to be small. Here, the sum runs over all particles in the lepton jet. An isolation cut is imposed on all reconstructed lepton jets by requiring that the scalar p T sum of all charged tracks, reconstructed by the inner detector inside a cone of radius 0.5 around the lepton jet axis, is below 3 GeV. Only charged tracks with p T > 400 MeV and with impact parameters |z 0 | < 10 mm and |d 0 | < 10 mm are considered in this procedure. An event is discarded if there are not exactly two lepton jets fulfilling |∆φ| > 1.0. The requirement of two lepton jets leads to the 6 different combinations 0-0, 0-1, 0-2, 1-1, 1-2, and 2-2, where e.g. 0-1 corresponds to an event reconstructed with one type-0 (muonic) and one type-1 (mixed) lepton jet. The sensitivity of the displaced analysis will depend strongly on the A decay mode. Therefore, we list in table III the most important A decay modes and their associated signatures. Check marks ( ) indicate in which region of the detector a particular decay has to happen in order to potentially be reconstructed as a displaced lepton jet. For decays that would fail our cuts, we give the reason why they are vetoed. For instance, due to the requirement of a low EM fraction, A → e + e − is only identified as a displaced lepton jet if the decay happens inside the hadronic calorimeter, while A → µ + µ − is also selected if the decay happens inside the electromagnetic calorimeter. Decays to any final states involving charged particles are vetoed if they happen in the inner detector. Regarding decays to mesons, it is important to note that π ± , K ± and K 0 L have lifetimes of order 10 −8 sec [64] and are thus stable on detector scales (i.e. they will be stopped in the hadronic calorimeter before they decay). Neutral pions, on the other hand, decay very fast. For K 0 S , the lifetime of 9×10 −11 sec may allow it to travel for several cm before decaying (mostly to π + π − or π 0 π 0 ). Therefore, A → π + π − , K + K − are accepted if they happen in the electromagnetic or hadronic calorimeter. The decay A → π + π − π 0 , on the other hand is vetoed if it happens in the electromagnetic calorimeter because the quasi-instantaneous decay of the π 0 to two photons  TABLE III: Illustration of where in the detector a specific A decay must happen in order to potentially be reconstructed as a lepton jet in our displaced lepton jet analysis (see sec. IV B and ref. [16]. For decays that will be vetoed, a reason for the veto is given. The type of lepton jet as which each decay mode is most likely to be reconstructed is given at the top of the table. deposits a large amount of energy, violating the cut on the EM fraction. The decay A → K 0 L K 0 S is successfully reconstructed as a displaced lepton jet in many, but not all, cases. It will be vetoed in the inner detector if both the A and the secondary K 0 S decay very fast (well within the 12.25 cm radius of the inner detector) in spite of their typically large boosts, and the K 0 S decay mode is to π + π − (branching ratio 69%). It will be vetoed by the cut on the EM fraction if the A decay happens in the inner detector or the electromagnetic calorimeter and the K 0 S decay mode is to π 0 π 0 (branching ratio 31%).
Note that, as can be read from table III, the decay mode A → µ + µ − is most likely to be reconstructed as a type-0 (muonic) lepton jet, while all other decay modes typically lead to type-2 (calorimeter) lepton jets. Type-1 (mixed) lepton jets are obtained only if several A are radiated from the same dark matter particle, with at least one of them decaying to µ + µ − and at least one of them decaying to a different final state. Therefore, the requirement for a type-1 (mixed) lepton jet are the most difficult to satisfy. The predicted signal event numbers in table. V confirm this.

C. Results
In tables IV and V we present the signal and background predictions as well as the observed event numbers for the two benchmark models defined in table I. Table IV is for the prompt lepton jet analysis from ref. [15], using 5 fb −1 of 7 TeV ATLAS data (see sec. IV A, while table V is for the displaced lepton jet search from ref. [16], using 20.3 fb −1 of 8 TeV ATLAS data (see sec. IV B). Note that, for the prompt search, we include only muonic lepton jets, which are easier to model without a full detector simulation. For the displaced search, we show event numbers with and without inclusion of the background-limited sample of events with two type-2 lepton jets. We see from tables IV and V that, as expected, benchmark point A is more easily detectable in the displaced search, while benchmark point B is probed more sensitively by the prompt search.
We also show our predictions for 100 fb −1 of 13 TeV data. In general, we expect better sensitivity at 13 TeV because of the larger integrated luminosity expected. Since an estimation of the multijet background at 13 TeV is problematic, we restrict the displaced lepton jet search at 13 TeV to type-0 (muonic) LJs, for which QCD backgrounds do not play a major role. The only relevant source of background in the 13 TeV displaced search are then cosmic rays, whose flux is independent of the collider center of mass energy as long as we use the same cuts as at √ s = 8 TeV. Our predictions for the sensitivity of a displaced lepton jet search at 13 TeV LHC are thus very conservative. For the prompt search at 13 TeV, we scale the background cross section from the 7 TeV analysis by a factor of 3.
We use the CLs method [65] to compare our predictions at √ s = 7 TeV and √ s = 8 TeV to the data and to constrain the production cross section for a DM pair, σ(pp → Z )BR(Z →χχ). We use the ATLAS background predictions [15,16] as the null hypothesis, and the sum of the ATLAS    [16]) at the benchmark parameter points from table I, compared to the background predictions and the observed event rates from ref. [16]. In the last two columns, the first error is the statistical uncertainty, while the second one is systematic. Our sensitivity study at √ s = 13 TeV includes only type 0-0 events (only muonic lepton jets) because a reliable extrapolation of the multijet background to 13 TeV is difficult.
background and our signal prediction as the test hypothesis. For √ s = 13 TeV, we compare our signal+background prediction to the background-only prediction to obtain the expected limit. We assume a systematic uncertainty of 30% on the signal normalization to account for errors in the signal and detector simulation. The results are presented in figs. 6 and 7 for the two benchmark points from table I, where in each panel one of the model parameters is varied while the others are kept constant. For the displaced search at √ s = 8 TeV, two different limits are calculated, one including all tagged events (solid red curves) and one excluding type 2-2 events, i.e. events with two type-2 lepton jets (solid black curves). Due to the larger background in the 2-2 sample, this second limit is stronger in certain parameter regimes than the limit including all lepton jets. For the displaced search at √ s = 13 TeV, (red dotted curves), we use only events with two type-0 lepton jets, as explained above. We also show in figs. 6 and 7 the predicted values for σ(pp → Z )BR(Z →χχ) as thin dotted curves.
Let us first discuss fig. 6, corresponding to benchmark point A. (We will comment on the differences compared to fig. 7 below.) Since for fixed m A , the dark photon decay length cτ depends only on the kinetic mixing parameter , figs. 6 (a) and (b) are equivalent. We see that the best limits from the prompt search are obtained at cτ 1 mm ( 10 −5 ), while for larger decay length, the cut on the impact parameter |d 0 | < 1 mm restricts the sensitivity. For the displaced search, the sensitivity is optimal around cτ = 10-100 mm, corresponding to few × 10 −6 . For shorter lifetimes (larger ) the dark photon decay vertices tend to be closer to the primary interaction  [15] (solid blue) and from the 8 TeV ATLAS search for displaced lepton jets [16] are shown. For the latter search, we show results including all lepton jet events (red solid) and excluding events with two type-2 lepton jets (black solid). The predicted sensitivity of similar analyses at √ s = 13 TeV is shown as blue/red dotted curves. The black dotted lines in each panel show the theoretically predicted production cross sections. point, so that the A decay products are more likely to leave tracks in the inner detector, thus failing the displaced lepton jet selection criteria. For too large lifetimes (smaller ), on the other hand, most dark photons will decay outside the ATLAS detector and not lead to a signal at all. The exact position of the most sensitive regime depends strongly on the A kinematics and its mass.
For the displaced search, notice also the slight shift of the most sensitive point in cτ and between the red solid, black solid and red dotted curves in figs. 6 (a) and (b). This is related mostly to the inclusion/exclusion of different types of lepton jets in the three cases. In particular, as explained above, we include only type-0 (muonic) lepton jets in the 13 TeV analysis. The decay A → µ + µ − can be successfully reconstructed for shorter A decay lengths than many other A decay modes (see table III) Consequently, the exclusion limit based on displaced type 0-0 events only (red dotted curves in fig. 6 (a) in (b)) is best at relatively small cτ (relatively large ), followed by the sample excluding type 2-2 events and by the sample including all types of lepton jets.
Similar arguments help us understand the exclusion limits as a function of the A mass in fig. 6  (c). The different peaks and dips in this plot reflect the dependence of the A decay branching ratios on m A illustrated in fig. 2, see also [64]. The most notable features arise due to the broad ρ 0 resonance around 780 MeV with the narrow ω resonance on top of it and the narrow φ resonance at 1 GeV. The ρ 0 resonance decays almost exclusively to π + π − pairs, so that around this resonance, the relative importance of all other decay channels is diminished. This implies that the prompt lepton jet search, which includes only muonic lepton jets, has a greatly reduced sensitivity around the ρ resonance. In the displaced search, most A → ρ 0 → π + π − decays are reconstructed as type 2 (calorimeter) lepton jets (see table III), therefore the sensitivity improves around the ρ 0 resonance if type 2-2 events are included in the analysis (red solid curve in fig. 6 (c)). If type 2-2 events are excluded (black solid and red dotted curves in fig. 6 (c)), the sensitivity decreases. The dominant decay of the ω resonance is to π + π − π 0 . This again decreases the sensitivity of the prompt search around the resonance. For the displaced search, we can read from table III that the decay A → ω → π + π − π 0 is more likely to be vetoed than A → π + π − , therefore the sensitivity is reduced at the ω resonance compared to the off-resonance region. The behavior at the φ resonance is similar to that at the ρ 0 resonance. The prompt search is insensitive because the φ mostly decays hadronically. For the displaced search, detection prospects for the dominant decay channel φ → K + K − (branching ratio 48.9%) as a type-2 (calorimeter) lepton jet are very similar to those of the ρ 0 resonance decaying to π + π − (see table III). The secondary decay φ → K 0 L K 0 S (branching ratio 34.2%) has a high probability of being reconstructed as a displaced type-2 lepton jet even if the A decay happens in the inner detector. For these reasons, the sensitivity at the φ resonance is further boosted in analyses including displaced type-2 lepton jets. For A masses above 1 GeV, the sensitivities of both the prompt and the displaced searches decrease rapidly. For the prompt search, the invariant mass requirement m µµ < 2 GeV restricts the realm of sensitivity to m A < 2 GeV. For the displaced search, the smaller A boost factor at larger m A implies that the A decay length becomes shorter and events are more likely to leave tracks in the inner detector and be vetoed. Moreover, the average number of radiated A decreases at larger m A . As discussed in sec. II, we do not consider the region 1.7 GeV < m A < 2.3 GeV, where numerous hadronic resonances appear but the quark model is not applicable yet. In fig.6 (c), we have smoothly interpolated through this region. As can be seen, the low-energy description in terms of hadron final states and the high-energy description in terms of quark final states match onto each other very well.
In fig. 6 (d), we find that, unsurprisingly, all limits improve as the dark sector gauge coupling α A increases because the probability for A radiation grows. Note, however, that in the very large α A region, our perturbative treatment of the dark photon shower is expected to break down.
production of the Z becomes impossible and the sensitivity decreases again. Finally, fig. 6 (f) shows that the sensitivity goes down when χ is made heavier because less boosted χ particles radiate less.
Let us now turn to a comparison of the two benchmark models from table I and compare figs. 6 and 7. The general features of the two figures are very similar, but in general benchmark point B is more easily detectable in prompt lepton jet searches due to the smaller cτ . Benchmark point B offers somewhat better sensitivity than benchmark point A also because the smaller values of m χ and m A imply that the average number of A radiated in each event is larger. This effect is especially pronounced in the displaced searches excluding type 2-2 events, which are limited by backgrounds. Note that in fig. 7 (c) the region with m A > 2m χ is not considered because it would lead to a very large branching ratio for A →χχ, making the A invisible to the detector.
To put our results in the context of other constraints on dark photons, we show in fig. 8 the dark photon parameter space spanned by and m A . For both of our benchmark points, we compare the limits we derived from the prompt and displaced ATLAS lepton jet searches (blue shaded/red shaded) and the predictions for 13 TeV (blue/red unshaded) to the exclusion limits from various low energy experiments. For the 8 TeV displaced analysis, we decide for each parameter space point individually whether or not to exclude type 2-2 events (events with two calorimeter lepton jets), depending on which analysis leads to the better expected limit.
We see that the parameter region probed by LHC lepton jet searches is complementary to the region probed by low energy searches. It is important to keep in mind, though, that low energy experiments probe any dark dark photon model, while our analysis is sensitive only to scenarios that, in addition, feature a light DM particle that can be pair-produced at the LHC. We see that benchmark point A is already excluded by the 8 TeV ATLAS search for displaced lepton jets, while benchmark point B can be probed by both the prompt and the displaced search at 13 TeV. The exclusion regions from the displaced searches move to smaller values of at larger m A because the A width increases and the lab frame decay length decreases with increasing m A due to the smaller boost factors. Both effects need to be compensated by a decrease in to avoid A decays in the inner detector. The impact of the various A decay channels on the sensitivity is again reflected by spikes, dips and kinks in the excluded regions. For instance, note that the prompt search at 7 TeV, which includes only muonic lepton jets, is insensitive close to the hadronically decaying ρ, ω and φ resonances. At 13 TeV, the gap is closed thanks to the much better statistics. The sharp edge at low A mass for the 7 TeV and 13 TeV regions corresponds to the opening of the decay channel A → µ + µ − at m A = 2m µ . Remember that, in our 7 TeV and 13 TeV analysis, we could only include muonic lepton jets, therefore, our analysis is insensitive below the muon threshold. A full experimental search would of course extend also into this region, like the 8 TeV ATLAS search for displaced lepton jets has done.

V. CONCLUSIONS
In summary, we have studied the possibility of revealing the properties of a dark sector of particle physics by using final state radiation from dark matter produced at the LHC. The characteristic experimental signature of this process is a pair of lepton jets. While our conclusions apply to a very large class of extended dark sector models, we have worked in the framework of a toy model where fermionic dark matter particles χ are charged under a new dark sector gauge group U (1) and coupled to the SM through a heavy mediator. When χ particles are pair-produced at the LHC via this heavy mediator, they can radiate several light U (1) gauge bosons A (dark photons) which subsequently decay to light SM particles (electrons, muons, mesons) through a small kinetic mixing with the SM photon. Due to the required smallness of the kinetic mixing, the A decay length can be macroscopic. We emphasize that our results for this toy model are easily generalized to any model with light dark sector particles charged under a new gauge interaction. We also remind the reader that, in order to account for all of the DM in the Universe, χ must be produced non-thermally, as for instance in asymmetric DM scenarios.
We have first developed two analytic treatments of dark photon radiation: in the first one, we use recursive expressions for the A and χ energy distributions, while in the second one, we compute the moments of these distributions fully analytically and then apply an inverse Mellin transform to obtain the distributions themselves. We have compared our analytic calculations with Monte Carlo simulations in Pythia, finding excellent agreement.
In the second part of the paper, we have extended these Monte Carlo simulations by a simplified description of the ATLAS detector that allows us to recast the existing ATLAS searches for prompt and displaced lepton jets into powerful limits on the DM pair production cross section and the dark photon parameters in radiating DM models. Our limits on theχχ production cross section range down to O(10 fb), depending on the A mass and lifetime, the U (1) gauge coupling, and the mass of χ. Regarding the A properties, we find that LHC searches for radiating DM can probe a region in m A -space that has been inaccessible to low energy dark photon searches so far, namely in the parameter range 10 −7 10 −3 for m A in the MeV-GeV range. Of course, these limits are subject to the condition that DM particles with large coupling to dark photons (α A ∼ O(0.1)) can be pair-produced at the LHC in significant numbers.
Looking into the future, we have also shown that LHC limits will significantly improve in Run 2 at 13 TeV center of mass energy. It is worth emphasizing here that our simplified 13 TeV analyses are still a far cry from what a full experimental search can achieve. Most importantly, with a more detailed detector simulation and data-driven background estimation methods, it will be possible to include not only muonic lepton jets, but also A decays to electrons and hadrons. To conclude, we have shown that lepton jets are an interesting and powerful tool to elucidate the dynamics of the dark matter sector. Together with searches for hadronic jets with non-standard properties, they will be essential in the hunt for dark matter at Run 2 of the LHC and might well be the first to find a positive signal.