PPPC 4 DM$\nu$: A Poor Particle Physicist Cookbook for Neutrinos from DM annihilations in the Sun

We provide ingredients and recipes for computing neutrino signals of TeV-scale Dark Matter annihilations in the Sun. For each annihilation channel and DM mass we present the energy spectra of neutrinos at production, including: state-of-the-art energy losses of primary particles in solar matter, secondary neutrinos, electroweak radiation. We then present the spectra after propagation to the Earth, including (vacuum and matter) flavor oscillations and interactions in solar matter. We also provide a numerical computation of the capture rate of DM particles in the Sun. These results are available in numerical form.


Introduction
The Dark Matter (DM) particles that constitute the halo of the Milky Way have a small but finite probability of scattering with a nucleus of a massive celestial body (e.g. the Sun or the Earth) if their orbit passes through it. If their velocity after the scattering is smaller than the escape velocity from that body, they become gravitationally bound and start orbiting around the body. Upon additional scatterings, they eventually sink into the center of the body and accumulate, building up a local DM overdensity concentrated in a relatively small volume. There they annihilate into Standard Model particles, giving origin to fluxes of energetic neutrinos [1]. These neutrinos are the only species that can emerge, after experiencing oscillations and interactions in the dense matter of the astrophysical body. The detection of high-energy neutrinos from the center of the Sun or the Earth, on top of the much lower energy neutrino flux due to nuclear fusion or radioactive processes, would arguably constitute one of the best proverbial smoking guns for DM, as there are no known astrophysical processes able to mimic it (except possibly for the flux of neutrinos produced in the atmosphere of the Sun, which however are expected to have a different spectral shape [2]).
While the above idea is rather old, recently it has gained more interest since large neutrino telescopes (ANTARES, ICECUBE and its planned extension PINGU, SUPERKAMIOKANDE and its proposed upgrade HYPERKAMIOKANDE, the future KM3NET) are reaching the sensitivity necessary to probe one of the most interesting portions of the parameter space, in competition with other DM search strategies.
Along the years there have been several computations of the ingredients necessary to predict the signal in a detector, in particular the DM capture rate [3] and the spectra of neutrinos [4,5,6,7,8]. In [5] a subset of the authors of this paper had computed the neutrino spectra: their production in solar matter and their propagation subject to interactions and flavor oscillations. Similar results, but in an event-based approach, had been subsequently obtained by WIMPSIM [6].
In this paper we improve the computation of the neutrino energy spectra with respect to previous computations and in particular to [5] in various ways: • we consider a full list of two-body SM annihilation channels; • we implement a better description of the energy losses of the primary particles in solar matter, including secondary neutrinos produced by scattered solar matter; • we extend to low energy the neutrino spectra in order to include the large flux of neutrinos generated by stopped µ, π, K (an effect pointed out in [7,8]); • we take into account EW bremsstrahlung; • we update the neutrino oscillation parameters adding the recently discovered θ 13 ≈ 8.8 • and the two possible mass hierarchies of neutrinos (normal or inverted).
Let us now review the main steps of the computation and present the results we provide. The final energy spectrum of the neutrino flux at the detector location is written as where d is the Sun-Earth distance, Γ ann is the total DM annihilation rate in the Sun, dN ν /dE ν is the energy spectrum of ν e,µ,τ andν e,µ,τ produced per DM annihilation after taking into account all effects. The computation proceeds as follows: 1. The total rate of DM annihilations Γ ann in the Sun is computed in section 2 in terms of the spin-independent and spin-dependent DM cross sections with nucleons.
2. In section 3 we compute the energy spectra of ν andν at production for DM annihilations into e + L e − L , e + R e − R , µ + L µ − L , µ + R µ − R , τ + L τ − L , τ + R τ − R , ν eνe , ν µνµ , ν τντ , qq, cc, bb, tt, γγ, gg, where q = u, d, s denotes a light quark, the subscripts L, R on leptons denote Left or Right-handed polarizations, the subscript L, T on massive vectors denote Longitudinal or T ransverse polarization. h denote the higgs boson, for which we assume a mass of 125 GeV.
Decays, hadronization, energy losses in matter and secondary neutrinos generated by matter particles scattered by DM decay products (an effect ignored in previous computations) are computed running GEANT4 [9] using the EU Baltic Grid facilities. As a check, the first three effects are also independently computed with PYTHIA [10], modified by us to include energy losses in matter.
3. As described in section 3.2, we correct the spectra at production including the dominant electroweak radiation effects (not included in MonteCarlo codes) enhanced by logarithms of M DM /M W , as described in [11]. EW corrections depend on the polarisation of the SM particles, which is why we need to specify it in eq. (2).
4. As described in section 4, we propagate the fluxes of ν andν at production, around the center of the Sun, to the Earth taking into account oscillations, matter effects, absorption and regeneration from collisions with matter using the neutrino density matrix formalism as in [5].
The final results are presented and briefly discussed in section 5. In section 6 we conclude.
The main numerical outputs of the computation are given on the PPPC 4 DM ID website ('DMν' section).

