Gravitino and Polonyi production in supergravity

We study production of gravitino and Polonyi particles in the minimal Starobinsky-Polonyi $\mathcal{N}=1$ supergravity with inflaton belonging to a massive vector supermultiplet. Our model has only one free parameter given by the scale of spontaneous SUSY breaking triggered by Polonyi chiral superfield. The vector supermultiplet generically enters the action non-minimally, via an arbitrary real function. This functon is chosen to generate the inflaton scalar potential of the Starobinsky model. Our supergravity model can be reformulated as an abelian supersymmetric gauge theory with the vector gauge superfield coupled to two (Higgs and Polonyi) chiral superfields interacting with supergravity, where the $U(1)$ gauge symmetry is spontaneously broken. We find that Polonyi and gravitino particles are efficiently produced during inflation. After inflation, perturbative decay of inflaton also produces Polonyi particles that rapidly decay into gravitinos. As a result, a coherent picture of inflation and dark matter emerges, where the abundance of produced gravitinos after inflation fits the CMB constraints as a Super Heavy Dark Matter (SHDM) candidate. Our scenario avoids the notorous gravitino and Polonyi problems with the Big Bang Nucleosynthesis (BBN) and DM overproduction.


Introduction
The Planck data [1,2,3] of the Cosmic Microwave Background (CMB) radiation favors slow-roll single-large-field chaotic inflation with an approximately flat plateau of the scalar potential, driven by single inflaton (scalar) field. The simplest geometrical realization of this description is provided by the Starobinsky model [4]. It strongly motivates us to connect this class of inflationary models to particle physics theory beyond the Standard Model (SM) of elementary particles. A reasonable way of theoretical realization of this program is via embedding of the inflationary models into N = 1 supergravity. It is also the first natural step towards unification of inflation with the Supersymmetric Grand Unified Theories (SGUTs) and string theory. Inflaton is expected to be mixed with other scalars, but this mixing has to be small. The inflationary model building in the supergravity literature is usually based on an assumption that inflaton belongs to a chiral (scalar) supermultiplet, see e.g., the reviews [5,6] for details. However, inflaton can also belong to a massive N = 1 vector multiplet instead of a chiral one. Since there is only one real scalar in a massive N = 1 vector multiplet, there is no need of stabilization of its (scalar) superpartners, and the η-problem does not exist because the scalar potential of a vector multiplet is given by the Dterm instead of the F -term. The minimal supergravity models with inflaton belonging to a massive vector multiplet were proposed in Refs. [7,8] by using the non-minimal self-coupling of a vector multiplet, paramaterized by an arbitrary real function [9]. These models can accommodate any desired values of the CMB observables (the scalar tilt n s and the tensor-to-scalar ratio r), because the corresponding single-field (inflaton) scalar potential is given by the derivative squared of that (arbitrary) real function. However, all models of Refs. [7,8] have the vanishing vacuum energy after inflation, i.e. the vanishing cosmological constant, and the vanishing vacuum expectation value (VEV) of the auxiliary fields, so that supersymmetry (SUSY) is restored after inflation and only a Minkowski vacuum is allowed. A simple extension of the models [7,8] was proposed in Refs. [10,11], where a Polonyi (chiral) superfield with a linear superpotential [12] was added to the action, leading to a spontaneous SUSY breaking and a de-Sitter vacuum after inflation.
A successful theoretical embedding of inflation into supergravity models is, clearly, a necessary but is not a sufficient condition. Even when these models are well compatible with the Planck constraints on the (r, n s ), they still may (and often, do) lead to incompatibility with the (Hot and Cold) Dark Matter (DM) abundance and the Big Bang Nucleosynthesis (BBN). A typical issue is known as the gravitino problem: in order to not ruin the BBN, gravitinos must not decay in the early thermal bath injecting hadrons or radiation during the BBN epoch [13,14,15,16]. As is well known, the BBN is very sensitive to initial conditions, while each extra hadron or radiation can radically jeopardize the BBN picture, leading to disastrous incompatibility with cosmological and astrophysical data. In addition, the so-called Polonyi problem was also pointed out in the literature: Polonyi particles decay can also jeopardize the success of the BBN [17,18,19,20,21,22]. Indeed, a generic supergravity model predicts a disastrous overproduction of gravitinos and/or Polonyi particles or neutralinos. Any specific predictions are model-dependent, because they are very sensitive to the mass spectrum and the parameter range under consideration. The mass pattern selects the leading production mechanism or channel: either thermal (WIMP-like) production or/and non-thermal production sourced by inflation and decays of other heavier particles. The last channel includes a possible (different) production mechanism due to evaporating Primordial Black Holes (PBH's) that may be formed in the early Universe [23,24,25]. Our minimal estimation of the probability of the mini-PBH's formation at the long dust-like preheating stage after inflation gives such a low value that the successive evaporation of the mini-PBH's doesn't lead to a significant contribution to gravitino production (Sec. 3). However, if inflation ends by a first order phase transition, the situation drastically changes, and copious production of mini-PBHs in bubble collisions [26,27] can lead to a huge gravitino overproduction, thus excluding the first order phase transition exit from inflation.
We consider the very specific, minimalistic and, hence, a bit oversimplified N = 1 supergravity model of inflation, with inflaton belonging to a massive vector multiplet. We demonstrate that our model avoids the overproduction and BBN problems, while it naturally accounts for the right amount of cold DM. We assume both Polonyi field, triggering a spontaneous SUSY breaking at high scales, and the massive gravitino, produced during inflation, to be super-heavy, and call it the Super-Heavy Gravitino Dark Matter (SHGDM) scenario. A production of super-heavy scalars during inflation was first studied in Refs. [28,29], whereas the gravitino production sourced by inflation was considered in Refs. [30,31,32], though without specifying a particular model. In this paper we apply the methods of Refs. [30,31,32] to the specific Starobinsky-Polonyi supergravity model proposed in Refs. [10,11]. The supersymmetric partners of known particles (beyond the ones present in the model) are assumed to be heavier than Polonyi and gravitino particles (in the context of High-scale SUSY), also in order to overcome several technicalities in our calculations. Some of the physical predictions of our model are (i) the Polonyi mass is a bit higher than two times of the gravitino mass, and (ii) the inflaton mass is slightly higher than two Polonyi masses. This implies that inflaton can decay into Polonyi particles that, in their turn, decay into a couple of gravitinos. We show that super-heavy gravitinos produced from inflation and Polonyi decays can fit the cold DM abundance. The parameter spaces of inflation and cold DM are thus linked to each other in a coherent unifying picture.
Our paper is organized as follows. In Sec. 2 we review our model. In Sec. 3 we consider the gravitino and Polonyi particle production mechanisms. Sec. 4 is our conclusion and outlook.

The model
In this Section we briefly review the inflationary model of Refs. [10,11]. We use the natural units with the reduced Planck mass M P l = 1. 1 The model has two chiral superfields (Φ, H) and a real vector superfield V , all coupled to supergravity, and having the Lagrangian in terms of the chiral scalar curvature superfield R, the chiral density superfield E and the superspace covariant spinor derivatives (D α , D • α ), a Kähler potential K = K(Φ, Φ) and a superpotential W(Φ), the abelian (chiral) superfield strength W α ≡ − 1 4 (DD − 8R)D α V , and a real function J = J(He 2gV H) with the coupling constant g.
The Lagrangian (1) is invariant under the supersymmetric U (1) gauge transformations whose gauge parameter Z itself is a chiral superfield.
The chiral superfield H can be gauged away via the gauge fixing of these transformations by a gauge condition H = 1. Then the Lagrangian (1) gets simplified to After eliminating the auxiliary fields and moving from the initial (Jordan) frame to the Einstein frame, the bosonic part of the Lagrangian (4) reads [10] with the scalar potential where we have used the supergravity (bosonic) field components defined by 3 2E| = e, DD(2E)| = 4eM , in terms of the vierbein determinant e ≡ dete a m , the spacetime scalar curvature R, and the minimal set of the supergravity auxiliary fields given by the complex scalar M and the real vector b m . The matter (bosonic) field components are defined by in terms of the physical fields (A, C, B m ), the auxiliary fields (F , X, D) and the vector field strength F mn = D m B n − D n B m . When the function J is linear with respect to its argument (i.e. in the case of the minimal coupling of the vector multiplet to supergravity), our results agree with the textbook [33]. In the absence of chiral matter, Φ = 0, our results also agree with Refs. [8,9]. 4 As is clear from Eq. (5), the absence of ghosts requires J (C) > 0, where the primes denote the differentiations with respect to the given argument. In this paper, we restrict ourselves to the Kähler potential and the superpotential of the Polonyi model [12]: with the parameters µ and β. Unlike Ref. [34], we do not impose the nilpotency condition Φ 2 = 0, in order to keep manifest (linear) supersymmetry of our original construction (1) and to avoid a concern about loosing unitary with the nilpotent superfields at high energies.
Then, on the one side, our model includes the single-field (C) inflationary model, whose D-type scalar potential is given by V (C) = g 2 2 (J ) 2 (8) in terms of an arbitrary function J(C), with the real inflaton field C belonging to a massive vector supermultiplet. On the other hand, the Minkowski vacuum conditions (after inflation) can be easily satisfied when J = 0, which implies [12] This solution describes a stable Minkowski vacuum with spontaneous SUSY breaking at arbitrary scale F = µ. The related gravitino mass is given by There is also a massive (Polonyi) scalar of mass M A = 2µe 2− √ 3 and a massless fermion in the physical spectrum.
As regards early Universe phenomenology, our specific model of Polonyi-Starobinsky (PS) supergravity has the following theoretically appealing features: • there is no need to "stabilize" the single-field inflationary trajectory against scalar superpartners of inflaton, because our inflaton is the only real scalar in a massive vector multiplet, • any values of CMB observables n s and r are possible by choosing the J-function, • a spontaneous SUSY breaking after inflation takes place at arbitrary scale µ, • there are only a few parameters relevant for inflation and SUSY breaking: the coupling constant g defining the inflaton mass, g ∼ m inf. , the coupling constant µ defining the scale of SUSY breaking, µ ∼ m 3/2 , and the parameter β in the constant term of the superpotential. Actually, the inflaton mass is constrained by CMB observations as m inf . ∼ O(10 −6 ), while β is fixed by the vacuum solution, so that we have only one free parameter µ defining the scale of SUSY breaking in our model (before studying reheating and phenomenology).
The (inflaton) scalar potential associated with the Starobinsky inflationary model of (R + R 2 ) gravity arises when [8] J that implies According to (5), a canonical inflaton field φ (with the canonical kinetic term) is related to the field C by the field redefinition C = exp 2/3φ .
Therefore, we arrive at the (Starobinsky) scalar potential The full action (1) of this PS supergravity in curved superspace can be transformed into a supergravity extension of the (R + R 2 ) gravity action by using the (inverse) duality procedure described in Ref. [8]. However, the dual supergravity model is described by a complicated higherderivative field theory that is inconvenient for studying particle production. Actually, there is also the F-type scalar potential in PS supergravity due to mixing of inflaton and Polonyi scalars, that leads to instability of the Starobinsky inflation described by the D-term alone. However, after adding a Fayet-Iliopoulos (FI) term [35] together with its supersymmetric completion [36] to the Lagrangian (1) and modifying the J-function above, the Starobinsky inflation can be restored, and the inflaton-Polonyi mixing can be suppressed [37]. The FI term does not affect the phenomenology discussed in this paper, as long as J is negative and close to zero [37], as we always assume.
Another nice feature of our model is that it can be rewritten as a supersymmetric (abelian and non-minimal) gauge theory coupled to supergravity in the presence of a Higgs superfield H, resulting in the super-Higgs effect with simultaneous spontaneous breaking of the gauge symmetry and supersymmetry. Indeed, the U (1) gauge symmetry of the original Lagrangian (1) allows us to choose a different (Wess-Zumino) supersymmertic gauge by "gauging away" the chiral and antichiral parts of the general superfield V via the appropriate choice of the superfield parameters Z and Z as Then the bosonic part of the Lagrangian in terms of the superfield components in Einstein frame, after elimination of the auxiliary fields and Weyl rescaling, reads [11] where h,h are the Higgs field and its conjugate.
The standard U (1) Higgs mechanism arises with the canonical function J = 1 2 he 2Vh , where we have chosen g = 1 for simplicity. As regards the Higgs sector, it leads to After changing the variables h andh as where ρ is the (real) Higgs boson, ν ≡ h = h is the Higgs VEV, and ζ is the Goldstone boson, the unitary gauge fixing of The same result is also achieved by considering the super-Higgs mechanism where, in order to get rid of the Goldstone mode, one uses the super-gauge transformations (2) and (3), and defines the relevant field components of Z and i(Z − Z) as Examining the lowest components of the transformation (2), one can easily see that the real part of Z| cancels the Goldstone mode of (17). Similarly, when applying the derivatives Dα and D α to (3) and taking their lowest components (then DαD α V | = σ m αα B m ), one finds that the vector field "eats up" the Goldstone mode as The Minkowski vacuum after inflation can be easily lifted to a de Sitter vacuum (Dark Energy) in our model by the simple modification of the Polonyi sector and its parameters as [11] It leads to a small positive cosmological constant and the superpotential VEV where a ≡ ( is the SUSY breaking vacuum solution to the Polonyi parameters in the absence of a cosmological constant (see Ref. [39] also).

Gravitino and Polonyi production
In this Section, we consider the gravitino (ψ µ ) and Polonyi (A) particles production during inflation and after it. We assume that all other (heavy) SUSY particles (not present in our model) have masses larger than those of Polonyi and gravitino (High-scale SUSY), with gravitino as the LSP (the lightest superpartner of known particles) and as the cold Dark Matter (SHGDM).
There are several competitive sources of particle production in our model. First, gravitino and Polonyi particles can be produced via Schwinger's effect (out of vacuum) sourced by inflation. Since the mass of a Polonyi particle is higher than two gravitino masses (Sec. 2), the former is unstable, and decays into two gravitinos, A → ψ 3/2 ψ 3/2 . Second, both gravitino and Polonyi can be produced by inflaton decays, during oscillations of the inflaton field around its minimum after inflation. A competition between the gravitino/Polonyi creation during inflation and their production by inflaton decays is known to be very sensitive to the mass hierarchy. It is, therefore, very instructive to study them in our model (Sec. 2) that is minimalistic and highly constrained.
Since the exact equations of motion in our model (Sec. 2) are very complicated, in this Section we take them only in the leading order with respect to the inverse Planck mass. Let us begin with Polonyi and gravitino production during inflation by ignoring for a moment their couplings to inflaton. The effective action of the Polonyi field in the FLRW background reads where the non-minimal coupling constant (of Polonyi field to gravity) is ζ = 1 in our case (Sec. 2), A is the Polonyi field, a is the FLRW scale factor, M A is the Polonyi mass, and R is the Ricci scalar.
A mode expansion of the Polonyi field in terms of the conformal time coordinate η reads where b, b † are the (standard) creation/annihilation operators, and the coefficient functions h, h + are properly normalized as It follows from Eqs. (24) and (25) that the equations of motion of the modes are given by and we have defined h = d 2 h/dη 2 with respect to the conformal time η. For our purposes, it is convenient to rescale Eq. (27) by some reference constants a(η * ) ≡ a * and H(η * ) = H * to be specified later, and rewrite it as in terms of the rescaled quantities Similarly, the gravitino field is governed by the massive Rarita-Schwinger action in terms of the gravitino field strength and the supercovariant derivative by using the notation γ µ 1 ...µn = γ [µ 1 ....γ µn] with unit weight of the antisymmetrization. The supergravity torsion is of the second order with respect to the inverse Planck mass, so that it can be ignored. The Γ ρ µν can be represented by the symmetric Christoffel symbols, but they are cancelled from the action (29).
The gravitino equation of motion now reads In the flat FLRW background, Eq. (32) becomes where we have A solution to Eq. (33) for the helicity 3/2 modes reads We find that the equations of motion for the 3/2-helicity gravitino modes have the same form as that of Eq. (27), namely, It is customary to write them down (similarly to the well known relation between Dirac and Klein-Gordon equations) as where, in our case, we have Equation (36) can be rescaled in the same way as Eq. (28).
Interactions of gravitino with matter fields can be described in terms of the effective gravitino mass M 3/2 that is a function of matter fields in the Rarita-Schwinger equation, with m 3/2 = M 3/2 . In our model with the matter given by inflaton and Polonyi scalars (Sec. 2), we find where α ≡ A and A = α +Ã, in terms of inflaton field φ and Polonyi scalarÃ with Ã = 0, and we have restored the dependence upon Planck mass.
In order to obtain the number density of produced particles, we perform a Bogoliubov transformation, This transformation is supposed to be done from the vacuum solution with the boundary condition η = η in , corresponding to the initial time of inflation, to the final time η = η f when particles are no longer created from inflation. Since we have a /a 2 1 and ba/k 1, we can take the extremes as η in = −∞ and η f = +∞ in the semiclassical approximation. Given such boundary conditions, the energy density of Polonyi particles produced during inflation is given by where we have used the standard notation Similar equations are valid for gravitino, with the power spectrum (ii) Our proposal about Polonyi particles produced during inflation reads where M A is the Polonyi mass, Ω R h 2 4.31 × 10 −5 is the fraction of critical energy density that is in radiation today, Ω A h 2 is fraction of the critical energy density of produced Polonyi fields (and a similar estimate for gravitino). Though we do not have a rigorous proof of Eq. (46), it can be argued by starting from where ρ R is the energy density of radiation, ρ A is the Polonyi energy density, T reh and T 0 are the temperature of the Universe at reheating time t reh and today t 0 , respectively. Assuming that Polonyi particles are produced after the de Sitter phase t e , when the transition to the coherent oscillation phase begins, the inflaton and Polonyi energy densities will be redshifted with almost the same rate. This scaling holds until the reheating stage finishes and the radiation dominated epoch begins. Assuming that most of the energy density is converted into radiation contribution with we obtain where H(t e ) is the Hubble parameter at a fixed time t = t e . Then Eq. (46) follows from Eq. (49).
(iii) According to Eq. (43), relating Hubble scale, Polonyi mass and the desiderata Polonyi energy-density, there is about 8th-orders-of-magnitude suppression of the energy-density. According to (i), the normalized power spectrum P A cannot provide such suppression with our values for M A and H inf . However, it comes from the dilution factor (ã) −3 = (a f /a i ) −3 in Eq. (43).
Our semi-analytical estimations for Eq. (43) indicate that almost all Polonyi particles are produced in an excursion of the inflaton field around φ e ≡ φ(t e ) with ∆φ 0.2. The value of the dilution factor can be estimated from where we have defined ∆Φ After integrating over the effective particle production region ∆φ, we find ∆Φ = 18.2, i.e.
(iv) To get the masses M A and m 3/2 ≡ m ψ in a different way, we have to add a few more assumptions about details of reheating. Since our SHGDM scenario is based on Starobinsky inflation, all cosmological parameters can be fixed modulo the e-foldings number N e that is between 50 and 60 for compatibility with CMB observations. This also allows us to estimate the error margin for the masses in question at about 20%.
Let us take N e = 55 as the best fitting (reference) point [6,60]. This leads to the following set of the cosmological (inflation) parameters [61]: n s = 0.964, r = 0.004, m inf = 3.2 · 10 13 GeV, H inf = πM P P g /2 = 1.4 · 10 14 GeV , (53) where M P = 2.44 · 10 18 GeV is the reduced Planck mass, and P g stands for tensor perturbations. In our SHGDM scenario, well below the inflaton mass scale, the low-energy theory is given by the Standard Model (SM) that has the effective number of d.o.f. as g * = 106.75. Then, it is reasonable to assume that all the SM particles were generated by perturbative inflaton decay, whose reheating temperature is well known in Starobinsky inflation [62,63,6], This value is also consistent with the successful leptogenesis scenario of Ref. [64] based on the seesaw type-I mechanism, that requires the reheating temperature to be higher than 1.4 × 10 9 GeV.
These masses are compatible with the correct abundance of cold DM, according to our numerical estimates in Fig. 1, that lends further support towards our conjecture in Eq. (46).
(v) The decay rate of Polonyi particle into two gravitinos is given by [40,41] This channel is a direct consequence of the gravitino mass generation mechanism from the nonvanishing Polonyi vacuum expectation value. In addition, it implies that Polonyi particles rapidly decay into gravitinos. Moreover, since we have m 3/2 = 7.7 · 10 12 GeV, it implies Γ H inf . This means that the decay time scale is much larger then the production time during inflation, i.e. τ A→ψ 3/2 ψ 3/2 τ inf lation . As a consequence, the decays of Polonyi particles into gravitinos are subleading and negligible during inflation. Therefore, gravitino and Polonyi particles are independently created during inflation. After the reheating, Polonyi particles will completely decay into gravitinos (see below). In particular, the Polonyi number density n S gets transformed into a contribution to the gravitino number density ∆n Ψ = 2n A . Since the Polonyi mass is about two times of the gravitino mass, the Polonyi energy density before its decays is completely converted into gravitinos, i.e. (∆Ω ψ )h 2 = Ω S h 2 .
In Fig. 1 we show a numerical simulation of the produced gravitino mass density as a function of the Polonyi mass. 5 The spectrum is composed of two contributions: (a) the Polonyi energy-density spectrum produced during inflation, converted into gravitinos after reheating; and (b) the energy-density spectrum of gravitinos produced during inflation. The first contribution largely dominates over the second one. Intriguingly, the Polonyi and gravitino masses inferred from inflation, reheating and leptogenesis bounds are well compatible with the correct amount of CDM.
Next, let us consider the gravitino and Polonyi production from inflaton decays. As regards gravitino, its coupling to inflaton arises from the Weyl rescaling of vierbein in the gravitino action. The gravitino kinetic term does not contribute because of conformal flatness of the FLRW universe, so that the only source of gravitino production is given by the gravitino mass term 6 where the G tot in our model is given by and the gravitino mass is by m 3/2 = e Gtot/2 . The mass term in the form (58) also shows the Polonyi and inflaton couplings to gravitino. The perturbative decay rate of inflaton φ into a pair of gravitino is given by [43] In our case, the factor G φ vanishes at the minimum of the inflaton scalar potential (14), because the inflaton VEV also vanishes. So, the perturbative production of gravitino from inflaton decays is suppressed. One may also wonder about a non-perturbative gravitino production from inflaton decays, as was studied e.g., in Refs. [44,45]. However, unlike the usual Yukawa couplings, the coupling of gravitino to inflaton in our model is given by the exponential factor in (58) that never vanishes. Hence, the gravitino production due to inflaton decays can be ignored in our model. The situation is different with the Polonyi production due to inflaton decays. Inflaton is heavier than Polonyi particle by the factor of two approximately, according to Eqs. (53) and (56). The perturbative decay rate of inflaton into a pair of Polonyi scalars is (see e.g., Ref. [6] for a review) One may expect that the non-perturbative pre-heating in this case can be significantly more efficient due to a broad parametric resonance [46] that is a rather generic feature for coupled scalars. Indeed, in our model, inflaton is mixed with Polonyi field via the scalar potential (6). An expansion of the scalar potential (6) with the J-finction given by Eq. (11), with respect to both scalar fields φ and A, gives rise to the quartic interacion term whose dimensionless coupling constant λ does not vanish. It implies that the broad parametric resonance can happen in our model along the lines of Ref. [46], while it could lead to the enhancement of the perturbative production of Polonyi particles up to the factor of O(10 5 ) -see. e.g., Ref. [47] for the example of numerical calculations. However, our direct calculation yields so that the coupling constant λ is of the order O(10 −7 ) or less in our model, and this value fully compensates any possible enhancement of the Polonyi production by the broad parametric resonance. In short, the Polonyi production from inflaton decays is just perturbative in our case. The Polonyi particles produced in this channel quickly decay into gravitinos, with the perturbative decay rate (57) implying that the gravitino production from inflaton decays is sub-leading with respect to Schwinger's effect shown in Fig. 1. The former channel would be kinematically suppressed if the Polonyi mass were higher than half of the inflaton mass, because then the inflaton two-particle decay into two Polonyi particles would be forbidden.
(vi) To the end of this Section, we note that the pre-heating stage with the dust-like equation of state p = 0, started at the end of inflation at t i ∼ 1/H and finished in the period of reheating at t f ∼ M P l /T 2 reh , is sufficiently long to provide the growth of density fluctuations and the formation of nonlinear structures of gravitationally bounded systems. In particular, the Primordial Black Holes (PBHs) may be formed at this stage. Later on, the PBHs can be evaporating and converting about 1/N (1g/M) of their masses to gravitinos, where N is the number of evaporated species with account of their statistical weight, and M is the PBH mass. Here we have taken into account that at M ≤ 1g the temperature of Hawking evaporation is T ev ≥ 10 13 GeV, so that the fraction of evaporated gravitino is determined by the ratio of their statistic weight to the statistic weight of all evaporated species. During the pre-heating stage at t i ≤ t ≤ t f , the gravitationally bounded systems are formed in the mass range where is the mass within the cosmological horizon at the beginning of pre-heating, and is the mass of the gravitationally bounded systems formed at the end of pre-heating. Here δ is the amplitude of density fluctuations. According to Ref. [23], the minimal estimation of the probability of the PBH formation at this stage is determined by consideration of their direct collapse in the black holes with the special symmetric and homogeneous gravitationally bound configurations [48], and is given by ∼ δ 13/2 [49]. For δ ≈ 10 −5 , the estimated fraction of the total density at the end of pre-heating stage β, corresponding to the PBH in the range (64) is of the order of Hence, the contribution of gravitino produced in the evaporation of these black holes is negligible at the beginning of the matter dominated stage with T md = 1 eV as where we have also used N ∼ 10 3 . A formation of black holes in the course of evolution of the gravitationally bounded systems formed during the pre-heating stage can significantly increase the value of β, though addressing the problem of evolution of the gravitationally bounded systems of scalar fields deserves a separate study. The situation drastically changes when inflation ends by a first order phase transition when bubble collisions lead to copious production of black holes at β ∼ 0.1 with the mass (65) [26,27]. Evaporation of these PBHs leads to a fraction of the total density by the end of pre-heating of the order 10 −4 in the form of gravitino, and it results in the gravitino dominated stage at T ∼ 10 −4 T reh ∼ 10 5 GeV. This huge gravitino overproduction excludes a possibility of the first order phase transition by the end of inflation. This is impossible in the single-field inflation, also by considering an extra axion-like field and extra moduli decoupled during the slow-roll epoch. The inflaton in our model has the characteristic Starobinsky potential that, as is well known, does not lead to any violent phase transition after inflation. The Polonyi field does not alter this situation. However, in a more general case, in which other scalar and pseudo-scalar fields enter the slow-roll dynamics, these extra fields can have scalar potentials ending in false minima during the reheating. In such case, the tunnelling from the false minima to the true minima can induce bubbles in the early Universe, catalising an efficient formation of PBHs. These issues deserve a more detailed investigation beyond the scope of our paper.

Conclusion and outlook
In this paper we studied the gravitino production in the context of Starobinsky-Polonyi N = 1 supergravity with the inflaton field belonging to a massive vector multiplet, and the mass hierarchy m inf > 2M A > 4m 3/2 close to the bounds (by the order of magnitude). On the one hand, we found the regions in the parameter space where the gravitino and Polonyi problems are avoided, and super heavy gravitinos can account for the correct amount of Cold Dark Matter (CDM). The dominating channel of gravitino production is due to decays of Polonyi particles, in turn, produced during inflation. On the other hand, we found that direct production of gravitinos during inflation is a subleading process that does not significantly change our estimates.
Intriguingly, our results imply that the parameter spaces of cold DM and inflation can be linked to each other, into a natural unifying picture. This emerging DM picture suggests a phenomenology in ultra high energy cosmic rays. For example, super heavy Polonyi particles can decay into the SM particles as secondaries in top-bottom decay processes. Cosmological high energy neutrinos from primary and secondary channels may be detected in IceCube and ANTARES experiments. A numerical investigation of these channels deserves further investigations, beyond the purposes of this paper.
Our scenario offers a link to the realistic Supersymmetric Grand Unified Theories (SGUTs) coupled to supergravity. The super-Higgs effect considered in Sec. 2 is associated with the U (1) gauge-invariant supersymmetric field theory. This U (1) can be naturally embedded into a SGUT with a non-simple gauge group. The Starobinsky inflationary scale, defined by either m inf or H inf is by three or two orders of magnitude lower, respectively, than the SGUT scale of 10 16 GeV. The SGUTs with the simple gauge group SU (5), SO(10) or E 6 are well motivated beyond the Standard Model. However, the SGUTs originating from the heterotic string compactifications on Calabi-Yau spaces usually come with one or more extra U (1) gauge factors as e.g., in the following gauge symmetry breaking patterns: E 6 → SO(10) × U (1), SO(10) → SU (5) × U (1), the "flipped" SU (5) × U X (1) and so on (see e.g., Ref. [50] for more). Alternatively, SGUTs can be obtained in the low energy limit of intersecting D-branes in type IIA or IIB closed string theories. Also in this context, extra U (1) factors in the gauge group are unavoidable. For instance, the "flipped" SU (5) × U X (1) from the intersecting D-brane models was studied in Refs. [51,52,53,54].
We propose to identify one of the extra U (1) gauge (vector) multiplets with the inflaton vector multiplet considered here. This picture would be very appealing because it unifies SGUT, inflation and DM. Moreover, the extra U (1) gauge factor in the SGUT gauge group may also stabilize proton and get rid of monopoles, domain walls and other topological defects [55].
Our scenario also allows us to accommodate a positive cosmological constant, i.e. to include dark energy (see the end of Sec. 2) too. Further physical applications of our supergravity model for inflation and DM to SGUTs and reheating are very sensitive to interactions between the supergravity sector and the SGUT fields. Demanding consistency of the full picture including SGUT, DM and inflation may lead to further constraints. For instance, a matter field must be weakly coupled to inflaton -less then 10 −3 -in order to preserve the almost flat plateau of the inflaton scalar potential. Among the other relevant issues, the Yukawa coupling of inflaton to a Right-Handed (RH) neutrino is very much connected to the leptogenesis issue. Inflaton can also decay into RH-neutrinos, in turn, decaying into SM visible particles. Of course, these remarks are very generic and have low predictive power, being highly model-dependent. But they motivate us for a possible derivation of our supergravity model from superstrings -see e.g., Ref. [50] for the previous attempts along these lines.
Another opportunity can be based on Refs. [56,57], by adapting the Pati-Salam model to become predictive in the neutrino mass sector and be accountable for leptogenesis in supergravity, as was suggested in Ref. [30]. Then one can generate a highly degenerate mass spectrum of RH neutrinos, close to 10 9 GeV, i.e. four orders smaller than the inflaton mass. In this approach an extra U (1) is necessary for consistency, while it can be related to the Higgs sector in supergravity. It is worth mentioning that our hidden sector includes only inflaton and Polonyi, and it may have to be extended. The scale invariance of the (single-field) Starobinsky inflation is already broken by the mixing of inflaton with Polonyi scalar, while its breaking is necessary for a formation of mini-PBHs during inflation. The physical consequences of the inflaton-Polonyi mixing demand a more detailed investigation, beyond the scope of this paper.
Our scenario can be reconsidered in the general framework of Split-SUSY and High-scale SUSY, by questioning its compatibility with the SM and the known Higgs mass value of 125.5 GeV in particular, e.g., along the lines of the comprehensive study in Ref. [58]. 7 Then the "unification help" from SUSY to GUT scenarios could be implemented in our model. In this case, several new decay channels are opened and new parameters enter. In particular, there is a scenario in which the Higgsino at 1 − 100 TeV scale is envisaged, with intriguing implications for future colliders. In such case, produced gravitinos can decay into Higgsinos that (in the form of neutralinos) could provide another candidate for Dark Matter. However, there also exist contributions from thermal production that may affect the Higgsino production. Actually, the upper bound on the scale of Split-SUSY, according to Fig. 3 of Ref. [58], is given by 10 8 GeV that already excludes compatibility of Split-SUSY with our scenario that requires a higher SUSY. In the case of High-scale SUSY, the upper bound in Fig. 3 of Ref. [58] is given by 10 12 GeV for a considerable part of the parameter space, so that this bound is again too low for our model. However, it is still possible to go beyond that bound in the case of High-scale SUSY, as is shown in Fig. 5 of Ref. [58], so that our SHGDM scenario is still allowed. A more detailed study of the compatibility of our scenario with the SM deserves further investigation, beyond the scope of this paper.

Appendix: the power spectrum of Polonyi and gravitino emissions
Below we provide some technical details about our calculation of the power spectrum of Polonyi and gravitino emissions, based on finding a numerical solution to Eq. (27) in the framework of adiabatic theory [67,68,69].
A change of variables from h k in Eq. (25) to W k as and plugging this into Eq. (25) yield where = ∂/∂η denotes the derivative with respect to the conformal time.
Applying the Bogoliubov transformation to h k leads to the following relation among W k and β: where ζ η k = W k η /2W η k .
The adiabatic approximation consists of considering the background metric to be slowly changing in time, so that the time variation can be treated by introducing a small parameter via the substitution ∂/∂η → ∂/∂η (this approximation can be reasonably applied during inflation), where the label (n) = (0), (2), (4), ... denotes the adiabatic order expansion. The numerical problem can be solved by defining the iterative map as follows: This map raises the adiabatic order as A few leading terms of the adiabatic order expansion can be straightforwardly calculated, with the following results: W (0) = w, W (2) = 3w 2 /8w 3 − w /4w 2 , W (4) = −k 1 (w ) 4 /w 7 + k 2 (w ) 2 w /w 6 − k 3 (w ) 2 /w 5 − k 4 w w /w 5 + k 5 w /w 4 , where k 1 = 297/128, k 2 = 99/32, k 3 = 13/32, k 4 = 5/8 and k 5 = 1/16. The j-th adiabatic order reads h (j) where the adiabatic vacuum state of the j-th order is defined by specifying the boundary conditions at a fixed value η * of η as A similar iterative procedure can be applied to gravitino also. The b µ modes can be decomposed into two vector-spinor fields, where u is a vector-spinor of spin 3/2, while h I,II has the structure similar to that of Eq. (69) in terms of other functions W ± . The equations of motion for W ± are very complicated, where we have used the notationσ = (I, −σ i ) and σ = (I, σ i ). One arrives at the following gravitino spectrum: with a similar equation for Polonyi fields. Plugging in the adiabatic solutions for β's and b's into the above equations leads to the final power spectra, though only numerically. Our numerical results can be compared with the semi-analytical results of Ref. [70], where the steepest descent method of integration was used, leading to the crude estimate where H i and R i are the Hubble parameter and the Ricci scalar during inflation. Inserting R i = 6(ä/a + H 2 i ) H 2 i , H i = 1.4 × 10 14 GeV and M A = 1.5 × 10 13 GeV into (82) yields |β k | 2 0.67, in basic agreement with our numerical results. It is worth noticing that, according to Ref. [70], Eq. (82) is subject to O(k/M A ) corrections.