Probing the BSM physics with CMB precision cosmology: an application to supersymmetry

The cosmic history before the BBN is highly determined by the physics that operates beyond the Standard Model (BSM) of particle physics and it is poorly constrained observationally. Ongoing and future precision measurements of the CMB observables can provide us with significant information about the pre-BBN era and hence possibly test the cosmological predictions of different BSM scenarios. Supersymmetry is a particularly motivated BSM theory and it is often the case that different superymmetry breaking schemes require different cosmic histories with specific reheating temperatures or low entropy production in order to be cosmologically viable. In this paper we quantify the effects of the possible alternative cosmic histories on the ns and r CMB observables assuming a generic non-thermal stage after cosmic inflation. We analyze TeV and especially multi-TeV super-symmetry breaking schemes assuming the neutralino and gravitino dark matter scenarios. We complement our analysis considering the Starobinsky R2 inflation model to exemplify the improved CMB predictions that a unified description of the early universe cosmic evolution yields. Our analysis underlines the importance of the CMB precision measurements that can be viewed, to some extend, as complementary to the laboratory experimental searches for supersymmetry or other BSM theories.


Introduction
The cosmic evolution before the Big Bang Nucleosynthesis (BBN) and after inflation is much unknown. To date there are no direct observational probes that can constrain this very early universe period, that can be called dark pre-BBN period. On the other hand, inflation that takes place at energy scales much higher than the BBN gives concrete predictions thanks to the presence of the quasi-de Sitter horizon. It is actually the dark pre-BBN cosmic phase that introduces an uncertainty at the inflationary predictions parametrized by the number of e-folds N * . This uncertainty could be minimized if the physics that operates beyond the Standard Model of particle physics (BSM) was known. Indeed, different BSM JHEP02(2018)118 scenarios often imply a different cosmic evolution in order to satisfy the BBN predictions and the observed dark matter abundance Ω DM h 2 = 0.12 [1,2].
The fact that the N * is modified by the details of the dark pre-BBN stage [3] motivate us to investigate this small but non-zero residual dependence of the inflationary predictions on the tentative BSM physics. In most of the inflationary models, a precise measurement of the spectral index n s (N * ) and tensor-to-scalar ratio r(N * ) value accounts for an indirect measure of the reheating temperature of the universe [4][5][6][7][8][9][10][11][12][13] and hence one could in principle examine the cosmology of theories beyond the Standard Model of particle physics as well as non-trivial extensions of the Einstein gravity [14]. From the inflation phenomenology point of view, for a given concrete BSM scenario a predictive inflationary model can be spotted on the (n s , r) plane, whereas from the particle physicist point of view, for a given predictable inflationary scenario the precise measurement of the (n s , r) observables is a measurement of the BSM effects on the cosmic evolution. In other words, we can say that the (n s , r) precision measurements provide us with a cosmic selection criterion for the assumed BSM physics. Planck collaboration has constrained the spectral tilt value of the curvature power spectrum and the tensor-to-scalar ratio at n s − 1 = −0.032 ± 0.006 at 1σ and r < 0.11 at 2σ respectively [1,2]. The current resolution of the temperature and polarizartion anisotropies of the CMB probes, although unprecedented, has not been powerful enough to support or exclude the different BSM physics schemes. There are promising prospects that the proposed next generation CMB experiments, such as the LiteBIRD [15], Core+ [16], CMB-S4 [17], PRISM [18], PIXIE [19], will improve significantly on this direction. The sensitivity forecasts for n s and r is of the order of 10 −3 and such a measurement will account for a substantial leap forward at the observational side.
We aim at this work to show how one can systematically extract non-trivial information about the BSM physics via the CMB precision measurements. We mostly focus on the supersymmetry since we consider it as a compelling BSM theory that remains elusive from the terrestrial colliders. A precise knowledge of the (n s , r) values can indicate us the duration of non-thermal phase after inflation and in this paper we use this information to examine whether different supersymmetry breaking schemes can fit in this picture of the early cosmic evolution.
From the experimental side, there is no signal that supports the supersymmetry hypothesis until today, see e.g. a recent analysis of searches at the LHC [20,21]. The absence of signals arouses increasing concern that supersymmetry does not fully solve the hierarchy problem suggesting that supersymmetry, if realized, may lay at energy scales much higher than the TeV scale. Multi TeV supersymmetry implies that the Large Hadron Collider (LHC) at CERN may find no BSM signal and the fiducial BSM physics scenarios will remain elusive for an unspecified long time. However from the telescopic observational side, the increasing sensitivity of the CMB probes has opened up a rich phenomenological window to the ultra high energy scales of cosmic inflation and indirectly to the dark pre-BBN period.
Definitely, the idea that the CMB studies may probe energy scales well above the TeV is not a new one. There are numerous of seminal works in the literature that examine the impact of BSM physics, and in particular supersymmetry, on the CMB power spectrum mainly either from the inflationary model building or from the dark matter perspective.

JHEP02(2018)118
However, successful inflation models can be consistently embedded into a supergravity framework often without any change in the inflationary dynamics since the inflationary trajectory may remain intact by the presence of additional supersymmetric fields that are efficiently stabilized. Moreover, it is often the case that studies of supersymmetric dark matter cosmology focus on the dark matter density parameter fitting, Ω DM h 2 = 0.12, neglecting other features of the scalar power spectrum.
The degeneracy between supersymmetric inflation models and with their nonsupersymmetric versions in terms of the n s (N ) and r(N ) observables can break due to the different post-inflationary evolution. The thermal evolution of a supersymmetric plasma is in general much different when supersymmetry is realized in nature [22]. Actually, the null LHC results push the sparticles mass bounds to larger values that spoil the nice predictions of the thermal dark matter scenario [23]. Therefore, assuming that the LSP is part of the dark matter in the universe the Ω LSP h 2 0.12 constraint reconciles only with particular radiation domination histories which may greatly differ to the simple scenario of a single and smooth radiation phase after the inflaton decay. An interesting point, that stimulates this work, is that the features of the radiation dominated phase depend on the details of the supersymmetry breaking patterns.
In order to extract information about the BSM supersymmetric scenarios from the (n s , r) precision measurements we utilize existing results on supersymmetric cosmology aiming at an analysis based on assumptions as minimal as possible. We consider that the MSSM plus the gravitino is the necessary minimal set-up that gives the most conservative results. We a priori consider the T rh and the supersymmetry breaking scale as unknown quantities. We estimate the neutralino and gravitino LSP abundances by scanning the sparticle mass parameter space. As a rule of thumb we adopt the classification of quasinatural, split and high scale supersymmetry when we scan the possible energy scales of supersymmetry breaking. As expected, see e.g. [24][25][26], we find that most of parameter space of supersymmetric theories yields an excessive dark matter abundance. Our perspective in this work is that the parameter space that yields an excessive dark matter abundance should not be faced as a cosmologically forbidden one but, on the contrary, as a parameter space that favours a different cosmic history for the very early universe. Namely, excessive LSP abundance implies either a low reheating temperature after inflation or low entropy production. Both cases have a non-trivial impact on (n s , r) observables, see e.g. [27] for a relevant analysis on non-thermal neutralino dark matter and [28] for a recent analysis on FIMP dark matter.
Departing from the minimal field content analysis, i.e. the MSSM, the overabundance problem in general deteriorates. Indeed, the dark matter abundance receives contributions from the perturbative and non-perturbative decay processes of the inflaton field [29] and from thermal scatterings, thermal and non-thermal decays of fields coming from the supersymmetry breaking sector such as the messengers. Extra fields can however decrease the DM abundance if they decay late and dominate the energy density of the early universe e.g. due to coherently oscillating scalars or scalars that cause thermal inflation. Such fields are rather common and well motivated in many BSM schemes such as supersymmetry; common examples are the moduli, supersymmetry breaking fields, the saxion, etc. Here

JHEP02(2018)118
we collectively label X any of this sort of scalars and explicitly refer to it as diluter, since what we actually measure on the CMB is the diluter impact on the expansion history. In our analysis, the diluter is the only field beyond the MSSM and gravitino that we consider. Finally, in order to perform a complete calculation of the spectral index value we consider the Starobinsky R 2 inflation model and we compare the R 2 inflation and R 2 supergravity inflation predictions by taking into account the effects of the post-inflationary phase.
Apparently one cannot exclude or verify supersymmetry by n s and r precision measurement, nevertheless one can indeed support the presence of BSM physics or, to put it differently, rule out the so-called BSM-desert hypothesis for a particular inflation model. This is a minimal but undoubtedly an exciting possibility given the fact that terrestrial colliders probe only a small part of the vast energy scales up to the Planck Mass, M Pl , and supersymmetry or any other BSM scale may lay anywhere in between. It is also exciting to note that the terrestrial experiments, such as colliders and direct detection experiments, are sensitive to low scale supersymmetry whereas the CMB observables are more sensitive to high scale supersymmetry. Hence precision cosmology can offer us complementary constraints to the parameter space of the supersymmetric theories. This prospect, though very challenging, is actually a feasible possibility.
The organization of the paper is the following. In section 2 we parametrize the uncertainty in the n s and r values coming from the unknown value of N * due to the dark pre-BBN era. We compute the shift in the spectral index and tensor-to-scalar ratio with respect to the dilution magnitude in a general BSM context. In section 3 we overview key results of neutralino, gravitino and briefly the axino cosmology regarding the LSP yield, that are necessary for the estimation of the dilution magnitude. In section 4 we analyze the implications of various supersymmetry breaking patterns to the early universe cosmology and examine the features of the possible alternative cosmic histories. In section 5 the Starobinsky R 2 inflation is used as a specific example to demonstrate a full computation of the spectral index and tensor-to-scalar ratio shift. A comparison between the theoretical predictions of the R 2 and supergravity R 2 inflation is also performed. In the last section we outline the main idea and the method proposed in this work and we comment on the future theoretical and observational prospects.

CMB observables and the post-inflationary evolution
It is convenient to expand the power spectra of the dimensionless curvature perturbation as where A s is the scalar amplitude and the powers of the expansion are the scalar spectral index n s , the running and the running of the n s . In general one can assume that the scale dependence of the spectral index to be given at leading order by the expression

