{'fraction_non_alphanumeric': 0.036121958546109945, 'fraction_numerical': 0.009942925803544608, 'mean_word_length': 4.295132813742644, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 24, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Weak Corrections are Relevant for Dark Matter Indirect Detection

The computation of the energy spectra of Standard Model particles originated from the annihilation/decay of dark matter particles is of primary importance in indirect searches of dark matter. We compute how the inclusion of electroweak corrections significantly alter such spectra when the mass M of dark matter particles is larger than the electroweak scale: soft electroweak gauge bosons are copiously radiated opening new channels in the final states which otherwise would be forbidden if such corrections are neglected. All stable particles are therefore present in the final spectrum, independently of the primary channel of dark matter annihilation/decay. Such corrections are model independent.


Introduction
There are overwhelming cosmological and astrophysical evidences that our universe contains a sizable amount of Dark Matter (DM), i.e. a component which clusters at small scales. While its abundance is known rather well in terms of the critical energy density, Ω DM h 2 = 0.110±0.005 [1], its nature is still a mistery. Various considerations point towards the possibility that DM is made of neutral particles. If DM is composed by particles whose mass and interactions are dictated by physics in the electroweak energy range, its abundance is likely to be fixed by the thermal freeze-out phenomenon within the standard Big-Bang theory. DM particles, if present in thermal abundances in the early universe, annihilate with one another so that a predictable number of them remain today. The relic density of these particles comes out to be: where σv ann is the (thermally-averaged) cross annihilation cross sections. A weak interaction strength provides the abundance in the right range. This numerical coincidence represents the main reason why it is generically believed that DM is made of weakly-interacting particles with a mass in the range (10 2 − 10 4 ) GeV. There are several ways to search for such DM candidates. If they are light enough, they might reveal themselves in particle colliders, such as the LHC, as missing energy in an event. In that case one knows that the particles live long enough to escape the detector, but it will still be unclear whether they are long-lived enough to be the DM [2]. Thus complementary experiments are needed. In direct detection experiments, the DM particles elastically scatter off of a nucleus in the detector, and a number of experimental signatures of the interaction can be detected [3]. In indirect searches DM annihilations or decays around the Milky Way can produce Standard Model (SM) particles that decay into e ± , p, p, γ and d , producing an excess in their cosmic ray fluxes. Present observations are approaching the sensitivity needed to probe the annihilation cross section suggested by cosmology, eq. (1). Furthermore, this topic recently attracted interest because the PAMELA experiment [4] observed an unexpected rise with energy of the e + /(e + + e − ) fraction in cosmic rays, suggesting the existence of a new positron component. The sharp rise might suggest that the new component may be visible also in the (e + + e − ) spectrum: although the peak hinted by previous ATIC data [5] is not confirmed, the FERMI [6] and HESS [7] observations still demonstrate a deviation from the naive power-law spectrum, indicating an excess compared to conventional background predictions of cosmic ray fluxes at the Earth. While the current excesses might be either due to a new astrophysical component, such as a nearby pulsar [8], or to some experimental problem, it could be produced by DM with a cross section a few orders of magnitude larger than in eq. (1), maybe thanks to a Sommerfeld enhancement [9,10].
In any case, it is undeniable that nowadays indirect search of DM is a fundamental topic in astroparticle physics, both from the theoretical and experimental point of view. Computing the energy spectra of the stable SM particles that are present in cosmic rays and might originate from DM annihilation/decay is therefore of primary importance.
The key point of this paper is to show that electroweak radiative corrections have a sizable impact on the energy spectra of SM particles originated from the annihilation/decay of DM particles with mass M somehow larger than the electroweak scale. The reason is in fact simple and should be familiar to readers working in collider physics: at energies much higher than the weak scale (in our case the mass M of the DM) soft electroweak gauge bosons are copiously radiated from highly energetic objects (in our case the initial products of the DM annihilation/decay). This emission is enhanced by ln M 2 /M 2 W when collinear divergences are present and ln 2 M 2 /M 2 W when both collinear and infrared divergences are present [11]. These logarithmically enhanced terms can be computed in a model-independent way through the well known partonic techniques based on the Collinear Approximation (CA). Our work will involve generalizing the partonic splitting functions to massive partons, because our 'partons' include the W, Z bosons.
Putting these technical details aside, what is important is that the emission of gauge bosons changes significantly some final energy spectra. Indeed, suppose that DM annihilates into a pair of leptons. The emitted gauge bosons give hadrons (resulting in ap flux) and mesons (giving a significant extra amount of photons via π 0 → γγ). The total energy gets distributed among a large number of lower energy particles, thus enhancing the signal in the lower energy region (say, (10 − 100) GeV), that is measured by present-day experiments, like PAMELA. This paper will be inevitably rather technical and therefore we have decided to defer as many as possible technicalities to the various Appendices. To diminish the burden, we present qualitative considerations in Section 2 and outline the quantitative computation in Section 3. The reader interested in the final results may jump directly to Section 4 where our findings are presented. Conclusions are presented in Section 5. In Appendix A we discuss EW evolution equations, in Appendix B we derive parton splitting functions for massive partons, and in Appendix C we list all splittings among SM particles, including the effect of the top Yukawa coupling.

Qualitative discussion
As mentioned in the Introduction, the presence of DM is probed indirectly by detecting the energy spectra of stable particles (p ± , e ± , ν, γ,d). At a first sight, since electroweak radiative corrections are expected to be small -weak interactions are weak, after all -they might seem to play no role in the DM indirect searches. At the typical weak scale of O(100) GeV, radiative corrections produce relative effects of O(0.1)%. For instance, this was the case for experiments that took place at the LEP collider. However at energies of the order of the TeV scale, like those probed at the LHC, things are different: electroweak radiative corrections can reach the O(30)% level [12] and they grow with energy, eventually calling for a resummation of higher order effects [13]. In a nutshell, what happens at energies much higher than the weak scale is that soft electroweak gauge bosons are copiously radiated from highly energetic objects that undergo a scattering with high invariant mass. This is much the same as photon, or gluon, radiation whenever the hard scale is such that the W, Z masses can be safely taken to be very small. Important differences with respect to unbroken gauge theories like QCD and QCD arise in the case of a (spontaneously) broken theory like the EW sector of the SM. It was found [14] that in hard processes with at least two relativistic non abelian charges, effective infrared divergencies that are manifest as double log corrections (α 2 ln 2 M 2 /M 2 W ) appear. They  Figure 1: DM annihilation/decay initially produces a hard positron-electron pair. The spectrum of the hard objects is altered by electroweak virtual corrections (green photon line) and real Z emission. The Z decays hadronically through a qq pair and produces a great number of much softer objects, among which an antiproton and two pions; the latter cascade decay to softer γs and leptons.
are not present in QED and QCD and this effect has been baptized "Bloch-Nordsieck Theorem Violation" [14]. We refer the reader to the relevant literature [14,15,16] for details. In the case at hand, since the initial DM particles are nonrelativistic, radiation related to the initial legs does not produce log-enhanced terms. Therefore, we only need to examine soft EW radiation related to the final state particles. The hard scale in the case we examine here is provided by the DM mass M > ∼ 1 TeV while the soft scale is the typical energy where the spectra of the final products of DM decay/annihilation are measured, E < ∼ 100 GeV. Even bearing in mind that weak interactions are not so weak at the TeV scale, one might wonder whether such "strong" electroweak effects are relevant for measurements with uncertainties very far from the precision reachable by ground-based experiments at colliders. In this context, and in view of our ignorance about the physics responsible for DM cross sections, it might seem that even a O(30)% relative effect should have a minor impact. This is by no means the case: including electroweak corrections has a huge impact on the measured energy spectra from DM decay/annihilation. There are two basic reasons for this rather surprising result.
• In the first place, since energy is conserved, but the total number of particles is not, because of electroweak radiation a small number of highly energetic particles is converted into a great number of low energy particles, thus enhancing the low energy ( < ∼ 100 GeV) part of the spectrum, which is the one of relevance from the experimental point of view.
• Secondly, and perhaps more importantly: since all SM particles are charged under the SU(2) L ⊗ U(1) Y group, including electroweak corrections opens new channels in the final states which otherwise would be forbidden if such corrections are neglected. In other words, since electroweak corrections link all SM particles, all stable particles will be present in the final spectrum, independently of the primary annihilation channel considered.
To illustrate these facts, consider for instance a heavy DM annihilation producing an electronpositron pair, see Fig. 1. Clearly, as long as one does not take into account weak interactions, only the leptonic channel is active and no antiproton is present in the final products. However, at very high energies there is a probability of order unity that the positron radiates a Z or a W . While the spectrum of the hard positron is not much altered by virtual and real radiative corrections (see [17]), the Z radiation opens the hadronic channel: for instance, antiprotons are produced in the Z decay. Moreover, also a large number of pions are produced, which in turn decay to photons (π 0 → γγ) and to low energy positrons (through the chain π + → µ + + X → e + + X). At every step, energy is degraded. Because of the large multiplicity in the final states, the total Z energy (already smaller than the hard M scale) is distributed among a large number of objects, thus greatly enhancing the signal in the (10 − 100) GeV region that is measured by present-day experiments, like PAMELA. The various processes of radiation are described by fragmentation functions D EW (x, µ 2 ) that evolve with the energy scale µ 2 according to a set of integro-differential electroweak equations [18]. When a value of virtuality of the order of the weak scale µ = M W is reached, the Z boson is on shell and decays. The subsequent QCD showering may be described with QCD traditional MonteCarlo (MC) generator tools, like PYTHIA.
At tree level, the spectra of hard objects emerging from DM annihilation are simply proportional to δ (1 − x), where x is the fraction of center-of-mass energy carried by a given particle. Once electroweak corrections are switched on and O(α 2 ) virtual and real corrections are calculated, the spectra D EW (x, µ 2 = M 2 W ) generically contains terms enhanced by log terms of the form α 2 ln 2 M 2 /M 2 W and α 2 ln M 2 /M 2 W . The presence of logarithmically enhanced terms is well-known in the literature both in the case of electroweak interactions [14,18] than for strong interactions through the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) equations [19] and is related to regions of phase space where the internal propagators become singular. In general the single log term is generated when two partons become collinear, while the double log arises when they are soft and collinear at the same time. The double log-enhanced contributions cancel in QED and QCD for physical scattering processes (Block-Nordsieck theorem [20]), while they are present in massive gauge theories [14].
The physical picture that arises is therefore the following: highly non relativistic DM particles annihilate, producing a particle-antiparticle pair belonging to the SM spectrum. This pair has a very high invariant mass, therefore it stars radiating photons and gluons, but also weak gauge bosons. The presence of collinear and/or infrared singularities allows to factorize leading logarithmic electroweak corrections with a probabilistic interpretation very similar to DGLAP equations, see Sec. 3. The exchange of virtual and emission of real electroweak bosons lead to the appearance in the final spectrum of all the stable SM particles, not only the ones initially emitted by the DM annihilation. Indeed, the higher is the mass of the DM, the more democratically distributed the final spectrum of DM particles is. Therefore, including electroweak corrections alters significantly the final spectrum of particles stemming from DM decay/annihilation and this has a large impact on indirect searches of DM.
Let us close this Section by recalling that, while in this paper we only consider DM annihilation/decay to two body final states, our approach is more general and model independent. Indeed, the only assumptions we make are that the physics up to the DM mass scale M is described by the SM and that the SM may be eventually extended by interactions that preserve SU(2) L ⊗ U(1) Y gauge invariance. While these assumptions exclude cases like the ones considered in Refs. [10] and [21] where gauge non invariant interactions where considered 1 , a large number of models can be examined with the techniques we describe here. For instance, let us extend the SM by adding a very heavy scalar S that interacts with the SM Higgs (H) and leptons (L, E), through an effective operator SLEH. Then, the dominant decay of the scalar is a three body decay, since the two body decay S → LE is suppressed by a relative factor M 2 W /M 2 . The framework described here applies as well, albeit with the additional complication that the three body decay with respect to which one factorizes electroweak interactions provides a distribution rather than a simple δ function. In this sense, provided assumptions specified above are fulfilled, our approach is completely model independent.