The DM annihilation rate
In this section we compute the DM solar annihilation rate, which sets the overall normalization of the expected neutrino flux. We essentially revisit the calculations in [3]. For definiteness, we focus on the case of the Sun, although we mention in passing the modifications needed for the Earth.
The number N of DM particles accumulated inside the Sun varies with time under the action of different competing processes: 1) DM particles are captured from the halo (N is thus increased by one unit) via the multiple scattering processes discussed above, 2) DM can pair annihilate (hence N decreases by two units) and 3) DM particles can be ejected (N decreases by one unit) by a hard scattering on a hot nucleus of the interior of the body, i.e. the inverse process with respect to capture. In formulae: where Γ capt is the capture rate, Γ ann is the DM annihilation rate and Γ evap is the DM evaporation (or ejection) rate. The evaporation process is important only for DM lighter than a few GeV 1 , so that we will neglect it for all practical purposes in the following. 1 This can be easily understood [3]: if one requires that the speed of a DM particle thermalized with the interior of the Sun v DM ∼ 2T /M DM be smaller than the escape velocity from the center of the Sun v esc 1387 km/s, one obtains the condition M DM 0.15 GeV. A more realistic computation takes into account that DM particles actually follow a (Maxwellian) velocity distribution and compares the time it takes to deplete the high-velocity tail (above the escape velocity) of such distribution with the age of the Sun, finding M DM 5 GeV. Yet more refined calculations find typically M DM 4 GeV [12].
The annihilation rate is proportional to N 2 : two DM particles annihilate (hence the square). It is given by where σv is the usual annihilation cross section averaged over the initial state 2 and n( x) is the number density of DM particles at position x inside the Sun, such that the total number of DM particles is N = d 3 x n( x). After capture, subsequent scatterings thermalize the DM particles to the solar temperature T , such that their density n( x) acquires the spherically symmetrical Boltzmann form where n 0 is the central DM number density and φ(r) = r 0 dr G N M (r)/r 2 is the Newtonian gravitational potential inside the Sun, written in terms of the solar mass M (r) enclosed within a sphere of radius r. Taking for simplicity the matter density in this volume to be constant and equal to the central density ρ , all the integrals can be explicitly evaluated . One finds that DM particles are concentrated around the center of Sun, Within this approximation, one obtains from eq. (4) Here ρ = 151 g/cm 3 and T = 15.5 10 6 K are the density and the temperature of matter around the center of the Sun. The same expression would hold for other astrophysical bodies, adapting these two quantities.
Neglecting Γ evap and solving eq. (3) with respect to time one finds where τ = 1/ Γ capt C ann is a time-scale set by the competing processes of capture and annihilation. At late times t τ one can approximate tanh(t/τ ) = 1. In the case of the Sun, the age of the body (∼4.5 Gyr) and the typical values of the parameters in τ indeed satisfy this condition (in the case of the Earth this is not generally the case). Therefore one attains the last equality of eq. (8). Physically, this means that the fast (compared to the age of the Sun) processes of capture and annihilation come to an equilibrium: any additional captured particle thermalizes and eventually is annihilated away.

