Enhanced anti-deuteron Dark Matter signal and the implications of PAMELA

We show that the jet structure of DM annihilation or decay products enhances the anti-deuterium production rate by orders of magnitude compared to the previous computations done assuming a spherically symmetric coalescence model. In particular, in the limit of heavy DM, M>>m_p, we get a constant rather than 1/M^2 suppressed anti-deuterium production rate. Therefore, a detectable anti-deuterium signal is compatible with the lack of an excess in the anti-proton PAMELA flux. Most importantly, cosmic anti-deuterium searches become sensitive to the annihilations or decays of heavy DM, suggesting to extend the experimental anti-deuterium searches above the O(1) GeV scale.


Introduction
The cosmological DM abundance is naturally produced in thermal freeze-out if Dark Matter (DM) has weak interactions and a TeV-scale mass M , that in appropriate models can be lowered down to the weak scale, 100 GeV. This scenario can be tested searching for DM annihilation (or decay) products in cosmic rays. In view of astrophysical backgrounds, a particularly sensitive signal is an excess in cosmic-ray anti-particles: positrons, antiprotonsp and anti-deuteriumd. According to the coalescence prescription [1], ad is formed when DM produces ap and an with momentum difference below p 0 ≈ 160 MeV. The standard formula for thed spectrum, obtained under the assumption of spherical symmetry of the events, in terms of the anti-nucleon (p plusn) energy spectrum per annihilation, dN N /dT , is [2,3,4,6,7] where the kinetic energies T = E − m are T p = T n = T d /2, m p = m n = m d /2 and kd = T 2 d + 2m d Td. Eq. (1) implies ad yield suppressed by 1/M 2 for large M . This result is qualitatively wrong. Increasing M just increases the boost of the primary DM annihilation products, giving rise, due to Lorentz symmetry, to an essentially constantd production rate with energy roughly proportional to M . The reason for this fundamental discrepancy is caused by the fact that the spherical approximation misses the jet structure of the DM annihilation products.
In this letter we show that for M m p the angular proximity of the produced p,n enhances thed yield, possibly by orders of magnitude. We critically compare the standard spherical approximation results with our Monte Carlo approach tod production, presenting thed energy spectra for the various DM annihilation or decay channels into W + W − , ZZ, qq, bb, tt, hh and comment on the astrophysicald background produced mostly in cosmic ray pp and pp collisions. We propagated in the Milky Way, studying the phenomenology and the prospects for DM producedd searches at AMS-2, in the light of the PAMELAp observations. We find that thed signal below 1 GeV is strongly enhanced increasing the chances ofd detection at AMS-2 even for the standard thermal DM annihilation cross section. This result is consistent with the lack ofp/p excess in PAMELA. Due to the qualitatively different large M behavior of the production rate, our result drastically enhance thed production at high energies. Therefore the cosmic rayd flux produced in heavy DM annihilations or decays exceeds the estimated background, and AMS-2 and futured experiments become sensitive to DM if they extend their sensitivity tod above 1 GeV.
2 Spherical-cow vs Monte Carlō d is formed viapn →dγ, that has a large cross-section due to the small binding energy ofd, which therefore has a spatially extended wave-function or equivalently a strongly peaked wavefunction ψ(∆k) ≡ d |pn in momentum space. Hered has momentum k d and energy E d and ∆k is the relativistically invariant relative momentum betweenp (with momentum k p and energy E p ) andn (with momentum k n and energy E n ): where we can neglect m p − m n = 1.29 MeV. The amplitude ford production in DM annihilations, DM DM →d, can be computed as giving the 'coalescence approximation' [1]: the probability | d |pn | 2 that ap and an coalesce to form ad is approximated as a narrow step function Θ(∆k − p 0 ), that drops from unity to zero if ∆k is larger than p 0 . Here p 0 is a constant (to be extracted from data later) that can be estimated as p 0 ∼ √ m d B d ∼ 60 MeV assuming thatd production happens until the relativep,n kinetic energy is smaller than the deuteron binding energy B d = 2.2 MeV. The totald yield therefore is In the non-relativistic limit k p,n m p,n and for small p 0 the region that satisfies ∆k < p 0 at fixed k n is a sphere in k p centered on k n with radius p 0 and volume 4πp 3 0 /3. In general the sphere gets dilatated along the direction k p ≈ k n by a relativistic Lorentz factor γ p ≈ γ n ≈ γ d . Multiplying eq. (4) times 1 = d 3 k d δ( k d − k p − k n ) we finally get, in the limit of small p 0 M/γ p,n , thed momentum distribution: where k p = k n = k d /2. 1 Eq. (5) is relativistically invariant as it contains the usual relativistic phase space