Quantitative computation
We now discuss in more technical terms the inclusion of EW gauge boson emission through the evolution equations. We start from a first principle definition of the energy spectrum for emitted particles and then we define the fragmentation functions as statistical objects describing the probability of a particle to be transformed into another with a certain momentum fraction. The full evolution equations for the fragmentation functions, containing EW and QCD interactions, are analyzed. We provide an expression that can be used to match the outcome of Monte Carlo codes adding EW corrections at leading order O(α 2 ). Our approach is similar in spirit to the one used, in a different context, in [23], with important differences.
The relevant quantity for indirect signals of DM is the energy spectrum dN f /dx of stable SM particles f = {e + , e − , γ, p,p, ν,ν,d} produced per DM decay/annihilation, where x = 2E f / √ s (0 ≤ x ≤ 1) is the fraction of center of mass energy carried by a stable particle f with energy 1 The analysis of the Infrared virtual corrections to gauge non invariant amplitudes, i.e. amplitudes proportionals to the higgs vev, has been recently performed in [22]. E f . For clarity we will sometimes specify the formulae assuming the case of non-relativistic DM annihilations, for which √ s = 2M such that x = E f /M ; it is immediate to obtain the corresponding formulae for DM decays, where √ s = M .
We assume that DM initially produces two primary back-to-back SM particles, and we consider all relevant cases: where q = u, d, s denotes a light quark; h is the Higgs boson; Left or Right are the possible fermion polarizations, and T ransverse or Longitudinal are the possible polarizations of massive vectors, that correspond to different EW interactions. Then, the spectrum can be written as: with a similar formula for the case of DM decay. The "X" in this equation reminds of the inclusivity already discussed in Section 2.
In each one of the possible cases I, MonteCarlo generators like PYTHIA allow to compute the inclusive spectrum dN MC I→f (M, x)/dx by generating events starting from the pair I of initial SM particles with back-to-back momentum and energy E = M , and letting the MC to simulate the subsequent particle-physics evolution, taking into account decays of SM particles and their hadronization, as well as QCD radiation and (partially) QED radiation.
Then, the spectra for a generic DM model that produce combinations of the two-body states I can be obtained combining the various channels: In some DM models, primary multi-body states can be important: one can obtain the final spectra without running a dedicated MC code by computing the model-dependent energy spectra D I (z) of each primary pair I (each one has energy E = zM with 0 ≤ z ≤ 1) and convoluting them with the basic MC spectra: Notice that we combine particle-antiparticle pairs because we assume that they have the same spectra, which is true whenever the cosmological DM abundance does not carry a CP asymmetry. Otherwise, hadronization can be significantly affected and dedicated MC runs would be necessary. The indices I, J = p +p denotes a primary particle p together with its antiparticlep, with the same energy spectrum. Factors of two are such that for complex particles one has dN p /dz = dNp/dz = D I , while for real particles (the Z, the γ, the Higgs h) one has dN DM→p /dz = 2 D I .