The DM capture rate
As a consequence of the attained equilibrium, the computation of the annihilation rate Γ ann crucially depends on the computation of the capture rate Γ capt , which acts as a bottleneck. The computation of the latter proceeds by summing the contribution to capture of all the individual shells of matter located at position r within the massive body. The result is [3] Its derivation is lengthy: we will only sketch it in the following by illustrating the individual pieces of this equation.
• ρ DM /M DM corresponds just to the local number density of DM particles at the location of the capturing body.
• The sum runs over all kinds of nuclei i with mass m i and number density n i (r), to be integrated over the volume of the Sun. We take the standard solar elemental abundances from [13]. The factor σ i is the low-energy DM cross section on nucleus i, assumed to be isotropic. In terms of the standard Spin Independent and Spin Dependent cross sections normalized on a single nucleon at zero momentum transfer 3 (denoted σ SI and σ SD ) we take σ H = σ SI + σ SD for Hydrogen and σ He 4 4 σ SI (m p + M DM ) 2 /(m He + M DM ) 2 for the spin-less Helium. Analogously, for heavier nuclei of mass number A one can take . This latter SD formula should be complemented by weighting factors due to the spin content in the different nuclei. In practice, only Hydrogen is relevant for our case, so that a more careful evaluation is not necessary.
• f (v) is the angular average of the DM velocity distribution with respect to the solar rest frame neglecting the gravitational attraction of the Sun. This quantity is connected to the f (v) with respect to the galactic rest frame as in terms of the solar velocity v ≈ 233 km/s (here c is the cosine of the angle between v and v ). Observations and numerical simulations suggest that f (v) can be parameterized as [14] f with normalisation ∞ 0 dv 4πv 2 f (v) = 1 and parameters 220 km/s < v 0 < 270 km/s, 450 km/s < v esc < 650 km/s, 1.5 < k < 3.5. (12) Here v esc is the escape velocity from the Galaxy (not to be confused with v esc , the escape velocity from the Sun) at the location of the solar system. For k = 0, f (v) reduces to a Maxwell-Boltzmann distribution with a sharp cut-off at v < v esc : . For k > 0 the cut-off gets smoothed. The resulting f (v) are normalised such that ∞ 0 dv 4πv 2 f (v) = 1. The actual choice of the DM velocity distribution turns out to have a rather limited impact, see e.g. [15] for a recent analysis.
• The gravity of the Sun is taken into account by v esc (r), which is the escape velocity (from the Sun) at radius r, such that v 2 + v 2 esc is the squared velocity that a DM particle acquires when arriving at r from a very distant point. One here neglects the effect of other bodies in the solar system, which presumably is a good approximation [16].
• The probability that a DM particle, with velocity v far from the Sun, gets captured when scattering on a nucleus located where the escape velocity is v esc is given in first approximation by where are the maximal and minimal fractional energy loss ∆E/E that a particle can suffer in the scattering process, provided that it is captured. I.e., the probability is computed as the ratio of the size of the interval in energy losses leading to capture (∆E min < ∆E < ∆E max ) relative to the whole possible interval (0 < ∆E < ∆E max ), assuming a flat distribution of the scattering cross section in energy. In general, however, one needs to introduce the form factors F i (∆E) that take into account the nuclear response as function of the momentum transfer. Explicitly, is the effective radius of a nucleus with mass number A i and mass m i ). The numerator of the 'ratio of sizes' becomes then an integral of the form factor over the energy loss ∆E: Eq. (13) and (14) mean that the fraction of scatterings that lead to capture is largest for nuclei with mass m i comparable to the DM mass M DM (∆ max is maximized) and for DM particles that are slow (small v) and in the central regions of the body (large v esc ).
Thereby the capture rate is proportional to 1/M 2 DM and can be approximated as where  We here parameterised f (0) = 1/π 1/2 v eff 0 3 , such that the parameter v eff 0 coincides with the parameter v 0 for a Maxwellian distribution with no cut-off and neglecting solar motion.
We also provide a numerical function which gives the annihilation rate Γ ann = Γ capt /2 in terms of M DM , the standard DM cross sections on a single nucleon σ SD and σ SI (as defined above) and the parameters of the velocity distribution v 0 , v esc , k.