JHEP02(2018)118
where N * is the number of e-folds remaining till the end of inflation after the moment the pivot scale k * exits the Hubble radius, N * ≡ t end t * Hdt = ln(a end /a * ). The N * is a critical quantity that determines the n s value. It carries the information of how much the observable k −1 * CMB scale has been stretched since the inflationary era. The uncertainty on the N * comes mainly from the post-accelaration stage and induces an uncertainty on the spectral index value given by the n s running that for the eq. (2.2) reads For ∆N ∼ 1 − 10 the ∆n s is of size O(1 − 10) , that is within the accuracy of the future observations. To explicitly estimate the N * value one relates the size of the scale k −1 * = (a * H * ) −1 , which exited the Hubble radius H −1 * during inflation, to the size of the present Hubble radius where the subscripts refer to the time of horizon crossing ( * ), the time inflation ends (end), the time BBN takes place (BBN), the radiation-matter equality (eq) and the present time (0). We defineÑ dark the number of e-folds from the end of inflation until the beginning of the BBNÑ wherew dark stands for the average value of the equation of state parameter during the dark pre-BBN period, andw dark = −1 has been assumed. We call this period dark due to the lack of observational evidences of the transition to the radiation dominated phase from the super-cooled conditions during inflation. Unless exotic forms of matter are assumed, such as thermal inflation or stiff fluid domination, we can estimate the maximum value of thẽ N dark to be around 56 forw dark = 0 and the minimum to be around 41 forw dark = 1/3. The observational uncertainty for temperatures T 1 MeV ∼ T BBN [30] implies an uncertainty at the e-folds of inflation about ∆N ∼ 15. We can split theÑ dark intõ whereÑ rh = ln(a rh /a end ) stands for the e-folds number of the postinflationary reheating period until the complete decay of the inflaton,Ñ rad the e-folds number of the radiation dominated era that preceded the BBN andÑ X stands for the e-folds number that take place during the domination of an arbitrary X field in the period after the decay of the inflaton and before BBN. After plugging in the value for the ratio a eq H eq /(a 0 H 0 ), the relation (2.4) is recast into [2]

JHEP02(2018)118
Utilizing the relation P R (k * ) = V * /(24π 2 * M 4 Pl ) = A s and after substituting numbers for the ratio k * /(a 0 H 0 ) we get We adopted the Planck collaboration pivot scale, k * = 0.002Mpc −1 and the measured value ln(10 10 A s ) = 3.089 [1,2]. We also introduced the ∆N dark factor to mark explicitly the uncertainty of the dark pre-BBN era on the N * value, We have split the ∆N dark into the contributions from the inflationary reheating, the Xdomination and the pre-BBN radiation domination period. It is ∆N rad = 0 sincew rad = 1/3 and where ρ dec X is the energy density of the thermal plasma right after the decay of the scalar X. In principle, for a concrete and predictable inflationary model thew rh and the reheating temperature after inflation can be estimated and hence the ∆N rh . The crucial quantity is the decay rate Γ inf of the inflaton which determines the reheating temperature. Assuming that the decay and the thermalization occur instantaneously at the time Γ −1 inf then the reheating temperature is found by equating (and omitting order one coefficients) Γ inf = H, The maximum temperature possible is achieved in the instant reheating scenario. Apparently when T rh = T max = ρ 1/4 end (30/π 2 g * rh ) 1/4 it is ∆N rh = 0. Note that the N * has a logarithmic dependence on g * rh , with g * rh being the effective number of relativistic species upon thermalization.
It is however well possible that after the inflaton decay the evolution of the universe could have been episodic with additional reheating events after inflation. Hence the cosmic thermal era could have started after the last reheating stage before primordial nucleosynthesis caused by other than the inflaton scalar field, for instance a modulus or a flaton [31] that we collectively label X. Here, we prefer to remain agnostic about the identity of X but we do utilize its property to cause efficient dilution and low entropy production. The X can dominate the energy density of the universe over radiation due to the slower redshift of its energy density stored. It is ρ X ∝ a −3 for a scalar condensate that coherently oscillates in a quadratic potential and ρ X ≈ constant for a scalar field with sufficiently flat potential that causes thermal inflation.

2.1
The shift in the scalar spectral index and tensor-to-scalar ratio due to late entropy production Let us now estimate the impact of the X domination era on the spectral index value. We call N (th) and n (th) s the thermal reference values, that is the e-folds number and the spectral

JHEP02(2018)118
index values respectively if there is no late entropy production after the inflaton decay, i.e. dilution effects. It is at leading order where, following eq. (2.8), At leading order the scalar tilt is generally given by the equation (2.2). Since precision is expected to increase in the future it is worthwhile to consider next-to-leading corrections. Due to the large number of inflationary models [32] there is no common form for the nextto-leading term [33]. A phenomenological way to parametrize it is based on the large N expansion 14) The parameters α and β(N ) are determined only after a particular inflation model is considered. In principle the parameter α can also be a slowly varying function of N [38].
In addition the expansion (2.14), for some inflation models, may involve parameters of the potential [33]. Here we assume that α is a constant and absorb possible complicated behaviors in the arbitrary β(N ) function. In section 5 we will explicitly estimate the shift in the spectral index for the Starobinsky R 2 inflation model where the parameters α and β(N ) have a particular form.
If ∆N X = 0, after Taylor expanding the n s N (th) − ∆N X , the spectral index n (th) s = n s (N (th) ) value is shifted by an amount ∆n s ≡ n s − n (th) s , . (2.16) The " denotes d/dN and β, β , β , β are estimated at N = N (th) . In the above expressions, given than ∆N X > 1 and ∆N X /N (th) < 1, terms of order O ∆N 4 X /N 6 and smaller have been neglected. We have also assumed that the terms in the parentheses in eq. (2.16) are roughly of order β. Otherwise, if β , β , β 1, the F β correction can be important, however such a behavior is not found in any of the known universality classes [35]. One can see that the next-to-leading correction β(N )/N 2 is at most of accuracy and for α∆N X > β the contribution to the spectral index shift is found to be subdominant with respect to the α-dependent terms.

JHEP02(2018)118
In order to specify the ∆N X , elements of the X scalar cosmic evolution have to be specified. When the scalar X coherently oscillates about the minimum of a effectively quadratic potential it isw X = 0. In such a case, at the cosmic time t dom X Γ −1 X the energy density of X is larger than that of the plasma and the universe enters a scalar dominated era that dilutes any pre-existing abundances of the relativistic degrees of freedom at the time of the X decay. The X field decays and reheats the universe with temperature T rh X ≡ T dec X . Considering instant decay of the scalar X, the dilution magnitude is estimated to be where we considered thatw X = 0. After plugging in the dilution magnitude we get The maximum value of the ∆N X ∼ 15 is achieved whenÑ rh → 0 andÑ rad → 0. This case corresponds to the maximum dilution scenario where the X field oscillations dominate the energy density of the universe right after the end of high scale inflation until the onset of BBN. The ∆N X = 0 case corresponds to an uninterrupted radiation phase following the post-inflationary reheating. If someone assumes the presence of X matter with exotic barotropic parameter the ∆N X limit values can be extended. Substituting ∆N X = 1 3 lnD X in the expansion (2.15) we obtain the shift in the spectral index, with accuracy |∆n s |/n s 1 , due to a post-inflationary dilution of the thermal plasma (th) . Notice that at leading order the (2.20) reads ∆n s = −α s ∆N X , where α s = (1−n (th) s ) 2 γ/α is the running of the spectral index at N (th) . We also mention that the three last terms in the brackets of the above equation can be neglected without significant cost in the accuracy. Plugging in the thermal reference value n the thermal plasma. We see that the ∆n s is negative which means that the spectrum tilt becomes more red when dilution of the radiation plasma takes place; this behavior is illustrated in figure 1. The precision of the expression (2.20) is sufficiently good, i.e one per mile, even for the extreme case ∆N X ∼ 15 or D X ∼ 10 20 .
Apart from scalar condensates, several BSM construction, mostly supersymmetric ones, predict the presence of singlets under the Standard Model symmetries that have a relatively flat potential. Such fields can realize the thermal inflation scenario and are generally called flatons [31]. Due to Yukawa interactions the flaton can be trapped at the origin of the field space by thermal effects. At some temperature that we denote T dom X the vacuum energy V 0 of the flaton dominates over the background radiation energy density and a period of thermal inflation starts. Thermal inflation ends at the temperature T 2 when the thermal trap has become too weak and the flaton field starts oscillating about its zero temperature minimum. The flaton finally decays at the temperature that we denote T dec X and we consider instant reheating. The dilution magnitude due to the flaton X domination is Respectively the ∆N X value for flaton domination is and each term can be written in a compact form ∆N X | FD = ∆N X | TI + ∆N X SC , where TI and SC stand for thermal inflation and scalar condensate respectively. The ratio of the relativistic degrees of freedom accounts for a small correction and one can see that it is actually ∆N X | FD ln D FD X /3. The dilution magnitude maximizes when the thermal inflation phase is followed by a scalar condensate domination phase, e.g. when the flaton field decays very slowly. The ∆N due to thermal inflation has an upper bound N TI ∼ 10 in order that the cosmological density perturbation remain intact. Nevertheless, the dilution magnitude can be many orders of magnitude larger than the dilution caused by a scalar condensate domination (2.17) and it can efficiently dilute any overabundant relic such as dark matter particles. Essentially, the shift in the spectral index due to thermal inflation, ∆n s | FD , is given again by the eq. (2.20) simply by replacing the lnD X with the lnD FD X , see Fig 1. It is remarkable that the shift in n s due to period of thermal inflation can resurrect ruled out inflationary models such as the minimal hybrid inflation in supergravity [36].
Finally, let us comment on the shift in the tensor-to-scalar ratio. The phenomenological parametrization of the scalar tilt n s = 1 − α/N implies that the first slow roll parameter = −Ḣ/H 2 writes [38] where A an integration constant coming from the differential equation +d ln /dN = α/N . At first order in slow roll we have r = 16 and the shift in the tensor-to-scalar ratio due to where r (th) = r(N (th) ). For ∆N X = lnD X /3, either due to a scalar condensate domination or thermal inflation, the relation ∆r = ∆r(D X ) is obtained. The scaling (2.23) depends on the potential that implements inflation. Different potentials yield different values for α and A. Moreover, accuracy of order ∆r ∼ 10 −4 requires to go beyond the approximate relation r = 16 and consider corrections at second order in slow roll. In section 5 we will explicitly estimate the ∆r for the Starobinsky R 2 inflation model with the next-to-leading order corrections taken into account. The general conclusion is that, according to eq. (2.24), a non-thermal phase withw X < 1/3 and durationÑ X = [(1 − 3w X )/4] −1 ∆N X increases the tensor-to-scalar ratio value. Summarizing, the duration of a non-thermal phase is encoded in the number of e-folds N between the moment a relevant mode exits the horizon and the end of inflation. If the radiation domination era, where w = 1/3, initiates at the moment of the complete inflaton decay and continues without break until the BBN epoch then the e-folds number, called here thermal e-folds number N (th) , can be explicitly determined by the dynamics and the full interactions of the inflaton field. If not, a non-thermal phase changes the aforementioned efolds by the amount ∆N X ∼ ln D X /3. A dilution of size D X = 20 translates into ∆N X ∼ 1 and a prolonged dilution e.g. of size D X = 10 13 into ∆N X ∼ 10. In order to estimate the shift in the spectral index and the tensor-to-scalar ratio one has to know the N (th) that is JHEP02(2018)118 given by the eq. (2.13). This is possible only after an inflationary model and the parameters describing reheating are chosen. Then from eq. (2.14) the n (th) s = n s (N (th) ) and the n s = n s (N (th) −∆N X ) can be estimated and hence the spectral index shift ∆n s , given by the eq. (2.15) or (2.20), is obtained. In figure 1 we illustrate the shift in the spectral index due to a non-thermal phase that is implemented after reheating and before BBN. In the left panel we considered the Starobinsky R 2 model that predicts T rh ∼ 10 9 GeV [37], and in the right panel a Starobinsky-like potential with non-gravitational interactions and a linear potential V ∝ φ both characterized by a fiducial reheating temperature T rh = 10 12 GeV. The knowledge of these inflaton features enables the explicit calculation of the n (th) s value, that corresponds to the red dots in the plots. A scalar condensate domination or thermal inflation shifts the spectral index value according to the formula (2.20) as illustrated in the figure 1.
From a more bottom-up approach, the postulation of a non-thermal phase during the pre-BBN era is not enough to determine the ∆n s and ∆r. Although a rough estimation of the spectral index shift can be done by the approximate expression (2.3) the result is far from accurate and cannot consistently constrain the early universe cosmic history. The best method is to choose an inflation model that is in accordance with a particular BSM description of the early universe (e.g. a supersymmetric, stringy or modified gravity framework) and estimate the ∆n s and ∆r according to the pre-BBN cosmology implied by the BSM theory at hand. Examples of BSM cosmic processes connected with the expansion history of the universe are the dark matter production and the baryogenesis processes. In the following we will consider the supersymmetric BSM scenario and determine features of the pre-BBN cosmology, such as possible non-thermal stages, that allow the accommodation of different supersymmetry breaking schemes assuming that the LSP is part of the dark matter in the universe. We will estimate the minimum dilution size dictated by the requirement Ω LSP h 2 ≤ 0.12 and determine the expected shift in the spectral index and tensor-to-scalar ratio when a particular inflation model, which in section 5 is the Starobinsky R 2 model, complements the description of the early universe evolution.