The spherical approximation
Previous computations proceed assuming spherical symmetry, d 3 k = 4πk 2 dk, and uncorrelatedp,n distributions: Writing the result in terms of the adimensional x i = T i /M (so that 0 ≤ x d,p,n < 1 and x p = x n = x d /2) one gets eq. (1) i.e.
which is explicitly suppressed by 1/M 2 for large DM mass M . This is qualitatively wrong. Consider for example the DM DM → W + W − annihilation mode. Increasing M increases the boost of each W , and thereby the boost of the antideuterons from W decay, but thed number stays fixed. Neglecting QED final state radiation (FSR), for M M W one should get a constant, M -independent function for dN d /dx d . Obviously the problem is in the 'spherical cow' approximation [8]. Due to the W ± boost the events are highly non spherical and SM particles are concentrated in two back-to-back jets, enhancing the probability of havingpn pairs with small momentum difference ∆k < p 0 . A similar argument applies to DM annihilations or decays into colored particles, such as qq. Hadronization leads to QCD jets, rather than to spherical events. Thereby the spherical approximation can grossly underestimate thed production.
Going to less relevant aspects that control order one factors, the analytic spherical approximation can also over-estimate thed yield, by neglecting anti-correlations between n andp or the fact that nod is obtained if only one anti-nucleon is present per event. As an example, we consider again the W + W − mode: within the spherical approximation ad can form coalescing ap from W − with an from W + , but this process is highly suppressed because the W + and W − go back to back.

The Monte Carlo approach
In order to take into account the jet structure of the events and the correlations between thep,n momenta we compute thed spectrum by searching event-by-event for then,p pair(s) which have relativistically invariant momentum difference ∆k smaller than p 0 . We verified that the spherical uncorrelated approximation of eq. (7) is reproduced if we first merge many events, and later coalescep withn without imposing that they come from the same event.
Various experiments extracted compatible values of p 0 from data aboutd production in hadronic and e + e − collisions. Presumably these studies adopted the 'spherical cow' approximation rather than performing a Monte Carlo computation. Giving the relatively low energies involved this should not make a large difference; anyhow we here prefer to directly extract p 0 from the ALEPH data [9]: one hadronic Z decay at rest gives rise to (5.9 ± 1.9) 10 −6d in the momentum range 0.62 GeV < k d < 1.03 GeV and angular range | cos θ| < 0.95. According to our Monte Carlo computation, this translates into p 0 = 162 ± 17 MeV. Should p 0 have a value different from the p 0 = 160 MeV adopted here for both the DM signal and the astrophysical background (as computed in [2,4]) thed energy spectra get rescaled roughly by an overall (p 0 /160 MeV) 3 factor.
We performed a Monte Carlo study by generating a huge number of events (up to 10 7 per DM DM annihilation, and up to 10 9 events when studying pp and pp collisions) with Pythia 8 [10], directly implementing the condition ∆k < p 0 = 160 MeV ford production. Such computing-power demanding results have been obtained using the EU Baltic Grid facilities [11]. Fig. 1 shows the total number ofd produced per DM annihilation as function of the DM mass for various annihilation modes, comparing our Monte Carlo result with the spherical approximation, which can under-estimate thed yield by various orders of magnitude. The same p 0 = 160 MeV is assumed in both cases. Fig. 2 shows our Monte Carlo results for thed spectra computed for three values of the annihilating DM mass M . The same spectra also hold for decaying DM, after replacing M ann → M dec /2. As we expected, the result has only a minor dependence on M and is thereby qualitatively different from the 'spherical-cow' approximation that would give a 1/M 2 suppression. There are three classes of qualitatively different cases: DM annihilations i) into W, Z, h (we assume a Higgs mass m h = 120 GeV); ii) into quarks q, b, t or iii) into leptons. The latter case gives nod. To compare the former two cases that gived, we focus on i) DM DM → W + W − and ii) DM DM → qq, and show thed spectra in fig. 3a  qq case thed spectrum at smaller x increases with M rather than being suppressed as 1/M 2 . This is due to QCD FSR that roughly scales as α 3 ln(M/m p ).
The Monte Carlo results differ both qualitatively and quantitatively from the previous studies of thed spectra in DM annihilations or decays. To draw conclusions about the detectability of the signal, we also need to study possible changes in the astrophysicald background, mainly generated by collisions of cosmic-ray p with energy E p on p at rest. In view of the kinematical threshold ford production (E p > ∼ 30 GeV) and of the energy spectrum of cosmic protons (roughly proportional to E −3 p ),d production is dominated by E p ∼ 60 GeV, in the range of validity of the parton model in Pythia. Our semiquantitative results for thed background suggest a reasonable agreement with the spectra of [2,4]. This is an expected result because the center of mass energy in cosmic pp collisions is small and, in this case, the uncorrelated spherical approximation is expected to work reasonably well. However this issue needs to be precisely investigated.
Some remarks are in order. First, we computed p, n allowing all other hadrons to decay despite that the life-time of some strange baryons, such as the Ξ = uss, is longer than the size of deuterium. This effect should already have been taken into account when extracting the value of p 0 from high-energy experimental data from its definition of eq. (5). Second, DM in general annihilates into various primary channels k. According to eq. (1) one should sum their contributions to thep,n spectra rather than to thed spectrum, getting ( k dN (k) /dx) 2 = k (dN (k) /dx) 2 . Our Monte Carlo result instead amounts to sum incoherently over all primary annihilation channels k as well as all secondary and tertiary contributions in the decay chain.

