Neutralino and gravitino dark matter with low reheating temperature

We examine a scenario in which the reheating temperature $T_R$ after inflation is so low that it is comparable to, or lower than, the freeze out temperature of ordinary WIMPs. In this case the dark matter relic abundance is reduced, thus relaxing the impact of the usually strong constraint coming from the requirement that the universe does not overclose. We first re-examine the dynamics of freezeout during reheating. Next we apply a Bayesian approach to study the parameter space of the MSSM with ten free parameters, the CMSSM and the singlino-dominated regions of the NMSSM. In each case we find dramatic departures from the usually considered regime of high $T_R$, with important implications for direct detection dark matter searches. In the MSSM we examine WIMP mass range up to ~5 TeV, and find regions of bino dark matter over the whole mass range, and of higgsino dark matter with mass over a similar range but starting from the ~1 TeV value of the standard high $T_R$ scenario. We show that the prospects for bino detection strongly depend on $T_R$, while the higgsino is for the most part detectable by future one-tonne detectors. The wino, which is excluded in the standard scenario, becomes allowed again if its mass is roughly above 3.5 TeV, and can be detectable. In the CMSSM, the bino and higgsino mass ranges become more constrained although detection prospects remain similar. In the Next-to-MSSM at low enough $T_R$ wide ranges of singlino-dominated parameter space of the model become cosmologically allowed. We also study the contribution to the DM relic density from direct and cascade decays of the inflaton. Finally, we consider the case of a gravitino as dark matter. We find strong bounds from overclosure and Big Bang Nucleosynthesis, and derive lower limits on $T_R$ which depend on the gravitino mass and on the nature of the lightest ordinary superpartner.


Introduction
In spite of persistent efforts of both experimenters and theorists, the Standard Model (SM) still reigns supreme as a correct phenomenological description of almost all data in particle physics. However, the existence of dark matter (DM) offers one of a few empirical hints pointing beyond the SM and suggesting that it has to be incorporated into a more fundamental theory. A well-motivated example of such a theory is the Minimal Supersymmetric Standard Model (MSSM) (for a review see, e.g., [1]), which -unlike the SM -offers a candidate for a DM particle. The most commonly discussed case, the lightest neutralino, which is a mixture of the fermionic superpartners of the gauge and Higgs bosons, represents a weakly interacting massive particle (WIMP) and is stable if it is the lightest supersymmetric particle (LSP). Its relic abundance is determined at so-called freeze-out, when the annihilations become inefficient due to a decrease in its number density in the expanding Universe and the production processes are already ineffective due to a drop in the temperature of the primordial plasma (the Lee-Weinberg scenario). The abundance of two other well-motivated DM candidates, a gravitino -a fermionic partner of a graviton -and an axino -a fermionic partner of an axion -(see, e.g., recent review [2]), if they are the LSP, is generated by scatterings of the primordial plasma particles and from out-of-equilibrium decays of the lightest ordinary supersymmetric particles (LOSP) which had previously undergone freeze-out.
The questions of the origin and the properties of dark matter remain among of the main driving forces of both experimental and theoretical research in physics beyond SM. The latter activity includes both performing increasingly accurate calculations of the DM detection rates and relic abundance, including a critical reappraisal of the conditions in which this abundance was determined. The importance of this twofold approach becomes obvious by noting that the evolution of the Universe has been empirically tracked back to temperatures as high as O(MeV), but to obtain an estimate for the DM abundance one typically needs to make bold extrapolations to much higher temperatures.
It is usually assumed that the early Universe underwent a period of cosmological inflation during which an accelerated expansion of the Universe was driven by the vacuum energy density of a scalar field -an inflaton. After inflation the large potential energy of the inflaton field was transformed into the kinetic energy of newly produced particles in thermal and chemical equilibrium. As a result of this process, dubbed reheating, the Universe entered a radiation-dominated (RD) phase, and its initial temperature T R is commonly called the reheating temperature. 2 Another commonly adopted assumption is that the scale of T R is much higher than the mass scale of the MSSM particles, which allows one to separate the dynamics of reheating from that of DM freeze-out. Although this assumption is convenient, there is no a priori reason that it has to hold in the early Universe. Intriguingly, a recent study [3] has found that in the most popular models of large-field inflation T R may be required to lie within one or two orders of magnitude from the electroweak scale if the value of the spectral index is to remain very close to its observationally determined central value.
In this paper we will explore the possibility that the reheating temperature is comparable to the temperature of freeze-out, and will investigate the ensuing implications for DM phenomenology relative to the standard case. A number of analyses along these lines have been performed before: in a generic case [4], as well as in the context of the Constrained MSSM (CMSSM) [5] and of more general supersymmetric models [6,7,8]. However, the discovery of the Higgs boson with mass m h ≃ 126 GeV [9,10], together with negative results of the ATLAS and CMS searches for supersymmetric particles with masses below ∼ 1 TeV point towards the soft SUSY breaking mass scale M SUSY at least an order of magnitude larger than the M Z scale. These results imply a significant shift in the standard paradigm for supersymmetric dark matter. Previously, from naturalness-based assumption of M SUSY ∼ < O(1 TeV) it followed that bino-like neutralino was considered as the most natural and attractive candidate for the WIMP [11] in the MSSM, the choice which was also most naturally realized in unified models [12,13] for comparable ranges of M SUSY . However, with increasing values of M SUSY the relic abundance typically exceeds the observationally determined value of Ωh 2 = 0.1199 ± 0.0027 [14] by orders of magnitude already for DM mass of a few hundred GeV, unless special mechanisms of resonant annihilations or so-called coannihilations are employed [15,16]. On the other, in the case of a higgsino-dominated neutralino, coannihilations are very effective [17], and the relic abundance remains too low until its mass increases to ≃ 1 TeV, which, intriguingly, is the scale implied by LHC limits on M SUSY and also by the Higgs boson mass of ≃ 126 GeV. Interestingly, just such a higgsino-like WIMP emerges in unified SUSY with a TeV scale of M SUSY [18]. For wino-like neutralino the cosmologically favored mass range is even higher, ≃ 3 TeV [19].
Here we will show that the problem of DM overabundance can be alleviated at low reheating temperatures. Hence this assumption will lead to the opening up of previously cosmologically disallowed regions in the WIMP parameter space. In particular, the wino can again become experimentally allowed, a multi-TeV higgsino can have a correct relic abundance, while the relic abundance of the singlino can be reduced to an acceptable level, which in the standard case is hard to achieve. Prospects for WIMP direct detection can also be significantly affected.
The plan of the paper is as follows. In Section 2, we briefly review the dynamics of freeze-out in order to set the stage and to understand the impact of a low reheating temperature on a cosmological evolution and on the relic WIMP abundance. In Section 3, we employ the Bayesian approach to investigate the parameter space of the MSSM, the CMSSM and the Next-to-MSSM (NMSSM), and will identify the regions that are phenomenologically acceptable, including producing the correct relic density of the neutralino DM at low reheating temperature. In Section 4 we discuss and quantify the additional non-thermal contribution to the DM relic density from direct and cascade decays of the inflaton field to DM species and show that it can increase the DM relic density up to the measured value in otherwise underabundant scenarios. In Section 5, we extend the analysis to include the gravitino and assume it to be the DM, taking into account bounds from the Big Bang Nucleosynthesis (BBN) that inevitably arise in the presence of a long-lived LOSP. We conclude in Section 6.