Supersymmetric dark matter cosmology
In the previous section we computed the shift in the spectral index and tensor-to-scalar ratio due to post-inflationary entropy production. The fact that the present universe contains dark matter with relic density Ω DM h 2 = 0.12 relates the amount of the dilution with the dark matter production. In this section, focusing on TeV and especially multi-TeV supersymmetric scenarios, we will overview the expected LSP yield. We will stress out that the dilution is generally required, hence the CMB inflationary observables should be non-trivially influenced by the post inflationary expansion history of the supersymmetric universe.
There are several fundamental theoretical reasons to believe that supersymmetry is a symmetry of nature. For the devotee of supersymmetry the central question is the scale that supersymmetry is realized. The direct superpartner LHC-limits for all colored sparticles exceed 1.5 TeV and suggest that we should depart from scenarios with natural supersymmetry paying the price of pushing the amount of tuning at the MSSM to less than JHEP02(2018)118 0.5 − 1 percent level. However, the absence of BSM signals in the LHC rules out only the electroweak scale supersymmetry and not supersymmetry in general.
BSM physics scenarios with unnatural supersymmetry are still very appealing. Gauge coupling unification, the presence of a stable dark matter particle, the possible baryogenesis processes and the stringy UV completion of the low energy theories do not link SUSY with the electroweak scale. Supersymmetry may appear at higher energy scales. In ref. [39], different supersymmetry breaking scenarios have been categorized according to the mass spectrum features into three representative cases: i) Quasi-natural supersymmetry, in which supersymmetric particles are heavier than the weak scale, but not too far from it, about in the 1 − 30 TeV range. ii) Split supersymmetry, in which only the scalar supersymmetric particles have masses of the order ofm, while gauginos and higgsinos are lighter, possibly with masses near the weak scale [40][41][42]. There are also the Mega-Split [43] or Mini-Split [44] scenarios. iii) Finally the High-Scale supersymmetry, see e.g [45,46] in which all supersymmetric particles have masses around a common scalem, unrelated to the weak scale. Them is constrained by the Higgs mass value according to the details, of each supersymmetry breaking scenario. Roughly in the Split supersymmetry the maximum value allowed form is 10 8 GeV when tan β is small, while in the High Scale supersymmetry them value can be up to 10 12 GeV.
For our analysis it is critical that the LSP is stable. The stability of the LSP dark matter is assured by the presence of a discrete symmetry of the supergravity Lagrangian, the R-parity. If the R-parity is violated then the cosmological constraint Ω DM h 2 = 0.12 is raised for the LSP. Although R-parity violating models have been actually constructed and have interesting phenomenological implications [47], there are strong arguments based on GUT models that support the R-parity conservation even when the scale of supersymmetry breaking is well above the electroweak scale [48]. These results motivate us to assume that the LSP lifetime is much larger than the age of the universe and thus the LSP is constituent of the dark matter.
Given the supersymmetry breaking scheme the stability of the LSP puts strong constraints on the thermal history of the universe. In the following subsections we overview the basic relevant cosmological aspects and results of the gravitino and neutralino LSP scenarios necessary for the goals of our analysis.

Gravitino dark matter
The gravitino is the supersymmetric partner of the graviton in supergravity and it can acquire a mass in the range of O(eV −m). The gravitino is naturally the LSP in gauge mediated supersymmetry breaking models (GMSB), see [49] for a review, and possibly it is the LSP in Split and High scale supersymmetry frameworks. The relic density of the gravitinos Ω 3/2 h 2 , which can be thermal or non-thermal, receives contributions from many sources.
Thermal gravitinos (freeze out): from scatterings (i) with the MSSM plasma, (ii) with the thermalized messenger fields. . Density and contour plot of the decadic logarithm of the required dilution for T rh = 10 9 GeV reheating temperature after inflation and gravitno the stable LSP. In the left panel degenerate spectrum for the sfermions and gauginos was considered, mf = mg, while in the right panel there is a split spectrum with mf = 10 3 mg. Thermal production of helicity ±3/2 and ±1/2 gravitinos from scatterings in the plasma, non-thermal production from decays of sfermions and the NLSP to helicity ±1/2 gravitinos have been taken into account. The contributions to the gravitino abundance have been conditionally added, i.e. in the parts of the contour that thermal equilibrium is achieved the total abundance is replaced by the thermal one. The magnitude of the logarithm of the required dilution is given by the contour numbers onto the density plot. Negative numbers correspond to underabundance, hence to no dilution contour area.
The less model independent estimation of the Ω 3/2 is achieved when only the MSSM sector is considered. The gravitino number density n 3/2 in the thermalized early universe evolves according to the Boltzmann equation [50]. A key quantity is the gravitino production rate, γ sc , in scatterings with thermalized Standard Model particles and sparticles The gravitinos obtain a thermal distribution via interactions with the MSSM for where mg 3 is the gluino mass evaluated at the reheating temperature, see eq. (3.4). If the reheating temperature is below the 3/2 . Furthermore, the heavier MSSM sparticles are unstable and will decay to gravitinos. The decay width into gravitinos is nearly the same for both gauginos and sfermions

JHEP02(2018)118
whereĩ =g,f . The total MSSM contribution to the gravitino yield is Y MSSM , and the relic density parameter reads The gravitino relic abundance sourced by the MSSM and messenger fields is illustrated in figure 2 and 3. In the case that the gravitino is the only sparticle with mass below the reheating temperature then the gravitino relic abundance is given by a much different expression with dependence Ω 3/2 ∝ m −3 3/2 T 7 rh [51,52]. Apart from particular cases, the above gravitino yield (3.3) cannot be final because we neglected sources beyond the MSSM. The supersymmetry breaking sector is a necessary ingredient for all the consistent supersymmetric BSM scenarios [49]. In general the extra fields only increase 1 the final Ω 3/2 , unless there is a late entropy production. For example, thermalized messengers fields generically equilibrate the gravitinos for broad range of values of the Yukawa coupling at the messenger superpotential, λ mess 10 −6 − 10 −5 [56], and the relic gravitino density parameter reads where the freeze out temperature is here equal to the messenger mass scale, T f.o. 3/2 ∼ M mess . Even if λ mess 1 the thermal scatterings of messengers contribute to gravitino relic density with Ω mess 3/2 h 2 ∼ 0.4 M mess /10 4 GeV GeV/m 3/2 (mg/TeV) 2 . In addition, the inflaton perturbative decay produces non-thermal gravitinos with rate [57,58] Also gravitinos are produced during the preheating stage via its non-perturbative decay of the inflaton [59][60][61][62][63][64][65], from the decay of the supersymmetry breaking field, see e.g. [66][67][68], or other moduli [69][70][71]. Therefore, the estimation of the gravitino relic abundance based solely on the MSSM sector gives a model independent albeit an underestimated and hence conservative value for the Ω 3/2 .
The Ω 3/2 result could decrease in the case that extra fields interrupt the thermal phase, e.g. due to the domination of a non-thermal scalar field that produces entropy at low temperatures. Thanks to the dilution the gravitino cosmologically problematic supersymmetric scenarios may become viable possibilities. The tentative low entropy production is caused by the scalar X that we do not identify and collectively call it diluter. We only assume that 1 It is though possible that the supersymmetry breaking sector leads to a suppressed Ω 3/2 , e.g due to R-symmetry restoration [53], or a high temperature decoupling of the messenger fields [26], or due to the dynamics of the sgoldstino field [54], or due to feeble couplings in the supersymmetry breaking sector [55]. Log D X , T rh 10 9 GeV , G LSP Figure 3. Density and contour plot of the decadic logarithm of the required dilution for gravitino LSP and reheating temperature T rh = 10 9 GeV. In the left panel the gravitinos are thermal (heavy gravitinos can thermalize due to the messenger sector). In the right panel the contribution of messengers plus MSSM is considered, assuming a small enough messenger coupling so that gravitinos do not thermalize by the interactions with messengers, but only due to the MSSM. The messengers scale is taken to be M mess = 10 8 GeV.
it interacts too weakly with the other fields, e.g. via gravitational interactions. Therefore the gravitino relic density parameter is the conditional sum where Ω MSSM 3/2 is the contributions of the MSSM (scatterings and decays), Ω mess 3/2 is the contribution of messengers (scatterings and decays), Ω inf 3/2 is the contribution of the infationary perturbative and non-perturbative decay and Ω SB 3/2 is the contribution of the supersymmetry breaking field. It is called conditional sum because the simple add of each contribution may result in an overestimate of the gravitino abundance. For example the presence of thermalized messengers modifies the gravitino production from the MSSM sector [56,72]. We mention that the sum (3.6) is not strictly exact: it is well possible that contributions from non-thermal decays, that take place below the T f 3/2 temperature, increase the Ω tot 3/2 value. Finally, the presence of a scalar X that produces low entropy modifies the result (3.6) as will be discussed in the section 4. In such a case the density parameter (3.6) value is renamed Ω < 3/2 in order to emphasize that it is sourced by processes taking place before the X decay.

Neutralino dark matter
The lightest neutralinoχ 0 is the most representative example of a WIMP dark matter and an appealing candidate thanks to its main merit: the thermal production mechanism.

