Dark Matter Spectra from the Electroweak to the Planck Scale

We compute the decay spectrum for dark matter (DM) with masses above the scale of electroweak symmetry breaking, all the way to the Planck scale. For an arbitrary hard process involving a decay to the unbroken standard model, we determine the prompt distribution of stable states including photons, neutrinos, positrons, and antiprotons. These spectra are a crucial ingredient in the search for DM via indirect detection at the highest energies as being probed in current and upcoming experiments including IceCube, HAWC, CTA, and LHAASO. Our approach improves considerably on existing methods. For example, we include all relevant electroweak interactions. The importance of these effects grow with DM mass, and by an EeV our spectra can differ by orders of magnitude from existing results.


Introduction
If the dark matter (DM) of our universe is a particle with a mass between the electroweak and Planck scales, then it could be discovered via indirect detection of standard model (SM) particles produced from its decay. Such decays can be initiated by an underlying hard process where the DM decays to two SM states, χ → XX. The SM states, injected with virtuality µ ∼ m χ , will shower and eventually hadronize, evolving down to on-shell stable particles such as photons, neutrinos, positrons, and anti-protons.  Figure 1. The prompt electron neutrino and photon spectrum resulting from the decay of a 2 EeV DM particle to ν eνe , as currently being searched for at IceCube [5]. Solid curves represent the results of this work, and predict orders of magnitude more flux at certain energies than the dashed results of Pythia 8.2, one of the only existing methods to generate spectra at these masses. In both cases energy conservation is satisfied: there is a considerable contribution to a δ-function at x = 1, associated with events where an initial W or Z was never emitted and thus no subsequent shower developed. Large disagreements are generically observed at these masses for electroweak dominated channels, while the agreement is better for colored initial SM states.
Calculation of the resulting prompt spectra is a central ingredient in testing the hypothesis of heavy DM. At present, a common approach is to simulate these events using Pythia [1][2][3], which accurately reproduces most of the relevant physics up to ∼TeV scales. Pythia is not, however, at present designed to operate well above these scales, for example it is missing interactions such as triple gauge couplings in the electroweak sector that can become increasingly important. In this paper, we propose an alternative approach, which in certain channels can produce spectra that differ significantly from existing results, as demonstrated in Fig. 1. We make our full results publicly available [4].
The remainder of this paper presents our approach to obtaining these. We begin by describing how the calculation of DM spectra can be mapped onto fragmentation functions (FFs), which can be evolved from the UV scale, µ ∼ m χ , down to the IR, µ ∼ 0. The computation can be performed in three stages: 1) Evolution from m χ to the weak scale q W ∼ 100 GeV; 2) Matching through q W ; and 3) Continued evolution down to 0. In the main text we outline the methods used at each step, providing full technical details in related appendices.

Framework
The flux of an observable particle S produced from DM decay depends centrally on the prompt spectrum, defined as 1 Here Γ is the inclusive decay rate of χ to S, Γ 0 = 1/τ is the inverse lifetime, and we use dimensionless variables x = 2E/m χ . If the decay is seeded by an underlying process χ → XX, for an arbitrary SM state X, the process begins with each particle at a virtuality scale m χ /2. The problem is then to determine the probability that X andX evolve to produce S carrying a fraction x of the initial energy. This process is described by a fragmentation function (FF) D b a (x; µ Q , µ 0 ), which determines the probability of an initial particle a at a scale µ Q evolving to produce a particle b at µ 0 carrying a momentum fraction x; in the absence of any evolution we would have D b a (x; µ Q , µ 0 ) = δ b a δ(1 − x). In this language, we can write the spectrum as 2 dN S dx = D S X (x; m χ /2, 0) + D S X (x; m χ /2, 0) .

(2.2)
At this stage, we have simply rephrased the problem. The power of Eq. (2.2) is that it allows us to bring to bear the considerable formalism of FFs to the calculation of DM spectra. In particular, the full evolution in virtuality can be decomposed into easier to compute segments, and then convolved together. For the present work we will exploit this result to break the calculation up as follows, The discussion is couched in the language of DM decay due to the Kamionkowski-Griest bound [54] providing a naive obstruction to DM annihilation at these masses. The bound can be evaded, see e.g. [55][56][57][58], and our results can be readily ported to annihilation with the simple identification m dec. χ = 2m ann.
χ . 2 Eq. (2.2) applies for a hard two-body decay. The formalism can be extended to (n > 2)-body decays, as described in the App. A. Figure 2. A cartoon of the three steps used to calculate the UV to IR evolution. An initial 2-body final state resulting from a hard interaction of χ → XX is DGLAP evolved down to q + W , just above the weak scale. At q + W , states with masses above the electroweak scale are integrated out. Finally, at q − W these results are matched onto Pythia which handles the subsequent evolution and hadronization effects.
The three pieces to be calculated are as follows. Firstly we evolve from the scale of the DM mass down to just above the weak scale, q + W , using the DGLAP equations [59][60][61] and in particular an implementation using all interactions in the unbroken SM, as well as a partial treatment of soft-coherence effects [62][63][64][65][66]. We next perform a matching by evolving across a parametrically small region through the weak scale, removing all particles with electroweak scale masses. Finally, these results are matched onto Pythia below q W , where it is used to calculate the subsequent showering, hadronization, and light particle decays in a regime where it has been extensively vetted. A simplified depiction of the full evolution is given in Fig. 2, and we next flesh out the details involved at each stage. associated with the leading collinear and collinear-soft divergences in the theory, both of which are described by the unregulated Altarelli-Parisi splitting functionsP (z). The evolution of the FFs under the 1 → 2 splitting interaction I encoded inP (z) is described by the DGLAP evolution equations, which take the schematic form (suppressing the fixed high scale) where α I is the corresponding SM coupling. The distribution of momenta amongst the particles, described by the FF D(x, µ), can only evolve in two ways: contributions are received by particles with momentum fractions greater than x splitting to exactly that value, and they are lost if a particle with fraction x splits at all. This describes the two terms in square brackets above, and they are associated with real and virtual emissions in the theory respectively. The evolution down to a desired scale is achieved by solving Eq.  [67]. To solve this problem, we use the results of Refs. [68][69][70], which included a treatment of polarization effects. These polarization effects are critical. Even if in the UV we start with an unpolarized initial state, the chiral nature of the electroweak interaction will generate a polarization, which can be considerable over the evolution scales we consider here. As the spectra of particles produced by decays of the electroweak states t, W , and Z depend on their polarizations, we will need to keep track of this through to the matching step described below.
The DGLAP evolution can effectively be described as a semi-classical shower that develops via consecutive 1 → 2 splittings, each occuring with probability determined by the appropriateP (z). Soft physics breaks this picture. As an example, consider a splitting W → dū, and then the subsequent emission of a gluon off either theū or d u d g Wū d g W If the gluon is emitted at a wide angle with respect to the dū pair, then in the soft limit, its wavelength will be such that it cannot resolve the individual quarks. As they came from an uncolored W , the gluon will now see 0 net color charge and cannot be emitted. In particular, there is destructive interference between the two diagrams, and consequently a large suppression of small-x states. The result is the well known angular-ordering effect, and in the context of QCD Monte Carlo simulations, it is understood how to include this effect [71,72]. As a modification to the full SM DGLAP equations, it is not. A complete treatment of the problem is beyond the scope of the present work, however, we introduce an identity in order to capture the largest effect of soft coherence: the reduction of real radiation at small-x. As we derive in the App. B, the equation that describes the leading effect of angular ordering, can be rewritten as a DGLAP equation, but evolved from a scale x × m χ /2, rather than m χ /2. This allows us to take our solution to the full form of Eq. (3.1), and augment them with the substitution m χ → x m χ . 4 This substitution allows for a simple inclusion of the soft physics, but it is not perfect. Soft-coherence not only reduces the real emission, it also increases the associated virtual noemission probability, and our result only accounts for the former. This deficiency manifests itself as a failure of momentum conservation generally at the level of ∼1 − 3%, although for particular states and masses it can be as large as 10%. Given the large impact of the effect on the spectra at small-x, we choose to accept this shortcoming, leaving the complete treatment as an open problem.