Including EW corrections
We now come back to the basic case of DM that annihilates or decays in one primary channel I and discuss how to achieve the goal of this paper: obtaining a set of basic functions dN I→f /dx that take into account EW radiation, replacing the functions dN MC I→f /dx computed via Monte-Carlo simulations. EW radiation is a model-independent subset of the higher order corrections discussed above, and gives rise to specific spectra of initial SM particles, such that its effect can be included in the primary basic spectra by a formula similar to eq. (5): where D EW I→J (z) is the EW I → J EW parton distribution: the J spectrum produced by initial I. Our normalization is such that we have the uniform normalization D EW I→J (z) = δ IJ δ(1 − z) at tree level for both real and complex particles. Some comments are in order: i) When including higher order effects, one must avoid overcounting and take into account that MC codes already include some particularly relevant higher-order effects: showering produced by strong (QCD) and electromagnetic (QED) interactions, up to details. 2 ii) For initial particles that do have strong interactions, eq. (6) misses the interplay between EW and QCD radiation; this limitation is not a problem because, as expected, in such cases EW radiation will turn out to be subdominant with respect to QCD radiation.
iii) For initial particles that do not have strong interactions, eq. (6) holds at leading order in the weak couplings: first they must do an EW splitting, and next one can add QCD splittings neglecting EW radiation. We emphasize an important different between e.g. a Z → qq splitting and the same Z → qq decay: in the decay the invariant mass of the qq pair is equal to the Z mass (such that Z → tt is forbidden by the heaviness of the top t); in the splitting the invariant mass can be much higher, and we approximate it as zM . This higher invariant mass strongly affects the subsequent QCD radiation from quarks, which is more abundant in the splitting case, leading to a higher multiplicity ofp and γ.
In Appendix A we give a detailed discussion of the interplay between EW and QCD radiation and of the level of approximation introduced by using eq. (6).

Computing the EW parton distributions
We define D EW I→J (z, µ 2 ) as the probability for a given parent particle I with virtuality of the order of µ to become a particle J with a fraction z of the parent particle's energy mediated by 2 We must include via eq. (6) only all those effects not included in MC codes. Existing MC codes have their own peculiarities, e.g. Phytia automatically includes γ radiation from charged particles but not from the W ± . Of course, an alternative more precise approach, that we do not purse, would be implementing the missing EW radiation effects into some existing MC code. EW interactions. At large virtuality, they take the tree level values: At low virtuality µ 2 ∼ M 2 W , they are the functions we need: . The evolution in the virtuality is described by integro-differential equations, that involve a set of kernels 3 P EW I→J (x, µ 2 ) that have been derived in [18]: Since we work at leading order in the EW couplings, eq. (8) with the boundary conditions of eq. (7) is solved by: Differently from QED and QCD, the EW kernels feature infrared singular terms proportional to ln µ 2 , so that the solutions (9) also include double logs beside the customary single logs of collinear origin: Our goal is to include the model-independent logarithmically enhanced terms. Electroweak radiation from the initial DM state is of course model-dependent: since DM is non-relativistic this effect only contributes to the non-enhanced terms D 0 , that we neglect. Notice however that for our purposes we need to include terms of the form (ln x)/x that are relevant in the region x → 0; this is discussed in detail in subsection 3.4.

Splitting functions
The leading-order parton distributions D I→J (z) can be computed by using the partonic splitting functions P summing over all possible SM splittings [18]; the relevant splitting functions are here collected in table 1.
A concrete simple example allows to clarify the procedure and the normalization factors: we consider DM producing an initial generic F ermion-antiF ermion pair with FF invariant mass √ s m F . We assume that F has charge q F under a generic vector V with mass M V √ s and gauge coupling α. (In the SM the vector could e.g. be a Z and the fermion a neutrino). The splitting process F → F V gives rise to: (11) 3 In this work we indicate with P (x, µ 2 ) the unintegrated kernels, while the splitting functions P (x), obtained by integrating in µ 2 , depend only on the energy fraction x; a list of the relevant splitting functions is given in Table 1.
splitting function: real and virtual  By replacing F → S one obtains the corresponding result for a pair for Scalars, and so on. The first term describes virtual corrections arising from one-loop diagrams, and the second term describes real emission corrections. We define: for any I and J; e.g. from P V →V in table 1 we have: While both virtual corrections, related to P vir I→J , and real corrections, related to P I→J , are listed in table 1, the simple relation (12) holds between them. The relationship between real and virtual contributions is dictated by the unitarity of the theory (see e.g. [17] for a more detailed discussion). Intuitively, this just amounts to say that when an F radiates, it disappears from the initial state.
One can verify that the partonic distributions D I→J satisfy a set of conservations laws (with corresponding identities for the splitting functions P ) • The conservation of splitting probability: where the factor 1/2 accounts for the fact that one particle splits in two.
• The conservation of total momentum: • The conservation of electrical charge: By combining these splitting functions with the appropriate electroweak couplings (including the top quark Yukawa coupling) we get the electroweak splittings among SM particles, described by the D I→J functions, explicitly listed in Appendix C. It is a simple exercise to verify that conservation laws are satisfied. For completeness, we also list the splittings involving photons (to be dropped if already included by MC codes), assuming for simplicity a photon with mass M W ; sending it to zero gives infrared divergences that can be regulated and dealt with using well known techniques that we do not need to discuss here.