JHEP02(2018)118
The thermal neutralino scenario however works best provided that the squark and slepton masses lie in the 50-100 GeV range [24,73] which has been excluded by collider searches. In addition the direct and indirect detection experiments shrink the parameter space of the neutralino with mass about the electroweak scale [23]. Specific neutralino types, such as the higgsino, see e.g. [74], or the annihilation mechanism can be invoked to match the Ωχ0h 2 to data, but in general a rather heavy neutralino cannot be a viable thermal relic.
The neutralinoχ 0 decouples from the thermal bath at a freeze-out temperature T f.o.
χ 0 the neutralinos reach thermal and chemical equilibrium and the relic density parameter is UV insensitive and depends on theχ 0 mass squared, Ω (th) When the sparticles masses lay well above the TeV scale the thermal neutralino scenario is disfavored and ususally non-thermal production scenarios are considered, e.gχ 0 production via the decay of heavy gravitinos.
The gravitinos, that are unstable, are produced via thermal scatterings, non-thermal decays of sfermions and possible decays of scalars beyond MSSM such as the inflaton [57][58][59][60][61][62][63][64][65] and the supersymmetry breaking field or other moduli [69][70][71]. Focusing on the MSSM sector the gravitinos dominate the universe either for large enough reheating temperature, T rh 5 × 10 14 (m 3/2 /10 5 GeV) 1/2 GeV, or large enough sfermion masses, mf 2 × 10 8 (m 3/2 /10 5 GeV) 5/6 [42]. The gravitinos decay when Γ 3/2 = H and the temperature after decay is T dec 3/2 = 6.8 Apparently, it has to be m 3/2 > 10 4 GeV to avoid BBN complications [75,76]. The 4 He abundance implies that it must be Y 3/2 10 −12 , for m 3/2 = 10−30 TeV and for smaller m 3/2 values the bound becomes much severer, see e.g. [57,75,76] for details. The gravitino decay populates the universe with neutralinos. Heavy enough gravitinos, m 3/2 10 7 GeV, decay promptly so that T dec 3/2 > T f.o. χ 0 and the neutralinos reach a thermal equilibrium. In the opposite case, the neutralinos produced by the graviton decay are out of chemical equilibrium and either have a yield Yχ0 ∼ Y 3/2 for a radiation dominated universe, or Yχ0 3T dec 3/2 /(4m 3/2 ) for a gravitino dominated early universe. Unless the reheating temperature is particularly high T rh > 10 14 GeV or the sfermions very massive, mf > 10 8 GeV the gravitinos do not dominate over the radiation, and the neutralino relic density parameter reads Thus one finds  . Density and contour plot of the decadic logarithm of the required dilution for neutralino stable LSP. In the left panel the neutralino abundance is the thermal one. In the right the reheating temperature is T rh = 10 12 GeV the sfermions are 10 3 times heavier than the gravitinos and hence the neutralino yield is dominated by the decay of gravitinos produced from sfermion decays; in the right bottom corner of the plot the neutralinos thermalize due to large T dec 3/2 . The gravitino mass is taken to be m 3/2 > 10 5 GeV to avoid BBN constraints.
where i runs up to N = 46 for sfermions heavier than the gravitino and T dec χ 0 were taken to be g * = 86.25. If gravitinos dominate the early universe that is D 3/2 1, then the relic density parameter of the non-thermally produced neutralinos is TeV (3.10) It is possible that the non-thermally produced neutralinos from the gravitino decay achieve a chemical equilibrium for nχ0 σχ0v > H(T dec 3/2 ). It is σχ0v ∝ 1/m 2 χ 0 and H(T dec 3/2 ) ∝ m 3 3/2 , hence for a wino-like neutralino at the TeV scale and m 3/2 > 10 5 GeV pair-annihilation can take place [77]. The neutralinos annihilate until their number density becomes n crit χ 0 ∼ 3H/ σχ0v and the relic density parameter is for this case, Ω χ 0 /T dec 3/2 , that is enhanced by the ratio (T f.o. χ 0 /T 3/2 ) compared to the thermal abundance. This is an appealing scenario, called annihilation scenario, because the critical value n crit χ 0 behaves as an attractor and determines the relic abundance of neutralino (mostly wino) LSP, making it independent of the primordial gravitino relic abundance [77]. Nevertheless it hardly works when one departs from the TeV scale neutralino. It is also much constrained from the indirect detection experiments.
In the section 4 the non-thermal scenario, that is often-called branching ratio scenario, where theχ 0 is produced non-thermally during the low entropy production caused by the diluter X field will be discussed, and in section 5 we will consider the production ofχ 0 from the supersymmetry breaking field aiming at a complete analysis.

Axino dark matter
In the sake of completeness of the basic LSP scenarios, we briefly comment here on the axino dark matter. In supersymmetry, the axion solution to the strong CP problem comes with an extra scalar, the saxion and a fermion, the axinoã. If the axino is the LSP it is a well motivated dark matter candidate [78,79]. It freezes out at high temperatures T f.o. a ∼ 10 11 GeV(f a /10 12 GeV) 2 , where f a the axion decay constant. At lower temperatures it can be produced from thermal scatterings and decays. In that case, for a radiation dominated universe, the axino relic density parameter is the sum of the contributions from thermal scatterings, the gravitino decay and the NLSP decays for T dec 3/2 below the NLSP freeze out temperature. We note that the two body decay of a squark to an axino is subdominant for gluino masses less than squark mass [80]. It is Ω MSSM(sc) a ∼ 2.8 × 10 8 (mã/GeV)Yã where Yã(KSVZ) ∼ 10 −7 (T rh /10 4 GeV)(10 11 GeV/f a ) 2 for the KSVZ axion model, see e.g [81], and Yã(DFSZ) ∼ 10 −5 (µ/TeV) 2 (10 11 GeV/f a ) 2 for the DFSZ axion model where µ the superpotential Higgs/Higgsino parameter, see e.g [82].
For axino mass not much smaller than the NLSP, the axino dark matter case is quite similar to the neutralino LSP. For mã TeV the axino dark matter is also cosmologically problematic since its relic density parameter generally violates the Ω DM h 2 = 0.12 bound, and the essential conclusion is that, in general, a special thermal history of the universe is required for the axino dark matter scenario as well. Remarkably in these models, the saxion can play the rôle of the diluter X for its condensate decay can produce late entropy that successfully decreases the LSP abundance [83], see also [84] for some recent results on the reheating temperature and the Ω LSP constraint.

Alternative cosmic histories and supersymmetry
The overview of the predicted relic density of supersymmetric dark matter in section 3 suggests that the observational value of Ω DM h 2 gets generally severely violated when the sparticle masses increase. For gravitino and neutralino LSP one can collectively write down a general scaling with respect to the mass parameters and temperature and where the exponents (α, β, γ, δ) and (α,β,γ,δ) are either positive or zero, depending on the dark matter production mechanism considered. The predicted supersymmetric dark matter overdensity for "unnatural" supersymmetry can be reconciled with the Ω DM h 2 bound if the reheating temperature is rather low or JHEP02(2018)118 late entropy production takes place. Remarkably, both solutions imply that an alternative cosmic history takes place if supersymmetry is a symmetry of nature. By the term alternative cosmic history we mean that the radiation domination phase after inflation was interrupted or delayed by a cosmic era, where a fluid X with barotropic parameter w X < 1/3 dominated the energy density of the early universe. As discussed in the introduction and in section 2 such a cosmic era impacts the observable values n s and r. In order to quantify this effect we consider in our analysis below different cosmic histories and different supersymmetry breaking schemes. We follow the base line framework of the benchmark supersymmetry breaking scenarios with either gravitino or neutralino LSP and degenerate or split mass spectrum. The scale of supersymmetry breaking, represented by the general sfermion mass m, is taken to be from the TeV scale up to the energy scale of the reheating temperature.

Low reheating temperature
The reheating temperature of the universe after inflation can be rather low if the inflaton decay rate, Γ inf , is small enough or if it is the result of the decay of a weakly coupled scalar unrelated to the inflaton. 2 In this case, the dark matter production due to processes sensitive to the maximum temperature gets suppressed.
We call low reheating temperature scenarios those with T rh 10 5 GeV. For gravitino LSP the yield from thermal scatterings decreases when the reheating temperature decreases, and the NLSP-decays to gravitinos account for the leading contribution to Ω 3/2 for m 3/2 ∼ m NLSP . On the other hand, for neutralino LSP the UV-sensitivity of the Ωχ0 to processes that take place at high temperatures is small. The neutralino abundance is IR-sensitive and it is mostly determined at the freeze out temperature T f.o. χ 0 . Both for gravitino and neutralino LSP, the observational bound Ω DM h 2 = 0.12 is generally violated for m LSP > O(TeV) andm < T rh . Apparently, in the MSSM the Ω DM h 2 = 0.12 can be satisfied for m LSP O(TeV) or for the particular case thatm ∼ T rh where the Boltzmann suppression may play a critical rôle. Note that the X domination cosmic phase is a decaying particle dominated phase, hence entropy is gradually produced for Γ X /H < 1, where Γ X the decay rate of the X particle. The maximum reheating temperature is greater that T dec X and this has implications for the relic LSP density [85]. If the LSPs reach chemical equilibrium before reheating, the relic LSP energy density, is roughly given by Ω LSP are the freeze-out temperatures for T rh m LSP and T rh m LSP respectively [85]. We have also called T rh the X decay temperature, T dec X . On the other hand, if T rh m LSP and the LSPs never reach a chemical equilibrium then the relic density has a dependence Ω LSP ∝ T 7 rh [86]. Finally, if the LSPs are produced non-thermally from the X decay and reach chemical equilibrium then the relic density reads Ω (th) LSP /T rh ), see [87] for a brief overview on the topic. Note that these scenarios that can reconcile heavy supersymmetry with the observational 2 Low reheating temperatures may be caused by the a scalar field X (or more than one scalar) with relatively long lifetime, ΓX Γ inf that dominated the energy density of the universe before the inflaton decay, e.g if the X is frozen during inflaton oscillations with the ρX ∼ m 2 X X 2 sufficiently large that sources some extra e-folds of X inflation. In this case, the reheating temperature at the expressions (4.1) and (4.2) is T rh = T dec X .

JHEP02(2018)118
bound Ω DM h 2 = 0.12 work mainly for the TeV neutralino dark matter scenario and share the common feature that T rh < T f.o. LSP . To this end, one draws the general conclusion that the gravitino or neutralino LSP relic density for "unnatural" supersymmetrym > O(TeV) and m LSP > O(TeV) requires the reheating temperature after inflation to be below or about the supersymmetry breaking scale, Otherwise the dark matter is overabundant. Let us mention that the very interesting scenario of EeV gravitino [51,52,88]