Dynamics of freeze-out
In this section, we review the dynamics of freeze-out for high and low reheating temperatures.

High reheating temperature
An evaluation of freeze-out at high reheating temperatures has by now become a standard textbook lore (see, e.g., [20]). One assumes that the Universe was initially in the RD phase and that the energy density of radiation with g * (T ) effective degrees of freedom was given by ρ R = (π 2 /30) g * (T ) T 4 , with the temperature T inversely proportional to the scale factor a, i.e., T ∼ a −1 . For some stable particle species which are pair-annihilated into radiation in equilibrium processes, the Boltzmann equations governing ρ R and the number density n of some relic species read: where H is the Hubble parameter, σv is a thermally averaged annihilation cross-section times velocity for the species, E is its average energy (which we approximate as m 2 i + 9T 2 ) and n eq is its equilibrium number density.
In the context of supersymmetric theories with the LSP being a DM candidate this description should, in principle, be generalized by considering a separate Boltzmann equation for each supersymmetric particle species that is heavier than the LSP. Owing to R parity, these states pair-and co-annihilate and their decay chains all end up with the LSP. Fortunately, it was shown in [16] that in this case the evolution of the Universe can still be effectively described by a system of equations (1), if one replaces the number density of a single particle species by n = i n i , where the index i runs over all the particle species, each with a number density n i , and σv is replaced by where n eq,i stands for the equilibrium number density of i-th particle species, n eq = i n eq,i and σ ij v ij stands for a thermally averaged (co)annihilation rate for ith and jth particle species (for a detailed discussion see, e.g., [16,21]). The effective average energy released in the (co)annihilations of relic species is given by This approach is sufficient for accurate determination of the DM abundance, since, due to aforementioned chain decays, already before freeze-out n becomes the number density of the single stable species, the LSP. Having justified using the formalism of a single particle species in the case of frameworks with many states, like the MSSM, we can now briefly describe the dynamics of freeze-out. Eqs. (1) can be approximately solved under assumption that σv = ( where m χ is the WIMP mass. Before freeze-out WIMPs undergo (co)annihilations but are also produced in inverse processes and remain in thermal equilibrium. These processes are efficient until (co)annihilation rate remains larger than the expansion rate of the Universe, i.e., n eq σv > H ∼ T 2 /M Pl , where M Pl is the Planck mass. The freeze-out temperature T fo below which this relation is no longer satisfied marks the onset of an era where the DM number density changes only due to the expansion of the Universe. The present DM density calculated from (1) is therefore given by where T 0 is the present temperature of the Universe, Ω R h 2 is the radiation relic density and x fo = m χ /T fo satisfies where g is the number of degrees of freedom of the DM.
The parameter x fo depends very weakly on the details of DM interactions. Therefore, where the subscript "fo" corresponds to the value at T fo and we used (4). This approximation remains valid also for DM relic density when σv fo is replaced by σv eff,fo .

Low reheating temperature
If the reheating temperature is comparable to the freeze-out temperature, WIMPs may freeze out before the inflaton field has fully decayed, i.e., when the energy density of the Universe is still dominated by the energy density ρ φ of the inflaton. Therefore, the system of Boltzmann equations (1) has to be extended to accommodate the decaying inflaton field [4]. At the beginning of the reheating period the temperature of the Universe rapidly increases from T ≈ 0 to some maximum value T max due to the inflaton decaying to radiation. 3 At this temperature -though radiation is still being effectively produced in inflaton decays -the effect of the additional dilution caused by the increased expansion of the Universe begins to dominate and the temperature starts to decrease with an increasing scale factor, scaling as T ∼ a −3/8 . In other words, the same drop in the temperature corresponds to a faster expansion of the Universe during the reheating period than in the RD epoch. The set of Boltzmann equations now reads: where Γ φ is the inflaton decay rate. The faster expansion during reheating is driven by the entropy production due to inflaton decays and it continues until the inflaton decays completely. One conventionally associates the end of the reheating period with the reheating temperature T R defined as the temperature of the Universe assuming that the inflaton decayed instantaneously, 4 at the time corresponding to Γ φ = H, The reheating temperature T R is a priori unrelated to the freeze-out temperature T fo defined in Section 2.1. In Figure 1 we illustrate, in the context of the MSSM, the temperature dependence of yield Y defined as both for T R ≫ T fo (high T R scenario) and T R < ∼ T fo (low T R scenario). The solid (dotted) curve represents the low T R (high T R ) scenario and supersymmetric mass spectra have been selected in such a way that both number densities reach their equilibrium values. Due to a faster expansion of the Universe, for low T R the freeze-out 3 The value of Tmax does not play a role in the determination of the DM relic abundance, since Ω DM h 2 is set mainly by the the rate of (co)annihilation processes near freeze-out. Other possible sources of DM are direct and cascade decays of the inflaton field and inelastic scatterings of the inflaton decay products [22,23]. We shall mention then only briefly at the end of our study, as they are model dependent and, moreover, in scenarios considered here the freeze-out temperature is very close to the reheating temperature, which, in principle, allows thermalization of DM. A recent discussion of these issues can be also found in [24]. 4 In the reheating scenarios considered here, at T R given by eq. (9) the Universe is typically still dominated by the inflaton field [4] and the radiation-dominated epoch actually starts at a somewhat lower temperature. occurs slightly earlier, with typical x fo = 10 − 25, than for high T R where it typically lies between 20 and 25. If the decay of the inflaton stopped at T fo , the DM abundance would be higher in the low T R scenario. However, a continuous entropy production keeps diluting it until the reheating temperature is reached. The end result is an overall reduction, 5 of the DM abundance relative to high T R scenarios [4]. Assuming again (4), an approximate DM abundance resulting from the set of Boltzmann equations (8) reads [4] Finally we obtain where, similarly to (7), the subscript "fo" corresponds to the value at T fo given by (12), which is slightly larger than the value of the freeze-out temperature obtained in the high T R scenario. Of course, in a full MSSM calculation one has to replace σv fo with σv eff,fo given by (2).

A comparison of the scenarios with a high and a low reheating temperatures
As shown in eqs. (7) and (13), the DM relic abundance in scenarios with high and low T R is determined by the value of σv eff at the respective freeze-out temperatures. Since the freeze-out temperatures are very similar in both cases, the following approximate relation holds: with (T fo /m χ ) 3 factored out since its value changes only in a narrow range. From (14) it immediately follows that in scenarios with low reheating temperatures, T R < T fo , the DM relic abundance is suppressed with respect to scenarios with high reheating temperatures. Since the latter case has been extensively studied and the DM relic density can be easily calculated for a given WIMP type and mass, it is useful to rephrase (14) in the following way. If Ω DM h 2 is fixed at the observed value of 0.12, a phenomenologically acceptable scenario is the one where the standard prediction for Ω DM h 2 (high T R ) is larger than the observed value by a factor of (m χ /T R ) 3 (T fo /m χ ) 3 . In other words, SUSY configurations which would be otherwise rejected as giving too large relic density become acceptable at low reheating temperatures. We shall explore this effect in Section 3 when scanning a parameter space of some specific SUSY models below. Although in practice eq. (14) is very useful for understanding the T R -dependence of Ω DM h 2 , it may also be slightly misleading, as it does not show a certain degree of correlation between T fo and Ω DM h 2 (high T R ). This correlation is easy to understand, since a large Ω DM h 2 (high T R ) results from a low (co)annihilation crosssection which, according to eqs. (6) and (12), drives T fo to higher values. An account of this effect is shown in Figure 2, which shows the relation between Ω DM h 2 (high T R ) and the true relic density Ω DM h 2 at some low T R for different values of m χ /T R . Obviously, in the high T R limit Ω DM h 2 approaches Ω DM h 2 (high T R ), while for values of m χ /T R of 20 and more we observe a stronger T fo dependence, as predicted by (14), which results in a slower increase of Ω DM h 2 with growing Ω DM h 2 (high T R ) and fixed m χ /T R . Of course, if the LOSP is the DM candidate, the phenomenologically relevant values of Ω DM h 2 belong to a narrow observed range. However, we shall see in Section 5 that for gravitino DM produced in LOSP decays even larger values of the LOSP relic density will become allowed.

Neutralino dark matter with low reheating temperatures
We will now apply the formalism presented in Section 2 to the MSSM with ten free parameters, to the CMSSM, and to the Next-to-Minimal Supersymmetric Standard Model (NMSSM) with a singlino-dominated DM.

The MSSM
In this subsection we will analyze the scenario with low reheating temperatures of the Universe in the context of the MSSM. Since a study of a completely general MSSM would be unmanageable, nor for that matter even necessary, we select a 10-parameter subset of the MSSM (p10MSSM) which exhibits all the features of the general model which are relevant for our discussion. The free parameters of the model and their ranges are given in Table 1. Our choice follows that of [26] (see discussion therein), except that we keep both the wino mass M 2 and the bino mass M 1 free in order to allow each of them to be DM. As we will see, the choice of ten free parameters will allow various accidental mass degeneracies which can contribute to coannihilations. Also, the ranges of parameters have been extended to obtain a wide range of Ω DM h 2 (high T R ) with m DM reaching up to 5 TeV.

Parameter
Range bino mass 0.  We scan the parameter space of p10MSSM following the Bayesian approach. The numerical analysis was performed using the BayesFITS package which engages Multinest [27] for sampling the parameter space of the model. Supersymmetric mass spectra were calculated with SOFTSUSY-3.4.0 [28], while B-physics related quantities with SuperIso v3.3 [29]. MicrOMEGAs v3.6.7 [30] was used to obtain Ω DM h 2 (high T R ) and DMproton spin-independent direct detection cross section σ SI p . The constraints imposed in scans are listed in Table 2. The LHC limits for supersymmetric particle masses were implemented following the methodology described in [26,37] [36] have been implemented as a hard cut.
calculated by solving numerically the set of Boltzmann equations (8), as outlined in [4]. In order to find the point where WIMPs freeze-out we adapted the method described, e.g., in [38] to the scenario with a low reheating temperature. Both σv eff and σv eff E eff as a function of temperature were obtained with appropriately modified MicrOMEGAs; we also checked that in the high T R limit we reproduced Ω DM h 2 obtained with the original version of this code.
The results of the scans -but without imposing the constraint on the DM relic abundance and direct detection rates -are shown in Figure 3, with lines of constant T R = 1, 10, 50, 100, 200 GeV superimposed along which Ω DM h 2 ≃ 0.12. The horizontal line corresponds to the correct DM relic density in the standard high T R scenario. Different colours denotes various compositions of the lightest neutralino: green, red and blue corresponds to the bino, higgsino and wino fraction larger than 95%. We will now describe the results for each of these three cases.
Bino DM. The region of bino DM covers most of the plane in Figure 3. In this case the relic density can vary by several orders of magnitude for a given m DM , since it is very sensitive to the details of the MSSM spectrum. Generically, bino annihilation rate is dominated by t-channel slepton exchange χχ → ll and for m B ≪ m l the bino relic density reads (see, e.g., [39,40]) where g * fo stands for the number of relativistic degrees of freedom at χ decoupling. By varying the bino and the slepton masses, one can obtain Ω B h 2 (high T R ) spanning a few orders of magnitude. The upper boundary of the allowed region in Figure 3 has no physical meaning -it simply corresponds to the maximum value of slepton masses in our scan which is ∼ 10 − 15 TeV.
It is well-known that the correct Ω B h 2 (high T R ) can be achieved for low m B typically thanks to coannihilations with the lighter stau or, for m A ≃ 2m B , to resonant annihilations through the s-channel exchange of the pseudoscalar Higgs boson (A-funnel region); however, Ω B h 2 (high T R ) ∼ 0.12 can also be obtained with the lighter Higgs boson h resonance [41], for bino-higgsino mixing or in the bulk region where the bino-dominated neutralino annihilates through t-channel exchange of sfermions (typically of sleptons as they are usually lighter than squarks -see a discussion of these regions in, e.g., [26]). Note that, for large bino mass, m B > 2 TeV, its relic density can still be reduced by the Higgs pseudoscalar exchange in the A-funnel region, but also through coannihilations owing to accidental bino-wino or bino-gluino mass degeneracies. This explains the presence of points with a very low bino relic density at large WIMP mass in Figure 3.
Higgsino DM. The results for the higgsino DM relic density agree well with other recent analyses (see, e.g., [26]). In Figure 3, Ω DM h 2 (high T R ) scales proportionally to m 2 DM , achieving the correct value at m DM ∼ 1 TeV. However, one can see that, for the whole range above that value one can obtain the observed value of the relic density provided T R is low enough, around 100 GeV.  Wino DM. Wino relic density is quite sensitive to a so-called Sommerfeld enhancement (SE) of the annihilation cross-section due to attractive Yukawa potentials induced by the electroweak gauge bosons [42] (see also, e.g., [43] for a recent and exhaustive discussion; we use enhancement factors from that reference in our numerical analysis). Incidentally, the SE is particularly important in the ∼ 2 − 3 TeV wino mass range, where the correct Ω W h 2 can be obtained for high T R . In our scan, the SE is responsible for a visible vertical broadening of the wino region around 2.5 TeV. When considering the wino as a DM candidate, one has to take into account that the SE is associated with enhanced rates of present-day wino annihilations giving rise to diffuse gamma ray background; therefore, stringent indirect detection bounds apply in this case. It has been shown [44,45,46] that the enhancement of indirect detection rates for m W 3.5 TeV is in conflict with current observational limits. On the other hand, wino DM with mass larger than 3.5 TeV generically has too large relic abundance, which excludes it as a DM candidate over the whole mass range in the standard high T R scenario.
For each of the three neutralino compositions discussed above, a suppression of the DM relic abundance at low T R leads to interesting, and often dramatic, consequences, allowing vast regions of the parameters space regarded as phenomenologically disallowed in the high T R limit. In the following we shall present a more detailed analysis of the parameter space of the MSSM with low T R .
Scenarios with a low reheating temperature allow choices of the MSSM parameters which at high T R would lead to too small DM annihilation rates and, as a consequence, too large relic density. Since small annihilation rates are usually associated with small direct detection rates, it is interesting to analyze the effect of the assumed low reheating temperature. We shall discuss here both the most recent constraints from the LUX experiment [36], as well as from expected future reach of the one-tonne extension of the Xenon experiment (Xenon1T) [47].
In Figure 4 we show -for fixed values of T R -the 2σ credible regions in the (m χ , σ SI p ) plane for the p10MSSM scans with the DM density constraint included. In the case of high reheating temperature (upper left panel) most points correspond to m χ 1.5 TeV: these are either bino-or higgsino-like neutralinos. Scenarios in which the neutralino is the bino with a few per cent higgsino admixture are typically characterised by enhanced σ SI p ; such points occupy the upper part of the bino DM (green) region and will be accessible to Xenon1T. An almost pure bino neutralino, instead, can have much lower direct detection cross-section and it often remains beyond the reach of current and future experiments. In the case of higgsino DM, a good fraction of points lie within the projected Xenon1T sensitivity. As we have discussed in Section 3.1, for higher m χ one needs specific mass patterns to obtain the correct relic density; as these are fine-tuned cases, one obtains fewer points for m χ 1.5 TeV than for lower DM mass values. The wino, which can have the correct relic density for m W ∼ 2 − 3 TeV, is not shown in the plot, since it is excluded by the indirect DM searches in this mass range [44,45,46].
As expected from Figure  virtually the same as for high T R . However, an important difference appears at m χ ∼ 3 − 4 TeV where one can obtain the desired value Ω χ h 2 ≃ 0.12 for the higgsino. In this region, the direct detection cross section σ SI p is high enough to allow testing the scenario by the Xenon1T experiment. We also note that, though Figure 3 suggests that for T R = 100 GeV one can have a higgsino-like DM with any mass in the scanned range, higgsino mass between 2 TeV and 2.5 TeV are disfavored because of Ω χ h 2 being often too large. As a result one observes a reduced number of higgsino-like points in this mass range.
For T R = 50 GeV (lower left panel), the low-T R relic density suppression is already effective for m χ ∼ 1 TeV and it is very strong for larger DM mass, making the higgsino strongly disfavoured. We can see just a few ∼ 1 TeV higgsino-like neutralinos characterised by Ω χ h 2 (high T R ) ∼ 0.2 − 0.4. On the other hand, for m χ > 1 TeV one can now easily obtain the correct relic density for a nearly pure bino without requiring any specific relation among soft supersymmetry (SUSY) breaking parameters. The region with m χ < 1 TeV now becomes less appealing, since it still requires some specific mass pattern to suppress the relic density, and we find only a few points there. As can be seen in Figure 4, only a fraction of the 2σ credible region lies above the Xenon1T expected reach in the range of ∼ 2 − 3 TeV mass.
For T R = 10 GeV (lower right panel), only points corresponding to m χ < 1.5 TeV are present in our scan. This feature does not have a physical origin, but it merely results from a finite, albeit generous, ranges of the superpartner masses which we have allowed; this limit can be seen in Figure 3. Since low-T R suppression is now very effective in the entire DM mass range, these points typically have large Ω DM h 2 (high T R ), hence low σ SI p and the experimental verification of such scenarios poses a challenge.
With the values of T R discussed so far we have not seen any acceptable points corresponding to wino DM. This can be easily understood by examining Figure 3 which shows that the wino DM with m χ > ∼ 3.5 TeV has the correct relic density for T R only between 100 and 200 GeV. The left panel of Figure 5 shows the reheating temperature for the points in the 2σ credible region in the p10MSSM corresponding to the wino with the correct Parameter Range common scalar mass 0.1 < m 0 < 10 common gaugino mass 0.1 < m 1/2 < 10 common trilinear coupling −15 < A 0 < 15 ratio of Higgs doublet VEVs 2 < tan β < 62 sign of µ parameter µ > 0 Table 3: The parameters of the CMSSM and their ranges used in our scan. All masses and trilinear couplings are given in TeV, unless indicated otherwise. Masses and trilinear coupling are given at the GUT scale. The nuisance parameters are the same as for the p10MSSM.
abundance -with and without the SE taken into account. With the SE neglected, the points form a narrow band with T R between 115 and 120 GeV. Since the SE leads to a suppression of Ω χ h 2 (high T R ), its inclusion allows one to obtain the measured DM relic density for slightly larger T R . The actual enhancement of the cross-section depends on the value of µ and can therefore vary for a given wino mass. Hence, including the SE one obtains Ω W h 2 ≃ 0.12 for a wider range of reheating temperatures 120 GeV T R 200 GeV. In the right panel of Figure 5 we show -for T R = 150 GeV -the 2σ credible region of the p10MSSM on (m χ , σ SI p ) plane. Regions with lower m χ corresponding to the bino or the higgsino are similar to the high T R case as expected. At larger m χ , a new region with the wino DM becomes allowed for mass of 3.5 TeV. It is not excluded by current limits from indirect detection experiments, but potentially can be in the future [46]. Although some of these points lie within the projected Xenon1T sensitivity reach, direct detection experiments will not constrain the scenario too strongly.
One may wonder whether the lighter Higgs boson mass, m h ≈ 126 GeV, constrains the low-T R bino DM scenarios in any significant way. The answer is negative: one obtains a sufficiently large m h by arranging large stop masses and/or a large left-right mixing in the stop sector, while the bino relic density depends mainly on bino and stau masses. Since in the p10MSSM discussed in this section the stop and bino/stau sectors are to a large degree independent, for all the points presented in Figure 4 the lighter Higgs boson mass comes out close to the experimentally measured value thanks to heavy squarks, well above the LHC limits for colored superpartners.

The CMSSM
We will now examine which features, if any, of the general MSSM with low reheating temperature will remain when we relate its many free parameters by the assumption of a grand unification. A prime example of this class of models is the Constrained MSSM [13], where unification conditions are imposed at the GUT scale. The parameters of CMSSM and their ranges are given in Table 3.
In Figure 6 we show 2σ credible regions of the (m 0 ,m 1/2 ) plane with high (left panel) and low T R = 10 GeV (right panel). In the high-T R scenario one can identify three well-known regions with low χ 2 (see, e.g., [37]) that correspond to the correct relic density of neutralino DM: from left to right, the stau coannihilation and the A-funnel regions, as well as the ∼ 1 TeV higgsino region. The focus-point/hyperbolic branch region is absent, since it has been excluded by the LUX limit on DM direct detection cross section for positive µ.
As we have seen in Figure 3, for T R = 10 GeV only the bino can produce the correct relic density. The lower left corner of the allowed region in (m 0 ,m 1/2 ) plane corresponds to stau coannihilation region, analogous to that obtained for high T R , (for such low WIMP mass values the suppression due to low T R is inefficient). For slightly higher values of the mass parameters, the suppression of the relic density by stau coannihilations is traded for low-T R suppression and we find acceptable points there. In that region, the bino relic density for a fixed T R and a fixed bino mass (or m 1/2 ) depends on many factors, in particular, on stau masses (which depend not only on m 0 , but also on tan β and A 0 ), as well as on the small but non-negligible higgsino fraction of the lightest neutralino.
Unlike in the general MSSM, in both the high-and low-T R regime the the Higgs boson mass and the DM  relic density depend in part on the same parameters of the model, so they are not completely independent. This is illustrated by the case of T R = 10 GeV. It is known that unless stop masses are in the few-TeV range, the Higgs boson mass must receive sizable contributions from large left-right mixing in the stop sector, possible for large tan β and/or large |A 0 |. However, in the CMSSM a large left-right mixing in the stop sector leads to a substantial left-right mixing in the stau sector, which in turn leads to a suppression of the mass of the lighter stau. For m 0 ∼ 2 − 3 TeV this results in the constraints tan β < 20 and A 0 < −5 TeV. Finally, for m 0 of a few TeV, the staus are so heavy that varying tan β or |A 0 | is not dangerous for the DM relic density, so the Higgs Parameter Range SH u H d coupling 0.001 < λ < 0.7 scalar cubic coupling 0.001 < κ < 0.7 soft scalar A-term −12 TeV < A κ < 12 TeV Table 4: Additional parameters in the p13NMSSM and their ranges given at the SUSY scale. The nuisance parameters are the same as for the p10MSSM.
boson mass measurement tends to push tan β to higher values. In Figure 7 the spin-independent direct detection cross section σ SI p is shown as a function of the neutralino mass for both the high T R scenario and for T R = 10 GeV. As it is already known [37], a significant part of the 2σ credible region in the high T R scenario can be tested in future one tonne experiments. On the other hand, in the T R = 10 GeV case prospects for DM discovery are much worse. Only a small fraction of the allowed region can be covered by Xenon1T; it is characterised by high m 0 , low m 1/2 and low |A 0 |, where, according to [48], µ can be suppressed by the negative m 2 Hu (SUSY) tending closer to zero, In this case the higgsino fraction of the bino-dominated DM goes up to even 5%.

The NMSSM
The superpotential of the MSSM contains a mass term µĤ uĤd with a mass parameter µ of the order of the soft SUSY breaking parameters. One therefore needs an explanation why µ should be much smaller than the other scales in the unbroken SUSY theory, such as the unification scale or the Planck scale. A simple and elegant solution to this 'µ-problem' consists in replacing a mass term with a Yukawa-like interaction between the chiral superfieldsĤ u ,Ĥ d and a chiral superfieldŜ which is a singlet of the SM gauge. Its scalar component can acquire a VEV, thereby generating an effective µ term (see [49] for a review). This framework is called the Next-to-Minimal Supersymmetric Standard Model (NMSSM). The fermionic component of the singlet multiplet, the singlino, carries no SU (3) or electric charges, so it can mix with the other four neutralinos. It is therefore possible that in the NMSSM a state which is mostly singlino-like is the lightest of the neutral, non-SM, R-parity-protected fermions and therefore a DM candidate.
The parameter space of the NMSSM contains three parameters absent in the MSSM. They come from new terms in the superpotential λŜĤ uĤd +1/3! κŜ 3 and from soft SUSY breaking potential 1/3! A κ S 3 (the coefficient A λ in the term A λ H u H d S of the soft SUSY breaking potential is then determined in terms of other parameters, including µ and m A ). We therefore extend the numerical analysis described in Section 3.1 to accommodate these three additional parameters; their ranges are given in Table 4. The spectrum and the decay widths are calculated with NMSSMTools 4.2.0 [50] and the high-T R relic density is obtained from an appropriately extended micrOMEGAs code [51]. Requiring perturbativity to hold up to the GUT scale requires λ, κ 0.7; this justifies our choice for the upper limit in the scan, but in practice there are always additional constraints. The condition that the singlino is lighter than the higgsino implies that κ λ/2 < 0.35 and the requirement that the DM is made up of an almost pure singlino (a parameter region which we focus on here) introduces an effective upper limit λ 0.1 for majority of points in the scan -this suppresses the respective off-diagonal entries in the neutralino mass matrix. As a result one typically finds κ < 0.05. The effective upper limit on A κ comes from positivity of the pseudoscalar mass matrix and it reads A κ 0.
The results of the scan projected onto the (m DM , Ω DM h 2 (high T R )) plane are shown in Figure 8. The range of the high-T R singlino relic density spans a few orders of magnitude, from 10 −2 to 10 7 . The largest values are ∼ 4 orders of magnitude larger than the largest values that we obtained for the bino LSP. This can be explained by the fact that a nearly pure singlino interacts very weakly; it annihilates mainly into scalar-pseudoscalar pairs (mainly H 2 A 1 ) with the associated couplings proportional to κ or λ. This dominant annihilation channel is characteristic of scan points with Ω DM h 2 (high T R ) > 10 5 . Smaller values of Ω DM h 2 (high T R ) require at least  . Solid black horizontal line corresponds to high T R limit. Shown scan points correspond to credibility levels of 95%; dark (light) brown triangles correspond singlino fraction > 99% (between 95% and 99%).
a partial mass degeneracy between the singlino and a heavier particle thus allowing coannihilations: these are mainly coannihilations with the bino for 10 3 < Ω DM h 2 (high T R ) < 10 5 and coannihilations with the higgsino for Ω DM h 2 (high T R ) ∼ 10 2 . We also find points with smaller values of Ω DM h 2 (high T R ), but they necessarily involve special mass patterns which permit coannihilations with either higgsino, wino, stau/sneutrino, stop or gluino.
The lines of constant Ω DM h 2 = 0.12 for different values of the reheating temperature T R shown in Figure 8 are the same as in Figure 3 and we arrive at a conclusion analogous to that of Section 3.1, namely that there exist vast regions of the parameter space of the NMSSM with an almost pure singlino DM which have been so far disregarded solely because of predicting too large a relic density; however, with sufficiently low T R the relic density can be suppressed enough to agree with the measured value and these regions become phenomenologically viable.

Direct and/or cascade decays of the inflaton field
We have so far made an implicit assumption that the inflaton field φ is very heavy and, therefore, that the direct and cascade decays of φ to DM species are negligible. It is, however, important to study the validity of this assumption for a range of inflaton mass, as inflaton decays can give an additional, non-thermal contribution to Ω χ h 2 . Our analysis follows here the model-independent approach used in [6,7].
Direct and cascade decays of the inflaton field to superpartners of SM particles correspond to an additional term in the Boltzmann equation (8) for n, which is now given by, 6 In the wino DM case we take indirect detection limits following [46]. For the reheating temperatures above thin dashed black lines the freeze-out of the DM particles occurs after the reheating period (i.e. in the RD epoch). The limit at ∼ 800 GeV comes from antiprotons and the one around 1.8 TeV from the absence of a γ-ray line feature towards the Galactic Center.
where b describes the average number of DM particles produced per inflaton decay described by the decay constant Γ φ and ρ φ denotes the inflaton energy density. We present our results in Figure 9 in the (m χ , T R ) plane in terms of the dimensionless quantity η = b · (100 TeV/m φ ) for higgsino (left panel) and wino (right panel) DM.
The relic density of DM in this case is a sum of the thermal and the non-thermal components. The thermal production with a low reheating temperature has been studied in Section 2 and shown to be an increasing function of T R . On the other hand, the magnitude of the non-thermal component may depend, for fixed η and m χ , on the reheating temperature in a non-monotonic way, as discussed in detail in [7]. When T R is sufficiently low, non-thermal production leads to Ω χ ∼ T R , while for larger reheating temperature DM relic density goes down with increasing T R . As a consequence, each curve corresponding to fixed relic density Ω χ h 2 = 0.12 and fixed η in Figure 9 is C-shaped. For the upper branch of each curve, corresponding to larger values of T R , the correct relic density is obtained for such values of m χ that freeze-out occurs only slightly earlier than the end of the reheating period. 7 As m χ increases required values of the T R become larger and finally reach the level at which freeze-out occurs after the reheating period, i.e., in the RD epoch, and therefore direct and cascade decays of the inflaton field play no role in determining Ω χ .
The additional, non-thermal contribution to the DM relic abundance can help reconcile with the measured value these regions of the MSSM parameter space for which Ω χ h 2 is otherwise too low even at high T R . Examples of such cases include the higgsino with mass below 1 TeV or wino with mass below 2 TeV, shown in Figure 9. For sufficiently large values of η, one can even generate too much DM from inflaton decays; this upper bound on η can be translated into a lower bound on the inflaton mass for which the direct production is negligible even for a branching ratio BR(φ → superpartners) ∼ O(1). In particular, for η < 10 −9 we obtain no significant non-thermal production of DM particles. This value corresponds to the inflaton mass of m φ = b · 10 13 GeV, which for typical values of b ∼ O(10 3 ) [52], points towards inflaton mass close to the unification scale. 8 5 Gravitino dark matter with low reheating temperature Many of the considerations presented in Section 3 can be applied to another theoretically motivated scenario where the DM is made up of the gravitino G, the supersymmetric partner of the graviton, assuming that G is lighter than all the superpartners of the SM particles.
Unlike the neutralino, for a sufficiently large mass the gravitino is not a thermal relic. Its abundance Ω G h 2 receives contributions from at least two sources: the thermal component Ω TP G h 2 is produced in scatterings and decays in the thermal plasma [53,54,55], while the nonthermal component Ω NTP G h 2 results from late decays of quasi-stable relic LOSPs after they freeze out [56,57]. Since Ω TP G h 2 is proportional to T R , for T R ≪ 10 6 GeV and m G > ∼ 1 GeV this component is much smaller than the measured value of the relic density, hence at low T R it is the nonthermal component of gravitino DM that is dominant, and the gravitino abundance can be related to the LOSP abundance by Long after they have frozen out, during or after BBN, the LOSPs decay into gravitinos and SM particles, thus initiating hadronic and electromagnetic cascades which can affect light element abundances (see e.g. [58,59]) and potentially lead to a violation of current observational limits.
Here we analyze the viability of the gravitino DM with low reheating temperatures, making use of the results of the scan described in Section 3.1 with an additional assumption that the gravitino is lighter than any of the superpartners of the SM particles and without requiring that the LOSP is neutral which allows the LOSP to be a neutralino (bino, wino or higgsino) or a slepton (a charged slepton or, with large enough splitting between right and left soft stau masses [60], a sneutrino). We follow Ref. [58] for the implementation of BBN constraints, which mainly depend on the LOSP mass m LOSP and abundance Ω LOSP h 2 , as well as on the LOSP hadronic branching ratio B h . We calculate Ω LOSP h 2 as described in Section 3.1 and for B h we use existing results for neutralinos [61], sneutrinos [62] and charged sleptons [63].
Typical results for m G = 10 GeV and 1 TeV obtained in the p10MSSM are given in Figure 10. We fix the gravitino abundance at the observed value, relate it to the LOSP relic density through (18) and then find the corresponding reheating temperature with the procedure described in Section 3. Similarly as in Section 3.1, we present the results in the (m LOSP , Ω LOSP h 2 (high T R )) plane. As one could expect from eq. (18), the line corresponding to the correct gravitino DM abundance in high-T R case is not horizontal, as it was the case for neutralino DM. Below the line the gravitino abundance is lower than the observed value and, in the absence of thermally produced component, 9 such points are not viable. We note that for the sneutrino LOSP it is mass degenerate with the lighter (left) stau, thus coannihilations do play an important role here. This typically makes Ω LOSP h 2 (high T R ) smaller for the sneutrino LOSP than for the (usually right) stau LOSP.
In the low T R regime, as long as m G 100 GeV, the bino as the LOSP is the only possibility for gravitino DM. In this case, however, Ω LOSP h 2 (high T R ) typically exceeds unity and B h ∼ 1; hence, in order to avoid bounds from the BBN one can simply require the LOSP lifetime to be < ∼ 0.1 s, which leads to [61] m LOSP > ∼ 1400 which is consistent with the results shown in the left panel of Figure 10. The interpretation of this bound is very simple: the LOSP number density is so large that the particle must decay before BBN in order not to affect its successful predictions; because of the lifetime-mass dependence, this places a stringent lower bound on the LOSP mass. While at low T R one can suppress the LOSP number density and alleviate BBN constraints, with a small m G /m LOSP in (18) this would lead to too small gravitino abundance. On the other hand, it follows from Figure 10 that a lower bound on m LOSP can be translated into a lower bound on T R . We show such bounds in Figure 11 as a function of the gravitino mass with and without efficient direct and cascade decays of the inflaton field to bino. As we argued in Section 3.1, the upper boundary of the points in Figure 10 corresponds to the maximum value of the stau mass, so the lower limits on T R with bino LOSP are presented for three maximum values of the stau mass: 5, 10 and 15 TeV.
A qualitatively different picture emerges when m G 100 GeV. The LOSP lifetime is then so large that the BBN bounds can only be evaded when B h is small and m LOSP 1 TeV with the number density reduced because of low T R . This is, however, only possible for the sneutrino and, very rarely, for the stau LOSP [70,72], as presented in the right panel of Figure 10 for m G = 1 TeV. 10 Hence also for m G 100 GeV we find a lower bound T R 150 GeV. This is true if direct and cascade decays of the inflaton field to the LOSP can be neglected; otherwise, the lower limit on T R becomes less severe, similarly to the bino LOSP case.
If one assumes gaugino mass unification at the GUT scale, then the lower limit on the chargino mass from collider searches, m χ ± 1 > 94 GeV [25], can be translated into a lower limit on the lightest neutralino mass m χ > 46 GeV. This in turns implies in our p10MSSM scan m G 0.1 MeV, where we assume soft scalar masses not to be greater than ∼ 15 TeV and T R low enough so that the gravitino is produced only in NTP. For much lighter gravitinos, in the keV mass range, the correct abundance can be obtained by thermal production for reheating temperature even of the order of a few hundred GeV (see, e.g., [73]).
It is important to note that the additional contribution to the LOSP relic density resulting from direct and/or cascade decays of the inflaton allows one to consider lower values of the reheating temperature in gravitino DM scenario. In such a case, the lower limit on T R becomes less severe, as it is illustrated in the right panel of Figure 11 for the bino LOSP; the same is true for the slepton LOSP.

Conclusions
Motivated by the observation that in scenarios with a low reheating temperature DM relic density is reduced with respect to the standard high-T R case, in this paper we studied the impact of assuming low T R on the phenomenologically favored regions of the (C)MSSM and the NMSSM with the singlino DM. We considered two distinct DM candidates: the LOSP and the gravitino.
In the case of the LOSP we found that, at low T R large regions of the parameter space open up which are normally considered excluded because of too large a relic density. With T R in the range 100 − 200 GeV, the DM can be the bino (coannihilating with staus), the heavy ( > ∼ 3.5 TeV) wino or the higgsino if it is not lighter than about 1 TeV. For T R = O(10) GeV the allowed regions of the parameter space mainly correspond to a bino-like neutralino in the bulk region, with a small fraction of solutions having a higgsino admixture of a few per cent. Similarly, in the singlino-dominated region of the NMSSM, when T R is less than about 200 GeV, large regions open up where at high T R the relic density can be very high.
If DM consists of nonthermally produced gravitinos only, then the relic abundance of LOSPs decaying into gravitinos must be greater than the observed dark matter abundance. Since the effect of low T R is to reduce the LOSP relic abundance, T R cannot be too low. In this case we obtain lower bounds on T R by combining the assumed generous upper bounds on the superpartner masses of a few TeV (which reflects our view that SUSY should not lead to a too severe hierarchy problem) and the BBN constraints. For bino (slepton) LOSP, we find the bound on T R of the order of 100 GeV for the gravitino mass in the range 0.1 − 10 (10 2 − 10 3 ) GeV. These limits are alleviated when significant direct and cascade decays of the inflaton field to the LOSP are present.
LR is also supported in part by a STFC consortium grant of Lancaster, Manchester, and Sheffield Universities. The use of the CIS computer cluster at the National Centre for Nuclear Research is gratefully acknowledged.