Splitting functions for massive partons
We see from table 1 that P vir I→J is explicitly given by linear or quadratic polynomials in ln s/M 2 V : the vector mass provides an infra-red regulator. This is unlike in the standard application of partonic techniques, where all partons are massless and some infra-red regularization is needed.
In this subsection we briefly describe how we generalized partonic functions to massive partons (such as the W, Z, t in the SM), and why particle masses modify the kinematics affecting the log-enhanced terms, as encoded in a non trivial universal function L.
Let us consider the splitting process i → f + f of a particle i into massive partons f, f : the corresponding splitting functions can be non zero only in the kinematically allowed ranges: where x = E f /E i , and similarly for x = 1 − x = E f /E i . The modification due to nonzero masses in the kinematical range is power suppressed and thereby mostly negligible; however 0.01 Thereby small parton masses must strongly affect such splitting functions, replacing the usual 1/x and 1/x divergent terms with new functions that vanish at the kinematical boundary.
We therefore need to correct the expressions for the splitting functions that are singular in the limits x → 0 and x → 0 (soft f and soft f ). These singularities arises when soft vectors are emitted, such that the correct expressions can be obtained using the eikonal approximation, which dictates the amplitude for emitting a soft vector. Since we stop at leading order, we can precisely define the splitting functions in terms of the energy spectra of the particles produced in three body scatterings: we just need to integrate such squared amplitude over the massive phase space in the region that produces the leading singularities. Broadly speaking, what happens is that the upper and lower limits of integration that appear in (9) are modified and become x-dependent; computations and checks are performed in Appendix B. The resulting splitting functions for massive partons are listed in table 1, and, for x → 0, they do not depend on spin, e.g. P S→V P F →V P V →V , as dictated by the eikonal amplitude. They contain the universal kinematical function L(x), that replaces the usual = ln s/M 2 V that holds for massive partons: This function indeed vanishes below the kinematical threshold (x < 2M V / √ s) and reduces to ln sx 2 /M 2 V well above it. This means that small parton masses, apart from providing kinematical   2 shows the splitting function between massive vectors (e.g. relevant for Z radiation from W ± ): it significantly differs from the massless splitting function (upper curve) even away from the kinematical boundaries, and it closely agrees with the full result of a 3-body computation in a specific model, which also includes not log-enhanced terms. These terms are subleading in the present example, where we here considered DM DM → W + T W − T annihilations with M = 3 TeV.

Results
Our main results are the energy spectra of all stable final SM particles f from any two-body DM non-relativistic annihilations or decays. We will make the code freely available in [24], and show here some examples. In table 2 we show the total number of e + ,p and γ, including EW radiation, for a few values of the DM mass M . Without EW radiation, increasing M would just boost all particles, such that dN/dx and consequently N does not depend on M . With EW radiation, the total number of particles increases with the DM mass: the typical effect is about a factor of 2 going from 0.3 to 3 TeV. Of course, such dependence on M is similar to the one already present in the quark channels (bottom row), due to QCD radiation.
Coming to the energy spectra, these new particles appear at low energy. In Fig. 3 we consider annihilations of DM with mass M = 3 TeV and compare the spectra with (continuous curves) and without (dashed curves) EW corrections.
• In the top row we consider DM annihilations into W + W − with transverse (left) or longitudinal (right) polarization. The final spectra do not significantly depend on the polarization of the W , but this is an accident: a T ransverse W splits into all quarks as demanded  : Comparison between spectra with (continuous lines) and without EW corrections (dashed). We show the following final states: e + (green),p (blue), γ (red), ν = (ν e + ν µ + ν τ )/3 (black). by their weak interactions, while a Longitudinal W (being the charged Goldstone in the Higgs doublet) mainly splits into t and b quarks, as demanded by the top quark Yukawa interaction. In both cases EW splittings increase the total number of final particles (e + ,p, γ, ν) increases by a factor of almost 2. The biggest effect is in thep spectra: without EW corrections they are strongly suppressed below E p < m p · (M/M W ) ∼ 100 GeV, because of the boost factor M/M W of the W . Adding EW corrections, a W splits into quarks with lower energy, producingp also at lower energy. This is important for interpretations of the PAMELA e + excess in terms of annihilations of very heavy DM particles, e.g. M = 10 TeV as predicted by Minimal Dark Matter [25]. EW corrections make more difficult to avoid unseen effects inp by having all the DM effects above 100 GeV where we do not have data yet. As illustrated in Fig. 4, in view of EW corrections there is some tension with data, even assuming the MIN model of diffusion of charged particles in the galaxy [26] (which gives the minimal amount ofp compatibly with cosmic rays data and theory).
• We do not plot our result for DM annihilations into the Higgs boson, because its mass and decay modes are not yet known. The Higgs channel is similarly affected by EW corrections as the W L and Z L channels (because they are the Goldstone components of the Higgs doublet), and the main effect is again the h → t splitting induced by the top Yukawa coupling.
• In the middle row of Fig. 3 we consider DM annihilations into charged leptons: e − L e + L (left) and µ − L µ + L (right). The spectra are significantly affected by EW corrections, because leptons split into W, Z bosons, finally producing quarks, and consequently a copious tail of e + , γ,p at lower energy. The main new qualitative feature is the appearance ofp from leptons, and in Fig. 4 (lower row) we show that DM annihilations into µ + L µ − L with M = 2 TeV (a scenario motivated by the PAMELA e + and FERMI e + + e − anomalies) also gives a possibly detectable excess ofp. We here assumed the favored MED model of propagation of charged particles in the Milky Way [26];p can be suppressed down to a negligible by considering the MIN model and/or DM annihilations into µ + R µ − R . Indeed EW effects are more significant for left-handed leptons which have SU(2) L interactions than for right-handed leptons which only have U(1) Y interactions.
• The bottom left panel of Fig. 3 shows that, in view of EW interactions, DM annihilations into neutrinos also induce a significant spectrum of γ and some e + ,p.
• The bottom right panel of Fig. 3 shows how EW corrections affect DM annihilations into γγ: an hypothetical γ line at high energy E γ = M ∼ few TeV must be accompanied by a comparable flux of γ with lower energy E γ ∼ (10 − 100) GeV, where we have more data. 4 • Finally, DM annihilations and decays into quarks are negligibly affected, so that we do not show our results.
The same results hold for DM decays. EW corrections are irrelevant in models where DM decays or annihilates into hypothetical 'dark' particles, lighter than the weak scale, that finally decay back to SM particles.
There is a striking difference between the two examples of e + fraction in Fig. 4: EW corrections are relevant in the upper case (W T channel) and negligible in the lower case (µ L channel). At first sight, this is surprising, because Fig. 3 shows that in both cases EW corrections induce significant low energy tails to the e + spectra at production. The qualitative difference is due to e + energy losses in the Milky Way that also gives a tail of e + at low energy: we explain how to understand which effect is dominant. The e + flux can be approximatively computed neglecting galactic diffusion and taking into account only energy losses, with the result where σ T is the Thompson cross section, u ⊕ is the energy density in radiation and magnetic fields and ρ ⊕ is the DM density, both at the location of the solar system. This means that the e + flux at a given energy E < M is proportional to the number of e + produced with energy between E and M . In the case of DM annihilations into leptonic channels, the tail of e ± at lower energy produced by EW radiation contains a number or e ± at most comparable to the amount of e ± already present at E ∼ M : therefore EW corrections negligibly affect the positron fraction. This is not the case for DM annihilations into W ± , where instead the tail at low energy is the dominant component of the total e ± number, such that EW corrections more significantly affect the e + fraction.