Late entropy production
If the inflaton decay reheated the universe it is generally expected that T rh 10 9 GeV, for inflaton mass m Φ ∼ 10 13 GeV. Such reheating temperatures mean that the dark matter generation processes that take place at high temperatures are critical for the determination of the dark matter abundance.
The abundance that is more sensitive to UV processes is that of the gravitino. In the MSSM framework, the leading contribution to Y 3/2 depends on the maximum temperature after reheating and the decays of the heavy thermalized sfermions. As illustrated in the figure 2 and 3, the Ω 3/2 h 2 /0.12 increases with them and m 3/2 . For T f.o.
3/2 < T rh the gravitino abundance is the thermal one, and in particular cases it may be enhanced due to late sfermion decays.
In the neutralino LSP case, theχ 0 thermal abundance freezes out at the temperature T f.o.
χ 0 and the Yχ0 receives extra contributions from the gravitino late decays. If the gravitino is heavy enough it can be T dec 3/2 > T f.o. χ 0 and the neutralinos from the gravitino decay equilibrate. Generally the neutralino abundance increases as m 2 χ 0 and the relic neutralino density is too large form and mχ0 beyond the TeV scale, see figure 4. Moreover TeV scale supersymmetry with neutralino LSP, although compatible with the Ω DM bound, is disfavoured for reheating temperatures T rh 10 8 GeV due to BBN constraints on the late decaying gravitino abundance [75,76].
One concludes that scenarios with reheating temperatures T rh >m and m LSP > O(TeV) are compatible only if late entropy production takes place. The above remarks are synopsised in the following conditions, wherem the sparticle mass scale. Hence, scenarios with high reheating temperature generally require an extra scalar field that causes dilution.

The diluter field X
In supersymmetric theories generically exist scalar fields with rather flat potentials and very weak or M Pl suppressed interactions. These kind of scalars, that are common in supergravity and superstring theories, are here collectively labeled X. The X domination, either due to its nearly constant potential energy or due to the energy stored in its oscillations about the vacuum, dilutes the LSP abundance D X times and supplements it with the contribution from the diluter decay where we labeled Ω < LSP the LSP abundance before the X decay. In order to specify the Ω LSP the system of the three interacting cosmic fluids of X, LSP and radiation has to be solved and we refer the reader to references [85,86,89] for detailed analytic results. For gravitino or axino LSP the above expression generally applies. For the neutralino LSP one should also check whether the conditions (i) T dec X < T f.o. χ 0 and (ii) nχ0 σv < H(T dec X ) hold. If not, then in the case (i) the neutralinos might reach a thermal equilibrium value Y (th) χ 0 . In the case (ii) pair annihilations take place until the neutralino yield reaches the value Y (th) χ 0 /T dec X ); this corresponds to the so-called annihilation scenario and works mostly for wino-like LSP with TeV mass scale. Let us mention here that the radiation produced from the decay of the X particles for the times Γ X /H < 1 can produce neutralinos even for T dec X < T f.o. χ 0 [85,89], which accounts for an extra contribution to Ω LSP that may be important in particular scenarios without, however, modifying the conclusions of the current analysis. Finally, the Ω X LSP depends on the branching ratio Br X LSP of the diluter into two LSPs (directly or via cascade decays) and the X decay temperature T dec X . The LSP yield from the X decay reads If the Y X LSP is subdominant the observed dark matter has to be produced by processes taking place at higher temperatures than T dec X and was appropriately diluted by the decay of the scalar X. On the other hand, if the dilution D X decreases the initial LSP abundance to negligible levels, then the LSP production from the X decay should fit the observed dark matter abundance. The constraint Ω LSP h 2 ≤ 0.12 implies which determines the dilution magnitude and consequently the shift in the spectral index (2.20). The D min X is referred as the required dilution throughout the text, necessary to give at most a critical density of LSP particles today.
The X decay is not free from constraints. It must decay before the BBN [75], not overproduce LSPs and not overproduce late decaying particles such as gravitinos. In the simple but quite unnatural case that the X is lighter than LSP then it is Br T rh 10 9 GeV E x c l u d e d Figure 5. The maximum possible dilution size, caused by a scalar X condensate, with respect to the LSP mass, for gravitino LSP (black, brown) and neutralino LSP (blue). We have made the conservative assumption m X m LSP that maximizes the diluter X lifetime. In the area above the lines it is Ω < LSP h 2 /D X > 0.12, hence it is an excluded parameter area. The solid and dashed lines correspond to c = 1 and c = 10 8 according to the parametrization (4.7). For a gravitationally decaying diluter (c=1), the thermal gravitino scenario (brown solid line) is excluded because the X spoils the BBN predictions. The plot demonstrates the decrease of the dilution efficiency for large supersymmetry breaking scale and T rh = 10 9 GeV. natural if m LSP < m X < 2m LSP since the channel X →GG orχ 0χ0 is forbidden due to kinematic constraints.
If the decay of the X produces LSPs or other late decaying particles the relevant branching ratios have to be considered. This is a model dependent issue and should be examined in the context of each model. In the next section we consider the supergravity R 2 inflation and we take into account the X decay rate and channels. Actually, the details of the X decay do not change any of the conclusions synopsised in the conditions (A) and (B). The minimum amount of dilution (4.6) is necessary regardless the diluter branching ratios, and this is a key point of this work.

The maximum possible dilution due to a scalar condensate
If the diluter mass is about or larger than the LSP mass, m X m LSP , then the dilution magnitude is correlated with the supersymmetry breaking scale. A late time entropy production takes place when the radiation dominated era gets interrupted by an X domiation era at T dom X < T rh , where T rh is the reheating temperature caused by the inflaton decay. For an oscillating scalar field the dilution magnitude is D X T dom X /T dec X . The decay rate of the X scalar can be parametrized as and the X decay temperature is T dec X (π 2 g * /90) −1/4 (Γ X M Pl ) 1/2 . For c ∼ 1 the X decays gravitationally and T dec X ∼ 4 MeV (M X /10 5 GeV) 3/2 . For c 1 non-gravitational decay channels exist; for example if the X field has Yukawa-like coupling y X to light degrees of JHEP02(2018)118 freedom then it is Γ X = y 2 X m X /8π. For the borderline case that T dom X = T rh and m X = m LSP the dilution magnitude due to an oscillating scalar field, D X , reaches a maximum value. Consequently, a minimum value for the Ω < LSP h 2 /D X exists which obviously must be below the observational value Ω DM h 2 = 0.12.
In particular, for gravitino LSP the lowest T dec where Ω < 3/2 = Ω MSSM(sc) 3/2 ,γ sc 1, see eq. (3.1), and the parameter c is explicitly written. Note that although the T rh is dropped out in the above relation it must be T rh > m 3/2 . The constraint (4.8) says that the abundance of gravitino LSPs produced from thermal scatterings in the plasma is possible to get diluted to observationally acceptable values by an oscillating scalar field that obtains mass from the supersymmetry breaking only if m 3/2 < 7 × 10 8 GeV . (4.9) The constraint becomes more severe ifγ sc 1, that is, if m 2 3/2 m 2 g < T 2 rh , or for a non-gravitational scalar X, c 1 or for m X m 3/2 . For thermalized gravitinos instead, the formula (3.4) applies and the maximum possible dilution magnitude gives the following constraint c 1/2 m 3/2 10 5 GeV 5/2 10 9 GeV T rh < where Ω < 3/2 = Ω eq 3/2 . We see from (4.10) that typical reheating temperatures T rh = 10 9 − 10 12 GeV imply a mass bound m 3/2 10 6 GeV for thermalized LSP gravitinos. Although such heavy gravitinos hardly get thermalized via interactions with the MSSM plasma, thermalized messengers can bring them to thermal equilibrium. Again here, the bound (4.10) becomes more severe for c 1 or for m X m 3/2 . When the neutralino is the LSP theχ 0 relic abundance is determined at the freeze out temperature that is T f.o.
χ 0 ∼ mχ0/20. If the decay temperature of the X field is below the T f.o.
χ 0 the neutralinos number density will get diluted. For thermally produced neutralinos the maximum possible dilution magnitude, for T dom X = T rh and m X ∼ mχ0, gives the constraint Thus, neutralino masses mχ0 > 10 7 GeV for T rh 10 9 GeV cannot be reconciled with cosmological scenarios where an oscillating scalar field dilutes the thermal plasma. If the leading contribution to Ω < χ 0 comes from the decays of gravitinos which are respectively produced from sfermion decays for mf > m 3/2 > mχ0 the expression (3.10) has to be used and another constraint for theχ 0 mass is obtained.

JHEP02(2018)118
This correlation among the dilution size due to scalar oscillations and the m LSP (or alternatively the supersymmetry breaking scale) is an extra and important constraint on these scenarios, see figure 5. The constraints on the LSP mass can be raised if thermal inflation takes place and then the dilution size is given by the expression (2.21). In such a case it is c 1 since the X is not a gravitationally decaying scalar due to the necessary presence of Yukawa couplings of X with the thermalized degrees of freedom, that regulate the decay rate. The gravitational diluter scenario is certainly a less model dependent and a more generic one.

A concrete example: the R + R 2 (super)gravity inflationary model
Inflation is the leading paradigm for explaining the origin of the primordial density perturbations that grew into the CMB anisotropies. If the early Universe is described by a typical model of inflation that naturally explains the statistical properties of the density and the expected tensor perturbations then the precision (n s , r) measurement can give us physical evidences for the radiation dominated era before the epoch of nucleosynthesis.
In this section we will apply the previously obtained results on the R 2 gravity and supergravity inflation models in order to perform a full estimation of the theoretically expected values for the n s and r observables. The Starobinsky R 2 inflation model is particularly motivated because it is placed in the center of the likelihood contour, nonetheless it is self evident that a similar analysis can be performed for any other inflation model. In the following, preliminaries of the R 2 gravity and supergravity inflation models will be reviewed. For the supergravity R 2 model new predictions for the (n s , r) observables will be derived, depending on the supersymetry breaking scheme, and the phenomenology of the two models will be compared.

The Starobinsky R 2 inflation
The Starobinsky model [90] is an f (R) gravity model described by the Lagrangian This theory is conformally equivalent to the Einstein gravity with a scalar field ϕ, the scalaron, minimally coupled to gravity From the CMB normalization [1] we get m 1.3×10 −5 M Pl . The inflationary predictions of the R 2 theory [91] at leading order are given by the following expressions of the primordial spectra and tensor-to-scalar-ratio r * = 16 * , Also, the tensor spectral tilt and running are respectively n t = −3/(2N 2 * ), dn t /d ln k −3/N 3 * .