The neutrino energy spectra at production
In this section we discuss the computation of the neutrino fluxes emerging from the DM annihilation process. Neutrinos are produced at several stages of the hadronic and leptonic cascades originating from the primary particles produced by annihilations, and these cascades develop within the dense matter of the astronomical body (the Sun). Such an environment has important consequences in determining the spectra. For example, some byproducts of the cascade can be absorbed by the surrounding matter; some others will be just 'slowed down'; others are negligibly affected. It is therefore important to model in the best possible way all these effects. We adopted and compared two strategies: 1) As described in section 3.1.1, we modified the public MonteCarlo code PYTHIA to include the effects of cascading energy losses with matter and absorption. The disadvantage is that in this way we cannot include the neutrinos emitted by matter particles scattered by products of DM annihilation.
2) As described in section 3.1.2, we run the GEANT4 code, dedicated to particle interactions with matter. The disadvantage is that GEANT4 is much more time-consuming than PYTHIA: to reach a reasonable statistics we employ the computational resources of the Baltic grid.
In section 3.1.3 we compare the two results, interpret them and decide to adopt the GEANT4 result. In section 3.2 we describe how we subsequently add electroweak bremsstrahlung corrections, which are not included in any Monte Carlo code.

PYTHIA
We here describe how we modified the PYTHIA public code to include the effect of cascading in matter in an event-by-event basis. We improved with respect to our previous computation [5] where we had only included averaged energy losses in matter. Moreover, we here use PYTHIA version 8.176 (version 6.2 was used in [5]). The process of particles with mass m and initial energy E 0 m decaying with life-time τ , while at the same time losing energy due to interactions with the surrounding matter, is described by a differential equation for their energy E and their number n: that is: the energy loss of particles in matter can by approximated as a constant term α plus a term β proportional to their energy, while their number n is just governed by τ (the factor m/E takes into account time dilatation). Therefore is solved by We thereby instruct PYTHIA to perform the decay of particles by first replacing their initial energy E 0 with a lower energy E randomly chosen according to the above distribution produced by energy losses, eq. (21), and then letting them decay. This equation takes different specific forms according to the particle considered, as we discuss now.

Dominant β: stopping of hadrons
The stopping of hadrons is well approximated by α = 0. In such a case eq. (21) simplifies to where E cr = m/βτ is a critical energy below which energy losses are negligible. The distribution of energies at decay is exponential in In the center of the Sun one has E cr ≈ 250 GeV for charmed hadrons and E cr ≈ 470 GeV for b-hadrons [4]. This improves over our previous computation [5] where we used a constant average energy at decay, given by Neutrons and negatively charged pions need however a peculiar treatment: once stopped, instead of decaying, they are absorbed by matter. Neutrons are absorbed by protons, forming deuteron. The π − are quickly Coulomb-captured by heavy nuclei in the dense and hot solar core environment, forming pionium which typicallys decay into a neutron and photons [17]. Since no neutrinos are produced in final states, to implement this phenomenology we simply modified our PYTHIA code in such a way that neutrons and pions are removed and not decayed.

Dominant α: stopping of charged leptons
The energy loss of a charged lepton can often be approximated as an energy-independent term (i.e. taking β = 0), which in the center of the Sun equals to The general solution of eq. (21) simplifies then to The exponent m/τ α equals 1/176000 for muons (which therefore lose almost all their energy before decaying) and equals 660 for taus (which therefore decay before losing a significant fraction of their initial energy). The average energy at decay employed in our previous computation is

GEANT
GEANT4 is a toolkit for the simulation of the passage of particles through matter [9], a very suitable tool for modeling cascades in the environment of the solar core. We follow the approach of [7,8]. Hadronization of the quarks and gluons produced by DM annihilations is performed by PYTHIA; the stable and metastable hadronization products (e, ν e,µ,τ , µ, K L0 , π, K, n, p, N D,He ) are injected into GEANT4 which adds the effect of particle/matter interactions. We model the matter around the solar core following the Solar Standard Model as discussed in Sec. 2, but with a simplified chemical composition (we just include 74.7 % by mass of H and 25.3 % of He). In our GEANT4 code, we consider a sphere with radius of 1 km, that is big enough for our purposes: we verified that only neutrinos from the cascades can reach its surface, while secondary products are contained. The passage of particles through matter is simulated using the QGSP BERT physics list, see [9] for further details. Notice that neutrinos have no matter interactions nor oscillations in the QGSP BERT physics model of GEANT4 (we will indeed discuss in detail how to deal with these effects in the subsequent neutrino propagation in matter in Sec. 4). It also shows (when available) the spectra previously computed in [5]. In very general terms, moving from low to high x = E/M DM , the new spectra are characterized by some very pronounced low energy humps or spikes (missing in [5]), an intermediate energy smooth shoulder and, for some channels, a high energy peak. These features are easily understood.