Conclusions
In this paper we computed the energy spectrum of SM particles stemming out from DM annihilation/decay. We have shown that EW corrections have a relevant impact on such spectra when the mass M of the DM particles is larger than the EW scale M W . Soft EW boson emission is enhanced in the collinear and infrared regime and this leads to ln 2 M 2 /M 2 W enhancement factors. The result of the inclusion of EW corrections is that all stable particles are present in the final spectrum, independently of the primary annihilation/decay channel. For instance, even if the DM particles annihilate/decay into light neutrinos, EW corrections cause the appearance of hadrons and photons in the very final spectrum. The inclusion of EW corrections is therefore an essential ingredient in order to have a physical picture of the correlated energy spectra of final stable particles. Our quantitative results may be inferred from Fig. 3 where the energy spectra dN/dE of stable SM particles are presented with and without EW corrections. These spectra are the necessary ingredients to predict the flux for indirect searches once the effect of diffusion and galactic energy loss are included. As a rule of thumb we may say that EW corrections are important in determining the final flux of stable particles whenever • the final flux of stable particles is dominated by the low energy tail of the dN/dE. One example is the case of DM annihilation/decay into gauge bosons and e ± final states (see fig. 3). This point becomes more relevant in the present experimental situation, where we mostly observe particles below 100 GeV, possibly much below the DM mass.
• the final flux of stable particles is absent when EW corrections are not taken into account. One example is the case of DM annihilation/decay into leptons and antiprotonsp final states (see fig. 4). This point is important also for neutrino fluxes from DM annihilations in the sun or in the earth, because all SM particles, even those that loose energy in matter before decaying into neutrinos, can radiate a W or a Z that promptly decays into neutrinos.
EW corrections may also significantly affect the fluxes of particles generated from heavy gravitino decays in supersymmetric theories. The computation of these energy fluxes is crucial in studying the dissociation of the light element abundances generated during a period primordial nucleosynthesis.
We expect that EW radiative corrections have a minor effect on the freeze-out cosmological DM abundance, because it is determined by just the total non-relativistic annihilation DM cross section. In the case of energy spectra instead, as explained in this paper, the low energy tails can be enhanced by orders of magnitude, while the high energy part of the spectrum is mildly depleted. The net effect on the total number of final particles typically is an enhancement by a factor of 2. We conclude that, when DM is around or above the TeV scale, one must take into account radiative EW corrections.
We computed EW corrections at leading order. Although we cannot give a sharp answer, we argued that it is not necessary to resum higher order EW corrections as long as DM is not too heavy: α 2 ln 2 (M 2 /M 2 W )/2π 1 at M 100 TeV. Indeed, higher order corrections are expected to give small effects for total cross sections at the TeV scale: a one loop effect of the order of 30% means that one expects higher order effects to be at the 1% level. It is more difficult to guess how low-energy tails might be affected without performing a dedicated computation, which could be done implementing EW corrections in a MonteCarlo: we provided splitting functions for massive partons and all ingredients.

Acknowledgments: D.C. during this work was partially supported by the EU FP6 Marie Curie
Research & Training Network "UniverseNet" (MRTN-CT-2006-035863). P.C. wishes to thank Luigi Lanzolla for many useful contacts and discussions, without which his work would have been much more difficult. We also thank G.F. Giudice for useful conversations. We used the high-statistics Pythia spectra computed on the Baltic Grid by Mario Kadastik, to appear in [24]; this work was supported by the ESF Grant 8499 and by the MTT8 project.

A Evolution Equations
At present the energy spectrum dN MC I→f /dx are computed with MC generators like PYTHIA by generating events starting from the pair I of initial SM particles with back-to-back momentum and energy E = M , and letting the MC to simulate the subsequent particle-physics evolution, taking into account decays of SM particles and their hadronization, as well as QED and QCD radiation. As it is evident the MC output are nothing else that the full QCD+QED fragmentation functions D QCD+QED I→f (x) that can be identified with dN MC I→f /dx. The fragmentation functions are related to a nonperturbative aspect of QCD, so that they cannot be precisely calculated by theoretical methods at this stage. The situation is similar to the determination of the PDFs, where high-energy experimental data are used for their determination instead of theoretical calculations. The µ 2 evolution for the fragmentation functions is calculated in the same way as the one for the PDFs by using the timelike DGLAP evolution equations. The splitting functions are the same in the Leading Order (LO) evolution of the PDFs; however, they are different in the Next-to-Leading Order (NLO). Explicit forms of the splitting functions are provided in [18]. The evolution equations are then essentially the same as the PDF case.
The correct determination of the energy spectrum dN I→f /dx of the final stable particle f , obtained through the particle-physical evolution of the initial pair I of SM particles with back-to-back momentum and energy E = M , needs the solution of a full set of evolution equations, including both strong and electroweak interactions. In previous works the QCD DGLAP formalism has been extended to EW interactions [18]. The analysis of mass singularities in a spontaneously broken gauge theory like the electroweak sector of the Standard Model has many interesting features. To begin with, initial states like electrons and protons carry nonabelian (isospin) charges; this feature causes the very existence of double logs, i.e. the lack of cancellations of virtual corrections with real emission in inclusive observables [14]. Secondly, initial states that are mass eigenstates are not necessarily gauge eigenstates; this causes some interesting mixing phenomena analyzed in [16].
On a general setting, the mathematical structure of the full system of EW and QCD evolution is provided by the following set of integral differential equations with kernels P QCD/EW for different gauge boson exchange (we omit for the moment the explicit QED evolution equations included in the EW part of the SM and the various index parametrizing the flavour and quantum numbers): 5 where we have made explicit the running range: for the EW interactions from M to M W where the gauge boson mass is freezing the running and for QCD corrections from M to Λ QCD where the nonperturbative effects generate an effective cutoff to gluon exchanges. In practice, from Λ QCD to M W the running is purely dictated by QCD (remember that QED interactions are not shown for simplicity, in this case the running scale for photons is until the electron mass m e ) while the EW-QCD interplay starts only above the M W scale. Note also that forgetting the underlying index structure can bring to the wrong conclusion that, simply due to the fact that α s > α 2 , the EW part is just a correction to the QCD dominant part. This would be a mistake as in many interesting annihilation channels (like all the leptonic ones or the W ± , Z or h) the QCD part is simple zero. 5 The precise index (flavour) structure is given in the Appendix C, the generic index structure is of the A numerical solution to the full (EW+QCD) problem is of course out of reach. Nevertheless, we can find some reasonable approximation taking advantage of the fact that we can simulate the pure QCD evolution also in the non perturbative regime with MC codes and that EW theory is in the perturbative regime. Technically specking, the evolution equations are Schrödinger-like equations with a time dependent Hamiltonian, where time is replaced by the µ 2 variable and the Hamiltonian by the P -kernel. A formal solution can be parametrized with the evolution operator: where P µ 2 is the µ 2 -ordering operator and I the identity in the flavour space. Due to the linearity of the Eqs. (A.1) 6 we can then formally write the full solution as: where we have separated the running from M to M W inside the evolution operator U (that can be perturbatively expanded as soon as α s ln M 2 /M 2 W and α 2 ln 2 M 2 /M 2 W are smaller than unity) and the purely QCD piece D QCD (x, M 2 W , Λ 2 QCD ) encoding also the non-perturbative low energy physics. In order to keep under control further simplifications we need also to know the matrix flavour structure. We display it under the form of a simplified four-dimensional space spanned by l=leptons, W = (W ± , Z, γ), q=quarks and g=gluons, for the EW and QCD kernels, P QCD and P EW : The above matrices do not commute; the EW and QCD sectors are connected through the channels W → q, q → W and q → q, furthermore the leptonic and hadronic sectors are connected through the mixed W → q, l, l, q → W channels. Reasonable approximate solutions are related both to the possibility to expand perturbatively the general solution (B.15), and to the outcome from MC generators which take automatically into account the full QCD plus QED evolution from M to Λ QCD scales (m e for QED). One way of proceeding is to define pure EW (D EW ) and QCD (D QCD ) fragmentation functions that evolve with their respective kernels, see eq. (A.4), for the energy range M 2 W < µ 2 < M 2 : whose formal solutions are: Then we introduce a new factorized EW ⊗ QCD fragmentation function: This is clearly not a solution of the true evolution equations (A.1) but can be a useful approximate solution. In order to relate the true solution D of eq. (A.3) with the new function D satisfying eq. (A.7)), we can use the fact that in the M 2 W < µ 2 < M 2 interval we are in perturbative regime also for the QCD side. Knowing that for two generic non-commuting operators A and B: (A.8) and identifying A = 2π P QCD , we can approximate the evolution operator U at any order in α s,2 . In particular at second order in O(α 2 s,2 , α s α 2 ) we have: where: (A.10) and the · · · stand for the fact that there is an infinite series of commutators with coefficients of order α m+1 2 α n+1 s with m + n ≥ 1. Starting from this expression we can write the perturbative relation between the exact solution D and the present outcome of the MC codes D QCD : The first order corrections to such a formula are obtained expanding also D EW at one loop: where we have explicitly shown the arguments x and µ of the matrix P EW . At order O(α 2 ln 2 M 2 /M 2 W , α s ln M 2 /M 2 W ) we find that the energy spectrum for the process I → f +X can be therefore written as in eq. (6): .13) where in the last passage, to be consistent, we have written EW in order to stress that only massive W ± and Z are included while QED is already encoded in the MC. Through this expression one can match the MC code with the first order EW corrections. Figure 5: Soft gauge boson real emission from spin-1 particle. It can be read as the sum of three scalar currents. Considering the process √ s → p 1 + k + p 2 we show explicitly only the bremsstrahlung contribution from the p 1 final leg.