JHEP02(2018)118
After the end of the inflationary expansion the inflaton is a homogeneous condensate of scalar gravitons. The scalaron universally interacts with other elementary particles only with gravitational strength and the inflaton perturbative decay process can be computed. The lifetime of the scalaron is rather long and ϕ decays after it has oscillated excessively many times about the minimum of its potential. The universe during scalaron oscillation phase evolves as a pressureless matter dominated phase and the effective value of the equation of state during reheating is to good approximation zero,w rh = 0, [37]. Thus the ∆N rh given by the expression (2.10) reads The energy density of the inflaton at the end of inflation is found to be ρ end Pl . The energy density at the end of reheating, ρ rh = (π 2 /30)g * rh T 4 rh , is determined by the reheating temperature T rh and the number of the degrees of freedom g * (T rh ) ≡ g * rh . In total, for the R 2 inflation the expression (2.12) is recast into The reheating temperature is estimated by equating Γ inf = H, where Γ inf ≡ Γ ϕ is the decay rate of the scalar graviton, Assuming only Standard Model degrees of freedom, at that energy scales it is g * rh = 106.75, thus T rh ∼ 10 9 GeV. For the R 2 we get for the first slow roll parameter * = (3/4)/N 2 * thus 1/4 ln * = −2.1 + 1/2 ln(54/N * ). In addition the R 2 plateau potential changes very slowly with the ϕ value and for N * = 45 − 60 it is 1/4 ln(V * /ρ end ) ≈ 0.2, hence In the above equation the logarithmic correction 1/2 ln(54/N * ) has been neglected because its value is less than 0.1 for relevant values of the N * . The thermal n (th) s value that the standard Starobinsky R 2 inflation model predicts at leading order is found when we substitute the thermal e-folds number N (th) = 54 into the eq. (5.3), that is n (th) s = 0.963. In terms of the e-folds number, the other two slow roll parameters for the Starobinsky model read η V −1/N and ξ V 1/N 2 . Since the corrections at second order in slow roll at the scalar tilt will not be negligible in the future it is crucial to go to order 1/N 2 . Also, going at next-to-leading order we could probe ∆N X ∼ 1 changes that could shed light on the pre-BBN cosmic history. For the Starobinsky model the expression (2.14) reads [33]

JHEP02(2018)118
Also, going to order 1/N 3 the tensor-to-scalar ratio and running read that is 2 larger than the leading order prediction. We also take at next-to-leading order r (th) R 2 = 0.0034 and α (th) Note that the r value is 17% smaller than the value obtained at leading order. Furthermore, going to accuracy level 1/N 3 the r = r(n s ) relation reads The eq. (5.12) was obtained from the expressions n s = n s ( V , η V , ξ V ) and r = r( V , η V ) written up to 1/N 3 order. In particular for the Starobinsky model it is If nature is successfully described by the Standard Model of particle physics and the R 2 inflation model then the ∆N X has to be zero and hence n s = n (th) s . Next we review and estimate the expected n s and r values for the R 2 supergravity inflation model.

The R + R 2 supergravity inflation
The embedding of the Starobinsky model of inflation in old-minimal supergravity in a superspace approach consists of reproducing the Lagrangian (5.1). This is achieved by the action [92][93][94][95][96] Modifications and further properties can be found in [97][98][99][100][101][102][103][104][105][106][107][108][109][110]. We mention that attention should be paid to the full couplings of the inflaton field that may yield a different reheating temeprature in each of these models since not all of them are pure supergravitational. The old-minimal supergravity multiplet contains the graviton (e a m ), the gravitino (G = ψ α m ), and a pair of auxiliary fields: the complex scalar M and the real vector b m . Lagrangian (5.13) when expanded to components yields R 2 terms and kinematic terms for the "auxiliary" fields M and b m . One may work directly with (5.13) but it is more convenient to turn to the dual description in terms of two chiral superfields: T and S and standard supergravity [92]. During inflation the universe undergoes a quasi de Sitter phase which implies that supersymmetry is broken, the mass of the sgoldstino S becomes large and it can be integrated out [111,112]. In this stage a non-linear realization of supersymmetry during inflation is possible [113][114][115][116]. The real component of T is not integrated out due to the non-linear realization and it is the only dynamic degree of freedom during inflation [93,94,96]. Eventually one finds the effective model (5.2).

JHEP02(2018)118
The inflationary predictions for the supergravity R 2 model are found to be identical to the non-supersymmetric Starobinsky R 2 predictions (5.3). In addition, the reheating phase is much similar and the inflaton decay rate roughly the same. Indeed, in the work of [117] the inflaton decay channels were identified and the branching ratios calculated. The total decay rate was parametrized as Γ sugra-inf = c m 3 Φ /M 2 Pl , where m Φ ≡ m inf and the reheating temperature was estimated to be (5.14) The fact that the reheating temperature is found to be about the same with that predicted in the non-supersymmetric R 2 model (5.6) means the supergravity and nonsupergravity versions of the R 2 inflation models are completely degenerate in terms of the inflationary predictions. However, the details of the expansion history of the universe after the decay of the inflaton should break the degeneracy between the supergravity-R 2 and gravity R 2 . We can directly apply the analysis and the results of the previous sections by minimally completing the supergravity R 2 sector with the MSSM and a basic supersymmetry breaking sector. Let us first examine the implications of the supergravity R 2 inflation to the abundances of superparticles.
The R 2 supergravity scenario can be distinguished in two basic cases: the ultra high scale supersymmetry breaking m 3/2 > m Φ and the sub-inflation supersymmetry breaking scale m Φ > m 3/2 case. The fist case is realized when the minimum of the inflationary potential breaks supersymmetry. Particularly in the model of [96], where a new class of Rsymmetry violating R + R 2 models was considered, it was found that it is possible inflation and supersymmetry breaking to originate from the supercurvature and obtain m 3/2 ∼ 2m Φ without invoking any matter superfields. The new properties of these models which distinguish them from the R-symmetric R 2 supergravity is that at the end of inflation the S field contribution becomes important. In such ultra high scale supersymmetry breaking scenarios the superparticles possibly play no rôle during the thermal evolution since the reheating temperature may not be sufficient to excite the superpartners of the Standard Model particles. Hence, the R 2 and supergravity R 2 models with m 3/2 > m Φ may be totally indistinguishable unless gauginos or some moduli fields are much lighter than the gravitino.
In the case that the inflaton field vacuum is supersymmetric an extra field is required to break supersymmetry, and the condition m 3/2 < m Φ is usually satisfied. The supersymmetry breaking spurion field called Z, i.e. the sgoldstino, although it can play the rôle of the diluter it overproduces LSPs and an extra scalar that we generically label X has to play the rôle of the diluter. 3 Assuming a simple supersymmetry breaking sector, with W SB = F Z + W 0 and K SB = |Z| 2 − |Z| 4 /Λ 2 , the gravitino yield due to the direct decay of the inflaton Φ is calculated to be [117] 3 Any late decaying scalar field, e.g stringy moduli, can be the diluter field.

JHEP02(2018)118
with branching ratio 16) The c is determined by the dominant decay channel, here the anomaly induced process [117]. In the case that the spurion field is heavier than the inflaton, m Φ m Z , and m 3/2 < m Φ the branching ratio maximizes, Br inf 3/2 (48πc ) −1 . Otherwise, the gravitino yield is calculated from the branching ratio (5.16) to be 17) For the supergravity R 2 inflation the above contribution to the gravitino abundance is small.
Apart from direct gravitino production from inflaton decays, gravitinos are produced via the decay of the Z field. The supersymmetry breaking field Z is produced as particles by the decay of inflaton with branching ratio (5.18) Considering the generic decay channel the Z decays dominantly into a pair of gravitinos when m 3/2 m Z < m Φ with the partial decay rate enhanced by the factor (m Z /m 3/2 ) 2 , Thus, the gravitino yield as a decay product of particle Z is found to be [117] where T rh the reheating temperature after the decay of the inflaton and Br Z 3/2 the branching ratio of the Z into a pair of gravitinos. In addition to the incoherent Z particles there are the coherent Z modes, produced by the inflationary de-Sitter phase, which may store a significant amount of energy. The precise VEV of Z is rather model dependent. TheG yield from the decay of the Z condensate can be computed if the initial amplitude of oscillations z 0 , the Z mass and couplings are known. Assuming a Z dominated universe it is and gravitinos and the LSPs are generally found to be overabundant. The initial value and zero temperature VEV of the scalar Z field are rather model dependent and it is possible

JHEP02(2018)118
that the Z scalar does not dominate the energy density of the universe. In the analysis of [117] the scalar Z is trapped near the origin during inflation. The zero temperature VEV, dictated by the Kähler, K SB , and the superpotential, W SB , is z = 2 √ 3(m 3/2 /m Z ) 2 M Pl . In the following we assume benchmark sparticle mass patterns and we estimate the corresponding shift in the spectral index and the tensor-to-scalar ratio in order the predicted dark matter density to be in accordance with observations. We assume gravitino and neutralino dark matter scenarios. We generally assume the presence of an extra scalar labeled X that dilutes the LSP abundance at the critical Ω LSP h 2 = 0.12 and sub-critical values. Particular hidden sector details concerning the X dynamics are left unspecified except for the requirement the diluter not to overproduce LSPs at the time of late entropy production. This is achieved is if the branching ratio to LSPs is very suppressed or m LSP < m X < 2m LSP .

The shift in the scalar spectral index and the tensor-to-scalar ratio for the Starobinsky R 2 inflation
The diluter X field dominates the energy density of the universe if T dom , for thermal inflation (5.22) the temperature that the energy density stored in the oscillating X field gets over the radiation energy density. The x 0 is the initial amplitude of the oscillations in a potential V (X) = m 2 X X 2 about the minimum and V 0 the vacuum energy of the flaton field X in the case of thermal inflation. The reheating temperature for Starobinsky and supergravity Starobinsky inflation is T rh ∼ 10 9 GeV and the decay temperature of the X field depends on its full interactions. The non-thermal X field domination induces a shift in the spectral index value n (th) s = 0.965 due to a change in the thermal e-folds number N (th) = 54. According to the formula (2.20) the size of the shift due to a non-thermal phase that lasts N X = [(1 − 3w X )/4] −1 ∆N X e-folds after the inflaton decay for the Starobinsky model is For a scalar condensate domination it is ∆N X = lnD X /3. The ∆n s depends on the dilution size D X plus a correctionĝ due to the change of the number of the effective degrees of freedom at the temperatures T dom X and T dec X . Keeping only the relevant terms and after analyzing lnD X = ln D X +ĝ the shift in the scalar tilt reads ∆n s (D X ,ĝ) = −2 × 10 −4 (ln D X +ĝ) 1 + 2 300 (ln D X +ĝ) , and T dec X , i.e. g * (T dom X ) = 10 g * (T dec X ), and the dashed line corresponds to no change, i.e. g * (T dom X ) = g * (T dec X ).

(5.26)
We have verified that the value r (th) + ∆r(D X ,ĝ), that the eq. (5.26) yields, agrees with 10 −4 precision with the value one gets from the relation r = r(n s ) given by the eq. (5.12) for n s = n (th) s + ∆n s (D X ,ĝ). Regarding the effective degrees of freedom, it isĝ O(1), hence the change in the number of degrees of freedom requires accuracy at the n s (r) measurement of the order of 10 −4 (10 −5 ) and one can safely neglect theĝ correction in the expressions (5.24) and (5.26) since the expected accuracy of the future CMB probes will be of the order of 10 −3 . Nevertheless observing that, in principle at least, one can additionally determine the number of the effective degrees of freedom at the thermal plasma from the (n s , r) precision measurement is certainly important and exciting, see figure 6.

The n s and r predictions for particular supersymmetry breaking examples
In this subsection we explore the impact on (n s , r) observables of the two base case dark matter scenarios of supersymmetry, the gravitino and the neutralino, when the initial