Cosmological fluxes
To compute thed flux in the solar system we consider four possible Milky Way DM density profiles ρ(r) [12]: (1 + r 2 /r 2 s )/(1 + r 2 /r 2 s ) isothermal, r s = 5 kpc (r /r)(1 + r /r s ) 2 /(1 + r/r s ) 2 NFW, r s = 20 kpc (r /r) 1. 16 (1 + r /r s ) 2 /(1 + r/r s ) 1.84 Moore, r s = 30 kpc exp(−2[(r/r s ) α − (r /r s ) α ]/α) Einasto, r s = 20 kpc, α = 0.17, keeping fixed the local DM density ρ(r = r ) = ρ = 0.3 GeV/ cm 3 . Concerning diffusion of chargedd in the galaxy, we approximate the diffusion region as a cylinder with height 2L centered on the galactic plane, a constant diffusion coefficient K = K 0 E δ and a constant convective wind directed outward perpendicularly to the galactic plane. We consider the min, med, max propagation models [13] forp,d, which are characterized by the following astrophysical parameters, Finally, one must take into account annihilations ofd on interstellar protons and Helium in the galactic plane (with a thickness of h = 0.1 kpc L) with rate Γ ann [2]. The solution to the diffusion equation for the energy spectrum of thed number density, f ,  acquires a simple factorized form in the "no-tertiaries" approximation that we adopt. Thē d flux in the galactic medium around the solar system can be written as fully analogous to the solution for thep flux in [14]. The function dNd/dT contains the particle physics input and was computed in the previous section. The function R d (T ) encodes the Milky Way astrophysics and is plotted in fig. 4b for various halo and propagation models. It roughly is some average containment time in the diffusion cylinder, and we verified thatd generated outside it provide a negligible extra contribution even in the min scenario, where most DM annihilations occur outside the diffusion cylinder: the probability of re-entering is sizable, but the probability of diffusing up to the solar system is small. Going from DM annihilations to DM decays with life-time τ one just needs to replace in eq. (11) σv ρ 2 /2M 2 with ρ /M τ ; we do not plot the corresponding R d (T ) functions for DM decay as they essentially coincide with the R d function for DM annihilations and the isothermal profile plotted in fig. 4b. Indeed, for all the considered DM profiles, DM decays close to the galactic center do not significantly contribute to thed flux at Earth, as for DM annihilations with the quasi-constant isothermal density profile. We notice that although R d (T ) is significantly uncertain (especially below a few GeV), the ratio with the corresponding astrophysical function R(T ) forp is essentially fixed, so that the non-observation of a DMp excess puts robust bounds on the possible DMd flux. Indeed thep flux has been observed below 100 GeV by PAMELA [15] and agrees with astrophysical expectations, which are believed to have an uncertainty of about ±20% [4,5]. Finally, we take into account the solar modulation effect, relevant only for nonrelativisticd: the solar wind decreases the kinetic energy T of charged cosmic rays such that the energy spectrum dΦd ⊕ /dT ⊕ ofd that reach the Earth with energy T ⊕ is approximatively related to their energy spectrum in the interstellar medium, dΦd/dT , as [18] The so called Fisk potential φ F parameterizes in this effective formalism the kinetic energy loss. We assume φ F = 0.5 GV i.e. eφ F = 0.5 GeV.