Results and comparisons
• The high energy peak occurs when DM annihilate into particles that directly decay into neutrinos. This is visible in fig. 2 as a slanted mountain top in the W + W − channel (panel (e)). A similar peak is of course present for the ZZ channel (not shown). For large DM masses, i.e. for large boosts of the primary W and Z, this feature smears into a smooth spectrum (see e.g. panel (a) of fig. 3).
• The smooth component of the spectrum arises from the neutrinos produced in the cascading event by primary and secondary particles (hadrons and leptons), that lose energy and rapidly decay.
• The low energy humps and spikes essentially arise from relatively long-lived particles that have been stopped in solar matter and then decay, essentially at rest. More precisely, recalling that n, π − , µ − and K − are mostly absorbed or captured by matter, the low energy neutrino peaks arise from the following processes: π + → µ + ν µ decays, which produce a monochromatic line in the ν µ spectra, at E ν = 29.8 MeV (visible in panel (d) of fig. 2). For numerical reasons we artificially broaden it.
µ + →ν µ ν e e + decays, which contribute to the ν e ,ν µ energy spectra producing the typical three body decay bump with an end-point at m µ /2 ≈ 53 MeV and peaked at panel (a)).
µ − → ν µνe e − decays, which similarly contribute to theν e , ν µ energy spectra, although the resulting bump is about two orders of magnitude less intense than the bump in ν e ,ν µ . Indeed π − are absorbed by matter before decaying into µ − .
-K − get absorbed, synthesised by nuclei into a Λ which decays into nucleons and pions. Their rare free decays negligibly affect the neutrino spectra.
-Decays of K 0 L into neutrinos are blocked by matter effects that break the quantum coherence between K 0 S and K 0 L , such that K 0 S are continuously regenerated and fastly decay hadronically.
While these general features are present both in the PYTHIA and the GEANT spectra, the results of the two MonteCarlos differ in the details.
• For hadronic channels (e.g. bb), PYTHIA systematically produces a softer high-energy spectrum than GEANT, as clearly visible in panels (a), (b) and (c) (and (d) too). The difference is small for small DM masses and increases for large DM masses. For the leptonic channels, on the other hand, the two results agree (see e.g. panel (f)). This can probably be ascribed to a difference in the way PYTHIA treats the energy losses of the energetic hadrons in the cascade.
• PYTHIA systematically underestimates the low energy humps and spikes, as it is clearly visible for instance in panel (f). 4 This is expected since, as we anticipated, PYTHIA includes only the neutrinos produced by the decays of pions, kaons and muons in the DM-originated shower. GEANT, instead, follows the fate of all the particles produced in the scattered matter, including therefore additional pions, kaons and muons that then decay into neutrinos.
Based on the considerations above, we adopt the GEANT spectra.