B Eikonal approximation and the improved splitting functions
The standard partonic approximation holds in QED and in QCD for the emission of soft massless gauge bosons (photons or gluons) from partons, showing the presence of universal logarithmical factors of collinear origin that multiply the usual splitting functions. This approach can't be naïvely applied to the electroweak case when a massive gauge boson, such as the W , is involved in the splitting process; considering for definiteness i → f + f and defying x ≡ E f /E i , in fact, the allowed kinematical range for this latter is: where particle masses act as cut-off for the soft singularities at x → 0, 1. These boundary regions in (B.1) are therefore extremely important and the standard partonic approximation have to be improved, introducing extra ln x and ln(1 − x) terms, well justified by the kinematical proprieties of the splitting process.
In this Appendix we derive our improved splitting functions for massive partons, following the logic outlined in section 3. In B.1 we use the eikonal approximation, that describes the amplitudes with soft gauge boson emission. We integrate it over the phase space using for it the Sudakov approximation in B.2 and using the exact expression in B.3: the Sudakov parametrization, commonly used in literature, don't respect the boundaries in (B.1). In B.4 we introduce, through an explicit example, the collinear approximation and its propriety of factorization. Finally, in B.5 we compare our results with those of a full full three body calculation (exact amplitude integrated over the exact phase space).

B.1 The eikonal amplitude
As well known, the spin of the emitting particles (scalar, fermion and vector) becomes irrelevant in the eikonal limit: for definiteness, and without losing generality we here consider the real emission of a particle with momentum k described through the three gauge boson vertex 3g.
Considering Fig. 5 and using the conservation of momenta in the splitting vertex p = p 1 + k we can write: where we have indicated with iM ρ the remaining part of the amplitude and taken for simplicity k 2 = 0. Roughly speaking using the eikonal approximation it is possible to neglect the soft momenta in the numerator in front of the hard one and in this example we can discuss the two opposite situation in which either p 1 or k are soft. Considering Ward identities the first term in (B.2) vanishes in the case of transverse gauge bosons, while for longitudinal degrees of freedom we refer the interested reader to [27] for an analysis of the electroweak symmetry breaking effects; therefore for our purpose we have: where: • the first term reconstructs the hard scattering amplitude iM ρ ε * ρ (p 1 ) and it survives when k is soft; • the second term reconstructs the hard scattering amplitude iM ρ ε * ρ (k) and it survives when p 1 is soft.
Squaring the amplitude (B.3), we sum over polarizations using the axial gauge [28]: and as a result the eikonal limit leads to the factorization of the process in the product of the hard cross section times an emission factor integrated over the allowed phase space of the soft particle: The two integrals can be obtained through the exchange p 1 ↔ k.

B.2 The Sudakov parametrization
We consider now the explicit evaluation of the eikonal integral in (B.5), and in order to perform this calculation we choose a convenient parametrization of the external momenta; fixing two basis vectors: the Sudakov parametrization consists in the following decomposition for the soft momentum p 1 : where, without losing generality, we have taken the spatial component of k ⊥ along the x direction. In order to highlight in a simple way the logarithmical behavior of the eikonal integral we choose to work considering as first step the massless case and approximating the two hard momenta k P and p 2 P , such that The eikonal integral takes the form: and logarithmical singularities clearly arise in two opposite kinematical regions: • x x: the soft gauge boson p 1 is emitted along the k direction; integrating over x using x = k 2 t /sx, the condition x x becomes an upper bound for the transverse momentum k 2 t sx 2 and therefore in terms of fragmentation functions in the p 1 soft limit x → 0 we obtain: • x x: integrating over x using the relation x = k 2 t /sx, the eikonal integral gives exactly the same logarithmical result previously discussed but in an opposite kinematical configuration since the soft gauge boson p 1 is emitted now along p 2 direction.
In order to clarify the consequences of the symmetry p 1 ↔ k between the two eikonal integrals, we can now discuss in more details the Sudakov parametrization in the region x x. First we generalize eq. (B.8) writing for k and p 2 : k = zP + zP + k ⊥ = (E(z + z), k t , 0, E(z − z)) , p 2 = yP = (yE, 0, 0, −yE), (B.12) and using the on-shell conditions in order to eliminate z, x writing z = k 2 t /sz, x = k 2 t /sx. The conservation of energy and spatial momentum gives the following relations between the kinematical variables x, y, z: and we can therefore generalize the result for k soft just considering the substitution p 1 → k =⇒ x → z = 1 − x, and as a consequence the kinematical end-points for the x variable in the Sudakov parametrization are: Comparing eq. (B.11) with: 15) where P (x, k 2 t ) is the usual unintegrated splitting function, we obtain its leading behavior in correspondence of the two kinematical limit in the Sudakov parametrization: (B.17) with:

B.2.1 Parton masses and the lower limit of integration
The upper bound on the integration over k 2 t is dictated by the kinematical proprieties of the collinear emission, and can be studied even working in the massless case; parton masses, contrarily, affect in a relevant way the lower bound of integration. To discuss this point, we need first to generalize the Sudakov parametrization in eq. (B.12): it's straightforward to verify that one can take into account of the on-shell conditions k 2 = m 2 k and p 2 1 = m 2 1 just redefining k 2 t → k 2 t + m 2 k for k and k 2 t → k 2 t + m 2 1 for p 1 . Then from the propagator of the collinear emission p → p 1 + k we have: .

(B.19)
Depending on which particles are massive, three different situations arise.
1. Only the emitted vector is massive. This happens e.g. in EW interactions of the W, Z vectors. Assuming m 2 k = M 2 W , m 2 p = m 2 1 = 0, eq. (B.19) becomes: where the latter passage holds in logarithmical accuracy. Integrating over k 2 t we have for the x → 1 singularity of the EW splitting function P F →F : (B.21) and therefore the lower x-dependence don't affect the soft limit x → 1.
2. Two particles are massive. This happens e.g. in the electromagnetic coupling of the W : m 2 p = m 2 k = M 2 W , m 2 1 = 0; eq. (B.19) becomes: Integrating over k 2 t we have for the x → 0, 1 singularities of the splitting function P V →γ : 23) and the lower x-dependence affects the x → 0 singularity for the soft photon.
3. All three particles are massive. This happens in the massive three gauge boson vertex, m 2 p = m 2 k = m 2 1 = M 2 W ; eq. (B.19) becomes: Integrating over k 2 t we have for the x → 0, 1 singularities of the splitting function P V →V : (B.25) and therefore the lower x-dependence don't affect the soft limits x → 0, 1.
• in the soft limit x → 0 from (B.29) it follows that z ≈ y ≈ 1, • since x 2 E 2 ≥ k 2 t and 0 ≤ x ≤ 1 it is possible to approximate xE 2 + k 2 t ≈ xE 2 .
As a consequence the scalar products in (B.31,B.32,B.33) are simplified: (B.37) and the eikonal integral, considering (B.30) and introducing the mass M W as a physical cutoff for the k 2 t → 0 singularity, reduces to: that shows the behavior: s. The symmetry propriety of the eikonal integral allow us to generalize this result for k soft just considering the substitution x → z ≈ 1 − x and therefore the kinematical end-points on the x variable for the exact parametrization are: As a conclusion we obtain the leading behavior of integrated splitting functions in correspondence of the two kinematical limit in the exact parametrization: with L(x) given in eq. (18).

B.4 The Collinear Approximation
The eikonal approximation allows to highlight the singular behavior of the improved splitting functions in the soft regions x → 0, 1, as shown e.g. in Eqs. (B.41,B.42). In order to extract the entire structure of the splitting functions and to show the factorization proprieties of our model independent approach, we need to go one step further, introducing the CA; we discuss now its main features, having in mind an illustrative explicit example. Following [17] we add to the SM Lagrangian a vector boson Z with mass M M W belonging to an extra U(1) gauge symmetry and singlet under SU(3) C ⊗ SU(2) L ⊗ U(1) Y .
In order to simplify our discussion, let us suppose that the Z couples only with left electron and neutrino At tree level we have two possible leptonic decay channel Γ ee ≡ Γ 2 (Z → e + L e − L ) and Γ νν ≡ Γ 2 (Z → ν L ν L ), related through total isospin conservation Γ νν = Γ ee ≡ Γ B . Considering the process Z → ν L ν L , labeling with P = (E, 0, 0, E) and P = (E, 0, 0, −E) the two back-to-back momenta of the two neutrinos (with E = M/2), and indicating with Γ 2 (Z → ν L ν L ) the corresponding two-body decay width for the amplitude squared we have at Born level: We calculate now the effect of adding one weak gauge boson emission, focusing on the three-body decay width Γ 3 related to the process Z (Q) → ν L (p 1 )ν L (p 2 )Z T (k) and using the CA.
The key point of this approximation is the following: in the high energy regime M M W the leading contributions to the three-body decay are produced by the region of phase space where the emitted boson is collinear either to the final fermion or to the final antifermion, and in this region the three-body decay width is factorized with respect to the two-body one. Introducing: it's possible to show in a simple way how the CA works both considering the factorization of the amplitude squared and of the phase space related to the final state.
Considering for definiteness the case in which the weak gauge boson k is emitted along p 1 direction, we depict in Fig. 6  Concerning the factorization of the amplitude squared we have, working for simplicity in the massless limit: where Υ ≡ if L γ µ µ (Q); therefore when the gauge boson is emitted along the p 1 direction p 1 · k → 0, and the first diagrams diverges while the second one is finite; as a consequence squaring the amplitude M 3 = M A + M B and using Feynman gauge divergences appear in the square |M A | 2 and in the interference term M A M * B + M * A M B . Using the Sudakov parametrization B.2 for the external momenta it's straightforward to verify that: now we use the only technical point of CA; writing: we simply observe that |M A | 2 diverges as k 2 t → 0; in the trace we can therefore approximate k ≈ (1 − x)P and p 2 ≈ P : all the terms excluded by these approximations, in fact, softens the divergence and can be neglected. As a consequence we have: (B.50) and the remaining trace is exactly the one obtained in eq. (B.44) related to the process with no gauge boson emission, showing our first step towards the factorization of the amplitude squared. Considering the interference term M A M * B + M * A M B it's possible to use the same trick writing: and as a result we obtain the factorization of the amplitude squared for the three-body decay in the CA respect to the two-body one: Referring to the eikonal approximation in eq. (B.3), we see that using CA it's possible to factorize the amplitude only respect to its squared. We can apply the CA also considering the three-body phase space of Γ 3 in eq. (B.45). Following [17] we therefore obtain a complete factorization: the three-body decay width can be expressed as the product of the two-body one times a collinear factor: Integrating over the final phase space we finally find: As commonly done in parton models, we interpret the factor multiplying the two-body decay width as the parton distribution for finding a "neutrino parton" in the neutrino: (B.55) where up to this point this expression take into account just the effects related to the real gauge boson emission. Adding a tree level delta term to account the process without electroweak emission and introducing virtual corrections at the same perturbative order of the real ones we obtain: where: A complete set of integrated splitting functions is collected in Table 1. eq. (B.56) represents a concrete one-loop example, obtained through a direct calculation, of the logarithmical structure of the electroweak fragmentation function D I→J . At this point one might be skeptical about the effective validity of the CA. In order to remove all doubt, we compare in B.5 our improved CA with the full result of a complete three-body calculation done in the context of the Minimal Dark Matter model [25], finding an excellent agreement. The one-loop fragmentation functions for the entire SM are discussed into Appendix C, in the more general context of the electroweak evolution equations [18]. At this point we are ready to evaluate in this example the energy spectrum of stable SM particles produced by the Z decay and for definiteness we consider the case of the neutrino. We get: where Γ tot = Γ ee + Γ νν = 2Γ B . From eq. (B.55) we have: and, in a similar way: with D e − L →ν L as in (C.2); finally we obtain the neutrino spectrum at perturbative order α 2 : (B.61)

