Gravitino LSP and leptogenesis after the first LHC results

Supersymmetric scenarios where the lightest superparticle (LSP) is the gravitino are an attractive alternative to the widely studied case of a neutralino LSP. A strong motivation for a gravitino LSP arises from the possibility of achieving higher reheating temperatures and thus potentially allow for thermal leptogenesis. The predictions for the primordial abundances of light elements in the presence of a late decaying next-to-LSP (NSLP) as well as the currently measured dark matter abundance allow us to probe the cosmological viability of such a scenario. Here we consider a gravitino-stau scenario. Utilizing a pMSSM scan we work out the implications of the 7 and 8 TeV LHC results as well as other experimental and theoretical constraints on the highest reheating temperatures that are cosmologically allowed. Our analysis shows that points with T_R>10^9 GeV survive only in a very particular corner of the SUSY parameter space. Those spectra feature a distinct signature at colliders that could be looked at in the upcoming LHC run.


Introduction
The phenomenology of supersymmetric scenarios both at colliders and in the early universe depends strongly on the nature of the lightest supersymmetric particle (LSP). The LSP is stable in the R-parity conserving case and thus is usually identified with the dark matter (DM) candidate, if supersymmetry (SUSY) is to explain this observation. In a neutralino LSP scenario with a gravitino mass of the order of the other sparticle masses, a cosmological problem appears once we want to explain the observed baryon asymmetry in the universe with the mechanism of thermal leptogenesis [1]. For this mechanism to work the universe has to be heated up to temperatures of T R 10 9 GeV [2,3] in the post-inflationary phase of reheating. On the other hand, from thermal scattering in the hot plasma gravitinos are produced [4,5] and the abundance of thermally produced gravitinos is proportional to T R [6,7,8]. Hence, a large T R leads to a large number density of gravitinos in the early universe. The Planck-suppressed couplings of the gravitino lead to a delayed decay into the LSP. These decays cause an additional energy release at or after the time of big bang nucleosynthesis (BBN) [4,9,10]. The abundances of light elements are very sensitive to such processes and thus from their precise determination strong bounds can be imposed on the abundance of late-decaying gravitinos [11]. These bounds clearly exclude a reheating temperature of T R 10 9 GeV. This problem is known as the gravitino problem [12].
One way of avoiding this problem is a gravitino LSP scenario. Indeed, the gravitino is a perfectly good DM candidate [13,14]. However, in this scenario the next-to-LSP (NLSP) usually becomes long-lived and might spoil successful BBN predictions [5]. In contrast to the former scenario, it is now the abundance (and the life-time) of the late-decaying NLSP which governs the phenomenological viability of the scenario in this concern. For an NLSP belonging to the sparticles of the MSSM-sharing the SM interactions-the abundance is determined by the thermal freeze-out (rather than the reheating temperature). The abundance of the NLSP depends upon the spectrum parameters of the model and could, in principle, be determined from measurements at colliders. One of the most promising NLSP candidates in this concern is a charged slepton leading to a rather clean signature at colliders [15,16]. In the upcoming high-energy run of the LHC such a scenario could reveal a rich phenomenology.
In a gravitino LSP scenario the maximally allowed reheating temperature can be constrained from the measured DM abundance. Since the abundance of thermally produced gravitinos is approximately inversely proportional to the gravitino mass [6,7,8], heavy gravitinos are favored from the requirement of large reheating temperatures whilst avoiding an over-closure of the universe. On the other hand, the gravitino mass governs the life-time of the NLSP. Since BBN bounds disfavor extremely large life-times, those bounds become more constraining for larger gravitino masses. This non-trivial interplay can be used to formulate upper bounds on the reheating temperature [5,8,17,18,19,20,21,22,23,24,25] on different levels of underlying assumptions.
In this work we consider a gravitino-stau scenario. We do not restrict ourselves to any constrained high-scale model but vary the SUSY parameters freely at the TeV-scale in the framework of the phenomenological Minimal Supersymmetric SM (pMSSM) [26]. Thereby we relax the particularly constraining [18,20,21] assumption of universal gaugino masses. Further, in this study we include the non-thermal production of gravitinos through the decay of the stau NLSP. This contribution can be very important for small mass differences between the stau and the gravitino and introduces a further dependence of the allowed values for the reheating temperature on the SUSY parameters. Consequently, low stau abundances are favored in two ways: by BBN constraints and by the desire for a small contribution of non-thermal gravitino production.
In [27] a survey for low stau abundances was performed in a Monte Carlo scan over a 17-dimensional pMSSM parameter space. In particular, the implications of a Higgs of around 125 GeV, constraints from direct SUSY searches, from MSSM Higgs searches, from flavor and precision observables and from charge or color breaking (CCB) minima on the phenomenological viability were highlighted. These results were obtained for a general super weakly interacting LSP. Here, we will specify the LSP to be the gravitino which allows us to apply constraints from cosmological observations and conclude on the allowed values for the reheating temperature. To this end we will extend the 17-dimensional parameter space introduced in [27] by the additional parameter of the gravitino mass. Requiring that the LSP abundance matches the measured DM density we will compute the corresponding reheating temperature by considering the thermal and non-thermal production of gravitinos. After computing the life-time and hadronic branching ratios of the stau we will utilize the BBN bounds presented in [28,29]. We will choose the conservative values for 6 Li/ 7 Li here. The analysis reveals the highest reheating temperatures that are consistent with bounds from BBN and other sensitive astrophysical observations, flavor and precision bounds, theoretical bounds from vacuum stability, bounds from direct SUSY searches at the 7 and 8 TeV LHC as well as bounds from the MSSM Higgs searches and the requirement of providing a Higgs around 125 GeV. Our analysis shows that points with large T R as required by leptogenesis only survive in a very particular corner of the SUSY parameter space. Those spectra feature a distinct signature at colliders [30] that could be looked at in the upcoming LHC run. In particular, it requires the triggering on very slowly moving heavy stable charged particles (HSCP) which is expected to be challenging in the high-luminosity run.
The paper is organized as follows. In section 2 we will review the relevant production mechanisms of gravitinos and discuss the underlying assumptions made for the non-thermal production in our setup. In section 3 we will describe the cosmological implications of a late decaying stau that are relevant for our analysis. The computational steps of the pMSSM parameter scan are introduced in section 4. In section 5 we present our results and discuss the implications for the upcoming high-energy LHC run. We will conclude in section 6.
Within this model the cold DM density has been measured with great precision [31]. Combining the Planck power spectrum data with the WMAP polarization measurements [32], BAO measurements [33,34,35,36,37] as well as ground based high multipole measurements performed by the Atacama Cosmology Telescope [38] and the South Pole Telescope [39] a best-fit value of was derived [31]. This value will be considered for the following analysis. There are two main production mechanisms for a gravitino which is not ultra-light and thus leads to a long-lived NLSP. On the one hand, this is the thermal production of gravitinos through inelastic scattering of particles participating in the thermal bath of the universe during the stage of reheating. On the other hand, it is the non-thermal production through decays of metastable supersymmetric particles into the gravitino. 1

Non-thermal production of gravitinos
In our setup the non-thermal production of gravitinos takes place via decays of the stau into the gravitino. Due to the assumed R-parity conservation each stau eventually decays into a gravitino. Hence, the number density of staus before their decay, n τ 1 , is equal to the number density of the gravitinos after all staus have decayed, n G , and thus However, this picture only remains true, if the decay of the stau takes place separated from the efficient annihilation of the staus into SM particles, i.e., if these annihilation processes do not compete with the decay. In order to quantify this requirement we consider the stau yield, Y = n τ 1 /s, where s is the entropy density. In figure 1 we show the evolution of the stau yield as a function of (decreasing) temperature T 0 and (increasing) time for a typical annihilation process 2 and for m τ 1 = 200 GeV and 2 TeV. We plot the relative deviation of the yield from its value for T 0 → 0 (if the stau were stable). This 1 Further sources of non-thermal production could arise from the decay of the inflation field. Since this contribution depends upon the actual model of inflation [40,41], we will not consider this contribution here. 2 We choose an annihilation process for which the thermally averaged annihilation cross section, σ eff v Møl , can be expanded in 1/x ≡ T /m τ 1 as where A is dimensionless, containing only numerical factors, mixing angles, couplings and mass ratios, see e.g. [42,43,44]. The first term in (3) often provides a good approximation [42]. The yield is then proportional to For a fixed x f = m τ 1 /T f , this expression uniquely determines the shape of the curves in figure 1 independent of the considered process.
Here, T f is the freeze-out temperature which is typically of the order T f ≃ m τ 1 /25 [43].
value is the quantity computed by micrOMEGAs [45] which will be used for our analysis. For cosmic times after 10 −4 sec the deviation is around or below one percent. Hence, for  Figure 1: Stau yield as a function of the temperature T0 = m τ 1 /x0 normalized to its value at T0 → 0 for the case of a typical annihilation process (see footnote 2).
We choose x f = 25 for this plot. The upper axis labeling denotes the corresponding cosmic time choosing g * (T ) according to the particle content of the SM [44]. By doing so we assume no additional relativistic degrees of freedom for temperatures T 10 GeV in our model.
significantly smaller life-times of the stau, decays take place while significant annihilation processes are still ongoing. With respect to the separated processes of annihilation and decay, this would lead to a higher gravitino abundance and would require incorporating the stau decay term in the Boltzmann equations. However, in this work we will focus on stau lifetimes larger than 10 −4 sec, first, because smaller life-times require gravitino masses which are far too small to achieve high reheating temperatures as desired for leptogenesis and thus are not of particular interest. Second, because BBN bounds that are subject to the investigation in this paper do not impose any restriction for lifetimes smaller than 10 −2 sec.

Thermal production of gravitinos
The relic abundance of thermally produced gravitinos, Ω th G , can be computed by solving the Boltzmann equation for the gravitino number density, where the collision term C G is determined by the thermal gravitino production rates. It has been computed to leading order in the involved gauge couplings considering the contribution from SUSY chromodynamics [7] and the full SM gauge group [8]. After the computation of C G , (5) can be solved analytically and yields [20] Ω th where g i and M i are the gauge coupling and the gaugino mass parameter, respectively, associated with the SM gauge groups U (1) Y , SU (2) L , SU (3) c and k i , ω i are corresponding numerical constants listed in table 1. The couplings and gaugino mass parameters are understood to be evaluated at the scale T R . Table 1: Assignments of the index i, the gauge coupling gi, the gaugino mass parameter Mi and the values of the associated constants ki and ωi to the SM gauge groups U (1)Y , SU (2)L, and SU (3)c. Taken from [20].

Implications of the stau decay
For a given MSSM parameter point all couplings of the gravitino to the MSSM particles are determined by the gravitino mass. We assume here that all heavier sparticles decay into the stau NLSP sufficiently fast so that direct decays of those sparticles into the gravitino are unimportant. The cosmological validity of a given parameter point then mainly depends on the yield, lifetime and the partial widths of the stau. For m τ 1 −m G > m τ the stau life-time, τ τ 1 , is dominated by the 2-body decay τ 1 → Gτ which can be computed from the relevant terms in the interaction lagrangian of a massive spin-3/2 gravitino [7,46,47], where ψ µ denotes the gravitino field and M Pl is the reduced Planck mass. For the general case of non-vanishing left-right mixing in the stau sector, τ 1 = cos θ τ τ L + sin θ τ τ R , we obtain the result The term proportional to sin 2θ τ (i.e., proportional to the amount of left-right mixing) can become significant for small mass splittings between the stau and the gravitino. It leads to a decrease or increase of the life-time depending on the sign of sin 2θ τ which corresponds to the sign of −X τ = −A τ + µ tan β (see, e.g., appendix B in [27]). This result reduces to the one given in [48] for the case of a purely left-or right-handed stau, θ τ = 0, π/2, π, and is analogous to the result found in [49] (published version) for the case of a stop NLSP. The scenario is subject to several bounds. The most important bounds come from BBN constraints. The particles that are emitted in the decay of the stau into the gravitino can induce hadronic and electromagnetic showers at cosmic times characterized by the life-time of the stau. The produced energetic hadrons and photons induce hadro-and photodissociation processes that potentially distort the predictions for the light element abundances of standard BBN [28,29,50,51,52]. Furthermore, staus may form bound states with the background nuclei potentially leading to a catalyzed overproduction of 6 Li [53,54]. For the application of the BBN bounds it is crucial to determine the hadronic branching fractions. The tau emitted in the 2-body decay of the stau, τ 1 → Gτ , has a hadronic branching fraction of roughly 65%. However, for cosmic times up to about 3 sec the interaction time of the tau is smaller than its life-time and the tau scatters off the background before decaying. This scattering leads to a purely electromagnetic energy release [55]. For later times the interaction time increases with decreasing temperature and hadronic decays of the tau become important. The mesons produced in the tau decays are unstable. In order to have a relevant effect on the BBN, the mesons have to scatter before their decay. This in turn only happens for cosmic times up to about 100 sec [50]. For later times BBN constraints are dominated by nucleons emitted in the stau decay. These nucleons stem mainly from the 4-body decays τ 1 → Gτ qq and τ 1 → Gν τ qq ′ with an invariant mass of the quark pair above the production threshold of the nucleon pair, m qq , m qq ′ 2 GeV [55].
If the life-time of the stau is very large, τ τ 1 10 12 sec, decays take place after the era of recombination and we can probe direct signatures of the stau decays in the measurements of the extra-galactic diffuse gamma ray background [56].
For even larger life-times much stronger bounds can be obtained from the searches for anomalously heavy hydrogen in deep sea water [57,58,59]. These measurements can be interpreted to provide a 95% C.L. bound on the yield of charged relics today, for the mass region m τ 1 ≤ 1600 GeV [59] and for the mass region 1600 GeV < m τ 1 ≤ 2000 GeV, where we chose an interpolated value between the ones given in [58] as an approximation. The limits translate into a maximal life-time, where t 0 is the age of the universe, t 0 = 4.354 × 10 17 sec [60], and Y is the stau yield before their decay. We will only consider parameter points that obey (11) in the following analysis. 3 Finally, we mention that one can also impose bounds on the life-time and abundance of late decaying particles from the observation of the CMB. The secondary particles produced in such a decay could affect the process of thermalization leading to a spectral distortion of the CMB away from a perfect black body spectrum [61,62,63,64,65]. However, the derivation and application of bounds from the CMB is beyond the scope of this work and is left for future investigations.