Adding electroweak bremsstrahlung
In the previous subsection we have presented the results of the detailed MonteCarlo computation of neutrino spectra from DM DM annihilations into pairs of SM particles. There are however higher order effects that can be relevant and are not included in such computations. In general, higher order corrections cannot be computed without having a full DM model: one needs to know the particles mediating the annihilation process in order to include all possible relevant diagrams. But some dominant higher order corrections can be computed in a modelindependent way: those enhanced by logarithms of ratios of particle masses, which describe bremsstrahlung. In terms of Feynman diagrams, such process corresponds to attaching soft or collinear particles to the SM final state particles. While the bremsstrahlung emission due to electromagnetic and strong interactions is automatically performed by MonteCarlo codes, the bremsstrahlung emission due to electroweak interactions is not included and is equally significant, if M DM M W .
We therefore include electroweak bremsstrahlung at leading order in the electroweak couplings by 'post-processing' the output of the MonteCarlo as described in [11]. At large DM masses, such bremsstrahlung corrections are enhanced by ln(M DM /M W ) logarithms, which become large for M DM M W . Such enhanced terms are model-independent: in our code we turn them on abruptly when M DM > ∼ M W . In a full DM model, these effects would actually appear in a smooth model-dependent way when increasing the DM mass. We instead neglect the finite non-logarithmic terms, that cannot be computed in a model-independent way.
In practice, we proceed as follows. In the previous section we computed dN MC J→ν /dx: the MonteCarlo spectra in x = E/M DM of ν produced by DM annihilations into a generic two-body state J. To include EW bremsstrahlung, we convolute such spectra with a set of electroweak splitting functions D EW I→J (z) (probability that I radiates a particle J with energy reduced by a factor z) as follows: The sum is over all EW splittings. The rationale behind this procedure is that splittings are kinematically different from decays and happen before decays, when particles have a large virtuality. The splitting functions are predicted by the SM and listed in [11]. For example, left-handed electrons can radiate neutrinos via the EW splitting e L → ν e W T described by the function where = ln 4M 2 DM /M 2 Z and From a phenomenological point of view, such effects can be relevant. Annihilation channels that would produce a low yield of neutrinos (such as DM DM → e + e − or DM DM → µ + µ − where muons are severely slowed down or absorbed by matter) are significantly affected by EW bremsstrahlung because they receive contributions from channels with large yield of neutrinos (such as DM DM → W + W − ). This case is shown in panel (b) of fig. 3. The other panels in fig. 3 show other cases in which the impact of EW effects is sizable. For a channel like W + W − (panel (a)), the radiation leads to a neutrino flux enhanced by about an order of magnitude at low energy. An 'extreme' case is reproduced in panel (c): for annihilations directly into tau neutrinos and looking at the flux of ν τ themselves, EW corrections add to the naïve monochromatic spectrum a broad low energy shoulder, as a consequence of the decay of the EW gauge bosons emitted by the primary ν τντ pair. Fig. 4 presents, for reference, an example of our final results for the neutrino spectra at production. We plot the spectra of ν (first column) andν (second column) for all the channels that we consider for a sample DM mass of 1 TeV. The spectra are normalized per one annihilation of two DM particles. The considered range of x = E/M DM covers from 10 −8 to 1. In the third column we zoom on the high energy part, relevant for neutrino detectors such as ICECUBE.

Spectra at production: results
All these spectra are provided in numerical form.

Neutrino propagation
Propagation of neutrinos from the interior of the Sun to the Earth is affected by flavor oscillations and (at energies above tens of GeV) by Neutral Current (NC) and Charged Current  (CC) interactions with solar matter, which give rise to absorption and (when a τ lepton is produced) to regeneration of neutrinos with lower energies. In the following we discuss each of these effects and summarize the formalism we employ.

Formalism
Coherent flavor oscillations and coherence-breaking interactions with matter simultaneously affect neutrino propagation. Following [4], the appropriate formalism that marries in a quantum-mechanically consistent way these two aspects, consists in studying the spatial evolution of the 3 × 3 matrix of densities of neutrinos, ρ(E ν ), and of anti-neutrinos,ρ(E ν ). The diagonal entries of the density matrix represent the population of the corresponding flavors, whereas the off-diagonal entries quantify the quantum superposition of flavors. 5 The matrices ρ(E ν ) andρ(E ν ) satisfy a coupled system of integro-differential equations in the distance r from the center of the Sun: with an analogous equation forρ.
• The second term in eq. (30) describes the absorption and re-emission due to NC scatterings ( ν ) N ↔ ( ν ) N * (where N is any nucleon in the Sun and with N * we denote its possible excited state after the collision), which remove a neutrino from the flux and re-inject it with a lower energy. So they contribute to the evolution equation as: where • The third term in eq. (30) describes Charged Current (CC) scatterings ( ν ) N → ± X of an initial neutrino ν with energy E ν , which remove the ν from the flux and produce a charged lepton and scattered hadrons X. They decay back into neutrinos ν and anti-neutrinosν with lower energy E ν : their energy distributions are described by the function f → (E ν , E ν ). When the initial neutrino is ν τ (ν τ ), the produced τ − (τ + ) decays promptly before losing energy, giving rise to energetic ν τ ,ν e ,ν µ (ν τ , ν e , ν µ ): this is the tau regeneration phenomenon [18]. When instead the initial neutrino is a ν e or ν µ we assume that the produced e, µ is totally absorbed and we neglect the corresponding low energy neutrinos. CC scatterings thereby affect the propagation of neutrinos with the term where Π is the projector on the flavor ν : e.g. Π e = diag (1, 0, 0). The matrices Γ CC ,Γ CC that describe the rates of CC interactions are given by In both NC and CC processes, we neglect low energy neutrinos that might emerge from the scattered hadrons and light leptons, i.e. the de-excitation of N * in NC and the decay of X and e, µ in CC. E.g. in particular we neglect the very low energy neutrinos with E ν ∼ m π,K,µ coming from the decay at rest of light hadrons/leptons. In order to include them, one should implement neutrino/matter interactions in dedicated codes such as GEANT, which currently do not include them. I.e. this would be analogous to the work we performed in 3.1.2, now with neutrinos as primary particles. We postpone this to possible future improvements. We estimate that such a neglected effect would only give a small enhancement in the final flux of neutrinos at very low energy.