B.5 Full computation in the Minimal Dark Matter model
In order to validate the eikonal and collinear approximations, we compare its results with a full computation performed in a specific predictive model. We consider the weak corrections to DM annihilations into W + W − as predicted by "Minimal Dark Matter" models, where DM only has electroweak interactions [25]. This generic situation is realized in the region of the MSSM parameter space where DM could be the neutral component of the fermionic wino weak triplet with a value of the mass M ∼ 3 TeV dictated by the cosmological DM relic density [25,29].
In the W spectra the first term is the leading order annihilation, the second one arises from 3-body processes with an additional Z, and the third term from processes with an additional γ. Eq. (B.62a) agrees with the result computed for fermionic DM in [30]. 8 We now compare the full result with its collinear approximation, where the same quantities in Eqs. (B.62a,B.62b) are described through electroweak fragmentation functions. In the previous 7 Incidentally we find that at this order the result is the same for both scalar and fermionic DM. Notice that this Taylor expansion in is only valid at x , and consequently, like the Sudakov approximation, does not correctly describe the kinematical boundaries. 8 Note that the diagram there called "QED Internal Bremsstrahlung" in our language is ordinary EW bremsstrahlung from the initial DM and it does not give any log-enhanced effect because DM is non-relativistic.
We also agree with [30] Figure 7: Left, Z from W radiation: Comparison between our full result in the Minimal Dark Matter model (continuous yellow line), with its limit for ≡ M W /M → 0 (blue dot dashed) and with our improved eikonal approximation (red dotted for the Sudakov parametrization and green dashed for the exact one). We show also the comparison with the naïve standard partonic approximation (black continuous line). Right, γ from W radiation: comparison between our full result (continuous red/blue line) with our improved splitting approximation in the exact parametrization (red/blue dashed) and the standard partonic one (red/blue dotted).
sections we derived the results treating the phase space either within the Sudakov parametrization [Eqs. (B.16,B.17)] or exactly [Eqs. (B.41,B.42)]. We therefore compare the full γ spectrum in eq. (B.62a) with: (B.63) and the full Z spectrum in eq. (B.62b) with: where L(x) is given in eq. (B.18) for the Sudakov parametrization and in eq. (18) for the exact phase space. Results are shown in Fig. 7, where we depict also the curve corresponding to the naïve partonic approximation, where the upper bound on µ 2 is chosen to be the typical hard scale of the problem s, obtaining a universal logarithmical factor ln s/M 2 W that multiplies the splitting function. At small values of x the Sudakov parametrization (red dotted line in the left panel) shows a bad behavior compared with the full calculation and don't respect the correct kinematical end-points: This is because in the Sudakov parametrization, considering the splitting process i → f + f , the variable x don't correspond exactly with the fraction of energy carried by the final particle f respect to its initial value. On the contrary the exact phase space (green dashed line) gives of course the correct kinematical boundary of the splitting process (B.40) and shows a correct agreement with the full calculation.

C One loop Electroweak Fragmentation Functions
In the following we collect one loop EW fragmentation functions obtained solving the EW evolution equations in [18] for the entire Standard Model particle spectrum. EW evolution equations have been already constructed in [18], exploiting the SU(2) L symmetry, and classifying the states looking to their total isospin quantum numbers; consequently one have to apply a projection technique (explained in details in [18]) in order to convert them to their QCD-like formulation, i.e. labeling splittings with particle names: D i→j is the single leg fragmentation functions that encode the probability for a single initial particle i to become a final particle j.
Since DM gives particle-antiparticle pairs, to reduce the combinatorics, we combine them into pairs of primary back-to-back SM particles I instead of a single particle i. A simple formula allows us to switch from the 'single leg' to the 'double leg' convention: where c = 1/2 for complex final particles (such as the W or the ν) while c = 1/4 for real ones (such as the Z or the γ).

C.1 Splitting of fermions
We start considering DM that produces two back-to-back fermions, ff with center-of-mass energy √ s, and compute the resulting partonic spectrum of other SM particles A, D f →A (z), with z ≡ E A /2 √ s. As the ff pair produces A andĀ with the same energy spectrum, we always average over particle and its anti-particle both in the initial and in the final state, even for real particles such as the Z. The initial fermion can be f = {e R , e L , ν L , u L , d L , u R , d R } and is identified by its T 3 = {−1/2, 0, 1/2}, its electric charge Q, its generation number. We define the usual coupling to the Z, g f = T 3 − s 2 W Q and to the photon, α em = s 2 W α 2 . We define the top Yukawa coupling y t = m t /v with v ≈ 174 GeV and α t = y 2 t /4π. Neglecting all other fermion masses we get: where f is the SU(2) L partner of f (e.g. f = e L for f = ν L and viceversa). The virtual term in the first line means that a fraction of the initial f with x = 1 disappears into f or other particles with x < 1.
All splittings of quarks are negligible with respect to the QCD splittings (not written). For the t and b quarks there are extra splittings (to be summed to the ones listed above due to gauge interactions) due to the top quark Yukawa interaction, which also have minor effects. Again these splittings depend on the quark polarization; for simplicity, we write them for average unpolarized t, b quarks:

C.2 Splitting of Higgses
The Higgs doublet H contains the physical Higgs h, as well as the goldstone components that describe the longitudinal polarizations W L and Z L of the SM massive vectors. Splittings induced by the top quark Yukawa coupling have a significant effect; we describe them without specifying the polarizations of the t, b quarks, which make a negligible difference.
For the physical Higgs h we have: The same expressions hold for an initial Z L . For the longitudinal W L we have: P vir,Yuk S→F , (C.5a) = δ(1 − x) 1 + α em 2π P vir S→S + α em 2π P S→S (x) + · · · , (C.5b) Eq. (C.5b) shows the QED corrections only, in case one needs to subtract them.

C.3 Splitting of vectors
For the transverse W ± we find: d, u; µ, ν µ , s, c; τ, ν τ , b, t}, (C.6d) and for the transverse Z we find: where we defined N c = 1 (3) if f is a lepton (a quark) doublet; N gen = 3, and: P vir SU(2) = 1 2 P vir V →V + 1 4 P vir V →S + N gen P vir V →F , (C.8a) The γ contributes to P W T →W T and can be excluded by dropping 1 = c 2 W + s 2 W → c 2 W in front of P V →V , both real and virtual. (Notice that PYTHIA does not include QED radiation from W ± ).
Finally for the photon we have: