Antideuterons from decaying gravitino dark matter

We study the possibility of improving the constraints on the lifetime of gravitino dark matter in scenarios with bilinear R-parity violation by estimating the amount of cosmic-ray antideuterons that can be produced in gravitino decays. Taking into account all different sources of theoretical uncertainties, we find that the margin of improvement beyond the limits already set by cosmic-ray antiproton data are quite narrow and unachievable for the next generation of experiments. However, we also identify more promising energy ranges for future experiments.


Introduction
Most current cosmic-ray detectors have the ability to determine the mass of the incoming particles -and in some cases also their charge. However, some isotopes in the cosmic-ray particle spectrum still remain to be observed -among those antideuterons. Donato, Fornengo and Salati were the first who suggested that antideuterons would be an interesting species for the search of dark matter (DM) via cosmic rays [1]. In the past 16 years since this pioneering work, various authors have followed up on this idea [2][3][4][5][6][7][8][9][10][11][12][13]. Indeed, antimatter atoms like antideuterium cannot exist in stars. Therefore, antideuterons -the corresponding nuclei -should not exist as astrophysical primary cosmic rays. 1 The only astrophysical process that is expected to produce it is the spallation of cosmic rays off the interstellar medium (ISM). These so-called secondary antideuterons, for kinematic reasons, cannot be produced at kinetic energies lower than few GeV per nucleon [14]. If, however, DM annihilations or decays lead to the production of antideuterons, their spectrum would not exhibit this kinematic threshold. This leads to the conclusion that the observed antideuteron flux from DM annihilation or decay in the Galactic halo could be orders of magnitude higher than the astrophysical background at very low energies [1].
This theoretical argument motivated the construction of dedicated instruments for lowenergy antideuteron searches like the General AntiParticle Spectrometer (GAPS) [15], which after a successful balloon-flight of a prototype in mid-2012 (pGAPS [16,17]) is expected to fly again in the future with its final design. Current multi-purpose cosmic-ray experiments, like AMS-02 and BESS, are also looking for this type of particles at somewhat higher energies.
Even though it has been understood that because of tertiary production [18] 2 and energy losses taking place during the cosmic-ray propagation [3] the difference between secondary JCAP07(2015)012 and primary fluxes may not be as high as originally expected, antideuterons are still considered one of the most interesting species for DM searches. More recently, this interest has been extended to other yet unobserved cosmic-ray species like antihelium nuclei from DM annihilations or decays [19,20]. The growing activity of the field expresses itself through the organisation of dedicated events such as the recent antideuteron workshop held at the University of California, Los Angeles. A summary paper on the status of the field has recently been completed [21].
In the present work, we investigate the case of gravitino DM within the framework of bilinear R-parity violation [22,23]. 3 In this type of scenario the gravitino would be long-lived enough to be the DM in the universe but would eventually decay and produceamong other species -antideuterons. A main motivation for this scenario is that it leads to a consistent cosmological scenario, explaining the baryon asymmetry of the universe via thermal leptogenesis and avoiding any cosmological gravitino problems [23]. For a more detailed introduction of the model, please consult our previous work, where we studied constraints on the gravitino lifetime from cosmic-ray antiprotons in the same theoretical framework [26].
The structure of this paper is as follows: In the next section, we briefly revise deuteron and antideuteron formation in the coalescence model, derive the coalescence momentum from collider data, and simulate antideuteron spectra for the gravitino decay channels. In section 3, after studying the cross sections at stake, we compute the propagated cosmic-ray antideuteron fluxes at Earth, taking into account various sources of uncertainties. In section 4, we discuss the prospects for the detection of gravitino decays via antideuterons, before we come to our conclusions. In an appendix to this paper, we present some cross section parametrisations relevant for antideuteron propagation.

The coalescence model
Deuterons and antideuterons may form in particle physics processes through final state interactions in a nucleon shower [27]. A usual approach to describe their formation is the phenomenological coalescence model [28,29]. In this prescription, independent of the details of the microscopic formation mechanism, a deuteron is formed when a proton and a neutron come sufficiently close in momentum space, i.e. when the absolute value of the difference of their four-momenta is below a threshold coalescence momentum: In order to determine the value of the coalescence momentum relevant for antideuteron production in gravitino decays we will make use of the antideuteron production rate measured by the ALEPH experiment at the LEP collider [30]. Although there are several measurements of antideuteron production in collider experiments -which in general do not lead to a common value for p 0 [10, 21] -we restrict to the ALEPH data since the gravitino decay channel ψ 3/2 → Zν produces antideuterons via Z boson fragmentation exactly as in the ALEPH measurement. Also the channels ψ 3/2 → W and ψ 3/2 → hν are expected to be more similar to the case of Z boson fragmentation than to the production of antideuterons in pp collisions [31][32][33], e − p collisions [34], or e + e − collisions on or nearby the Upsilon resonance [35,36].

JCAP07(2015)012
A common approximation for deuteron formation is that the distributions of neutrons and protons in momentum space are spherically symmetric and uncorrelated [27]. This leads to an energy distribution of deuterons that can be directly calculated from the product of the individual spectra of neutrons and protons [4,6]: where T p = T n = T d /2 are the kinetic energies of protons, neutrons and deuterons. Therefore, this approximation is called factorised coalescence.
As pointed out in [6], this approximation is qualitatively wrong and significantly underestimates the deuteron yield in high-energetic processes. This is due to the fact that the distributions of neutrons and protons are actually neither spherically symmetric nor uncorrelated. For instance, in the decay of a DM particle into a Z boson and a neutrino the probability for the formation of a deuteron in the fragmentation of the Z boson should be independent of the DM mass. This is due to the fact that the Z boson fragmentation process is always the same as viewed from the Z boson rest frame. However, the factorised coalescence approximation gives a lower yield of deuterons for larger DM masses, since the protons and neutrons are distributed over a larger phase space for higher injection energies. Another qualitatively wrong behaviour of this approximation is the possibility that protons and neutrons from distinct DM decays form a deuteron. This is due to the fact that the spectra of protons and neutrons are simply multiplied, while in principle the coalescence condition on the four-momenta of protons and neutrons should be applied on an event-by-event basis.
This can be achieved in a Monte Carlo simulation of the decay process [6]. For instance, using an event generator like Pythia 6.4 [37] one can simulate the hadronisation of massive gauge and Higgs bosons and then apply, event by event, the coalescence condition on the protons and neutrons. 4 This method leads to plausible results, e.g. the deuteron yield in the DM decay to final states including W , Z or Higgs bosons is independent of the DM mass. However, this strategy also requires a lot of computing time to generate smooth spectra as only one deuteron or antideuteron is produced in O(10 4 ) fragmentation processes. Moreover, the authors of [38] point out that baryon pair distributions are not used to tune event generators and therefore their predictions for deuteron production are not without uncertainty. In fact, the authors of [9] find a discrepancy of up to a factor 2-4 when comparing deuteron production between the Monte Carlo generators Herwig++ [39,40] and Pythia 8 [41].
Due to the deuteron's binding energy of 2.224 MeV [42], the coalescence of free protons and neutrons is forbidden since energy and momentum cannot be conserved at the same time if no further particles are involved in the process. However, in a hadronic shower with many particles, it is easily conceivable that the overall process conserves energy and momentum and thus allows for deuteron coalescence. 5 When simulating deuteron coalescence in a Monte 4 In fact, there are some combinatoric ambiguities in the deuteron coalescence in Pythia. Namely, it is possible that a proton fulfils the coalescence criteria with two different neutrons or vice versa. We checked numerically that the number of events where this is the case is very low. The probability increases with p0, but even for p0 = 250 MeV less than 0.1 % of the deuteron events are affected. Therefore, we conclude this not to be a significant effect. In our analysis we always combine the first pair of protons and neutrons fulfilling the coalescence criterion according to the particle order in the event record. 5 While low-energy deuteron formation proceeds via the process p n → d γ, this is not necessarily the case in high-energy nuclear collisions. In their pioneering work on deuteron formation, Butler and Pearson considered an explicit interaction with the optical potential of the target nucleus to take care of energy and momentum -3 -

JCAP07(2015)012
Carlo event generator, the four-momenta of protons and neutrons are given explicitly and we have no way of realistically treating the multi-body process leading to deuteron formation. The most common approach in previous works on antideuteron signals from DM annihilation or decay is thus to consider momentum conservation and to determine the deuteron kinetic energy from the proton and neutron momenta and the deuteron mass of 1.8756 GeV [42]: We will also follow this approach since it turns out that the uncertainty on the spectra introduced by this ambiguity is marginal for our discussion. 6 More recently, it was realised that the condition |p p − p n | < p 0 is not not enough to guarantee a physically meaningful coalescence prescription in Monte Carlo event generators [10]. Some of the antiprotons and antineutrons emerging in DM annihilations or decays are generated by the decay of metastable mother particles like the baryonsΛ 0 , and Ω − b , and the mesons D ± s , B 0,± and B 0 s , which have lifetimes of O(10 −10 -10 −13 ) s [37,42]. By contrast, the nuclear interactions leading to antideuteron formation take place on femtometer length scales, corresponding to lifetimes of O(10 −23 ) s for relativistic particles. Therefore, it is necessary to include a condition on the spatial separation of protons and neutrons as well.
To accommodate this, Ibarra and Wild excluded weakly decaying baryons with lifetimes τ > 1 mm/c from the decay chain [10]. 7 Although this condition removes the most relevant metastable mother particles, namelyΛ 0 andΣ ∓ , several other metastable particles producing antiprotons and antineutrons have lifetimes just below 1 mm/c. Fornengo et al. introduced the condition ∆r < 2 fm on the spatial separation [11]. In Pythia 6.4, this treatment is basically equivalent to considering only antiprotons and antineutrons with a production time t = 0 since only non-zero lifetimes that are relevant for displaced vertices in collider experiments are stored in the event table. This treatment removes all antiprotons and antineutrons coming from metastable mother particles and this is the treatment we will adopt for our analysis below.

Determination of the coalescence momentum
Using Pythia 6.4, we simulated 10 9 events of the process e + e − → Z to determine p 0 from ALEPH data. The authors of [30] found the number of antideuterons produced per hadronic Z boson decay in the process e + e − → Z →d + X at the Z pole to be 8 Rd = (5.9 ± 1.8 ± 0.5) × 10 −6 (2.4) conservation [27]. However, this model was later found to be in conflict with experimental data [43,44]. The factorised coalescence model does not touch the issue at all since it does not go into the details of the microscopic formation mechanism [29]. In the 1980s it was argued that no explicit interaction with a third body is needed since deuterons are formed in a tiny space-time region and thus the uncertainty principle applies [43], or that the dominant process for deuteron formation in a hadronic shower involves a proton and a neutron that are slightly off-shell [44,45]. Nonetheless, other authors still stressed the necessity of an explicit interaction with a third body in the nuclear shower [46,47]. The authors of [48] recently discussed an alternative deuteron formation model that makes use of measured cross sections for deuteron production processes with photon or pion emission. 6 We checked how alternative approaches affect the resulting deuteron spectra. Starting from energy conservation, i.e. determining the deuteron kinetic energy from the proton and neutron energies and the deuteron mass to be T d = Ep + En − m d , the resulting spectra differ only notably at deuteron momenta below the coalescence momentum. 7 The exact criterion is not stated in the paper; S. Wild, private communication (2014). 8 The first error gives the statistical error, while the second error corresponds to systematic errors.

Coalescence Model
Coalescence Table 1. Coalescence momenta for antideuteron production determined from the ALEPH data. We compare the factorised coalescence approach using antiproton and antineutron spectra generated with Pythia to the event-by-event coalescence using either only the momentum condition or additionally excluding antiprotons and antineutrons produced in decays of long-lived mother particles.
within the momentum range 0.62 GeV < pd < 1.03 GeV and in the angular range | cos θ| < 0.95. In the left panel of figure 1 we present the simulated antideuteron yield as a function of the coalescence momentum. We show the total antideuteron yield per Z decay event as well as the yield per hadronic Z decay applying the momentum and angular cuts of the ALEPH analysis. The hadronic branching ratio of the Z boson is 69.91 %. Leptonic Z decays do not contribute at all to antideuteron production [42]. We overlay the ALEPH result, adding statistical and systematic errors in quadrature. The best-fit values for the coalescence momentum and the 1-σ ranges can be basically read off this plot and are summarised in table 1 along with the value following from the factorised coalescence model, see eq. (2.2). In the right plot of figure 1 we demonstrate that in all cases the dependence of the antideuteron yield on the coalescence momentum is in very good agreement with the expected p 3 0 behaviour. In the left panel of figure 2, we present the antideuteron spectrum in Z decays as a function of the momentum for the three different treatments of coalescence we discussed above. The spectrum derived from the factorised coalescence approximation is completely different from the spectra based on event-by-event coalescence. By contrast, the exclusion of antiprotons and antineutrons from metastable mother particles only slightly changes the shape of the high-energy part of the spectrum. The single ALEPH data point cannot give any constraints on the shape of the spectrum. Therefore, all coalescence treatments fit the Antideuteron Momentum p GeV Antideuteron Spectrum p dN dp Gravitino Mass m3 2 GeV Deuteron Antideuteron Yield Figure 2. Left: Simulated antideuteron spectra from the process e + e − → Z →d+X compared to the measurement of the ALEPH experiment at LEP. The spectrum obtained from factorised coalescence is shown in green. The spectra derived from Pythia are shown in blue (only ∆p < p 0 ) and red (∆p < p 0 and t = 0). The two Monte Carlo spectra are very similar except for the high-energy part. The spectrum derived from uncorrelated antiprotons and antineutrons exhibits a very different shape. Right: Antideuteron yields in the decay channels Zν, W and hν as a function of m 3/2 . The yields obtained from the Pythia simulation are independent of the gravitino mass and very similar to each other. By contrast, the factorised coalescence prescription leads to unphysical behaviour, with yields falling off like m −2 3/2 . data equally well and our choice in favour of the Monte Carlo simulation with ∆p < p 0 and t = 0 is only based on the theoretical considerations discussed above. For the rest of the paper we will stick to the event-by-event coalescence using p 0 = 203 +20 −25 MeV and excluding antiprotons and antineutrons produced at displaced vertices.
Let us conclude this section with a comparison to other p 0 determinations in the literature. Our p 0 values differ somewhat from the values found by the authors of [6,10] and [11]. Kadastik et al. used Pythia 8 and found p 0 = 162 ± 17 MeV for the ∆p < p 0 condition; roughly 10 MeV below our result. Ibarra et al. also used Pythia 8 and found a coalescence momentum of p 0 = 192 ± 30 MeV when removing long-lived mother particles; again 10 MeV below our result. Fornengo et al. used Pythia 6.4 and found p 0 = 195 ± 22 MeV when excluding metastable mothers. For the ∆p < p 0 condition they found p 0 = 180 ± 18 MeV and for the factorised coalescence approximation they found p 0 = 160 ± 19 MeV. 9 Putze finds p 0 = 202 MeV when using Pythia 6.4, a value similar to ours, and p 0 = 195 MeV when using Pythia 8, compatible with Ibarra et al. [49]. As mentioned before, differences between different event generators are not completely unexpected. However, in some cases there are also different results for the same tool. These differences are a bit worrisome since the Monte Carlo statistical uncertainty is only of the order of few MeV. A potential source of systematic uncertainty is the choice of the Monte Carlo tune, but to our knowledge all of the studies listed above use the default tune of Pythia 6.4/8.

Antideuteron spectra from gravitino decay
In models with bilinear R-parity violation, the gravitino has four main decay channels: ψ 3/2 → γν, Zν, W and hν [50,51]. As for the case of antiprotons [26], only the channels containing a massive gauge or Higgs boson in the final state are relevant for antideuterons.  Table 2. Multiplicities of deuterons + antideuterons from gravitino decays after the event-by-event coalescence process simulated with Pythia 6.4. The central values are given for p 0 = 203 MeV (t = 0) and p 0 = 173 MeV (p < p 0 ), respectively. The quoted errors of roughly 30% correspond to the 1-σ uncertainty in the determination of p 0 . Kinetic Energy x 2 T m 3 2 Deuteron Antideuteron Spectrum x dN dx Kinetic Energy x 2 T m 3 2 Deuteron Antideuteron Spectrum x dN dx  For the generation of the antideuteron spectra from gravitino decay we simulated 10 9 events with Pythia 6.4 for each of the decay channels and for a set of gravitino masses of roughly equal distance on a logarithmic scale: m 3/2 = 85 GeV, 100 GeV, 150 GeV, 200 GeV, 300 GeV, 500 GeV, 1 TeV, 2 TeV, 3 TeV, 5 TeV, and 10 TeV. Using the same strategy as in [26], we started the Pythia simulation with a resonance decay into two particles, Z boson and neutrino, W boson and charged lepton, and Higgs boson and neutrino, respectively. In this way the Z, W and Higgs bosons are treated as decaying isotropically in their rest frames.

JCAP07(2015)012
The antideuteron yields in these decay channels are expected to be independent of the gravitino mass since the antideuterons are produced in the fragmentation of an on-shell gauge or Higgs boson. Larger gravitino masses should thus only lead to boosted spectra, while the underlying physical process remains unchanged. 10 Our simulation confirms this expectation and shows that the factorised coalescence prescription leads to unphysical results, see the right panel of figure 2. The antideuteron yields per gravitino decay are summarised in table 2. 11 The resulting spectra are presented in figure 3 for the central value of p 0 = 203 MeV. 12 Although these spectra by eye appear to have the same shape as the antiproton spectra from gravitino decay (see [26]), rescaled by a factor of O(10 −4 ), there is no simple scaling relation among them. 10 Note that we did not take into account weak corrections that would lead to a slight increase of the yields with increasing gravitino mass, see for instance [52]. 11 Note that these results differ from those reported in [53] due to an erroneous treatment of the Pythia routine in the earlier study. 12 The decay spectra are available in tabulated form at http://www.desy.de/~mgrefe/files.html.

JCAP07(2015)012
We will now move on to the discussion of the antideuteron flux observed at Earth. The flux expected from gravitino decay is simply a linear combination of the fluxes for the individual decay channels: The branching ratios for the different decay channels depend on the choice of the supersymmetry parameters and are discussed in detail in [26,53]. Note that they do not depend on the amount of R-parity violation but only on the mass of the gravitino and the mass hierarchy of the neutralino sector of the supersymmetric particle spectrum. It is hence convenient to label the different cases by the name of the next-to-lightest supersymmetric particle (NLSP).
In this work we use the values presented in figure 3 and table 3 of [26]. For illustration, we restrict to an example case where the NLSP is Bino-like.

Antideuteron flux from gravitino decay
As explained in the introduction, antideuterons allow for an almost background-free search for an exotic DM component in certain parameter ranges. In fact, no cosmic-ray antideuterons have been observed so far and there exists only an upper limit on the antideuteron flux from the BESS experiment [54]. In addition, BESS-Polar II looked for antideuterons using more than ten times more cosmic-ray data than the previous BESS analysis. No candidate antideuterons were observed and a flux limit is expected to be published in the near future [55,56]. In addition, several experiments are currently taking data or will start operating within the next years to improve this situation: The AMS-02 experiment operating on the International Space Station is expected to greatly improve on the current sensitivity to antideuteron fluxes [57]. 13 Moreover, the General AntiParticle Spectrometer (GAPS) is expected to perform several balloon flights, starting with a first Antarctic campaign in the austral summer 2019/2020 [21]. A prototype flight of the GAPS experiment was successfully carried out in June 2012 [17]. Several studies on antideuteron fluxes from DM annihilations or decays can be found in the literature [1][2][3][4][5]. These early studies, however, employ the factorised coalescence approximation for antideuteron formation, which, as discussed in section 2, is in general insufficient to describe the actual production rate. Only more recent studies employ the Monte Carlo approach [6][7][8][9][10][11][12][13]. As discussed in section 2, in this work we employ decay spectra obtained by this latter method.
A relativistic antideuteron, formed by the coalescence of an antineutron and an antiproton within the hadronic shower of a gravitino DM decay in the Milky Way halo, will then propagate through the ISM and might eventually arrive at a detector at Earth. As for all cosmic rays, the propagation is described by a diffusion equation of the cosmic-ray phase-space density ψ: For a description of the individual terms we refer to the appendix of [26] and references therein.

Cross sections
The spallation term Γ spal = vd (n H σ ann dp + n He σ ann dHe ) deserves a separate discussion. This term corresponds to the interaction of cosmic-ray antideuterons of velocity vd with the ISM, which is mainly composed of hydrogen and helium. We have considered n H = 0.9 cm −3 and n He = 0.1 cm −3 [58]. Contributions of heavier elements have been neglected. The annihilating inelastic cross section σ ann dp = σ inel dp − σ non-ann dp is practically given by the inelastic cross section since the non-annihilating inelastic cross section is very small due to the small binding energy of antideuterons [3,4,9,18]. The inelastic cross section is the difference of the total and the elastic cross sections, σ inel dp = σ tot dp − σ el dp . A common issue in cosmic-ray antideuteron analyses is that there are no experimental data on the cross section of inelastic antideuteron-proton scattering. A typical assumption is hence to rely on the hypothesis that σ inel dp (Td/n) = 2 × σ inel pp (Tp = Td/n), where Td/n is the antideuteron kinetic energy per nucleon, see [10] using a direct parametrisation of σ inel pp data by Tan and Ng [59]. 14 A similar path was followed by [3], calculating the cross section as σ inel pd (pp) = 2 × (σ tot pp (pp) − σ el pp (pp)). Making use of the reasonable assumption of symmetry under charge conjugation, one can also use the total antiproton-deuteron cross section σ tot pd (pp), for which data exists [42]. In this case, the inelastic antiproton-deuteron cross section was calculated as σ inel pd (pp) = σ tot pd (pp) − 2 × σ el pp (pp) [4,9,11]. However, although the assumptions presented above give the correct order of magnitude of σ inel dp , they are not entirely supported by experimental evidence. In fact, 2 × σ tot pp is larger than σ tot pd by roughly 10%. This is theoretically expected due to Glauber screening [60] and in nucleon-nucleus collisions one rather expects a geometric scaling σ pA A 2/3 × σ pp than a linear scaling with nucleon number A [61][62][63].
In this work, we will thus use an inelastic antiproton-deuteron cross section based on parametrisations of available data. Indeed, in many previous works it was incorrectly stated that there were no data on the inelastic and elastic antiproton-deuteron processes. In figure 4, we compare our parametrisation of σ inel pd with other parametrisations used in the literature. We observe that the approach of rescaling σ inel pp by a factor of two overshoots the available low-energy data. The parametrisation by Tan and Ng [59] is only based on σ inel pp data from 50 MeV to 100 GeV in antiproton kinetic energy and clearly leads to an incorrect high-energy behaviour. The approach of Donato et al. [3], based on the difference of σ tot pp and σ el pp data that extend to much higher energies, leads to a more reasonable high-energy behaviour. Dal et al. [9], 15 using parametrisations of σ tot pd and σ el pp , find a better agreement with low-energy data. Our parametrisation gives a comparable result, but makes use of better motivated functional forms for the parametrisations and explicitly takes into account available lowenergy data for the antiproton-deuteron process. See appendix A for a detailed derivation of the parametrisation used in this work. For σ non-ann pd we use the parametrisation of Dal et al. [9]. Our result is certainly not yet satisfactory, but we think that it is definitely an improvement compared to the treatment of cross sections in earlier antideuteron studies. We hope that our work serves to stimulate further discussion in this area.
In addition to the antideuteron-proton cross section, we need the cross section for annihilating antideuteron-Helium scattering. Since no experimental data are available for this 14 Some earlier studies used the assumption σ ann dp (Td/n) = 2 × σ ann pp (Tp) [5,7]. This relation, however, underestimates σ ann dp since the annihilating antiproton-proton cross section falls quickly with rising energy, in contrast to the antideuteron-proton process, see appendix A. 15 The used parametrisations are not given in the paper; L. Dal  process, we have rescaled the antideuteron-proton cross sections by the geometric factor 4 2/3 , hoping that it suffices to give a reasonable estimate of the true cross section.
Concerning the production of tertiaries, for which the differential cross section of the non-annihilating inelastic antideuteron-proton process is required, unfortunately no data are available. The usual method is thus to take the integrated cross section of the chargeconjugate process, σ non-ann pd , and to multiply it by an appropriate energy distribution of the outgoing deuteron. In the limiting fragmentation hypothesis, this is simply 1/T , where T stands for the kinetic energy of the incoming deuteron (see for instance [18,59]). However, one can also follow the authors of [3,18] and -inspired by the pp process -use the functional form suggested by Anderson et al. [65]: where dσ Anderson dp ≡ 2π π 0 d 2 σ(pp → pX) dp dΩ sin θ dθ and In the latter expression, the energy and momentum of the incoming proton are labelled E and p, respectively; p t stands for the transverse momentum of the outgoing proton in units of GeV. The boost from the laboratory frame to the centre-of-mass frame is ruled by the Lorentz coefficients β and γ. We follow this latter method in this work.

Propagated fluxes
As for the case of antiprotons, no primary antideuterons are expected from astrophysical objects and the dominant background for DM searches are secondary antideuterons created in spallation processes of cosmic-ray protons and helium nuclei impinging on the ISM, i.e. hydrogen and helium gas. In order to estimate the background for the signal from gravitino decay, we employ here the same calculation of the astrophysical secondary antideuteron flux as used in [3]. As one can see in figure 5, various processes have to be taken into account for the calculation of the antiproton and antideuteron fluxes. The relevant processes for antiprotons are presented in table 3 and those for antideuterons are presented in table 4. For the antiproton processes, we have used the cross section parametrisations given in [66]. For the calculation of secondary antideuterons, we have used the antiproton cross sections along with the factorised coalescence prescription for antideuteron formation in the final state, see section 2. 16 Since both primaries and secondaries produce tertiaries, i.e. cosmic rays produced by the non-annihilating inelastic scattering of high-energy antideuterons on the ISM, we have not shown this component separately but rather added it directly to our estimates of the primaries and secondaries, respectively. However, we stress that this process, as well as convection and diffusive reacceleration, have important impact on the low-energy part of the computed fluxes and should not be neglected. Note also that in the case of antideuterons the secondaries coming from the spallation of primary antiprotons, i.e. antiprotons of DM origin, on the ISM are to be considered as signal and not background. This component is much lower than the primary component and the background when the gravitino mass is low, but for masses higher than 10 TeV, it is of the same order of magnitude as the primary component and hence should not be neglected. Figure 5 displays these different components for a gravitino mass of 100 GeV.
In order to be able to compare the expected antideuteron fluxes to experimental results, one has to take into account the effect of solar modulation [68]. Ideally, one should try to fully model this effect, for instance with a simulation like HelioProp [69] as done in [11]. 16 Alternatively, one could also use the Monte Carlo approach for calculating secondaries [67].  The primary antiproton and antideuteron components come from gravitino decay; for the case of antideuterons also from the interaction of primary antiprotons with the ISM. The tertiary component corresponds to a redistribution of high-energy cosmic rays to lower energies due to non-annihilating inelastic scattering off the ISM. Since every component creates tertiaries, they have been incorporated directly and are not displayed as separate components.
However, this requires knowing precisely the status of the solar environment at the time of data taking. In the absence of data and because we do not know yet the time of the data taking, we think it is untimely to go through such a precise modelling. For figure 5, we hence satisfy ourselves with the so-called Fisk approximation [70]: assuming a Fisk potential of φ F = 500 MV as an example since this values allows to have a good agreement with the most recent PAMELA antiproton data [71]. In the subsequent figures, where we do not show antiproton data, we will only consider interstellar fluxes since we do not know what the Fisk potential will be when data will finally be taken.
In figure 6, we present the interstellar antideuteron spectrum from gravitino DM decays and compare it to the expected astrophysical background and the flux limit obtained by the BESS experiment [54]. In addition, we present the projected sensitivity regions of the BESS-Polar II [56], AMS-02 [57] and GAPS [17] 17 experiments. In the left panel, we show the uncertainty band due to the lack of precise knowledge of the propagation parameters as constrained by measurements of the boron-to-carbon ratio in cosmic rays [72]. The coloured bands correspond to a full scan over the allowed parameter space. As an illustration we also 17 The GAPS sensitivity assumed in this work corresponds to three Antarctic long duration balloon flights with a total duration of 105 days [21].   Table 5. Parameters of cosmic-ray propagation models that correspond, respectively, to the best fit of cosmic-ray boron-to-carbon data (MED) as well as the minimal (MIN) or maximal (MAX) antiproton flux compatible with cosmic-ray boron-to-carbon data. Figures taken from [73].
display the fluxes obtained with three benchmark models often used in the literature (see table 5). We use a common decay lifetime of τ 3/2 = 10 28 s for illustration. Note that, as for the case of antiprotons, the MIN/MED/MAX benchmark models do not size the full extent of the uncertainty band. This shows again the importance of performing scans over the full propagation parameter space allowed by other data rather than checking only a few cases. In the right panel of figure 6, we display the same fluxes but setting the decay lifetime to the minimum values allowed by the constraints obtained using the antiproton measurements (see [26]). Note that the coloured bands correspond to the extreme cases obtained for antiprotons and are not the full uncertainty band. Indeed, the limiting cases correspond to situations where the antiproton flux is maximised at a given energy bin. This does not mean that the antideuteron flux is maximal over the whole energy range.
As one can see from figure 7, if the Fisk approximation describes correctly the influence of solar modulation on antideuterons, then this should not affect our conclusions dramatically. Indeed, unlike many other cosmic-ray species, antideuteron fluxes are relatively flat below 10 GeV and a shift of the antideuteron energy does not affect the flux very strongly, at least within a reasonable range of the Fisk potential.  . Same as the right panel of figure 6 for the MED propagation parameters. The only difference is that the fluxes displayed are now corrected for solar modulation (φ F = 300 MV and φ F = 700 MV). As one can see, because the fluxes are relatively flat below 10 GeV, solar modulation has little impact on the expected fluxes.

Discussion of the detection prospects
A striking feature of figure 6 is that for masses as low as 100 GeV the gravitino decay signal can be of the same order as the astrophysical background below a few GeV, even for lifetimes as large as 10 28 s, a value not yet excluded by gamma-ray and antiproton observations (see for instance [74]). 18 In this respect, it would also be interesting to see what antideuteron flux could be expected for even lower gravitino masses. It could thus be worthwhile to study this region in a future work using the spectra obtained from gravitino three-body decays [53].
But also for larger gravitino masses the antideuteron signal could be at the same order as the background. This was not observed in earlier studies as the signal for large DM masses is artificially suppressed in the factorised coalescence prescription. Therefore, for decaying DM candidates there is in principle also the possibility of observing an exotic component in the higher-energetic part of the spectrum, where currently no experiments are planned. Note, however, that both background and signal are extremely low and quite challenging for experimentalists as this would mean improving sensitivity by at least four orders of magnitude in flux, but also to reach much higher energies.
When taking into account the constraints derived from antiproton observations [26], we find that the remaining parameter space for having a gravitino decay signal significantly higher than the astrophysical background becomes quite small but does not vanish completely. Since the coalescence process still suffers large theoretical uncertainties, one cannot exclude that all the fluxes are in fact larger (or smaller) than what we assume here. The 1-σ uncertainty in p 0 could lead to an increase of roughly 30% in the antideuteron flux from JCAP07(2015)012 gravitino decay, irrespective of the antiproton constraints. In addition, making use of Monte Carlo methods to estimate the production of secondary antideuterons instead of the factorised coalescence model used here, tends to predict a slightly lower background level at low energies [67]. This would increase the signal-to-background ratio, independently of any other constraints.
Still, it seems clear that the current generation of experiments will not be able to observe any antideuteron events. Only an improvement of the flux sensitivity by two orders of magnitude should at least allow for a detection of astrophysical antideuterons -or even those coming from gravitino decay. The main hopes seem to reside either in the highest energy range (above ∼ 50 GeV) or in the lowest one (below ∼ 1 GeV). Note, however, that the latter is affected by solar modulation, a phenomenon that to date is not fully under modelling control.
This clearly challenges antideuterons as the golden channel it has long thought to be. Indeed the antiproton channel has become extremely constraining thanks to the PAMELA data [71]. A forthcoming release of AMS-02 antiproton data could make these constraints a bit more stringent, especially at higher energies.

Conclusion
In this work, we have studied the potential of detecting cosmic-ray antideuterons produced in the decay of gravitino dark matter within a framework of bilinear R-parity violation. This work was a natural sequel of a related work concerning cosmic-ray antiprotons. After discussing the deuteron formation in hadronic showers and calculating antideuteron spectra from gravitino decay, we have determined the gravitino decay signal at Earth. We have also assessed the uncertainties affecting the expectations for the astrophysical background and the signal. We have shown that there is some room left for a discovery of gravitino decays through antideuterons, however not within the sensitivity of the current and planned generation of experiments.
We have also shown that -once the antiproton constraints are taken into accountthe remaining parameter space for a detection of gravitino dark matter with bilinear R-parity violation via antideuterons is quite small. On the other hand, since not much progress is expected in the background-limited antiproton channel in the coming years, the antideuteron channel could still serve to put stronger constraints on the strength of the R-parity violation in the future. If the detection technology were to improve considerably, also the high-energy regime (above ∼ 50 GeV) would become interesting for the search of gravitino dark matter.

A Cross sections
In this appendix, we briefly present the parametrisations for the antiproton-proton and antiproton-deuteron cross sections we used to estimate the inelastic scattering cross sections relevant for the spallation term of antideuteron propagation in the Milky Way.
Antiproton-proton cross sections. For the total cross section of antiproton-proton scattering, a large number of data points ranging from roughly 200 MeV to 2 PeV in antiproton momentum (in the rest frame of the proton target) exist in the literature [42]. A useful parametrisation is given by [78,79]: where s = 2 m 2 p + 2 m p m 2 p + p 2 p is the centre-of-mass energy of thepp system and σ tot asmpt (s) = a 0 + a 2 ln 2 ( In the definition of R 0 , the factor 1 mb = 2.568 GeV −2 simply accounts for unit conversion. The remaining parameters in these expressions were determined from a fit to experimental data in [78]: With these parameters, eq. (A.1) gives a relatively good fit to the available data. Including also high-energy proton-proton data between 10 TeV and 2 EeV in proton momentumwhere the antiproton-proton and proton-proton cross sections should be equivalent -we get a goodness of fit of χ 2 /dof = 2821/460. 19 Given that many early experimental works, especially from the 1960s and 1970s, only quote statistical errors without any assessment of systematic uncertainties, the fit is quite acceptable (see also the remark on errors in the introduction of [64]). For the elastic cross section of antiproton-proton scattering, there is also a lot of data ranging from roughly 200 MeV to 2 PeV in antiproton momentum [42]. Using typical basis JCAP07(2015)012 functions for cross sections (see [64]), we fit the data using the piecewise ansatz: 20 and requiring differentiability at pp = 5 GeV. We find the parameters giving a goodness of fit of χ 2 /dof = 364.6/167. 21 Besides data on elastic antiproton-proton scattering we included data on elastic proton-proton scattering with incoming momentum between 100 GeV and 34 PeV [42].
For the inelastic antiproton-proton cross section there is considerably less data than in the previous cases, ranging only from 300 MeV to 175 GeV [64]. This cross section is, by definition, the difference between the total and the elastic cross section and we just check the consistency with available data. Including also data on inelastic proton-proton scattering between 1 and 2 TeV, we find a goodness of fit of χ 2 /dof = 404/43. Given that only three of the 43 data points come with an estimate of systematic uncertainties, we think this is an acceptable result. 22 An overview of the antiproton-proton cross section data is presented in figure 8. In addition to the total, elastic and inelastic cross sections, we present data on annihilating inelastic antiproton-proton scattering ranging from 240 MeV to 22 GeV [64]. These data can be parametrised by the function Antiproton-deuteron cross sections. The total cross section for antiproton-deuteron scattering is almost twice as large as the total antiproton-proton cross section. The exact result depends on the size of the elastic and inelastic Glauber shadow or screening corrections [60,81,82]: σ tot pd = 2 σ tot pp − δ el − δ inel . If these correction terms are known, also the elastic and inelastic cross section components can be easily determined individually: Arkhipov describes antiproton-deuteron scattering and the corresponding correction terms in [82]. Unfortunately, his parametrisations do not describe very well the available data. Therefore, we will follow an alternative route. 20 The authors of [80] present a parametrisation of the elastic antiproton-proton cross section based on the methods of [78,79]. However, it appears that their parameters are faulty since we were not able to reproduce their result. 21 We have also checked that the cross section data for the elastic antiproton-neutron [42] and antineutronproton [64] scattering processes are compatible with the antiproton-proton result. 22 We have also checked that the inelastic antiproton-neutron scattering data [64] are compatible with the antiproton-proton result.   Figure 8. Left: Total, elastic, inelastic and annihilating inelastic cross sections of thepp process. Data points are taken from the Particle Data Group [42] and the Landolt-Börnstein compilation [64]. The lines are the parametrisations discussed in the text. Right: Same as left panel but for the cross sections of thepd process. For the total cross section of antiproton-deuteron scattering, there is only data from roughly 300 MeV to 280 GeV in antiproton momentum [42]. For the high-energy part we thus assume that the cross section matches the antiproton-proton cross section, rescaled with a suitable factor. We find that eq. (A.1), rescaled with a factor of 1.85, gives a reasonably good fit to the antiproton-deuteron data above 6 GeV. Using the same basis functions as for the elastic antiproton-proton cross section, we then fit the data using the piecewise ansatz: giving a goodness of fit of χ 2 /dof = 452.9/159. For inelastic antiproton-deuteron scattering, there are data from roughly 300 MeV to 2 GeV [64]. As for the case of antiproton-proton scattering, we will determine the inelastic cross section by subtracting the elastic from the total cross section. Lacking a physically motivated ansatz for the elastic antiproton-deuteron cross section, we rescale eq. (A.2) by a suitable factor to match the data points. We find that σ inel pd (pp) = σ tot pd (pp) − 1.75 × σ el pp (pp) (A.5) gives a goodness of fit of χ 2 /dof = 15.8/13. We still have to check if our approach for the elastic cross section matches the available data ranging from roughly 300 MeV to 5 GeV [64]. As can be seen from the right panel of figure 8, the shape of the rescaled elastic cross section does not entirely follow the data points. The goodness of fit thus is a very poor χ 2 /dof = 535/9. 23 Considering only the low-energy data would result in an acceptable χ 2 /dof = 46/9. 23 Note that the agreement for the elastic cross section would be even much worse if we had determined the inelastic antiproton-deuteron cross section by rescaling the inelastic antiproton-proton cross section.