Weak Matching
We now take our DGLAP evolved, soft-coherence corrected, FFs and evolve them across the electroweak threshold. Formally we evolve across q ± W = q W (1± ), with 1, a parametrically small separation of scales, ensuring this step cannot generate large logs from the evolution. Instead, the point of this step is a matching from the unbroken to broken SM, where we integrate out the electroweak mass states t, W , Z, and h. For h we let Pythia handle the decay. For t, W , and Z, we instead need to account for the fact that our evolution at the first step can generate a significantly polarized spectrum for each state. In order to ensure this physics persists into our final results, we decay each of these states analytically, accounting for the polarization.
For our purposes, the details of the polarized decays are sufficiently described by the differential spectra obtained from the tree level diagrams, such as depicted below.
The calculations are straightforward, although for the top slightly involved, and so we postpone the full details to the appendices. In each case, we obtain an analytic or simple parameterization for D f W 0,± (x; q + W , q − W ), where f represents the fermions that can be produced from a W decay, and the equivalent functions for Z 0,± and t ± , where 0, ± correspond to the longitudinal and ± transverse polarizations, respectively. Full details are given in App. C. 4 Technically this substitution is only appropriate for the single logarithmic terms associated with the isosinglet evolution. As such we need to factor out the non-cancelling electroweak double logs from this change of variables. See App. B for details.
For states without m ∼ q W , other than the contribution they receive from decaying particles, the threshold is uneventful. We remove polarization at this stage by, for example, combining f L + f R = f , because its effects in evolution below the electroweak scale are small and Pythia does not take account of them. This step also represents the first appearance of electroweak masses, which are neglected in the initial evolution. For m χ ∼ q W , we will accordingly underestimate the phase space suppression of electroweak states, and therefore do not quote results for m χ < TeV. There are many options available at these scales, including Pythia or PPPC4DMID [73,74]. Another possibility would be to match our results to fixedorder matrix elements without double counting, in the way proposed in Ref. [75]. 5

Low Scale Evolution with Pythia
Our final step is to take each particle that resulted from the weak matching, and continue its evolution down to lower scales using Pythia, which will include the remaining showering, soft coherence, light particle decays, and importantly the non-perturbative hadronization. For each state, starting at a scale q W , the spectra of stable states, S ∈ {γ, e ± , p ± , ν e,µ,τ ,ν e,µ,τ }, are determined. When convolved with the earlier steps the full DM result is obtained.
For our purposes there is one deficiency with a full Pythia treatment at this step. Pythia models the final state radiation (FSR) emission of photons off charged particles, f → f γ, only down to an isolation or p T cut. Photons that are highly collinear with the parent charged particle cannot be separated in the environment of a collider like the LHC. If they travel over galactic or cosmological distances, they certainly can however. As such, we want to include the photons that Pythia deliberately excludes at small-x. To do so, we turn photon FSR off in Pythia and instead include it analytically to first order using the appropriate result given by, Here e N is the charge of the emitter, and m N its mass. If N indexes a quark, we take m q = max(m q , Λ QCD ), as the evolution will stop at the higher of the two scales. This treatment means we neglect processes like subsequent splitting of the FSR photons into charged fermions but as these are formally higher order and not enhanced by large logs, the above treatment is sufficient for our purposes.
In addition, we also need to account for the fact that turning off photon FSR in Pythia means the charged particles that would have emitted photons now have a higher momentum fraction. In other words, including only Eq. (5.1) is inconsistent with momentum conservation. However, it is straightforward to derive the required modifications to the charged particle FFs, simply from momentum conservation, and these are applied. The full expressions are provided in App. D.  Figure 3. The photon spectrum resulting from χ → bb, shown across the entire mass range considered in this work. A curve is shown for every decade in mass between a TeV and the Planck scale.

Discussion
Combining these three FFs, with the aid of Eqs. (2.2,2.3), we obtain the desired prompt DM spectra. 6 An example is shown in Fig. 1. The difference between our results and those of Pythia is driven by the lack of the full electroweak interactions in the latter. In Pythia, Z emission from the hard neutrino leads to an endpoint contribution to the neutrino spectrum, whilst the low energy bumps are associated with QCD states. However, a full electroweak shower cannot form in Pythia as it can in our results, and by m χ ∼ EeV, the shower can involve a large number of electroweak bosons that produce significantly more emission across the allowable energy fractions. For the same reason, in other channels where electroweak showers contribute significantly, such as χ → e − e + , large disagreements are also observed. Smaller differences occur for channels where QCD showers dominate. Yet in all cases, for m χ well above the electroweak scale, running Pythia will eventually become impractical given the size of the generated showers. Our approach encounters no such obstacle, as exemplified in Fig. 3 where we show the photon spectrum from χ → bb for DM masses all the way to the Planck scale. The spectrum of hard photons is seen to increase with m χ , whilst the lower energy peak associated predominantly with neutral pion decays progressively softens. These are only a few representative examples. In total we compute 638 FFs D S X : X indexes the 58 states in the unbroken SM, and S ∈ {γ, e ± , p ± , ν e,µ,τ ,ν e,µ,τ }. The full spectra are available [4]. 7 An exploration of the impact these results have on searches for DM with neutrino final states is provided in Ref. [38].
Our results are not the final word on heavy DM spectra. There are a number of directions in which our results can be systematically improved, including: • A rigorous treatment of the relevant soft physics as part of a full next-to-leadinglogarithmic (NLL) calculation.
• Matching to fixed-order electroweak matrix elements using the method of Ref. [75].
• Inclusion of estimated theoretical uncertainties associated with both the DGLAP and parton shower evolution, for the latter see Refs. [80][81][82].
• Polarized decays of the µ and τ .
• We assume no additional thresholds are crossed between m χ and q W . Results going beyond this have been considered for supersymmetric QCD [83][84][85] and also the minimal supersymmetric SM [86,87].
At present, our spectra carry quantifiable theoretical uncertainties. Our calculation is performed with double logarithmic accuracy, and therefore single logs are not resummed. Accordingly, our results carry a global uncertainty of order exp[α I L]/(1 + α I L) with L = ln(m χ /q W ) and I = 2, 3. At an EeV this corresponds to an O(10%) error, growing to O(20%) by the Planck scale. For x 10 −3 there are even larger uncertainties associated with an incomplete treatment of the soft physics associated with coherent emission of gauge bosons. We can obtain an uncertainty estimate by comparing our default procedure of implementing the m χ → x m χ substitution at the high scale to two alternative approaches: 1. an underestimate of soft-coherence, where no accounting for the effect is performed; and 2. an overestimate where we apply the substitution to the FFs combined across all three scales, which will double count the soft-coherence already present in Pythia. Doing so, we conclude that our treatment can produce O(1) errors at small-x values for an EeV, although the effect at a given x decreases with mass. A fuller discussion of uncertainties is given in App. E.
Even with these uncertainties, our spectra represent a manifest improvement over existing treatments: we include effects that are demonstrably important, and our formalism extends all SM states to arbitrarily high masses. Combining these results with astrophysical probes of the high energy universe, the heavy DM hypothesis will be put to the test in the coming years. In the event of an excess, we may finally begin to unravel the particle nature of DM.

Appendix
As mentioned, we have postponed many of the technical aspects of our calculation to the appendices, and we will organize that discussion as follows. In App. A we begin by providing an expanded discussion of the connection between DM spectra and FFs, and a comparison between the approach to the problem presented in this work with ideas discussed previously. In Apps. B, C, and D we provide an unabridged discussion of the three steps of our evolution. Afterwards, in App. E we provide an estimate of the theoretical uncertainties associated with our spectra. Appendix F contains additional results highlighting the physics inherent in our spectra, and finally App. G outlines the details of the public spectra and code.

A Dark Matter Spectra and Fragmentation Functions
This section expands upon the connection between DM spectra and FFs. In particular, we detail the steps between Eq. (2.1) and Eq. (2.2), as well as outlining the corresponding result for DM annihilation. In doing so, it will become clear what steps are needed to modify our approach when the hard interactions is (n > 2)-body.
To begin with, let us establish our conventions for indirect detection. We follow Ref. [88], and refer there for additional details. The DM differential energy flux into an observable state S for decay (dec.) annihilation (ann.), is given by 8 In both cases, dΦ/dE gives the number of S states, per detector area, per observation time, per energy interval, carrying units [particles/cm 2 /s/TeV]. The D and J-factors dictate the variation in flux across the celestial sphere for decay and annihilation, controlled by the DM density or density squared, respectively. The focus of the present work is the result of fundamental DM-SM interaction, and this enters into the particle physics (pp) factor. For a single annihilation or decay channel (for multiple channels simply sum over each weighted by the appropriate branching fractions), this can be written as These expressions expose the prompt spectrum per decay or annihilation, which represented our starting point in Eq. (2.1). We note that this expression neglects propagation effects such as oscillation, redshifting, and other interactions which will transform the spectrum from the point of production to detection. For details of how these effects can be incorporated, see e.g. [76,78,79]. Although Eq. (A.2) is how the particle physics factor is usually represented, the spectra should really be thought of as emerging from the differential decay width or annihilation cross-section, so that instead we have The spectra can then be defined as from which Eq. (A.2) follows, with the identification of the lifetime and cross sections appearing there as taking the tree-level or Born values, indicated by the 0 subscript. As we will consider spectra over a wide range of masses, it is convenient to move to dimensionless variables. Specializing to the case of decay from now on, we take x = 2E/m χ , and then the above becomes exactly Eq. (2.1). From now we will discuss the case of decay exclusively, although results for annihilation follow by rescaling m χ → 2m χ .
To calculate the decay spectrum we need to determine dΓ/dx. The decay will be initiated by a hard process occurring at a scale µ ∼ χ, where the DM decays to a set of SM particles each denoted by I. We can then approximate each of these SM particles evolving separately down in virtuality to a set of stable observable SM particles, labelled by S, which is exactly described by a FF as described in the main body. Corrections to this picture are described in the next subsection. For the moment, however, the spectrum is given by the convolution of the hard process with each FF, In the particular case where the hard interaction is a simple two body decay χ → XX, then so that the convolution is trivial, and exactly as in Eq. (2.2). If the initial hard process is more complicated than a two body decay, then it can be incorporated with a simple modification. In detail, the appropriate generalisation of Eq. (A.6) should be substituted into Eq. (A.5), and combined with the FFs we provide [4]. We emphasize that the individual FFs are made available, and can be used directly. As an additional example of how the FFs could be utilized, one signature of the formation of DM bound state formation could be the emission of a photon carrying away the binding energy, see e.g. [90][91][92]. For heavy DM, this photon can be sufficiently energetic that higher order processes become relevant, and a the photon should be replaced by the appropriate FF, thereby modifying the observable signature. Further details are provided in App. G.

A.1 Comparison to Existing Approaches
At present, experimental collaborations are primarily using either Pythia or PPPC4DMID to determine the DM decay spectra at masses above the electroweak scale. Nonetheless, a number of approaches to calculating these results have been presented in the literature.
Having just expanded upon the basic underpinnings of our approach, here we summarise several alternatives proposed in the literature, highlighting where we differ. Our first comparison to existing results appeared already in Fig. 1. As emphasized there, the dramatic difference observed between Pythia and our work is driven by the absence of the full electroweak interactions in the former; in particular, Pythia includes W and Z emission off fermions [94], but not the electroweak triple gauge couplings W W Z and W W γ. These terms are central in the development of electroweak showers at higher energies, and the spectra of particles they produce. Even at lower masses these differences can be important. This is shown on the top left of Fig. 4: electroweak effects have a visible impact for the hard photon spectrum. In particular, we see that Pythia does not predict a hard photon contribution, which arises from polarization effects it does not include. There we also compare our results to PPPC4DMID, 9 which augmented an earlier version of Pythia with leading order electroweak 9 The version of PPPC4DMID we used was downloaded in August 2019.   Fig. 1 for the photon spectrum, but for a DM mass of m χ = 2 TeV. Even at this lower mass, the additional electroweak effects in our calculation are producing clear differences with Pythia. We also demonstrate that our spectra differ from those in PPPC4DMID, a point further discussed in the text. The disagreement is not generic, however, and instead restricted to scenarios where the initial decay involves neutrinos. We show a single example in the top right, showing significantly improved agreement, and in general for these other channels the agreement is within the uncertainty we have ascribed to our results. On the bottom we show photon spectra for GUT scale DM (m χ = 10 16 GeV) decaying to light quarks, q = (u + d + s)/3. Here our spectrum is compared to results obtained using pure QCD DGLAP of a low energy fragmentation function, taken from Fig. 1 in Ref. [93]. Even for this hadronic channel, there is a clear difference in the hard photon contribution associated with electroweak evolution effects.
corrections, and still clear differences are observed. There are several possible sources for this disagreement, although at present we cannot isolate its exact origin. Primarily, at TeV-scales, finite electroweak masses become important, and while these were not included in our results (see App. E), they were in PPPC4DMID. However, we note that the simple expectation would be that including electroweak masses, would suppress the development of electroweak showers and the emission they generate. This is not realized in Fig. 4, where instead the peak of the PPPC4DMID spectrum is roughly a factor of two larger. We cannot resolve this difference without augmenting our calculation with fixed order corrections, which is beyond the scope of the present work. For this reason, at TeV-scale masses we do not claim our result is the more accurate, and instead suggest that the difference between the two results can be used to estimate the theoretical uncertainty on the spectra. The situation will ultimately be resolved by the addition of fixed-order corrections to our calculation. In spite of the pronounced differences visible in this figure, for non-neutrino channels we find agreement, within the theoretical uncertainties on our spectra, with Pythia and PPPC4DMID at mass scales below 10 TeV, a single example of which is shown on the top right of Fig. 4. In the bottom of Fig. 4, we contrast our results to an approach similar in spirit to that taken in the present work. The idea is to take a FF measured experimentally at low energies (rather than relying on Pythia), and DGLAP evolve this to the desired DM mass. This approach was put forward in Refs. [83][84][85]95], where the DGLAP evolution used either pure QCD or its supersymmetric analogue. 10 These methods have been taken up more recently by the authors of Refs. [43,93,[96][97][98][99][100] in order to set strong constraints on heavy DM. For hadronic channels, those works disregard electroweak interactions in the evolution, whereas our formalism includes these effects for all initial states. A comparison for light quarks is shown in Fig. 4 to a spectrum taken from Ref. [93]. The results are qualitatively similar, although with clear differences. In particular, the pure QCD evolution misses the hardest photons that result from the electroweak shower, and further there are differences at lowerx likely associated with our inclusion of soft-coherence. A related approach was taken in Ref. [101], where the authors used a hybrid approach of evolving the hadron FFs with QCD DGLAP, and decaying the particles with Pythia. As that work did not include electroweak effects, the differences to our results are similar to those shown on the bottom of Fig. 4.

B Details of the High-Scale Evolution
Here we expand upon the high-scale evolution from µ ∼ m χ to q W , describing both the full DGLAP calculation, and our treatment of soft coherence.

B.1 Review of DGLAP Evolution in the Unbroken Standard Model
For the DM masses considered in the present work, the starting point for our evolution is a scale far above electroweak symmetry breaking, µ q W , where the SM can be accurately described by an unbroken SU(3) × SU(2) L × U(1) Y gauge theory. One can therefore treat the electroweak gauge bosons, as well as all fermions as massless degrees of freedom. DGLAP evolution in this theory can therefore be used to evolve fragmentation functions gives the distribution at the scale µ of the momentum fraction x for particle species k in a shower initiated by a parton i (labeled by both type and helicity) produced in a hard process at momentum scale Q, and d k i (x; Q, µ) denotes the corresponding momentum weighted fragmentation function. The DGLAP equations take the standard form where the two terms in brackets correspond to the virtual and real contributions. This result can be viewed as the more complete version of the schematic form presented in Eq. (3.1). The details of how to solve these evolution equations was presented in Ref. [69]. We repeat only the salient points and refer the reader to the original work for details. The sum over I in Eq. (B.2) runs over the different interactions in the SM, and we denote by I = 1, 2, 3 the pure U(1) Y , SU(2) L and SU(3) gauge interactions, and by I = Y the Yukawa interactions. Besides these contributions to the evolution, there is also a mixed interaction, denoted by I = M . This originates from an interference contributions, where the particle i originates from a U(1) Y gauge boson B in the amplitude and the W 3 of the SU(2) L gauge group in the complex conjugate amplitude (or vice versa). The coupling of the mixed interaction is therefore proportional to The maximum cutoff on z in the integration of the real radiation is dependent upon both the splitting and interaction type. We define where V is the set of vector bosons. This prescription ensures that an infrared cutoff m V , of the order of the electroweak scale, is applied when a B or W boson is emitted. To evolve in the full (unbroken) SM, one needs to differentiate the two chiralities of the fermions, the two transverse polarizations of the gauge bosons, the mixed B/W 3 state, and include all 4 components of the complex Higgs field (instead of the longitudinal polarizations of the heavy gauge bosons). The complete set of states required is summarized in Table 1 Here i labels the states required at the high scale Q ∼ m χ . As we work in the full unbroken standard model, i can take 58 values: 42 fermions (2 × 3 × {e L , e R , ν L , u L , u R , d L , d R } for particle/antiparticle and three generations), 2 gluons (helicities), 4 charged electroweak bosons (two helicities and two charges), 6 neutral electroweak bosons (two helicities for each of B, W 3 , and the mixture, collectively labeled X 0 ), and the four degrees of freedom of the SU(2) L Higgs doublet. For k at the weak scale, the counting is similar although slightly rearranged. We now have 38 states as we only distinguish the helicity of the electroweak states (as this information is used in the weak matching), so we now have 26 fermions sensitivity as well. Double logarithmic contributions arise from the limit where radiated particle are simultaneously soft and collinear relative to the particle they were emitted from. In the strong interaction, these simultaneously soft and collinear contributions cancel between the virtual and real terms in the DGLAP equations. This occurs as an arbitrarily soft emission of a gluon cannot be observed experimentally, so the divergence associated with this emission must cancel against the virtual contribution. This is different from the case of the soft emission of a W boson, which can always be observed through the change of flavor (or SU(2) quantum numbers) of the emitting particle. Thus, as long as a process is sensitive to the SU(2) quantum numbers of the external states, soft radiation of W bosons from these particles leads to an incomplete cancellation of the soft and collinear divergences, which gives rise to double logarithms. The form of the soft boson cutoff in Eq. (B.4) ensures that these double logarithms have the correct coefficients [67,68].
As discussed in Ref. [69], the set of evolution equations can be decoupled to some degree by switching to a basis of well-defined isospin T and CP. The definition of all FFs in this new basis was given in Ref. [69]. To provide examples, for left-handed fermions one can write in a basis d TCP while for the SU(2) bosons we have Using this isospin basis allows to isolate the double logarithmic dependence. In particular, for isosinglets (with T = 0) there cannot be any double logarithms generated, since the emission of an isosinglet cannot change the isospin of the emitting particle. In general, the double logarithmic term is given by a Sudakov factor, which depends on the total isospin, and takes the form One can then show that the rescaled FF has only standard single logarithmic evolution. As mentioned in the main body, rather than perform the DGLAP evolution from µ ∼ Q = m χ /2 down to µ ∼ q W , we instead start at the electroweak scale and evolve upwards. 11 In detail we start with q W = 100 GeV as the starting point for our evolution. For quarks and leptons (k = f ), assuming that the helicity of the fragmentation product is not detected, we take as input setting all other initial FFs to zero. The only exception to this is the top quark, where in order to correctly account for its decay at the electroweak scale, we evolve t L and t R separately. Similarly, as neutrinos (k = ν) have no right-handed states, the non-zero initial conditions becomes For fragmentation into a gauge boson V we keep track of helicity (as we did for the top), because for the electroweak scale vectors we will use this information in the matching. Accordingly, the non-zero initial conditions are Ultimately, the DGLAP evolution will provide expressions for Several example outputs at this stage are shown in Fig. 5 (including soft coherence), where we have chosen to highlight the generation of polarization effects through the evolution. In particular, due to the chiral nature of the SM, the chiral fermion states and boson polarizations do not evolve identically. See Ref. [70] for an extended discussion of this effect. Because we are neglecting any SU(2) breaking effects in this evolution, one expects a break down in this description for Q ∼ q W . We discuss this in more detail in App. E. To continue the evolution below q W we will first need to remove particles with masses m ∼ q W , and we will describe how to do so in App. C. Before doing so, however, we next outline how to incorporate a partial treatment of soft coherence effects into the DGLAP evolution.

B.2 Incorporating the Soft-Coherence of Real Radiation
At low energy fractions the simple DGLAP evolution described above is not reliable: there are large double logarithms of x that need to be resummed. The physical origin of these logarithms is the soft-coherence or angular-ordering effect described in the main body. In this section we will demonstrate that taking the output of our DGLAP evolution and applying the substitution m χ → x m χ incorporates a partial treatment of these effects. This treatment is manifestly incomplete: in particular no accounting for the virtual effects of soft-coherence will be included. One manifestation of this shortcoming will be incomplete momentum sums, which we discuss in App. F. A full treatment of these effects is left as an open problem.
To begin with, consider the soft coherence of gluons in the context of QCD, as reviewed in Ref. [66]. The resummation of the gluon-to-gluon FF at small x is given in leading-logarithmic approximation (LLA) by Eq. (5.21) of [66] as 12 This result provides a FF with soft-coherence included, at least as far as it impacts real radiation. As we will now show, this FF can be determined as the solution to the DGLAP equation, but in q = xQ rather than Q.
In order to expose this, first we introduce xD g (x; Q) = xD g (x; q) ≡ d g (x; q) -we will find thatD g satisfies the unmodified DGLAP equations. Introducing this notation to rewrite Eq. (B.13), and the substitution Q = q/x, after differentiating with respect to ln q 2 , we find (B.14) We can rewrite this using the identity which holds for n > 0. Accordingly, (B. 16) In the second step we used Eq. (B.13), and in the final step we used the fact that in the small-x limit P gg (x) = 2C A /x. Accordingly, d g (x; q) satisfies the DGLAP equation to LLA. Therefore to obtain the correct small-x gluon FF we should solve the DGLAP equations for evolution in q and then set q = xQ where Q is the hard process scale. At large x, and for those terms in the evolution equations that do not give rise to extra small-x logarithms (i.e. those without a 1/z singularity in the splitting function), this gives rise to unenhanced NLO contributions, which are in any event beyond the precision of our treatment. Therefore the same procedure can be applied to include the small-x suppression of all the QCD FFs. Note for x < q W /Q the substitution samples D g (x; Q) for Q < q W . As we begin our evolution at q W , these results simply vanish, and thus at the high scale there is an artificial cut in the distribution at small-x. This effect is washed out after convolution with weak matching and Pythia, but is unphysical and a manifestation of our incomplete treatment of soft coherence. The same large double logarithms of x arise whenever a sequence of emissions with 1/z singularities in the splitting functions can occur, which is also the case for W boson fragmentation in the unbroken SM. We expect a similar small-x behavior of the SU(2) evolution equations, with α 3 replaced by the appropriate gauge coupling α 2 . There is, however, an additional complication in this case, namely the double-logarithmic evolution of FFs with non-zero weak isospin. As discussed above, these FFs are suppressed by the Sudakov factor in Eq. (B.7), due to a mismatch of real and virtual contributions, which is not relevant to small x. As such, the substitution q = xQ does not apply to the isospin suppression factor, and the formula for small-x resummation becomes We apply this prescription to the FFs above the EW breaking scale, D h h (x; q > q + W ), before the matching corrections discussed below. The boundary conditions for D andD at q + W are identical, as they are proportional to δ(1 − x). We do not apply the prescription to the FFs at lower scales q < q W provided by Pythia, since these already take coherent emission into account and are tuned to experimental data at such scales.
The above expands upon the m χ → x m χ substitution mentioned in the main text. In practice, this effect generically suppresses small-x contributions to FFs, as the states that the suppressed bosons would have split into are also removed. An example is given in Fig. 6 for the example of a photon evolving to a tau neutrino.

C Matching at the Electroweak Scale
Evolution through the electroweak scale is handled analytically. As already emphasized, the chiral nature of the SM ensures that quite generically the spectra of states resulting from the high scale evolution will be significantly polarized. Through a matching procedure outlined in this section, we ensure the polarization information of the electroweak states is not discarded. In detail, we will compute FFs where q ± W = q W (1 ± ), i.e. we evolve the states through a threshold of differential width at the electroweak scale. We choose 1 in order to ensure no large logarithms associated with the evolution can be generated at this step.
For states that do not have electroweak masses, the threshold is uneventful. For fermions, we combine chiralities at this stage. For example, we only track the evolution of e below q W , whereas e L and e R are evolved separately at the high scale. For a fermion f we take 13 Similarly, the helicities of the photon and gluon are combined into a single unpolarized state. Note this procedure washes out any remaining information about the polarization: we explicitly assume that the experiments are searching for unpolarized states. Even then, this assumption is associated with an imprecision in that muon and tau decays do depend on helicity. Extending our formalism to evolve the polarization down to the scale of leptonic decays is left to future work. For the Higgs, Z, W , and top, q W marks the end of their evolution. In the following subsections we outline how their momentum is redistributed amongst their decay products. For each, the strategy is as follows. We begin by calculating the differential energy and, where relevant, angular spectrum of the decay products for all polarizations or spins in the rest frame of the electroweak state. At this stage we account for all electroweak scale masses, although particles that will continue their evolution below q W are left massless. Then, in order to match these results onto the high scale distributions where all particles were treated as massless, we perform an infinite boost, i.e. boosting the electroweak states to an energy E q W . At this stage angular differences in the rest frame are transformed into energy differences in the boosted frame, indicating why this information was retained. The distribution of the energy fractions amongst the boosted decay products then gives us exactly D(x). We now implement this procedure case by case.

C.1 Higgs Decays
As a scalar, the Higgs carries no polarization information: there are no initial polarizations or spins to account for. For this reason, we can extract the D(x) simply by generating the spectrum of boosted Higgs decay products in Pythia. To do so, we generate e + e − → HH events at √ s = 200 TeV. 14 We forbid initial and final state showering of any kind, and turn off hadronization. The only states we allow to decay are the W and Z. The energy distribution of leptons, neutrino, quarks, gluons, and photons associated with each Higgs is collected, and used to form

C.2 Z Decays
We begin by computing the differential decay spectrum of f in Z → ff , working in the Z rest frame. Here, and throughout, we work in unitary gauge. In order to establish our conventions, the coupling between f and Z is determined by Here c L = I 3 W − Qs 2 W and c R = −Qs 2 W , with I 3 W weak isospin and Q electric charge, whereas c W and s W are the cosine and sine of the Weinberg angle, respectively. In the broken phase of the SM, for each generation, ν and u carry I 3 W = 1/2, whilst e and d have I 3 W = −1/2. The charges are, of course, Q = 0, −1, 2/3, −1/3 for ν, e, u, and d. For the moment, we will perform the calculation for the case of a single f with arbitrary couplings, and then we can ultimately weight this result by the appropriate branching fractions for the specific states in the SM. Continuing with our conventions, we define our coordinates such that the ff are produced in the x-z plane, with f produced at an angle θ from the z-axis. When we eventually boost the Z, we will do so in the z direction, and accordingly in the rest frame we choose our polarization vectors as follows, In terms of these quantities, we can now compute Γ(Z ±0 → ff ), and thereby determine the following distribution of angles, (1 − cos 2 θ) .

(C.6)
Observe that each expression is a normalized probability distribution for cos θ, and thus we will label each of these expressions as p Z ±0 (cos θ), denoting the distribution the f emission angle is drawn from in the Z rest frame, for a given polarization. Note, these distributions also contain the information on the emission angles off . Definingθ to be the anglef makes with the z-axis, we havep Z ±0 (cosθ) = p Z ±0 (− cos θ).
14 At this energy, the minimum energy fraction for two-body decay product is x ∼ m 2 H /4E 2 H 10 −6 , and therefore below the smallest values considered in this work.
Armed with these results, we will now determine the energy distribution of f andf in the boosted Z frame. We wish to boost the Z to energy E Z , propagating along the +z direction, which we can achieve through the boost After this boost, the energy of f , which was m Z /2 in the rest frame, is now where in the last step we implemented the large boost approximation. In this limit, the energy fraction carried by f is given by , with cos θ drawn from the appropriate distribution given in Eq. (C.6). We then determine the distribution of x f by the following change of variables, where Θ is the Heaviside step-function. Similarly, In terms of these expressions, we can calculate the expected energy fractions for an arbitrary p(cos θ) as follows, and an identical calculation yields xf = (1 − cos θ )/2. Accordingly, independent of the exact form of p(cos θ), we have x f + xf = 1, consistent with momentum conservation. This also indicates that the momentum weighted fragmentation functions should be associated with x p(x), and so for the case of a single fermion we have where again, the explicit p Z ±0 (cos θ) are given in Eq. (C.6). In detail, (C.13) The equivalent results forf follow by taking c L ↔ c R in each case. Observe that for each polarization, we have explicit momentum conservation independent of the values of c L and c R , We can extend the result immediately to all the states in the SM by weighting these distributions by the appropriate Z branching fractions, and in all cases inserting the appropriate values of c L and c R .

C.3 W Decays
Having determined in detail the spectrum for the Z decay products, the result for W decays follows almost immediately. As the W couples only to left-handed fermions, we simply take Eq. (C.13) (and the analogous result forf ) and set c L = 1 and c R = 0. Again weighting the different final states by the appropriate branching fractions, this specifies all relevant fragmentation functions.

C.4 Top Decays
As the top decay is a three-body process, the distribution of energy fractions amongst the decay products is a more involved calculation than for the states already considered. We take |V tb | = 1, considering only t → bW , and to be explicit let us take the decay W → ff , so that the full process is The translation of the final results to all relevant W -decay modes is straightforward, and each channel can then be weighted by the appropriate W branching fraction. Once again, we will begin the calculation in the top rest-frame. We measure the spin of the top along the z-axis, denoting spin up and down by t ± . Ultimately we will boost the top in the z direction, so that in the infinite boost limit, spin up and down will be associated with positive and negative helicity, or right and left handed chirality. In that sense, while the sign of t ± will represent spin in the first part of the calculation, it will later be translated directly to helicity.
Consider first the spectrum of b-quarks produced from t ± → W b (ignoring any contribution from W decays at this stage). The intermediate W is produced on-shell, and thus the calculation proceeds similarly to that of the Z. A key difference, however, is that here and throughout we will retain both the top and W mass, defining the ratio ρ = m 2 W /m 2 t ∼ 0.2. All other fermions, including the b-quark, are left massless. Further, we need to account for the three possible polarizations of the W . Due to the chiral nature of the weak interaction, the resulting b-quark must be left handed. This implies in the limit where m b = 0, only a longitudinal and negative helicity W can participate in the process. The ratio of the branching fractions is so that the longitudinal polarized W s are produced twice as often. We can use this to gain a rough intuition for the resulting spectrum of b-quarks. For a longitudinal W , the left-handed b will be preferentially emitted in the opposite direction to the spin of the top. After boosting, we then expect a softer spectrum for t + than t − . Performing the calculation, we find (C.16) The step-function enforces the condition x b ∈ [0, 1 − ρ] ∼ [0, 0.8]: the finite m W ensures the b cannot carry away all the energy. The distribution of energy fractions is presented in Fig. 7. We now turn to the distribution of W decay products. As the W is produced on shell, we can start in the W rest frame, where the energy spectrum is a two-body δ-function. Yet, there will be a non-trivial angular dependence associated with the W polarization, which will be converted into energy differences when we boost to the top rest frame, and again when we boost the top itself, at which point we need to account for the angular distribution of the W as a function of boson polarization and fermion spin. Further, we must account for the fact that as the W is an intermediate state, we will have interference between the different polarizations.
We begin by establishing our coordinates. In the W rest frame, we represent p f with a general lightlike four-vector The equivalent expression for pf can be obtained by sending θ → π − θ and φ → π + φ. In the top rest frame, we choose coordinates such that the b and W are emitted in the x-z, with the W at an angle θ from the z-axis, i.e.
where E W = m t (1 + ρ)/2 and p W = m t (1 − ρ)/2. In terms of this, to transform p f to the top rest frame, we boost in the W direction by an amount after which Eq. (C.17) becomes (C.20) In terms of these coordinates, the fully differential spectrum in the top rest frame is For both spins this is a normalized probability distribution for the three relevant angles, which we denote p t ± (φ, cos θ, cos θ ). We can now use this to determine the energy fraction carried by f in the boosted top frame. To determine this, we boost Eq. (C.20) in the −z direction by an amount γ = E t /m t and β ≈ 1, to obtain Note that x f ∈ [0, 1]. At this stage, we know the energy fraction as a function of the angles, and the distribution from which the angles are drawn, so we can formally write down the relevant FF,  Table 2. Parameters of the fitting function Eq. (C.24) used to provide an adequate description of the W -decay products resulting from polarized top decays. The spectra themselves are depicted in Fig. 7.
The equivalent expression forf follows by sending θ → π − θ and φ → π + φ in the argument of the δ-function. These expressions can be readily computed numerically, and are depicted in Fig. 7. To facilitate rapid evaluation, we further determined piecewise polynomial fitting functions. Recall that x b < 1 − ρ. Momentum conservation then requires x f + xf > ρ, and thus both distributions display a discontinuous derivative at x = ρ. As such, we determine fitting functions of the form For each spectrum, p is fixed by the asymptotics as x → 1, whereas N and N were chosen such that {a n , b n }, determined by a least squares fit, provided a satisfactory description. Explicit values are provided in Table 2.
We do not need to repeat any calculations in order to obtain the equivalent spectra for anti-top decays, instead we obtain the result by a CP transformation. In detail, CP flips the helicity off all states, and also interchanges particles and antiparticles. Recalling that helicity and chirality are identified for a massless particle, but opposite for anti-particles, we have Consequently, the spectrum ofb is given directly by Eq. (C.16). The distribution forf is given now by Eq. (C. 23), and that for f can be obtained by taking the same equation, but with θ → π − θ and φ → π + φ in the δ-function argument.

D Low Scale Evolution with Pythia
The evolution of our FFs below q W is computed with Pythia v8.235 [1][2][3]. In this section we expand upon this, outlining the options used in running Pythia, and also the modifications we made for our purposes such as an improved treatment of FSR, and how we incorporate the proton mass. Firstly, let us outline the basic details of how we used the program. In order to calculate d S X (x; q W , 0) appearing in Eq. (2.3), we simulate events with a hard interaction e + e − → XX at √ s = 2q W , with initial state radiation switched off so that the e + e − operates as an energy injection, starting the X at µ ∼ q W . We then determine d S X (x; q W , 0) from the spectrum of stable S particles in the hemisphere of the initial X. This is all handled using a modified version of the example script main07. All long lived particles that may not necessarily decay on collider scales are forced to decay, in particular muons, pions, kaons, and neutrons. We leave off all electroweak radiation effects, as these were in the higher steps of our evolution. This means that neutrinos do not evolve at all at this lower stage, although their spectra can still receive contributions from various decays. Finally, a number of details of photon emission are modified; the choices and motivations here are discussed next.

D.1 Improved Treatment of FSR for DM
Below the weak scale, charged fermions will contribute to the photon spectrum via FSR, in the form χ → ff γ. The expression for this contribution is known analytically, as given in Eq. (5.1). We can arrive at this result by, for example, calculating the photon spectrum from a three-body decay of a vector, V , via the process V → ff γ. Giving the vector a mass m V = 2q W , and defining = m 2 f /m 2 V , we can extract the leading order result from the analogous QCD calculation in Ref. [102], where e f is the charge of the fermion, and α = α EM . In the limit q W m f , we have 1 and we can expand the above result to arrive at twice the result in Eq. (5.1) (the factor of two arises as this is the spectrum of f +f ). 15 The first correction to Eq. (5.1) is O( ), and so this is an excellent approximation for all remaining SM fermions below q W .
With the analytic expression in hand, we can compare this to the output of Pythia for an example state, f = e. The result of doing so is shown in Fig. 8. At lower x values Pythia is under-predicting the number of photons. The origin of this mismatch is that Pythia imposes a p T cutoff on the QED shower: photons that would carry a p T < 1 keV by default are not produced in the shower (this is controlled by TimeShower:pTminChgL). In order to validate this interpretation, note that the analogue of Eq. (5.1) in the presence of a p T cut is, after expanding in (p cut T ) 2 /q 2 W 1. This expression is also plotted in Fig. 8, and we see the Pythia result trending towards it at low energy fractions. As discussed, given concerns about the accuracy of our results at small-x, we do not produce spectra below x = 10 −6 . Further, for most channels, due to showering into other states that can produce photons through hadronic channels, FSR is often a subdominant contribution to the photon spectrum at low-x, particularly as we move above q W and the shower can develop. Nevertheless, at low masses certain states such as right handed electrons still receive a large contribution from FSR and thus we outline a procedure we used to correct this issue. We note this treatment may also be useful for studies of DM at lower masses also.
Our approach is to turn off FSR in Pythia and instead deal with it analytically. Showering for all fermions is switched off, although we leave on γ → ff , as this is not subject to the issue described above. The full photon spectrum for a given initial state can then be determined using Eq. (5.1) and, so that the low scale FSR is handled analytically instead of with Pythia. Adding Eq. (D.3) alone is inconsistent with momentum conservation. The charged particles that would have emitted photons also need to have their distributions corrected to ensure the momentum is correctly subtracted. In particular, there will be both real and virtual corrections to the charged fermions' FFs. For real emission, if we have a process where the high scale evolution produced a charged fermion, f , from an initial particle X, then we need to account for the fact the f would have lost momentum to radiated photons via f → f γ in the evolution below q W , which we now exclude. In particular, if the fermion has an initial momentum fraction z and the photon carries away a fraction 1 − u, then the fermion FF receives a correction, lost by the sum of all charged fermions. Looking at one in particular, we have Summing over f , this is exactly the negative of Eq. (D.7) as claimed, demonstrating using both terms restores momentum conservation. In summary, we implement the corrected FSR treatment with a combination of Eq. (D.3) and Eq. (D.6). The latter needs to be computed carefully, as [D γ f (1 − z; q W , 0)] FSR diverges as z → 1, although that divergence is regulated by the vanishing of the two terms in square brackets in the same limit.

D.2 Incorporating the Proton Mass
The finite proton mass provides a physical cutoff in the energy fraction the p/p final states can carry of x ≥ 2m p /m χ . In our calculation at the high scale, however, all states are treated as massless, allowing for in principle arbitrarily small energy fractions. At the low-scale, however, the evolution in Pythia explicitly includes m p = 0, and therefore will have a sharp cut-off at x = m p /q W ∼ 10 −2 . In summary, the finite mass of the stable hadrons we are interested in is not treated consistently throughout the calculation. By default, this inconsistency will manifest in two ways. Firstly, as at the high scale we allow arbitrarily small energy -or, as our particles are massless, equivalently momentum -fractions, after convolving this with Pythia, by default we will have a non-zero spectrum for x < 2m p /m χ . Secondly, as the proton mass is only treated in one part of our calculation, a feature at x ∼ 10 −2 can generically be expected to appear. This occurs because, if the high scale spectrum has a sharp feature near x = 1, a residual bump at 10 −2 can remain after convolution with Pythia.
Both of these problems are on display in Fig. 9, focusing on the comparison between the unweighted spectrum and Pythia for the moment. On the right, we see a clear feature at 10 −2 . This results from the fact that for a relatively light DM mass of m χ = 2 TeV, the b-quarks will not always undergo significant evolution by q W , and there will still be a large contribution to the fragmentation functions near x = 1. When convolved with Pythia, the result is a clear bump. Fixing this particular problem in detail is left to future work, however we note the feature rapidly becomes less pronounced as we increase in mass, and for states that do not have a large colored component near x = 1, as is the case for neutrinos viewed on the left, the issue is far less apparent.
The second issue, of a non-zero spectrum for unphysical energy fractions, we will resolve. Our default output at this stage is the "unweighted" dotted distributions, exhibiting this exact behaviour. We will reweight these results as follows. To begin with, we treat all x values as momentum fractions (these are interchangeable at the high scale as discussed). So for a given x and DM mass m χ , a proton would have a momentum p = x p m χ /2, and hence an energy fraction Importantly, as x p → 0, x E ≥ 2m p /m χ , as required. This allows us to transform to the relevant energy fractions, and further the spectrum as The above change of variables from x p → x E will ensure that the spectrum is cut-off at the correct x value, however, it does not ensure the spectrum exhibits the correct asymptotics as x p → 2m p /m χ . In general we expect there will be power corrections to our result in this limit, which manifest here in the form of a threshold factor (x p /x E ) k , so that By comparing our results near the threshold with Pythia, we found good agreement for k = 3. The combined result of the rescalings is depicted as the solid orange curves in Fig. 9.

E Estimating the Accuracy of our Results
As discussed in the main text, there is considerable scope for systematic improvement of the results presented in this work. In this section we will present a more quantitative discussion of this point, outlining the formal accuracy of our calculation and estimating the size of neglected terms. As we will show, at high DM masses our spectra for x ∈ [10 −3 , 1] are accurate to O(10%). For x 10 −3 , the uncertainty increases (particularly due to an incomplete treatment of soft coherence) to the O(1) or even order-of-magnitude level. Nevertheless, in light of the pronounced differences with existing results as shown in Figs. 1 and 4, our results represent a significant improvement. Obtaining a full NLL calculation with theoretical uncertainty bands and O(5%) errors at the level that has been achieved for specific DM spectra, see e.g. [104], represents a clear target for future work. Our treatment of fragmentation function evolution above the electroweak scale q W ∼ 100 GeV is according to the leading-order DGLAP evolution equations, with neglect of all particle masses. The evolution code uses Heun's method for direct solution of the integro-differential equations, with typical precision of a few parts per mille. However, the evolution equations only treat terms enhanced by logarithms of m χ /q W and therefore results are only reliable well above the electroweak scale. While mass corrections are expected to be of order q 2 W /m 2 χ , other terms that are not logarithmically enhanced can be important for TeV-scale DM masses. These could be included in future by matching to fixed-order matrix elements, as discussed in Ref. [75]. Meanwhile, for DM decay or annihilation at or below the TeV scale, a direct parton shower simulation taking account of masses, such as Pythia, would be more reliable for the majority of final states. The exception would be final states where the spectra is dominated by electroweak showers, in particular neutrinos as highlighted in the left of Fig. 4.
For all except the SU(2) L electroweak interaction, terms resummed by leading-logarithmic (LL) DGLAP evolution from the electroweak scale to scale ∼ m χ are of order α n I L n where L = ln(m χ /q W ). However, our calculation fails to resum all αL terms related to soft emissions. Nevertheless, αL remains less than 1 to the Planck scale, as shown in Fig. 10. In detail α 3 L ∼ 0.75 and α 2 L ∼ 0.80, taking into account the running of both couplings.
In the case of SU(2) L , fragmentation functions in general have non-isosinglet contributions of order α n 2 L 2n . These, together with a class of terms down to order α n 2 L n+1 , sum up to yield leading-logarithmic Sudakov factors, of the form exp[−Lg 1 (α 2 L)] where g 1 is a known function. The missing NLL terms are of order α 2 L, just as in QCD. As shown in Fig. 10, these terms are smaller for masses up to m χ ∼ GUT scale, and at higher masses remain comparable to α 3 L and crucially less than unity. Failure to fully resum these terms induces an error of size O(10%) up to the EeV scale. A full NLL calculation would suppress these effects by an additional α I ∼ 0.1, and reduce these errors further to the 1% scale. Clearly the αL 2 terms are rapidly larger than one, and thus must be resummed for a reliable result for most of the mass range considered in this work, and our calculation achieves this aim. The αL terms, which we do not fully resum, remain less than unity across the entire mass range. In both cases we account for the running couplings.
All the above relates to evolution of fragmentation functions due to collinear enhancements in the relevant matrix elements, giving rise to large logarithms of the energy scale ratio m χ /q W . At small values of the energy fraction x there are also large (double) logarithms of x that need to be resummed. As discussed above, these have the effect of strongly suppressing fragmentation at small x, essentially due to destructive interference between different amplitudes involving soft gauge bosons. These soft coherence effects are well understood at the leading-logarithmic level in QCD, but less so at NLL in QCD and even at LL in the electroweak sector. Our treatment takes account of soft coherence at LL level in QCD, with a plausible extension of the same effects to the full SM. However, it is difficult to estimate the quantitative uncertainty of this procedure. A qualitative estimate of the uncertainties of this procedure is provided in Fig. 11. The figure depicts spectra arising from χ → bb for three different treatments of the small-x physics. Firstly, we show the result if you use pure DGLAP evolution with no treatment of soft-coherence. This is certainly an overestimate of the true soft multiplicity. On the other extreme, we show a result where we apply our m χ → xm χ correction to the final FF, i.e. the substitution is applied to the left-hand side of Eq. (2.3), rather than just to the high-scale or DGLAP FF, as we do by default. This results in most likely an overly aggressive suppression. Pythia already has a partial accounting for soft-coherence, so applying the substitution globally we are double counting the effect at this stage. Further, as mentioned in App. B.2 this substitution should not be applied to the electroweak double logs, however we cannot factor out their contribution to the final FF. Finally, we show our default procedure where the correction is only applied to the high-scale evolution. We see that our approach sits in between the two extreme alternatives. For x 10 −3 the differences are not that pronounced, but by x ∼ 10 −6 there is now more than an order of magnitude separation between the approaches. This can be taken as a rough estimate for the uncertainty we have at these scales, although Fig. 11 also demonstrates that by the Planck scale the differences are reduced to a factor of ∼ 2 (we emphasize that the difference between the three approaches shown in the two plots is driven entirely by the mass scale). The results also clarify that in order to reduce these uncertainties, a serious study of coherence effects in the whole SM is required.
Turning to the matching at the electroweak scale, the calculation was performed at leading order. Next-to-leading QCD corrections, of order α 3 ∼ 10%, could be included to improve precision, although this would require care to avoid double counting in the subsequent Pythia shower. Electroweak corrections would contribute at the few percent level. Finally, below the electroweak scale, the Pythia parton shower and hadronization generator generally agrees with collider data at such energies at the 10% level. This would be difficult to improve significantly without major advances in event generator technology In summary, our predictions at high DM masses and moderate to high x values are subject to uncertainties at the few times 10% level, which could be reduced somewhat by inclusion of higher-order corrections at the evolution and matching stages. Uncertainties increase markedly at lower energy fractions, due to a lack of precise understanding of soft coherence effects in the full SM.

F Additional Results
Having outlined the details of our calculation and discussed their accuracy, we now present a number of additional outputs from our formalism. Firstly, we present a number of additional spectra in the same spirit as Fig. 3, highlighting additional physics inherent in our results. Afterwards we present a non-trivial cross-check on our results, demonstrating the extent to which momentum is conserved as it is repartitioned among the various states through our evolution. In a similar vein, we will then show how the momentum is distributed between the various states for selected processes, demonstrating how this varies as the mass is increased. F.1 Additional Spectra Figure 12 and Fig. 13 furnish additional examples of the output of the formalism introduced in this paper. In all cases, we show the spectrum of stable SM final states for DM masses between 1 TeV and m P ∼ 10 19 GeV, with a spectrum shown at each decade in mass. Quite generically, we see that when plotted in dimensionless variables, the spectra vary most rapidly for scales around the TeV scale, often slowing as the mass approaches the Planck scale. Just above the electroweak scale, new channels become kinematically available through the emission of states with m ∼ q W , and as the evolution heads into the full unbroken SM. This is clear in the case of the muon neutrino spectrum from χ → e + e − shown in Fig. 12, which can only arise from electroweak boson or hadronic decays, both of which are primarily accessed through electroweak states. This can be contrasted with the positron spectrum resulting from χ → HH. The rich decay pattern of the SM Higgs already involves many SM states, explaining the lack of significant evolution in the spectrum. Nevertheless, the hardest emissions near x = 1 do evolve considerably, growing rapidly as multiple electroweak emissions become available, and then softening again as the size of electroweak showers develops. In all cases, the evolution eventually slows down as the various channels become well mixed. This is perhaps unsurprising, given the connections between the DGLAP equation and diffusion. We reiterate that our calculation assumes that there are no new-physics thresholds crossed between q W and m P . If in fact there were, rapid variations could again be observed as momenta is redistributed amongst the newly available channels. For results in this direction, see e.g. [83][84][85][86][87].
As emphasized a number of times already, the chiral nature of the SM plays a central role in the high-scale evolution. This point is further emphasized in Fig. 13, where we depict the electron neutrino spectrum resulting from the three polarizations of the massive Z-boson. In general the softest emissions, dominated by QCD hadron decays, is comparable between all three states. The hard emissions, however, differ dramatically, particularly when considering the spectrum of a purely chiral state.

F.2 Confirming Momentum Sums
Conservation of momentum implies that the momentum weighted FFs must satisfy the following consistency condition, , this equation is satisfied trivially. However, it must also remain true as the momentum is repartitioned amongst different states through the evolution in virtuality, and the result becomes an important check on the evolution. In this section we discuss how well our results satisfy Eq. (F.1), taking various values of Q, and µ 0 ∼ 0 appropriate for the end of our evolution.
The results are provided in Table 3. The table shows a select set of values for a and Q in Eq. (F.1). In particular, it is clear that near the electroweak thresholds, the deviation from perfect momentum conservation is O(5−20%), whereas by the highest scales considered in this work, they have shrunk to O(2−4%). These results represent an irrefutable uncertainty in our results. Nevertheless, although there are a many steps in our calculation where momentum is redistributed, each of which could contribute to the errors shown in the table, the uncertainty is almost exclusively due to a single source: our procedure for incorporating soft coherence described in App. B.2. Indeed, when we do not implement the corrections described in that section, we find momentum conservation is obeyed to better than 1% in all cases, limited by the numerical precision used in our calculation.
Our procedure for implementing soft-coherence removes the real-radiation associated with destructive color-interference. However, we do not account for the associated virtual corrections: the suppression of soft emission also increases the probability for a state to not emit  and thereby retain a larger fraction of its momentum. Accordingly, the presence of an offset in Table 3 is unsurprising, although also representative of a clear target for improving the treatment of color-coherence in our results.

F.3 Momentum Distributions Amongst Final States
In addition to checking the overall conservation of momentum, we can also consider how that momentum is redistributed amongst the stable SM final states. In Fig. 14 we show exactly that, plotting the momentum fractions for all final states, for two example spectra, χ → ν eνe and χ → bb. In both cases, particles represent the momentum carried by the state and its conjugate (where applicable), so p represents the momentum carried by p +p. For the neutrino initial state, we see considerable variation as a function of mass. The difference between the e and ν e FFs, being an isovector quantity, evolves double logarithmically according to Eq. (B.7), so that their fractions become almost equal at high enough masses. At high enough masses, the fractions are identical for e and ν e , as a result of their connection in the unbroken SM. For the representative hadronic channel, the distributions are highly stable as a function of mass. Of course, even if the total momentum fraction deposited into a given SM state is constant, how that momentum is divided between individual states can still evolve considerably, as demonstrated in Fig. 3. The tight correlation between the photon and muon neutrino can be understood as follows. Strong isospin implies charged and  Figure 14. Momentum fraction carried by all stable SM final states considered in this work for two example decays, χ → ν eνe (left) and χ → bb (right). In both cases we show how the fractions evolve as a function of DM mass. In the legend, particle labels are a proxy for the contribution of both particles and antiparticles, so ν µ labels the momentum fraction carried by both ν µ andν µ .
neutral pions will be produced in these decays at a ratio of two to one. The decays of π 0 will produce two hard photons, each carrying a large momentum fraction, whereas the π ± decays will produce only a single hard muon neutrino (and a softer one which will not generally carry a large momentum fraction).

G Details of the Public Code
The spectra generated in this work are publicly available at github.com/nickrodd/HDMSpectra. Examples of how to generate spectra for arbitrary initial states or even individual FFs is provided. For several cases, the code can also be used to extract the coefficient of δ(1 − x) in the spectrum, an example output is shown in Fig. 15. Further, the repository contains the details of how to reproduce many of the figures in this work. In this section we outline several additional details of how those results were computed, but for details of how to use them we refer to the repository. We emphasize once more that all spectra provided in the repository are prompt: no propagation effects are included.
With the exception of the weak matching, all details of our calculation are computed numerically. For the high scale evolution, we solve the DGLAP equations using the procedure outlined in [69,70]. We note that in this stage of the calculation we made use of LHAPDF [105]. At the low scale, we evolve our results using Pythia. In each case, we determine the FFs d(x) as ln(x) spaced histograms, and then perform the convolution using the approach in App. G.1.
At the end of the procedure, we have a collection of 616 FFs evaluated at a set of Q values between 500 and 10 19 GeV. We then implement a reduction algorithm to reduce this to a minimal set of points necessary for retaining the details of the spectra at the level of accuracy of our calculation for 500 GeV < Q < 10 19 GeV and 10 −6 < x < 1. In this reduction of points we ensure that all spectra are unchanged to within 1% in the region 10 −4 < x < 0.99, while we allow for larger deviations in regions where the precision of our calculation is expected to be worse. To be precise, the accuracy as a function of x we use is We then discard as many points as possible, while maintaining this accuracy, in order to compress the output dataset. This data is then packaged into a single file which is accessed via the public code.

G.1 Computing the Convolution of Binned Fragmentation Function
A central tool in the present work was the use of the convolution expression satisfied by FFs, This result was exploited to simplify our calculation into steps as shown in Eq. (2.3), and appears frequently in our calculation of the DM spectra. A key ingredient in the final result is the spectra obtained from Pythia, which are inherently binned. It is convenient to have a form of Eq. (G.2) appropriate for binned FFs. Such a result is presented in this section. Before doing so, let us briefly provide some intuition for Eq. (G.2). The result quantifies that the probability a particle a at a scale µ 1 produces a particle c at µ 3 carrying momentum fraction x, is given by the combination of the probability of a → b at an intermediate scale µ 2 , and then b → c, summing over the allowed states and momenta for b. One may worry about only keeping track of a single state across the threshold µ 2 , rather than all details of the evolution. For example, if b is quark or gluon, it will be color connected to other objects, and this information is lost across the threshold. Here we invoke Amati-Veneziano preconfinement [106]: as long as the various scales are sufficiently separated, what we color connect b to is asymptotically irrelevant. As our non-matched evolution always satisfy µ 1 /µ 2 1, we are always in this regime. For this reason, when simulating our Pythia results to compute the low scale FFs, we always initiate these at √ s = 2q W for states XX, with the X andX color connected to form a singlet. There will be corrections to this picture at higher order, although this is sufficient for the level of accuracy of the present work. In practice we evaluate the FFs over a wide dynamic range, and thus it is convenient to change variables to l x = ln x, yielding i.e. the high-scale values are described byd, and the low scale by d.
To determine the convolution for a given l x , it is convenient to define three additional quantities. The first is simply the bin width l n+1 − l n = ∆. The next two specify the location of l x . We introduce an integer m defining the bin l x falls in, in detail l m < l x < l m+1 . From here, the fractional bin width to the point l x , is defined as δ = (l x − l m )/∆ ∈ [0, 1]. In terms of these auxiliary quantities, we can then evaluate Eq. (G.3), which can be readily evaluated numerically. Note if m = N the summation within square brackets above expression vanishes. As a simple check, imagine µ 2 = µ 1 , so thatd b a,n = δ b a δ N n . Then the above becomes, which is an appropriately weighted sum.