Transition probabilities
We numerically solve the full evolution equation (eq. 30) starting from the initial condition dictated by the spatial distribution of DM annihilations inside the Sun, proportional to n(r) 2 given by eq. (6). In practice, we numerically computed the full transition probabilities   P ± (ν (E ) → ν i (E)) from the Sun to the Earth, with = {e, µ, τ,ē,μ,τ }, i = {1, 2, 3,1,2,3} and E ≤ E in the two cases of normal (P + ) and of inverted (P − ) neutrino mass hierarchy.
The transition probabilities incorporate all propagation effects. In order to have a qualitative understanding of their behavior, however, we plot them only in the limit E = E that excludes the neutrinos produced by regeneration (at E < E ). Fig. 5 and 6 show the transition probabilities P (ν → ν ) between flavour eigenstates of neutrinos (continuous curves) and anti-neutrinos (dashed) from the Sun to the surface of the Earth (before a possible crossing of the Earth, which will be considered in the next subsection). Fig. 5 shows the probabilities computed for θ 13 = 0. In fig. 6, θ 13 is restored to its actual value θ 13 = 8.8 • and therefore we need to distinguish the hierarchy: the upper row refers to a normal spectrum of neutrinos (m 1 m 2 m 3 ), while the lower row to an inverted spectrum (m 1 m 2 ≈ m 3 ). Considering for example P (ν e → ν e ): it goes from 1 − 1 2 sin 2 2θ 12 at E MeV (averaged vacuum oscillations) to sin 2 θ 12 at larger energies (adiabatic MSW resonance for θ 12 ). If neutrinos have normal hierarchy, at E ν 100 MeV also θ 13 is enhanced by an adiabatic MSW resonance, and P (ν e → ν e ) drops. 7 Many effects happen when E ν ∼ 10 GeV: MSW resonances cease to be adiabatic, and the solar oscillation wave-length becomes comparable to the size of the Sun. These two effect cause an increase of P (ν e → ν e ) towards its vacuum oscillation value. However this increase gets stopped by neutrino absorption due to interactions with solar matter, which causes all probabilities to drop to zero in the limit of large energy. The other oscillation probabilities can be similarly understood.
Our propagation results agree, in the limit cases, with known results of the past, and in particular with [5] and the works that made use of it. We verified that the unknown neutrino CP-violating phase that enters into oscillations has negligible effects. Furthermore, such results do not depend on the Majorana or Dirac nature of neutrinos.

Earth crossing
The spectra of neutrinos and anti-neutrinos at the Earth are finally computed (in terms of mass eigenstates ν i ) by convoluting the transition probabilities with the spectra at production as dN ± The final step consists in taking into account the oscillations in the matter of the Earth. If neutrinos do not cross the Earth, the energy spectra for the neutrino flavour eigenstates are simply given by If instead neutrinos cross the Earth with zenith angle ϑ (cos ϑ = −1 corresponds to the maximal vertical crossing, and cos ϑ = 0 corresponds to the minimal horizontal crossing), the neutrino fluxes at detection are given by where the oscillation probabilities P earth are readily computed adopting the standard Earth density model. We neglect neutrino absorption within the Earth (they would be relevant only for energies above ∼ 10 TeV and neutrinos with those energy essentially do not emerge from the Sun, as discussed above). Some Earth oscillation probabilities are plotted in fig. 7 for illustration; for large and small neutrino energy they approach the limiting values |V i | 2 .