Fig. 5 compares our Monte
Carlo results for thed flux with the spherical approximation. 2 The shading indicates the enhancement. We here assumed the NFW profile, MED propagation. and the DM annihilation cross section σv = σv cosmo ≡ 3 10 −26 cm 3 /sec that reproduces the cosmological DM abundance via thermal freeze-out.  Figure 6: Thep/p ratio, for the same DM models described in the caption of fig. 7, showing that they are compatible with the PAMELAp data (red points). Shading indicates the expected astrophysical background.  Rather than relying on theoretical assumptions, in order to explore the maximald flux from DM compatible with present data, we assume σv = max(1, M/300 GeV) 2 · σv cosmo .
Indeed, fig. 6 shows that these assumptions give ap flux compatible with (and comparable to) the PAMELAp/p data. As discussed in the previous section, thep/d ratio is negligibly affected by astrophysical uncertainties. Furthermore, the assumed cross section is about one order of magnitude below what is needed to explain the PAMELA [19] e + excess and is compatible with the bounds from galactic γ and ν observations [20] as well as with the diffused γ-ray constraints [21]. Then, the upper row of fig. 7 shows our Monte Carlo results for thed flux, T · dΦd/dT, while the lower row shows the corresponding much lowerd flux obtained in the spherical approximation. The caption describes all various assumptions.
Comparison of these two results shows that the signal is enhanced in our Monte Carlō d computation: almost an order of magnitude for smalld kinetic energies (T ∼ 1 GeV) or lighter DM (M ∼ 100 GeV) and orders of magnitude at higher energies or for heavier DM. Such enhancement does not depend on the assumed value of σv Before our calculation it was believed that only the sub-GeV energy region is suitable for searches of a DM-induced d signal. Our result implies that heavy DM, as suggested by PAMELA and FERMI data, also induces detectabled signal at high energies. Therefore our result has important implications on the strategy of DM searches using thed signal.
The red line in fig. 7a roughly shows the Minimal Dark Matter [14] prediction: DM with M ≈ 10 TeV that makes Sommerfeld-enhanced annihilations into W + W − , giving rise top and consequently tod at energies (per nucleon) above m p M/M W , not yet explored by PAMELA.
The PAMELA [19], FERMI [22] and HESS [23] e ± excesses suggest a DM interpretation in terms of multi-TeV DM that annihilates dominantly into leptons with a Sommerfeld-enhanced cross section [17,24]. An interesting class of models with these properties is obtained by assuming that DM annihilates into a new vector with mass m < 2m p , that subsequently can only decay into the lighter e, µ, π [25]. We notice that this condition is not strictly necessary neither for the Sommerfeld enhancement nor for compatibility with PAMELAp data: indeed if m > ∼ 2m p one would obtainp with energy larger than m p M/m, where M/m is the boost factor of the new vector. This boost is large enough not to give an unseenp excess below 100 GeV (the energy range explored by PAMELA so far) even if m is several tens of GeV, as in [17]. Similarly to the Minimal Dark Matter case, these models would give a flux ofd above 100 GeV. Our enhancedd signal should also be used to re-evaluate prospects of discovering supersymmetric Dark Matter candidates, which often annihilate into the W + W − or bb modes we considered.

Conclusions
We computed thed flux at Earth produced by DM annihilations or decays in the Milky Way using an event-by-event Monte Carlo technique run on the Grid, improving on previous computations that assumed spherically symmetric events and obtained a 1/M 2 suppression of thed yield for heavy DM masses M . Due to the jet structure of high energy events implied by relativity no such suppression is present, and thed signal is strongly enhanced: by orders of magnitude ford energies above 10 GeV or DM masses above 1 TeV, as illustrated in fig. 7. Thed astrophysical background seems not to be significantly affected, being dominantly generated by low-energy cosmic ray collisions. While thep andd fluxes suffer from significant astrophysical uncertainties, their ratio is robustly predicted. Thereby the non-observation of ap excess in PAMELA data implies an upper bound on thed DM flux. In the light of our enhancedd fluxes, we find that ad DM signal is still possible. For example, heavy DM models [14,25] that can account for the PAMELA e ± excess can lead top andd excesses above 100 GeV/nucleon.
Most importantly, our result implies that the experiments searching for cosmic rayd become sensitive to M ≥ TeV mass DM, provided that the DM annihilation cross section is larger than what naively suggested by thermal freeze-out. Therefore it is important to extend future searches ford above the GeV energy range. For the moment, the AMS-2 experiment is expected to achieve a very energy-dependent efficiency tod detection, so that AMS-2 would have a sensitivity to ad flux down to 5 10 −7 /(m 2 sec sr GeV/nuc) in the energy ranges 0.2 GeV/nuc < T < 1 GeV/nuc (where time-of-flight is enough to discriminated fromp) and 2 GeV/nuc < T < 4 GeV/nuc (where the magnetic spectrometer is needed) [26]. According to previousd DM computations based on the spherically symmetric approximation, only the lower energy range was promising for DM searches. We have shown that the DM signal can manifest itself also at higher energies, where it is less affected by astrophysical uncertainties.