Scan over the 17-dimensional pMSSM
In this work we employ the Monte Carlo scan performed in [27]. In this subsection we will briefly summarize the computational steps and the constraints imposed on the parameter space. For further details we refer to [27]. We scanned over the 17-dimensional pMSSM parameter space with the following input parameters and scan ranges: For the third generation sfermions the spectrum parameters were chosen as input parameters. For simplicity we set m Q 1,2 = m u 1,2 = m d 1,2 . We imposed several hard constraints on the parameter space. The lighter stau was taken to be the NLSP, hence we only accepted points where Further, we required that at least one of the neutral CP -even Higgses, m h , m H , can be identified with the recently discovered Higgs boson at the LHC [66,67], We generated the sparticle spectrum with SuSpect 2.41 [68]. For the third generation sfermions we used tree-level relations in order to translate the chosen input parameters into soft parameters that feed into the spectrum generator. The Higgs sector was recalculated using FeynHiggs 2.9.2 [69]. We computed the stau yield with mi-crOMEGAs 2.4.5 [45].
We imposed several experimental and theoretical constraints on the parameter space. Lower bounds on the sparticle masses were derived from searches for heavy stable charged particles (HSCP) at the LHC. To this end and in order to discuss the perspective for a future discovery at the LHC, we determined all relevant cross sections for a center-ofmass energy of 7, 8 and 14 TeV. We computed the direct stau production via s-channel Higgses h, H with Whizard 2.1.1 [70]. The cross sections for all other contributions were estimated via a fast interpolation method using grids computed with Prospino 2.1 [71,72,73,74] as well as grids from the program package NLLfast [75,76,77,78]. The cross section upper limits were estimated from a reinterpretation of the HSCP searches for the 7 and 8 TeV runs reported by CMS [79]. For spectra with mass-degenerate staus and colored sparticles the respective R-hadron searches were taken into account. The decay widths and branching ratios were computed with SDecay [80,81] and Whizard 2.1.1 [70].
The point density was adjusted to the expected variation of the yield. In co-annihilation regions and regions around resonances or thresholds proportionally more points were accumulated (see [27] for details). We use a set of 10 6 pMSSM scan points 4 obeying the hard constraints (13) and (14).

Extension of the pMSSM parameter scan
We will now extend the 17-dimensional pMSSM scan described in [27] incorporating the gravitino LSP. For each point of the 17-dimensional pMSSM parameter space we perform the following computational steps. First, we determine the possible mass range for the gravitino under the following restrictions depending on the stau mass, the stau mixing angle and the yield of the given parameter point. The resulting life-time of the stau is required to be greater than 10 −4 sec-motivated by the arguments given in section 2.1-and smaller than the upper bound from (11). From (8) this imposes a lower and upper bound on the gravitino mass. Furthermore, the non-thermal contribution to the gravitino abundance (2) should not exceed the measured DM abundance (see below for further details). This requirement imposes an additional upper limit on the gravitino mass which can be both either more or less restrictive than the upper bound from (11). Second, for a given point we randomly generate 10 values for m G in the required interval. Since the interval spans over several orders of magnitude we use logarithmic priors here. The following steps are then performed for each of the 10 gravitino mass points.
We computed the non-thermal contribution to the gravitino abundance from the stau yield with (2). By demanding that the resulting total gravitino abundance matches the measured DM abundance, Ω non-th G h 2 + Ω th G h 2 = Ω CDM h 2 , we compute the required abundance of thermally produced gravitinos 5 , Ω th G h 2 . For Ω CDM h 2 , we chose the best-fit value (1). 6 From (6) we compute the reheating temperature, T R , that provides Ω th G h 2 for the given parameter point. Since M i and g i have to be evaluated at the scale T R , these quantities are functions of T R and the equation has to be solved iteratively. However, we achieved a fast convergence within 2 to 4 iterations to a more than sufficient accuracy. For the evaluation of g i and M i we take into account the one-loop running and the fact that see e.g. [95]. In (15), b i are the MSSM coefficients of the 1-loop renormalization group equations, (b 1 , b 2 , b 3 ) = (11, 1, −3) and Q in is the input scale, which we choose to be the electroweak scale here. 7 For the interpretation of BBN bounds and bounds from diffuse gamma ray observations we compute the life-time, (8), and the hadronic branching ratios, B h , of the stau. For τ τ 1 100 sec the relevant contributions to B h stem from 4-body decays, 5 Note, that the result (6) was obtained using hard thermal loop resummation [94] which requires weak couplings. Hence, the result might not be reliable for small reheating temperatures TR 10 6 GeV [20]. 6 The 68% confidence interval for the ΩCDMh 2 [31] is much smaller than the expected precision of the computations performed here. Therefore, we refrain from varying the ΩCDMh 2 within the confidence interval by a Monte Carlo method. The effect of such a treatment would be marginal. 7 We tolerate a slight overestimation of the couplings gi(TR) that could arise from the fact that the running with the MSSM coefficients starts below the precise mass scale of the corresponding SUSY particles. The effect on the final results is, however, expected to be marginal.
where Γ tot is the total width, which we approximate by the 2-body decay, Γ ( τ 1 → Gτ ) being the dominant decay mode. The partial widths Γ ( τ 1 → Gτ qq) and Γ ( τ 1 → Gν τ qq ′ ) include the decays into all kinematically accessible quark-anti-quark pairs. However, the contributions from diagrams containing top quarks in the final state are found to be negligible for all situations relevant here. We perform the computation of B h with the spin-3/2 extension of HELAS [96] implemented in MadGraph [97]. This program package supports the computation of arbitrary tree-level amplitudes with external gravitinos interacting with MSSM particles. In order to save computing time we determine the hadronic branching ratios in two steps on an increasing level of accuracy.
In the first step we conservatively estimate B h on the basis of a precomputed grid. To this end we computed B h as a function of the stau life-time for various choices of the stau masses and use an interpolation routine to obtain the values for arbitrary masses. For the computation of the grid we ignored left-right mixing effects and considered a purely right-handed stau taking into account diagrams with Z/γ-exchange only. Equally, we set the masses of all sparticles heavier than the stau to 3m τ 1 . This way diagrams involving EWinos (and squarks) are suppressed and do not contribute. Those diagrams can potentially increase the hadronic branching ratios. As an example, in the case of a right-handed stau with m τ 1 = 500 GeV and m G = 100 GeV we found a maximal enhancement of B h for almost mass-degenerate squarks of the first two generations and the bino-like neutralino, m q ≃ m χ 0 1 ≃ 510 GeV, by a factor of three. The branching ratios computed in this way are in rough agreement with results found earlier [55,98]. 8 In the second step, for each point that passes the bounds described in section 4.1 as well as the BBN bounds described further below (under the assumption of the conservatively estimated B h ) we recompute the hadronic branching ratios with MadGraph from the full spectrum. To this end we consider all diagrams of the processes τ 1 → Gτ qq and τ 1 → Gν τ qq ′ containing an intermediate vector boson, an intermediate light or heavy Higgs (for the process τ 1 → Gτ bb) as well as all diagrams containing an intermediate lightest neutralino or chargino. For a large fraction of scan points the contribution from τ 1 → Gν τ qq ′ -mediated via W ± -or χ ± -exchange-is found to be the most important. It can exceed the contribution from τ 1 → Gτ qq (q = d, u, s, c) by up to an order of magnitude. The contribution from τ 1 → Gτ bb is less important in our scan and we found Γ ( τ 1 → Gτ bb)/Γ ( τ 1 → Gτ qq) ≃ 3 at most, where q = d, u, s, c again. This contribution can potentially be enhanced from a Higgs exchange in the presence of large stau-Higgs couplings. As argued above for all computations we impose the lower cut on the invariant mass of the quark pairs m qq , m qq ′ > 2 GeV.
For life-times τ τ 1 100 sec the interactions of the mesons produced in the decays of the tau can become important. We estimate the corresponding hadronic branching ratio by using the results given in [98].
We apply the constraints from BBN derived in [28,29]. This analysis takes into account effects from proton-neutron interconversion, hadro-and photodissociation as well as all currently known bound-states effects. The constraints are based on the following observationally determined limits on the light element abundances: Here a conservative choice was made concerning the value of 6 Li/ 7 Li. As the BBN bounds derived in these references are given in terms of the life-time of the relic, its mass and its hadronic branching ratio, we do not compute the hadronic energy release nor simulate the hadronization of primary partons here. Rather we directly apply the computed values for τ τ , B h to the bounds given in [28,29]. These bounds are given for two masses of the relic m X = 100 GeV, 1 TeV and for (at least) six values for B h as a function of the life-time of the relic τ X . For life-times below 10 7 sec, where the hadronic energy release is important, the maximal yield which is compatible with the bounds, Y max , almost scales like B −1 h and m −1 X . Therefore we apply a linear interpolation (and extrapolation for masses above 1 TeV) in log(B h ) and log(m X ) between the corresponding values of Y max for a given lifetime. We take the bounds for 10 2 sec < τ X < 10 9 sec from [28] (erratum from 2009). As the bounds in [28] are only given for this interval, for life-times 10 −2 sec < τ X < 10 2 sec and 10 9 sec < τ X < 10 12 sec we estimate the constraints by using the results of [29], where we ignored the curves for B h > 0.01 in the latter interval. The constraints in this analysis are, however, derived for a neutral relic. As stated in [28], for large B h -typically achieved for very small life-times-the constraints on charged and neutral particles are almost identical. This is why we expect the analysis to apply for the former interval. For life-times in the latter interval, effects of photodissociation are the most relevant effects from the decaying staus. We expect the corresponding constraint to be similar to the bounds on decaying neutral relics for B h > 0.01, which is indeed the case for life-times 10 8 sec < τ X < 10 9 sec for which the constraints are given in both analyses.
For very large life-times τ τ 1 > 10 12 sec we consider bounds derived from the observation of diffuse gamma ray emissions [99]. We apply the relic density bounds for 2-body radiative decays derived in [56]. These bounds become restrictive only for life-times of τ τ 1 5 × 10 12 sec which corresponds to a mass splitting m τ 1 − m G 10 GeV in the considered scan region for m τ 1 . Consequently, the electromagnetic energy release in the stau decay is relatively small. We estimate the electromagnetic injection energy times photon branching ratio by where the pre-factor 0.3 conservatively takes into account the energy taken away by neutrinos emitted in the tau decays [100]. In the most relevant region 10 13 sec τ τ 1 10 15 sec the constraints on Y E inj B γ grow almost linear in E inj for small E inj , i.e., the displayed curves for E inj = 25 GeV, 50 GeV and 100 GeV are almost identical for these life-times. Assuming a linearity down to even smaller E inj , we apply the limits for the smallest value for the injection energy given, E inj = 25 GeV. 9

Results and discussion
The left panel of figure 2 shows the domains of the contributions to thermal gravitino production associated with the different gauge couplings. In blue, green and yellow we plotted points where the SU ( The right panel of figure 2 shows the ratio between the non-thermal and the thermal production of gravitinos. For small m G the non-thermal contribution is unimportant and the band spanned by the resulting reheating temperature grows linearly with the gravitino mass. Once the gravitino mass approaches the mass of the other superpartners we encounter two effects. First, according to (6), the linear growth of T R turns into a decrease when approaching small mass splittings between the gravitino and the gaugino masses. This effect causes the points with the highest T R to lie around gravitino masses of a few hundred GeV. The maximal T R reached by the generated points in our scan depends on the lower limits of the scan ranges for the gaugino masses, in particular for M 3 . 10 Here, having chosen M 3 > 1 TeV, it reaches T R ≃ 4 × 10 9 GeV in accordance with the conservative limits found in [25]. As a second effect, once the gravitino approaches the stau mass non-thermal contributions become important. Depending on the stau yield of a considered point the required reheating temperature is pushed down by a more or less significant amount. The points that still lie within the linearly rising band when m G approaches m τ 1 tend to be those with rather small yields. However, we found points with yields Y 10 −13 for T R 10 9 GeV. For these points the non-thermal contribution to the gravitino production is of the same order of magnitude as the thermal contribution and cannot be neglected.
In figure 3 we show the effect of the bounds imposed on the (17 + 1)-dimensional parameter space in the τ τ 1 -Y τ 1 plane and in the m G -T R plane. The blue and yellow points are rejected by the HSCP searches and by the additional bounds from flavor and precision observables, HiggsBounds and CCB bounds, respectively, as they have been described in section 4.1. The red points are rejected by the BBN bounds or the bounds from the diffuse gamma ray spectrum. The left panel of figure 3 reveals the effect of the BBN bounds on our parameter space. The border-line between the green and red points falls down relatively rapidly for life-times above 1000 sec according to the stronger bounds from hadrodissociation processes as well as bound-state effects. For life-times above 10 6 sec photodissociation processes become most restrictive. As a consequence we do not find allowed points with τ τ 1 > 10 7 sec in our scan. However, the point density starts to dilute for τ τ 1 > 10 7 sec as a consequence of our logarithmic prior in the scan over the gravitino mass (rather than over the stau life-time). Further, we do not encounter any point which is allowed by all other constraints but lies close to the bound on the yield imposed by the diffuse gamma ray spectrum. The spot of red points in the region Y 10 −12 and τ τ 1 10 2 sec stems from the energy release of mesons originating from tau decays.
Note that the BBN constraints from [28] show almost no dependents on the hadronic branching ratios for τ τ 1 > 10 5 sec and for the typically achieved hadronic branching ratios in this region that are well below B h = 10 −2 . Hence, the BBN constraints are not sensitive to the precise computation of B h in this region.
The right panel of figure 3 shows the parameter points in the m G -T R plane. The search for HSCP at the 7 and 8 TeV LHC imposes very restrictive limits on the gluino and wino masses, e.g., conservatively m g 1.2 TeV, M 2 800 GeV [27]. 11 These limits exclude all points with a reheating temperature above T R ≃ 2.3×10 9 GeV (cf. blue versus yellow points). Bounds from flavor and precision observables, MSSM Higgs searches and CCB vacua further reduce the parameter space leaving a maximal reheating temperature of slightly below 2 × 10 9 GeV (cf. yellow versus red points). The application of BBN bounds has the most significant effect in the region of large Ω non-th The analysis reveals the existence of points which provide reheating temperatures T R > 10 9 GeV and are consistent with all discussed bounds and with a Higgs mass of around 125 GeV. All these points share very distinct features. First, these points feature a heavy gravitino, 300 GeV < m G < 1.4 TeV, resulting in a relatively large stau life-time, 10 4 sec < τ τ 1 < 10 7 sec. It is interesting to note that the upper bound on the life-time (coming from BBN bounds) still causes a separation of the stau and gravitino masses of at least 200 GeV in our scan. Second, all points lie within the resonance region where m A ≃ 2m τ 1 . In this region exceptionally small stau yields can be achieved due to annihilation via a resonant s-channel heavy Higgs. For most points (88 points) the dominant annihilation process is resonant stau-pair annihilation [102]. 12 For three points effects of co-annihilation are important: we found that one and two points feature resonant stop and EWino co-annihilation [27] as the dominant annihilation process, respectively. Note that EWino co-annihilation via a resonant heavy Higgs requires no particularly large Higgs-sfermion couplings. Thus, the viability of these points does not depend upon constraints from CCB vacua.
Third, for most points the yield is smaller than 10 −14 . However, we encountered a few points with 10 −14 < Y < 3 × 10 −14 . In order to compensate for the slightly larger contribution of non-thermal gravitino production, those points were driven into a region of small gaugino masses and thus very small mass splittings between the stau and the gauginos, M 2 /m τ 1 < 1.2, M 1 /m τ 1 < 1.3 and M 3 /m τ 1 < 1.5. This strong tendency for small gaugino masses is in fact relaxed for Y < 10 −14 . Still, we found no points with M 2 > 2.1m τ 1 (cf. right panel of figure 4), M 1 > 3.1m τ 1 or M 3 > 3.7m τ 1 for T R > 10 9 GeV. The fact that (at the low scale) M 1 and M 3 are less constrained than M 2 is due to the smaller coupling in the former case and due to the slower running up to the scale T R in the latter case. The tendency for small stau-gaugino mass splittings is in fact the result of two effects. On the one hand, according to (6), the gravitino mass that maximizes the reheating temperature for a given Ω th G grows with increasing gaugino masses. On the other hand, the preference for smaller stau life-times from BBN bounds favors larger mass splittings between the stau and the gravitino. As a consequence the strong bounds on m g and M 2 also lift up the stau masses for points with T R > 10 9 GeV in our scan, which we found to lie above m τ 1 ≃ 800 GeV (see left panel of figure 4).
Finally, we want to comment on the prospects of studying these scenarios at the upcoming long-term run of the LHC. Figure 5 shows the full SUSY cross section of the points that have passed all bounds discussed above. The points that are closest to the exclusion limit from the HSCP search at 7 and 8 TeV typically provide a SUSY cross section at the 14 TeV LHC run of σ SUSY 14 TeV ≃ 100 fb, corresponding to the red points in figure 5. Since the cross section can have a strong dependence on sectors that are rather decoupled from the physics that constrain the reheating temperature-like the masses of the first generation squarks-the variation of the point color is relatively uncorrelated. However, we see that the uppermost stripe of the allowed band in the left panel does not contain points with very small cross sections due to the generically lighter gauginos for larger reheating temperatures. Many points in our scan with T R > 10 9 GeV provide cross sections around 1 fb or higher.
Since the points with T R > 10 9 GeV all feature the resonant configuration m A ≃ 2m τ 1 , at the LHC the direct stau production via a resonant heavy Higgs in the s-channel will be an important production mechanism [30]. For this process the production near threshold is significantly enhanced and the velocity distribution of the staus peaks at rather low values β 0.4 [30]. Such a signature is expected to be challenging for the current trigger settings at ATLAS and CMS and may require an extended buffering of the tracker data as pointed out in [104]. Further, providing rather slow staus, a noticeable amount of staus might be trapped inside the detector and eventually decay into the gravitino and a tau. Potentially this enables the determination of the stau life-time [105,106,107]. This is particularly interesting regarding the fact that a possible determination of the gravitino mass from the detection of the tau requires the tau energy, to deviate significantly from m τ 1 /2, i.e., m 2 G / ≪ m 2 τ 1 [108]. In the right panel of figure 5 we show the allowed points in the plane spanned by 1 − m 2 G /m 2 τ 1 and T R . Points with large T R tend to have values 1 − m 2 G /m 2 τ 1 that deviate significantly from one. Therefore, the prospects of testing supergravity by the simultaneous measurement of m τ 1 , τ τ 1 and m G [48,109]-allowing the verification of (8)-are significantly better in these scenarios, featuring large gravitino masses, than in scenarios with smaller gravitino masses and therefore smaller T R .

Conclusions
We worked out the interplay between constraints on the SUSY parameter space and the highest possible reheating temperatures in a gravitino-stau scenario. We performed a Monte Carlo scan over a (17 + 1)-dimensional parameter space. By demanding that the gravitino abundance matches the measured DM abundance we computed the required reheating temperature for each scan point taking into account the thermal and non-thermal production of gravitinos. Both quantities depend non-trivially on the MSSM spectrum parameters. We derived the cosmological viability from the application of bounds from BBN and the diffuse gamma ray spectrum. According to the strong constraints imposed for large stau life-times, τ τ 1 10 7 sec, from photodissociation processes causing an overproduction of 3 He, we do not encounter allowed points with stau life-times larger than 10 7 sec.
We found valid points with a reheating temperature high enough to allow for thermal leptogenesis, T R 10 9 GeV. These points are consistent with BBN bounds, flavor and precision bounds, theoretical bounds from vacuum stability, bounds from the HSCP searches at the 7 and 8 TeV LHC as well as bounds from the MSSM Higgs searches and the requirement of providing a Higgs around 125 GeV. All these points lie in the resonant region, m A ≃ 2m τ 1 . In this region annihilation dominantly takes place via the exchange of an s-channel heavy Higgs. For most of these points stau-pair annihilation is the dominant channel. However, we also found points where pair-annihilation of co-annihilating stops or EWinos is dominant. Most of the points with T R 10 9 GeV have exceptionally low stau yields 10 −16 < Y < 10 −14 . Further, the separation in the mass between the stau and the gauginos tends to be small especially for points with larger yields. This tendency is most pronounced for M 2 . This is due to the fact that the abundance of thermally produced gravitinos is approximately proportional to g 2 i M 2 i evaluated at the scale T R . Compared to M 2 the slower running of M 3 up to the scale T R over-compensates the effect of the larger coupling for the strong interaction.
For most of the points with T R > 10 9 GeV the dominant production mode at the 14 TeV LHC will be the production of EWinos or gluinos being relatively close in mass to the stau. However, due to the resonant configuration, m A ≃ 2m τ 1 , resonant stau production via the s-channel heavy Higgs will be an important contribution. This leads to the signature of extremely slowly moving heavy stable charged sparticles. For such a signature one would greatly benefit from an extended buffering of the tracker data in the LHC detectors increasing the trigger efficiencies for staus that arrive largely delayed in the muon chambers. Further, the signature can lead to a large amount of staus that are stopped in the detectors. This could provide the intriguing possibility of measuring the stau life-time. Moreover, especially for a heavy gravitino as required in order to obtain a high reheating temperature the determination of the gravitino mass might be possible from the measurement of the energy of the tau that is produced in the decay of the stopped stau. The combination of a variety of bounds on the low-scale SUSY parameters has pointed us to a very interesting corner in parameter space that should be looked at in the upcoming LHC run.