Final result
Recapping the results of the previous sections: we provide the energy spectra at detection of neutrino flavor eigenstates dN ± ν /dE ν produced by one DM annihilation in the Sun, after taking into account all effects that neutrinos experience during their journey. We provide two sets of spectra: dN + ν i /dE ν corresponding to neutrinos with normal mass hierarchy, and dN − ν i /dE ν corresponding to neutrinos with inverted mass hierarchy. Fig. 8 presents, for reference, an example of our final results for the neutrino spectra at detection, analogously to fig. 4. The spectra are, as always, normalized per one annihilation of two DM particles.
In fig. 9 we present a more detailed comparison of the effect of propagation, for a few selected masses and channels. One sees, for instance:  • The effect of flavor vacuum oscillations: for an annihilation into τ + τ − the flux of electron and muon neutrinos is greatly enhanced and the corresponding flux of tau neutrinos is depleted; for an annihilation intobb, the opposite happens since ν e,µ mostly emerge from the b channel.
• The effect of solar matter absorption: moving towards higher masses, the spectra are significantly degraded in energy; the case of the Z spectrum (peaked at production) is the most apparent. For even larger m DM all spectra approach a limit, 'bell-shaped' exponential spectrum dictated by the maximum energy to which the Sun is transparent [5].
• The effect of Earth crossing oscillations: the wiggles at around 1 to 10 GeV.
All these spectra are provided in numerical form.
The neutrino fluxes can then be converted into fluxes of detectable particles (up-going muons, through-going muons, showers). Various experiments searched for DM neutrinos from the Sun, producing bounds: SUPERKAMIOKANDE [21], ICECUBE [22], ANTARES [23] and BAKSAN [24]. These collaborations report results as bounds on DM annihilations assuming some DM annihilation channel and rates of events (such as through-going-muons) with cuts optimised assuming specific energy spectra. Unfortunately, however, they do not presently report experimental bounds on the neutrino fluxes from the Sun (this kind of analysis could be done assuming monochromatic neutrino spectra at different energies). 8 Therefore we cannot compare our improved neutrino fluxes with experimental data and we cannot at present derive improved bounds on all the DM annihilation channels.

Conclusions
In this paper we have reconsidered the computation of the high energy neutrino fluxes from the annihilation of DM particles captured in the Sun. With respect to other DM indirect detection search strategies, such a signature is particularly interesting and timely given that: i) It is unique in the sense that it cannot be mimicked by known astrophysical processes; it is therefore virtually background-free at the source, except for solar corona neutrinos [2] (the background for detection consists of atmospheric neutrinos). ii) The propagation of neutrinos from the Sun to the Earth is under control, in particular now that the neutrino oscillation parameters are (almost) all measured and with good precision. iii) The current and upcoming generation of neutrino telescopes are reaching the sensitivity necessary to probe one of the most interesting portions of the parameter space, in competition with other DM search strategies.
We have reviewed the capture of DM particles in the Sun, considering the relevant elements in solar matter and different assumptions for the halo velocity distribution.
We have computed the neutrino spectra from the annihilation into all two-body SM annihilation channels, implementing a detailed description of the energy losses of primary particles in solar matter, including secondary neutrinos produced by the scattered matter and low energy ones from the decay at rest of light hadrons/leptons. We performed this computation with the GEANT MonteCarlo, after checking against PYTHIA. We added on top of this the effect of ElectroWeak radiation, which is relevant for DM masses above some hundreds of GeV. We span the range of masses from 5 GeV to 100 TeV.
We then perform the propagation of neutrinos from the center of the Sun to the detector, taking into account (vacuum and matter) oscillations and interaction with solar matter, as well as Earth crossing. We adopt up-to-date oscillation parameters and consider normal or inverted neutrino mass hierarchy.
The main numerical outputs of the computation are given on the PPPC 4 DM ID website ('DMν' section). We provide: a numerical function for Γ ann , the neutrino spectra at production, including EW radiation effects, the neutrino spectra at detection, including all propagation effects.
These results can be readily used to derive predictions for the experimental observables.