JHEP02(2018)118
conditions for the hot Big Bang are set by the supergravity Starobinsky inflation. We consider both thermal and non-thermal dark matter production from the hot plasma and scalar decays. We examine different and illustrative supersymmetry breaking schemes and we quantify how the expected values for the inflationary observables change due to a nonthermal post-reheating phase dictated by the universal constraint Ω LSP h 2 ≤ 0.12. We mention that this analysis, that probes cosmologically a BSM scheme, can be applied to any other inflationary model after the appropriate adjustments regarding the reheating phase, the reheating temperature and the inflaton field branching ratios.
Example I: gravitino dark matter. The gravitino is the LSP if the supersymmetry breaking is mediated more efficiently to the MSSM than to the supergravity sector. The standard paradigm is the gauge mediation scenario [49]. In such a scenario the supersymmetry breaking Z field decays dominantly into MSSM fields with non-gravitational interactions. Following realistic models [66][67][68], it is the imaginary part of the Z field that decays last and the dominant channel is onto a pair of gauginos, in particular binos, with the decay temperature given by The LSP gravitinos are produced from thermal scatterings and decays in the plasma and from the non-thermal decay of the Z scalar field and the inflaton. The inflaton contribution to the gravitino abundance is given by eq. (5.17) and in general is found to be subleading in R 2 inflation. The decay rate of the Z scalar to gravitinos is given by eq. (5.19). If the Z decay produces late entropy then the gravitinos from the Z decay, with branching ratio Br Z 3/2 will be part of the dark matter in the universe with yield Y Z 3/2 ∼ (3/2)Br Z 3/2 T dec Z /m Z and relic density parameter [67] Ω Z where mg the mass of the bino. We mention that it is also possible that the spurion field does not dominate the energy density due to thermal effects [56,118]. Before proceeding with the survey of particular examples, let us mention that the gravitino relic density parameter violates the observational bound unless the sparticles lay in the TeV and sub-TeV scale. Another scalar field X is required to dilute the thermally produced gravitinos and the energy stored in the oscillations of the supersymmetry breaking field, in case of Z domination. In order the precise dilution size to be determined the knowledge of the m Z , m 3/2 and the MSSM mass pattern is necessary.
Let us now consider four benchmark mass patterns for the supersymmetry breaking sector plus the MSSM, with different sizes of supersymmetry breaking scale. We also consider the presence of messenger fields and the diluter X field necessary to decrease the LSP relic density and which dominantly decays to visible sector fields and not to gravitinos.

JHEP02(2018)118
1. m 3/2 10 2 GeV, mf ∼ mg ∼ m Z 10 4 GeV and M mess 10 8 GeV. The messengers get thermalized since M mess < T rh and the scalar spurion field Z follows the finite temperature minimum without sizable oscillations and hence does not dominate the energy density [56,118]. We also assume that the messenger coupling is small enough, λ mess 1, so that the gravitinos do not get thermalized. The gravitinos produced from scatterings of thermalized messengers would have a relic density parameter Ω < 3/2 h 2 ∼ 10 4 , see below eq. (3.4). The Ω 3/2 h 2 ≤ 0.12 bound implies that the thermally produced gravitinos are sufficiently diluted if D X 10 4 . This dilution can be caused by scalar condensate X with D X T dom X /T dec X . The shift in the spectral index and tensor-to-scalar ratio are respectively |∆n s | 2 × 10 −3 and ∆r 4 × 10 −4 .
2. m 3/2 10 3 GeV, mf ∼ mg ∼ m Z > m 3/2 and M mess < T rh . Messengers get thermalized and Z does not dominate the energy density of the universe. The gravitinos obtain a thermal equilibrium abundance due to interactions with the thermalized messengers [56] and their relic density would be Ω < 3/2 h 2 ∼ 10 10 , see eq. (3.4). The Ω 3/2 h 2 ≤ 0.12 bound implies that the thermally produced gravitinos are sufficiently diluted if D X 10 10 . The diluter can be either a flaton field that causes thermal inflation or a scalar condensate. In the later case the X field dominates the energy density of the universe shortly after the reheating in order such a dilution size to be realized. The shift in the spectral index and tensor-to-scalar ratio are respectively |∆n s | 5 × 10 −3 and ∆r 10 × 10 −4 .
3. m 3/2 10 4 GeV, mg 10 5 GeV , mf ∼ m Z 10 6 GeV and M mess > T rh . The Z field does not receive thermal corrections because the messengers are not thermalized. The Z scalar oscillations generally have a large enough amplitude and Z does dominate the energy density of the universe. Equations (5.27) and (5.28) say that the spurion Z decays at T dec Z 1 GeV and produces non-thermally gravitinos that exceed about 10 6.5 times the observational bound. In order the Z condensate to get diluted the X field has to be a flaton and cause thermal inflation. In this case, the shift in the spectral index and tensor-to-scalar ratio are respectively |∆n s | 3 × 10 −3 and ∆r 7 × 10 −4 . 4. m 3/2 = few GeV, mg ∼ mf ∼ m Z = few TeV. There are scenarios in the literature that reconcile gravitino cosmology with high reheating temperatures [26,53,54,56] and generally assume non-minimal features for the hidden sector. For example when the messengers masses lay in the range M mess 10 6 GeV and the goldstino does not reside in a single chiral superfield [56], or when the messenger coupling is controlled by the VEV of another field [26] it is possible that gravitinos have the right abundance. These supersymmetry breaking schemes do not require dilution and predict ∆n s = 0 and ∆r = 0. We mention that these scenarios, in their original versions, work better when supersymmetry is broken about the TeV scale. Features of these scenarios are currently tested by the LHC experiments.  Table 1. The n s and r prediction for gravitino LSP and a gauge mediation scheme for the R 2 supergravity model. In the cases # 1, 2 and 4 the gravitinos are produced from thermal scatterings of messengers and MSSM fields while in the case # 3 from the non-thermal decay of the supersymmetry breaking Z field. In cases # 1, 2 and 3 dilution is required to decrease the LSP abundance below the observational bound. In the case # 4 non-minimal hidden sector features have been assumed. The masses are in GeV units.
Example II: neutralino dark matter. For gravity or anomaly mediation of supersymmetry breaking the gravitino mass is naturally heavier than the neutralinos. The gravitino decay populates the universe with neutralinos. Here we assume the gravitino mass to be above 10 5 GeV not to spoil BBN predictions at the time of decay. The gravitinos are produced non-thermally by the decay of the inflaton, see eq. (5.17), which generally accounts for a subleading contribution in the framework of R 2 supergravity inflation, and by the decay of the supersymmetry breaking scalar field Z. Contrary to the GMSB case the Z scalar oscillations are not thermally damped and generally the Z produces late entropy if displaced from the zero temperature minimum. The temperature that the Z field decays is estimated by considering the various partial decay rates. The dominant decay channel is into a pair of gravitinos, when m Z m 3/2 , and the total decay rate yields the decay temperature If the Z field oscillations dilute the thermal plasma then the gravitinos coming from the Z decay are the leading source of dark matter neutralinos at the gravitino decay temperature T dec 3/2 . The neutralinos are generally found to be overabundant when supersymmetry breaks at energies beyond the TeV scale and dilution is required. Hence we assume the presence of a diluter field X that decreases the LSP relic density via late entropy production. We mention that according to the general constraint (4.11) the neutralinos with mass mχ0 > 10 7 GeV are impossible to get diluted by the oscillations of the X scalar and thermal inflation is required.
Let us now consider benchmark mass patterns for the supersymmetry breaking sector plus the MSSM, characterized mainly by split and quasi-natural sparticle mass spectrum.

mχ0
10 3 GeV, m 3/2 ∼ mf 10 6 GeV, m Z 10 7 GeV. Here we assume the annihilation scenario where the neutralino has an annihilation cross section few orders of magnitude higher that the conventional value. The universe is generally dominated by the Z scalar that decays to gravitinos at the temperature T dec Z ∼ 12 GeV. In turn, the gravitinos produced from the Z decay dominate the energy density and decay at T dec 3/2 ∼ 0.2 GeV producing neutralinos that annihilate rapidly and acquire a relic density Ω (ann) The resulting LSP abundance can fit the observed value and here the rôle of the diluter is played by the Z field and the gravitinos. We note that if mχ0 > TeV then a diluter scalar X is required. It is D 10 2 but the dilution size due to Z oscillations can be many orders of magnitude larger. This scenario is currently tested and constrained by the LHC and indirect detection experiments [27]. This minimum value of the dilution magnitude yields |∆n s | 1 × 10 −3 and ∆r 2 × 10 −4 .

mχ0
10 3 GeV, m 3/2 ∼ mf 10 8 GeV, m Z 10 9 GeV. In this example we assume that the T dec 3/2 > T f.o. χ 0 and neutralinos thermalize after the decay of gravitinos. A Z dominated early universe becomes in turn gravitino dominated at T dec Z ∼ 10 4 GeV. For TeV and sub-TeV scale neutralinos the observational bound Ω DM h 2 ∼ 0.12 can be satisfied, see eq. (3.9). Here again the rôle of the diluter is played by the Z field and the gravitinos and it is D 10 2 , but it can many orders of magnitude larger. This scenario is currently tested by LHC and direct detection experiments. This dilution magnitude induces a shift in the spectral index and tensor-to-scalar ratio respectively at least of size |∆n s | 1×10 −3 and ∆r 2×10 −4 .

mχ0
10 5 GeV, m 3/2 ∼ mf 10 7 GeV, m Z 10 8 GeV. In this example the neutralinos are produced from the gravitino decay and they are out of chemical and kinetic equilibrium. The LSP relic density, given by eq. (3.10), is Ω < χ 0 h 2 ∼ 10 8 . The LSP abundance has to be decreased eight orders of magnitude down and this is possible only if the gravitinos and the Z scalar condensate are sufficiently diluted. Thermal inflation is required with D X 10 8 . This dilution magnitude induces a shift in the spectral index and tensor-to-scalar ratio respectively at least of size |∆n s | 4 × 10 −3 and ∆r 8 × 10 −4 . This is a phenomenologically viable example not constrained by terrestrial experiments.

mχ0
10 3 GeV, m 3/2 ∼ mf ∼ m Z > 10 5 GeV. As a last example we consider the conventional thermal WIMP scenario assuming that the Z scalar field is not displaced from the zero temperature minimum and never dominates the energy density of the universe, hence no non-thermal phase is required (although a non-thermal JHEP02(2018)118 χ 0 is not ruled out in general). In this scenario it is ∆n s = 0 and ∆r = 0. This supersymmetry breaking scheme is currently tested at LHC, direct and indirect detection experiments.
The above benchmark examples for the neutralino dark matter scenario are synopsized in the table 2 and figure 8.
Let us finally note that the LSP particles produced from the gravitino and X decay are warmer than the LSPs produced from thermal scatterings and this changes the free streaming length of the LSP dark matter, which has the effect of potentially washing out small scale cosmological perturbations, see e.g. [120,121]. This is a very interesting possibility that could provide further constraints to these scenarios, though the mass scales and lifetimes considered here yield free streaming lengths that are not in conflict with the Lyman-α forest observations [122].

Distinguishing the R 2 and the R 2 supergravity inflationary models
The supersymmetric and non-supersymmetric R 2 inflation models predict the same reheating temperature, T rh ∼ 10 9 GeV and the same expressions for the n s = n s (N ) and r = r(N ). However, the degeneracy between the two models that appears during the accelerating and the reheating stage breaks after the inflaton decay. 4 In the case of supergravity R 2 inflation, ifm < T rh , sparticles will be constituents of the thermalized plasma of the reheated universe. In addition to thermal processes, the presence of the supergravitational inflaton and the supersymmetry breaking field produce a significant number of gravitino particles after inflation, as the expressions (5.17), (5.20) and (5.21) make manifest.
The BBN and the Ω DM h 2 = 0.12 constraints imply that the thermal cosmic history is influenced by the change of the supersymmetry breaking pattern. The LSP is found to be overabundant in the greatest part of the MSSM parameter space and it receives further contributions when the supersymmetry breaking field is taken into account [117] for both gravitino and neutralino LSP scenarios. Therefore, R 2 supergravity inflation is compatible with the cosmological observations only if the thermal history of the universe is not perpetual from T rh ∼ 10 9 GeV until T BBN ∼ 1 MeV. A non-thermal phase that dilutes the supersymmetric thermal relics and potentially supplements the universe with dark matter particles can fully reconcile the R 2 supergravity inflation model with observations. The required dilution generally increases with increasing the sparticle masses. Henceforth, we conclude that the degeneracy breaking of the inflationary predictions between the R 2 and supergravity R 2 models depends on the energy scale and the pattern of supersymmetry breaking, see figure 7 and 8.
The fact that the R 2 supergravity automatically alters the details of the thermal history and possibly the expansion history of the universe compared to the simple R 2 case, where sparticles and supersymmetry breaking fields are absent, allows the discrimination between the two inflationary models. Considering only the MSSM degrees of freedom as the less model dependent and conservative analysis, the conditions (A) and (B) of section 4, when JHEP02(2018)118 Figure 7. This is a compound plot consisting of 3D graph and a density-contour plot. The 3D graph shows the decadic logarithm of the required dilution magnitude as a function of the gravitino LSP and gaugino-sfermion masses with mf = mg for a reheating temperature T rh = 10 9 GeV. The dilution is calculated by requiring the Ω 3/2 h 2 not to exceed the observational bounds. The densitycontour plot demonstrates the change in the n s value, magnified 1000 times on the contour labels, for inflationary models that predict a reheating temperature T rh = 10 9 GeV. The information that one extracts from this graph is that supersymmetric models (e.g. quasi-natural, split, high scale) can be compatible with the CMB data only for particular values for the scalar tilt n s . true, imply that the supergravity R 2 inflationary model predicts Tensor to scalar ratio Figure 8. Constraints on the (n s , r) contour plane from Planck-2015 in the pink, and the schematic illustration of 2σ forecast constraints from a future CMB probe with sensitivity δn s ∼ 10 −3 and δr ∼ 10 −3 depicted with the dotted and dashed ellipsis. The R 2 model is targeted with a fiducial value of r ∼ 4×10 −3 . The red asterisks correspond to the predictions of the four benchmark models (#1, 2, 3, 4) with gravitino LSP and the green asterisks to the four benchmark models (#1, 2, 3, 4) with neutralino LSP, as explained in the text and tables 1 and 2 respectively. If the future CMB experimental probes select the area inside the dashed ellipsis then either the R 2 or the SUGRA-R 2 inflation model is selected plus a roughly continuous thermal phase with reheating temperature, T rh ∼ 10 9 GeV. The selection of the dashed ellipsis area will exclude a large class of supersymmetry models that predict a too large LSP abundance for that reheating temperature. On the contrary, if the dotted ellipsis area is selected then the duration of the thermal phase before the BBN is much limited and extra scalar particles should be present above the TeV scale, hence supporting the SUGRA-R 2 model rather than the R 2 inflation model plus "desert".
Let us also mention that a similar result to (5.30) can be obtained if one simply assumes the presence of extra scalars that dominate the energy density of the early Universe. For example, for a supergravity R 2 inflation model, a gravitational modulated reheating was assumed [123] and non-Gaussianity was additionally predicted. In the present work the postulation of a non-thermal phase has been motivated by the general requirement to fit the universal constraint Ω LSP h 2 ≤ 0.12 that in turn implies the result (5.30). Last but not least, we emphasize again that the precise measurement of the (n s , r) cannot "prove" or disprove the existence of supersymmetry. It will only indicate the presence of extra scalar degrees of freedom, that supersymmetry or any other BSM scenario will be challenged to explain. . This graph demonstrates the analysis followed in this work to cosmologically probe BSM scenarios.

Discussion and conclusions
The cosmic energy window from about 1 MeV up to the inflationary energy scale is shuttered to the current observational probes and the corresponding timescale can be reasonably called a dark early universe cosmic era. Any understanding of the cosmic processes that take place before the BBN will provide us with critical insights into the microphysics that operates at that energy scales. One significant prospect to contemplate the early dark cosmic era is through the precision measurement of the CMB observables n s and r. The essential fact is that the (n s , r) are not strictly scale invariant, hence important information about the background expansion rate and the reheating temperature of the universe can be obtained.
The inflationary paradigm can be used as a concrete and compelling framework for the theoretical determination of the n s and r values. However, our ignorance about the reheating process and the subsequent evolution of the universe, encoded in the dependence on N * , is rather strong and will become significant as the accuracy on the observations are expected to be further improved the next decades. An inflationary prediction that is independent of N * is the contour line r = r(n s ) which can distinguish different inflation models. Furthermore, if inflation is followed by a continuous thermal phase then a concrete inflation model predicts a specific number for the number of e-folds between the moment the relevant modes exit the horizon and the end of inflation, hence predicts a specific spot on the r = r(n s ) line that corresponds to what we called thermal values for the e-folds number, scalar tilt and tensor-to-scalar ratio, N (th) , n (th) s and r (th) respectively. Motivated by the advertised sensitivity of the future CMB probes in this paper we quantified the effect of a generic primordial non-thermal phase on the spectral index value (2.20). The n s value is possible to have been shifted by the amount ∆n s /n s ∼ O(1 − 6) from the expected thermal value, n (th) s , due to a scalar condensate or a flaton JHEP02(2018)118 field domination. The observation of non zero ∆n s and ∆r along a contour line r = r(n s ) is an indirect observation of a non-thermal phase and connects cosmology to microphysics since it has to be attributed to a BSM scalar field domination.
Moving a step further we applied our general results to study the observational consequences on the CMB of a supersymmetric universe. Supersymmetry is one of the most motivated theories that is extensively used to describe the very early universe evolution.
Although it lacks any experimental support, it provides an appealing framework that consistently accommodates high energy processes such as inflation and dark matter production. Actually, the fact that the LHC probes only a small part of the vast energy scales up to the Planck mass, while supersymmetry may lay anywhere in between, strongly motivates the systematic cosmological examination of supersymmetric scenarios. Supersymmetry can be cosmologically manifest if supersymmetric degrees of freedom get thermally excited or produced non-thermally during the dark early cosmic era.
The most direct cosmological implication of supersymmetry is that the LSP expected to be stable and hence contributes to the dark matter density. The LSP abundance is the key quantity that we estimate in different classes of supersymmetry breaking schemes and examine how it can be cosmologically reconciled with the observational value Ω DM h 2 = 0.12. We find that a non-thermal phase or low reheating temperatures are generally required if supersymmetry UV completes the Standard Model of particle physics. We quantified the effect of the different expansion histories on the (n s , r) and we broadly related it with the different supersymmetry breaking schemes. In this paper we mostly focused on ultra-TeV scale supersymmetry since low scale supersymmetry models with thermal WIMPs are in growing conflict with collider data and direct detection experiments. In our analysis we have not assumed that the LSP accounts for the bulk dark matter component in the universe. If it is actually Ω LSP h 2 0.12 then the expected change in the (n s , r) values due to a non-thermal stage becomes greater.
A complete understanding of the pre-BBN thermal phase and the CMB observables requires the knowledge of the initial condition for the thermal Big Bang, which are successfully provided by the inflationary theory. In this work we suggested a unified study of inflation and the subsequent reheating stage. Actually it is often the case that supergravity inflationary models are degenerate, in terms of the inflationary observables, with their the non-supersymmetric versions. However, the supersymmetric degrees of freedom can be excited either thermally or non-thermally after the end of the inflationary phase. For the sake of completeness we considered in this paper the R 2 supergravity inflation and we performed a theoretical estimation of the (n s , r) observables. Our findings point out that the ultra-TeV scale supersymmetry leaves a more clear cosmological imprint on the CMB observables. This fact is particularly exciting because high scale supersymmetric scenarios can be cosmologically falsified while the low mass range supersymmetric scenarios are directly tested at the terrestrial colliders.
Undoubtedly any non-trivial cosmological information about the BSM physics is of major importance. Certainly the results of this cosmological analysis, illustrated in the graph of figure 9, cannot discover or disprove supersymmetry. The only concrete cosmological information that we get from the n s and r observables concerns the expansion rate of the JHEP02(2018)118 very early universe. The identity of the matter content that controls the cosmic expansion rate cannot be revealed and it is only subject to interpretations. Nonetheless, if the (n s , r) deviate from their thermal values then new physics exists in high energies. In this paper we focused on supersymmetry, though any BSM scenario can be analyzed accordingly. In the event of detection of primordial gravitational waves, that is an observation of r = 0 together with possible features of the tensor power spectrum, then the selection of a particular inflation model is possible. In such a case our analysis has the power to rule out the BSM desert scenario and indicate possible features of candidate BSM theories, as the figure 8 illustrates.
From the theoretical side a more complete analysis should also take into account baryognesis scenarios and the details of thermalization process. The generation of the matter-antimatter asymmetry in the universe, seems to have a critical dependence on the temperature, as e.g. the thermal leptogenesis scenario [124] suggests. Moreover the understanding of several distinct stages in the reheating process that leads to thermalization of the universe in a radiation dominated phase at some reheating temperature T rh is necessary in order a more accurate value for the equation of state parameter w rh and the reheating e-folds numberÑ rh to be estimated, see eq. (2.10). A thorough understanding of the reheating process can also bring out new observables that can further constrain the reheating temperature of the universe, see e.g. [125] for a review. We should mention here that the oscillatory epoch and the reheating process of the R 2 inflation model is well understood, a fact that makes the results obtained in section 5 reliable [37].
From the observational side, future CMB primary anisotropy measurements should play a decisive rôle in probing the pre-BBN cosmic era. Complementary observational programs, such as the direct observation of tensor perturbations, should contribute significantly to this endeavor as well. Information on the thermal history after inflation is imprinted in the gravitational wave spectrum in the frequencies corresponding to the reheating energy scales, which can be probed by future space-based laser interferometers such as DECIGO [126]. Presumably, the synergy of different cosmological surveys will enable a leap forward in precision cosmology giving us, at the same time, access to the physics that operates beyond the Standard Model of particle physics, at energy scales much higher than can be obtained at CERN.