Reheating in runaway inflation models via the evaporation of mini primordial black holes

We investigate the cosmology of mini Primordial Black Holes (PBHs) produced by large density perturbations that collapse during a stiff fluid domination phase. Such a phase can be realized by a runaway-inflaton model that crosses an inflection point or a sharp feature at the last stage of inflation. Mini PBHs evaporate promptly and reheat the early universe. In addition, we examine two notable implications of this scenario: the possible presence of PBH evaporation remnants in galaxies and a non-zero residual potential energy density for the runaway inflaton that might play the role of the dark energy. We specify the parameter space that this scenario can be realized and we find that a transit PBH domination phase is necessary due to gravitational wave (GW) constraints. A distinct prediction of the scenario is a compound GW signal that might be probed by current and future experiments. We also demonstrate our results employing an explicit inflation model.


Introduction
The detection of gravitational waves (GWs) from black hole mergers has confirmed that black holes are undoubtedly present in our universe with a notable abundance and a quite extended mass range. Black holes with masses from about 3 solar masses (M ) up to 160M have been reported [1]. The black hole catalog is extended by the electromagnetic observations of supermassive black holes in the center of galaxies with masses that reach tens of billions of solar masses. Apparently massive black holes are easier to detect. Nevertheless, it is natural to ask whether black holes spread across the mass spectrum towards smaller values as well.
From simple theoretical considerations, the minimum black hole mass is of the order of the Planck mass m Pl = G −1/2 ∼ 10 −5 g. Upper bound does not exist; the entire universe energy density could in principle collapse into a black hole. What is of special interest in the small mass limit, in sharp contrast to the large mass limit, is that the predicted Hawking radiation is manifest [2,3]. The emitted thermal radiation has a temperature T BH ∝ 1/M BH and black holes with mass M BH less than 10 15 grams have a lifetime less than the age of the universe and are not expected to wander around in space. Moreover, black holes with such a small mass, far below the Chandrasekhar limit [4] of 1.4M , cannot be produced by stellar collapse processes. However, small mass black holes can be produced from the collapse of overdense regions in the highly dense early universe [5,6].
Black holes of small mass are therefore relevant for the early universe cosmology. They are produced in fractions of the first second of the cosmic evolution and evaporate promptly. Indeed, black holes with mass less than 10 9 grams evaporate before big bang nucleosynthesis (hereafter referred as BBN) and leave no apparent observational signature behind [7,8]. The fact that mini black holes transform nearly their entire mass into thermal radiation is a mechanism that contributes to the entropy observed in our universe. Actually, mini black holes can reheat the universe if no other radiation sources are efficient enough. This is the central topic of this paper.
In the conventional inflationary scenario, reheating is understood as the decay of the inflaton field itself into entropic degrees of freedom. This is achieved either via parametric resonance or by single particle decays [9,10]. However, a stage of an oscillatory phase for the inflaton field is necessary and this is not always the case. There are models where the inflaton potential V (φ) gradually decreases after the end of the inflationary phase or even does not have a minimum at all. These are the so-called runaway models. This class of post-inflationary models has also been dubbed "NO-models" due to their non-oscillatory behavior [11]. Reheating in these models seems problematic because it is generally too feeble. The main source of radiation is the gravitational particle production, [12] which, in its original version, fails to reheat the universe successfully. Reheating can become efficient if the inflaton sector is augmented with a tailor made coupling with an extra scalar field that depletes the inflaton energy density [11,13] via instant preheating. A third post-inflationary reheating mechanism is via the primordial black hole production and evaporation, which has been first discussed in Ref. [14,15] Inflationary models without a minimum are predicted in many theories beyond the standard model of particle physics. Some notable early examples can be found in stringy set ups, supergravity, braneworlds as well as in many phenomenological models, see [16][17][18][19][20][21][22][23], just to mention a few. An additional phenomenological motivation to study these models is the unification they offer in describing the two accelerating stages of our universe: inflation and dark energy. In fact, runaway models have been widely introduced to describe quintessence fields. This unified framework has been called quintessential inflation [24,25], and recently notable progress on this direction has taken place [26][27][28][29][30].
In this work, we consider runaway inflationary models and implement the reheating of the universe via the Hawking evaporation of primordial black holes produced by the collapse of overdense regions of the inflaton field itself. Actually, a kination phase due to a stiff energy component, w > 1/3, can be realized only at sufficiently early times of the cosmic evolution, where w is the equation of state. A simplistic comparison of the different red-shifts a −3(1+w) of the various energy components shows that the scenario of a stiff fluid domination at the early universe and a cosmological constant domination at the late universe is a plausible guess. Furthermore, since the inflaton potential energy density does not vanish, we contemplate the scenario that the dark energy density of the universe is the present-day energy density of the inflaton. We employ a particular inflation model, built in the framework of α-attractors [31,32]. Notably, the fact that the generated PBHs are ultra light implies that the (n s , r) inflationary predictions lay in the sweet spot region of Planck and BICEP [33,34]. This scenario, which is remarkably economic in terms of ingredients, has been introduced and analyzed in Ref. [35].
We specify the conditions under which this scenario can be realized extending the work of Ref. [35]. We take into consideration the impact of primary and secondary GWs on the BBN and CMB observables and we constrain the (M PBH , β) parameter space. We find out that a transit PBH domination phase, taking place between kination and radiation domination, is necessary. Interestingly enough, the allowed parameter space can become even more narrow by near future GW probes. Last, but not least, an additional striking implication of the very same scenario is that the evaporation of the mini black holes might leave a stable remnant behind [15,[36][37][38][39][40][41][42][43][44][45][46][47][48][49][50]. We notice that remnants must have a mass in a particular range and in correlation with the PBH mass in order to have a cosmologically significant abundance.
The paper is organized as follows. In the next section we discuss the process of PBH formation during a kination era. In Section 3 we motivate the formation of ultra-light PBHs and we review the basic elements of PBH evaporation and PBH remnants. In Section 4, we study the evolution of the energy densities of the inflaton, PBHs and remnants and we specify the relevant for our discussion cosmological quantities. In Section 5, we examine the primary and secondary GWs and their impact on BBN and CMB. That section contains the main result of this paper, since we derive the main constraints on the PBH evaporation reheating of the runaway inflation models. In Section 6, we introduce an explicit runaway inflation model that features an inflection point, and in Section 7 we discuss the quintessential inflation scenario and study the post-inflationary evolution of the inflaton field. Finally, in Section 8 we conclude. The alternative to this scenario mechanisms, gravitational reheating and instant (p)reheating, are reviewed in Appendix A.

PBH Formation during a Kination Era
Let us assume that at the early moment t form a fraction β of the energy density of the universe collapses and forms PBHs. The mass density of the PBHs is at the moment of formation, where ρ tot = 3H 2 M 2 Pl with M Pl = m Pl / √ 8π, the reduced Planck mass. The γ parameter is the fraction of the collapsing mass that finds itself inside the black hole. If the collapsing region is the Hubble horizon, the mass of the black hole is Pl t, w the equation of state of the background fluid and t the cosmic time. The formation probability is usually rather small, β 1, and the background energy density ρ bck = (1 − γβ)ρ tot is approximately equal to ρ tot . The PBHs are non-relativistic matter and their number density, scales like a −3 , while the background energy density scales as ρ ∝ a −3(1+w) . If PBHs have an extended mass spectrum one writes n PBH (M ) = dN PBH (M )/d ln M , where N PBH is the integrated PBH number density up to mass M . Let us assume that the bulk energy density is in the form of stiff fluid, realized by the domination of the kinetic energy of a scalar field. Such a fluid has a barotropic parameter w 1 and the phase is called kination (KIN). A non-oscillatory inflaton field can give rise to a kination phase. During kination the scale factor goes like a(t) ∝ t 1/3 , the Hubble horizon mass is M hor = (3/2)m 2 Pl t, and the scaling of the background energy density redshifts like ρ KIN ∝ a −6 .
Let us now turn to the parameter β, which is of central importance for the PBH cosmology. It is the fraction of ρ tot that collapses into black holes and it can also be interpreted as the probability of such an event to happen. Assuming Gaussian statistics, the black hole formation probability for a spherically symmetric region is [51] β(M ) = δc dδ 1 that is approximately equal to β √ 2π σ/δ c e − δ 2 c 2σ 2 for σ δ max , where δ max the upper value of the integration interval. During kination pressure is maximal and we expect the overdense regions to be spherically symmetric collapsing nearly immediately after horizon entry. The PBH abundance has an exponential sensitivity to the variance of the perturbations σ(M ) and to the threshold value δ c . In the comoving gauge Ref. [52] finds that δ c has the following dependence on w, see Figure 1. For the conventional radiation cosmology, w = 1/3, it is δ c = 0.414. In our case of kination dominated early universe the PBH formation occurs when the density perturbation exceeds the threshold, Notice that the expression (4) is different from the result of Carr [6]. In Carr's analysis an overdensity collapses if its size at the maximum expansion is smaller than the particle horizon R hor and larger than Jeans length R J ∼ √ wR hor . The former is comparable with the curvature scale and for a pure kination phase, where the sound speed is equal to one, Jeans length and curvature scale are equal. According to Carr's condition, it is δ c w and the pressure gradient force seems to rule out any collapse of a density perturbation during a pure kination phase.
However, Carr's result has been refined in Ref. [52] where it was pointed out that collapse is realized if the sound crossing time over the radius is longer than the free fall time in the interval from the maximum expansion to complete collapse. In their analysis the Jeans length has been identified with [52] where a max coincides with the Hubble radius of the background Friedmann universe in the uniform density slice. This expression for the Jeans length implies that the threshold value of primordial black hole formation is given by Equation (4). It is interesting that the analytically found δ c value reaches its maximum for w ∼ 0.4 and decreases as w increases. The analytic δ c decrease for large w values is caused by the shortening of the dynamical time of the collapse due to the contribution of the pressure to the source of gravity. We note that the decrease of the threshold value for large w values (w 0.4) has not been observed in numerical simulations completed for 0.01 ≤ w ≤ 0.6 in Ref. [53], and the analytic expression (4) might underestimate the threshold value. Nonetheless, what is important for us is that the threshold value for large w appears to be less than the maximum allowed, (3 + 3w)/(5 + 3w), and a collapse during kination domination can be in principle realized. In the right panel the increase of β, and consequently the PBH number density, with the hardening of the EOS is shown. The β(w) value has been normalized by the β(w = 1/3) ≡ β rad value. The effect becomes stronger for smaller σ values.
Taking the analytic result (4) at face value, we see that for w 0.4 we have a quite similar behavior to the case where the equation of state softens by taking values w < 1/3 as, e.g., during the QCD phase transition where an enhancement of the PBH formation probability is expected [54]. In our case the EOS becomes harder due to the stiff fluid domination and there is a δ c (1/3) − δ c (1) = 3.8% decrease in the analytic value for the threshold. Smaller values of w but still close to one also give an enhanced PBH formation probability, δ c (0.9) = 0.385, δ c (0.8) 0.395. Thus, a smaller amplitude for the density perturbations is required during kination regime, see Figure 1. We recall that the power spectrum of the comoving curvature perturbation is where σ 2 ∼ P δ is the variance of the density perturbations in a window of k. For kination domination it is σ 2 ∼ (1/4)P R and approximately we obtain P R ∼ 4δ 2 c /W 0 [(2πβ 2 ) −1 ], where W 0 is the principal branch of Lambert W function.

PBH Evaporation and PBH Remnants
The generation of PBHs via the collapse of large inhomogeneities at the time of reentry require a particular shape for the spectrum of the primordial curvature perturbations, P R (k). The spectrum has to be nearly flat at scales observed at the CMB sky, k ∼ 0.05 Mpc −1 , with value P R ∼ 2 × 10 −9 . On the other hand, PBHs form only if the power spectrum is amplified roughly by seven orders of magnitude at the wavenumber k peak and this amplification must not, by any means, spoil the CMB predictions. For the class of inflationary models that predict a spectral index value n s ∼ 1 − 2/N * and tensor-to-scalar ratio r ∼ 12α/N 2 * , such as α-attractors, the Planck 2018 measured value (in the 68% confidence region), n s = 0.964 ± 0.004 (8) constrains significantly the position of the P R peak on the k-space. In Ref. [35] it was pointed out that this class of inflationary models generate PBHs and predict a spectral index value inside the 68% CL region of the Planck 2018, n s 0.96 (9) only if PBHs are ultra-light with mass M PBH 10 5 grams [35]. This rough bound comes from (9) and the relation ∆N ∼ 2/(1 − n s ), where ∆N are the efolds of inflation that take place in between the moments of horizon exit of the scales k = 0.05 Mpc −1 and k peak . M PBH admits larger values for n s > 0.955 which is the lower bound of the 95% CL region. This issue has been also discussed recently in [55], where they found agreement with the CMB data for M PBH 10 8 g in the framework of α-attractors.
The fact that ultra-light PBHs seem to be favored by the CMB measurements of the spectral index n s and the running of the spectral index α s prompt us to investigate inflationary models that generate mini PBHs. The striking characteristic of the mini PBH scenario is first, the Hawking radiation emitted by the evaporating PBHs and second, the remnants that might be left behind at the end of the evaporation process.

Evaporation
Hawking has shown that black holes radiate thermally with a temperature [2,3] assuming the Schwarzschild solution, that is a black hole without charge or angular momentum. Black holes are expected to evaporate on a time scale t evap ∼ G 2 M 3 BH /( c 4 ). In the second line of Equation (10) we took k B = c = = 1 and we will use this convention in the following.
The mass-loss rate of an evaporating black hole is where g H (T BH ) is a spin-weighted number of degrees of freedom of emitted particle species and it takes the value 108 for the Standard Model particles for hot enough black holes, T BH > 100 GeV or M BH < 10 11 g; in the limit of cool black holes it is g H (T BH ) ∼ 7. Integrating the mass-loss rate over time, the time dependent mass of the BH, for t f t evap , and the BH lifetime, are found [56]. For t evap equal to the present age of the universe, 13.8 billion years, one finds a critical black hole mass M BH, cr 10 15 (g H (T BH )/108) 1/3 ∼ 5 × 10 14 g. Apparently, light BH with mass M BH M BH, cr are irrelevant for late time cosmology but might play a crucial role in the early cosmic evolution, see Figure 2.
PBHs with mass M PBH ∼ M BH, cr are a source of gamma rays today and can constitute only a small fraction of the dark matter in the galaxies. Lighter PBHs generate an entropy at the evaporation that might alter big bang observables such as BBN. In particular, for lifetimes τ = 10 2 -10 7 s, which correspond to M PBH = 10 10 -10 12 g, hadrodissociation processes become important and the debris deuterons and nonthermally produced 6 Li constrain β(M PBH ); for lifetimes τ = 10 7 -10 12 s, which correspond to M PBH = 10 12 -10 13 g, photodissociation processes overproduce 3 He and D and put strong constraints on β(M PBH ) [57][58][59][60]. In addition, the heat produced by PBHs evaporation after the era of recombination damps small-scale CMB anisotropies and in the mass range 2.5 × 10 13 g M PBH 2.4 × 10 14 g only a very suppressed PBH abundance is allowed [56,61,62]. The conclusion is that PBHs with mass in the range 10 9 -10 17 g are significantly constrained [7]. In the smaller mass limit the tight constraints on PBHs abundance are raised.
Finally, let us add some interesting remarks. First, if the temperature of a PBH is initially smaller than the background cosmic temperature, accretion effects should be taken into account, but this practically should not modify PBH lifetime. Accretion decreases the temperature of a PBH but this decrease is negligible for small enough PBH mass, whereas the cosmic temperature falls fast due to the expansion. Second, the Hawking temperature T BH is expected to reach a maximal value at some point and afterwards decrease. In this last and most uncertain stage of the PBH evaporation the rate dM/dT BH should turn into positive.

Remnants from the Evaporation of PBHs
Apparently the evaporation process of BHs is uncertain at masses of the order of the Planck scale. The semiclassically derived Hawking temperature T BH diverges for vanishing M BH and there is a logical guess that Hawking radiation halts somewhere near the Planck scale, leaving behind a stable black hole remnant 1 with mass M rem . The PBH remnant mass can be written in terms of the Planck mass where K is a factor that parameterizes our ignorance of the physics that operates at the relevant energy scales. Energy conservation [36], extra spatial dimensions [37,63], higher order corrections to the action of general relativity [38], the information loss paradox [39] might prevent complete evaporation. Different theories predict stable black hole relics of different mass. K may be of order one, with remnant masses characterized by the fundamental scale of gravity, m Pl = G −2 , but other values for K are also admitted. It can be K > 1 if BH have quantum hair [41], or K < 1 in treatments where a generalized uncertainty principle is applied [40]. In our analysis we let K be a free parameter and we remain agnostic about the fundamental physics that might prevent black holes from complete evaporation.
Equation (11) written in the compact formṀ BH = −A M −2 BH yields the time dependence of the black hole mass for t ≤ t evap and A a dimensionful parameter dependent on T BH , A = πGg H (T BH )M 4 Pl /480, where G is the gray-body factor. The decay of the BH mass implies that we can define a decay rate, Γ BH , according to the formulaρ ≡ −Γ BH ρ. We attain Γ BH = −Ṁ BH /M BH and, contrary to the conventional perturbative particle decays, here it is the black hole mass that decays not the number density. The time dependent decay rate is where c K = K 3 480 √ 8π/(πGg H (T BH )). As t → t evap the final explosive phase of black hole evaporation takes place.
At the moment right after the PBH evaporation the energy of a PBH remnant is E rem (t + evap ) = p 2 rem + M 2 rem where p rem = M rem v rem / 1 − v 2 rem is the momentum obtained by the remnant after the completion of the evaporation, with v rem the recoil velocity of the remnant. The order of magnitude of the recoil velocity can be determined only after particular assumptions about the final stage of BH evaporation are made.
As discussed above, it is expected -but not firmly confirmed-that the result for the Hawking temperature, Equation (10), stops being valid in the vicinity of Planck densities. Let us assume that there is a maximal temperature T max for the black hole during the course of evaporation and afterwards a cooling down occurs during the emission process. A maximum temperature implies the presence of a black hole remnant of mass M rem , that might or might not be a horizonless object. Accordingly, the surface temperature of the remnant, T rem , might be zero or in equilibrium with the background cosmic temperature. In both cases, the temperature of the final state is well below the Planck energy scale, whereas it is plausible to assume that T max is not far from the Planck scale. Hence we take . From the moment where T BH = T max until the moment of the remnant formation a number of N f quanta has been emitted. The typical energy of these quanta is E ∼ T max ∼ ∆M/N f where ∆M = M BH (T max ) − M rem is the mass radiated away in the final stage of the evaporation. The number N f of quanta depends on the model assumed for the the description of the final stage. Ref. [50] discusses these issues and gives a number N f ≤ O(10 2 ). The momentum of each quantum emitted is about ∆M/N f and the evaporating black hole performs a random walk in the momentum space leaving behind a remnant with average momentum [50]

Remnants CDM Scenario
The remnants, though initially close to relativistic, are cold dark matter today. The momentum of remnants redshifts as a −1 , and being kinetically decoupled, they free-stream a distance where subscript "0" denotes present values and "eq" equality era. It has to be k −1 fs < 0.1 Mpc in order that the smallest structures observed in the matter power spectrum are not washed out by the free stream of remnants. This constraint translates into an upper bound on t evap , which, in any case, precedes the BBN epoch. As noted in [64], the remnant dark matter scenario is a CDM scenario in agreement with the cosmological observations.
Next we overview the cosmological evolution of mini PBHs and their remnants.

Evolution of Mini PBHs and Their Remnants during a Kination Era
In this section, we will discuss elements of the cosmological evolution of mini PBHs and their remnants with the background energy density dominated by a stiff fluid. Initially, we will examine the exact system of equations that describes the evolution of the different cosmic fluids. Afterwards, we will pursue an approximate analytic description that will show the characteristic mass, time and temperature scales of our scenario and the parameter space that a viable cosmology is realized.

Energy Densities
We assume a universe initially dominated by the kinetic energy of the inflaton scalar field ϕ.
Its energy density has a maximal redshift a −6 . During the kination phase the production of PBHs takes place with β 1, that is the initial amount of PBHs is tiny compared to the bulk energy density. The PBHs evolve as pressureless matter a −3 until the cosmic time becomes comparable to the evaporation time, t evap . Assuming that the evaporation of PBHs results in relativistic particles mostly, leaving also behind a remnant mass, we have a universe with energy density partitioned among three fluids: scalar field, evaporating PBHs and radiation in the first phase, and scalar field, radiation and PBH remnants in the second phase. This dynamical system is described by the Friedmann together with the continuity equations, where dN = d(ln a) = Hdt is the differential of the e-fold number. Figure 3 depicts the evolution of this system of fluids with a transit PBH domination phase and for arbitrary initial conditions. For t t evap the PBH energy density scales such as e −3N and that of radiation increases as N . The scalar field rolls down its steep potential without decaying, but redshifting as e −6N . The produced radiation will dominate the energy density reheating finally the universe. In the late universe the dark energy of the runaway scalar field ϕ might be the cause of the accelerated expansion observed today. PBHs and radiation, ρ ϕ , ρ PBH and ρ rad , respectively, normalized by the initial critical energy density ρ crit, 0 , are shown for a transit PBH domination phase. The runaway scalar field is modeled as a stiff fluid, the PBHs initially evolve as a pressureless fluid but they promptly evaporate producing radiation. The truncation of the PBH evolution is dictated by Equation (16); for comparison the effect of a constant in time decay rate is also depicted. The radiation initially evolves as N and after the completion of the evaporation as e −4N . In the right panel we depict the β values that a PBH domination phase is realized/avoided for a background dominated by stiff fluid. We also contour plot five values of the reheating temperature due to PBH evaporation.
After the PBH evaporation, the remnant of each PBH is left behind. Their energy density evolves as It is understood that in Equation (22), ρ PBH is replaced by the remnant energy density. The PBH remnants will be dark matter constituents of the galaxies today.
The exact evolution can be traced by solving the above system of equations. In the following, for clarity and simplicity, we will make the approximation of instantaneous evaporation for PBHs. This way, analytic expressions for the energy densities and the constraints, as well as threshold conditions, can be obtained.

Reheating the Universe via PBH Evaporation
Runaway inflation scenarios in their simplest version cannot reheat the universe. A modelindependent source of radiation comes from the de Sitter horizon that is characterized by a Hawking temperature. If radiation domination is realized by this process, gravitational reheating is achieved [12,65], but it is rather inefficient. A completion with specially designed interactions is necessary in order for the thermalized early universe to be realized, see Appendix A. In our scenario, the Hawking radiation from mini PBHs automatically reheats the universe. The special characteristic is the need of a "feature", such as an inflection point or a step, which enhances the curvature power spectrum.
Let us start initially assuming a general EOS w for the background. In the approximation of instantaneous evaporation, the moment right before evaporation, t − evap , the energy density of the PBHs over that of the scalar field is where we considered a(t) ∝ t 2 3(1+w) and we plugged in the horizon mass. The coefficient g =g(g * , t − evap ) is equal to one unless the universe is radiation dominated; in the latter case it is ρ rad ∝ g * T 4 and T ∝ g −1/3 s a −1 where g s are the entropic degrees of freedom and g * the thermal ones. Substituting the mass inside the curvature radius at the evaporation moment of the PBH M hor (t evap ) = 3M 3 PBH (1 + w)/4m 2 Pl , a threshold β(M PBH ) value is found that specifies whether or not the universe has become PBH dominated.
Let us assume that the formation of black holes and the subsequent evaporation takes place during an early cosmic era of stiff fluid domination, called kination dominated phase (KIN). In our context it is the inflaton field ϕ that rolls down a runaway potential that gives rise to kination phase. According to Equation (26) and for w = 1 the energy density of the PBHs at the moment right before evaporation is Neglecting the PBH mass loss due to the evaporation process, the assumption of a kination phase is valid roughly for γ 2 β < m 2 Pl /M 2 PBH , otherwise a PBH dominated universe is realized. Plugging in benchmark values, the PBH domination is avoided for β < β thresh where, In Figure 3, we depict the bound above and highlight the parameter space that a kination or a PBH domination phase is realized for PBHs with mass M PBH < 10 9 g.

Reheating without PBH Domination
After PBH evaporation the background energy density is partitioned between the runaway inflaton, ρ ϕ , the entropy produced by the PBH evaporation, ρ rad and remnants ρ rem . Radiation is about M PBH /M rem times larger than ρ rem (t evap ). Assuming fast thermalization of the evaporation products, the radiation redshifts like ρ rad ∝ g * g −4/3 s a −4 and sooner or later it will dominate the runaway scalar field that redshifts as ρ ϕ ∝ a −6 . At some moment the radiation dominates the energy budget of the universe and when ρ ϕ (t) = ρ rad (t) we define the reheating moment t = t rh . Equivalently, a reheating temperature of the universe is defined that reads [35] T rh ≈ 10 −2 GeV β 10 −28 It is notable that the reheating temperature depends only on the β value.

Reheating after PBH Domination
Let us now consider the case β > β thresh where PBHs dominate the energy density of the universe before evaporation. At the time t − evap we can approximate ρ PBH ≈ 3H 2 M 2 Pl and at the moment of evaporation nearly the entire ρ PBH transforms into radiation with density π 2 g * T 4 rh /30 and H = 2/(3t). The reheating temperature of a PBH dominated universe is Comparing with Equation (27), we see that T rh ∝ M

PBH Remnants Abundance
The moment right after the evaporation, which we label t + evap , the energy density of the PBHs has either vanished or a remnant is left behind. In the latter case the evaporation halts somewhere close (or not) to the Planck scale and the energy density contained in the form of remnants is (M rem /M PBH )ρ PBH (t − evap ), plus their kinetic energy. The factor (M rem /M PBH ) is much smaller than one and nearly the entire energy density of the initial population of PBHs turns into radiation apart from a tiny amount. The PBH remnants have a number density, where ρ tot = ρ ϕ + ρ rad + ρ rem is the sum of energy densities for the runaway inflaton, radiation and remnants.

Remnants Abundance without PBH Domination
Let us estimate the remnants abundance ignoring any interaction between the different fluids and approximating again the redshift of each fluid with the expression ρ ∝ a −3(1+w) for constant w. Until the moment t rh the energy density of the PBH remnants increases relatively to the runaway scalar field as ρ rem /ρ ϕ ∝ a 3 and afterwards, that radiation dominates, it increases as wherec S = 2 1/4 (g(T rh )/g(T eq )) −1/4 Ω Matter /Ω DM and M eq ∼ 6 × 10 50 g is the horizon mass at the moment of radiation-matter equality. For times t < t rh the Hubble radius mass increases like M hor ∝ a 3 and given that M hor (t evap ) = (3/2)M 3 PBH /m 2 Pl the above expression can be simplified, (31) For the benchmark values β ∼ 10 −10 and M PBH ∼ 10 5 g the remnants are found overabundant, except for masses much less than the Planck mass. In Figure 4, we depict the allowed parameter space for the kination scenario with reheating via PBH evaporation when a remnant is left behind. Since Ω rem /Ω DM ≤ 1, a large part of the (M PBH , β) space is ruled out because remnants overclose the universe. Smaller β values are in accordance with the observed dark matter values for M rem ∼ m Pl . However, such β values are rather small for a successful reheating. We will further comment on the allowed β values in the following.

Remnants Abundance with PBH Domination
If the fraction of the universe energy density that collapses into PBHs is above β thresh then a PBH domination phase is realized. At the moment of PBH evaporation the remnants contain a small amount of the energy density, where T rh is given by Equation (28). Taking into account that ρ rad ∝ g * g −4/3 s a −4 , ρ rem ∝ a −3 and Equation (13) that gives the time of evaporation, we find the remnants abundance today, where we took the entropic and thermal degrees of freedom equal. After normalizing with benchmark values we rewrite, We note that the expression (33) can also be derived from Equation (31) after plugging in the expression for β thresh .
The PBH remnants constitute a significant part of the cosmic dark matter only if M BBH ∼ 10 6 (M rem /m Pl ) 2/5 . Since it must be M PBH < 10 9 g, the remnants cannot have a mass larger than 10 7 m Pl . Otherwise, the remnants have a small or negligible contribution to the total dark matter abundance.

Minimum Masses
Right after the end of inflation the Hubble horizon mass has a minimum value 4πρ end H −3 end /3 and steadily grows. We consider PBHs that form from the collapse of strong overdensities that reenter the horizon (different scenarios have been also considered, see, e.g., [66,67]). Hence, PBHs cannot have mass smaller than the value, where r * is the tensor-to-scalar ratio at CMB scales and, in the right hand side we used that H end < H * ≈ (π 2 A R * r * /2) 1/2 M Pl . We also normalized r * with the most recent bound from BICEP collaboration [34]. Plugging in the minimum PBH mass in Equation (34) (or in Equation (31) and assuming maximal β value so as to maximize Ω rem ) we find a limiting minimum mass for the remnants, M rem,min ∼ 10 4 γ 5/2 0.036 r *
The implication of this cosmological bound is that if a particle/object with mass less than M rem,min ∼ 10 4 GeV is found, then this particle/object is unlikely to be a remnant from PBH evaporation.
Recalling that the PBH mass has to be M PBH < 10 9 g, we conclude that, together with the lower bound (36), the remnants must have a mass in the range in order to have a non-negligible cosmic abundance.

Inflaton Residual Energy Density
BBN is a process with distinct observables that are sensitive to modifications of the background dynamics. Energy contributions beyond radiation must be sufficiently suppressed to prevent a too fast expansion during BBN, A kination regime has two important implications. First, ρ ϕ is non zero during BBN and second, the power spectrum of primordial GWs is blue shifted.
Parametrizing the extra energy components by an equivalent number of additional 2 neutrinos 3 the expansion rate has to be less than [70] H H rad 2 where H rad is the Hubble parameter for a universe filled only with the Standard Model radiation fluid, which is made of photons, electrons, neutrinos and their antiparticles. In our scenario, part of the total energy density 3H 2 M 2 Pl during BBN is reserved by the inflaton field and Equation (39) is equivalent to an upper bound on ρ ϕ /ρ rad at BBN epoch. For a field ϕ that fast rolls, ρ ϕ =φ 2 /2 a bound on the reheating temperature T rh O(10) MeV has been mentioned [71]. We note that the value 0.234 that normalizes Equation (39) comes from the CMB data [72]. ∆N ν,eff = 3.28 − 3.046 is the difference between the cosmologically measured value and the SM prediction for the effective number of neutrinos after e + e − annihilation [73,74]. Apart from BBN, CMB anisotropy measurements also constrain the effective number of neutrinos [75]. Increasing the radiation density (or ∆N ν,eff ) at the time of recombination increases the angular scale of the photon diffusion length on the CMB sky and this reduces power in the damping tail of the temperature spectrum. Therefore, at the time of recombination it has to be ρ GW /ρ γ ≤ (7/8)∆N ν,eff .
We now turn to gravitational radiation that gives the stringent constraints to our scenario.

Gravitational Waves and BBN/CMB Constraints
There are two sources of primordial GWs that we are going to discuss. First, the quantum tensor perturbations generated during inflation, called inflationary or primary GWs [76]. Second, the classical GWs generated by density perturbation, called induced or secondary GWs [77,78] and for a recent review see [79]. There are also the Hawking radiated gravitons [8,[80][81][82], as well as other aspects (see, e.g., [83]), whose implications are going to be examined in a separate work.

Reheating after a Kination Era
Kination era follows inflation in our scenario. If β < β thesh , Equation (26), kination era ends by the Hawking radiation that dominates the energy budget of the universe.

BBN/CMB Constraints on GWs from Inflation and a Kination Domination Phase
A stringent constraint comes from the gravitational wave energy density, which is enhanced in the GHz region during the kination regime [84][85][86][87]. The spectrum of GWs is described in terms of the fraction of their energy density per logarithmic wavenumber interval ρ −1 crit dρ GW (η, k)/d ln k. In terms of the tensor power spectrum it is written as, where H ≡ a −1 da/dη = 2η −1 /(1 + 3w) is the conformal Hubble parameter and the overline denotes oscillation average. Assuming that the modes with wavenumber k enter the horizon during radiation domination at the temperature T k , the present-day value is where η c denotes the conformal time that the generation of the primordial tensor modes has been completed. We take the present day thermal and entropy degrees of freedom to be g * (T 0 ) = 3.36 and g s (T 0 ) = 3.93, respectively, and equal at the high temperature T k . The present value of the radiation energy density parameter is h 2 Ω rad (t 0 ) ≈ 4.18 × 10 −5 . Detailed dependence on the degrees can be found in [88]. Primordial tensors generated by quantum fluctuations during inflation are where H inf is the Hubble parameter during inflation and n t is the tensor spectral tilt. On the other hand, the power spectrum of curvature perturbations is, Assuming that a scale k enters the horizon in the early universe either during the early kination phase or during radiation domination that follows, and utilizing the ratio P h /P R = 16 ≡ r, the energy density parameter of the inflationary GWs in the scales of interest is where we introduced the parameter α k ≡ H 2 k /H 2 * to include the decrease of the Hubble parameter during inflation. In the above equation we took into account that gravitational radiation increases like a 2 during a kination phase and that the scale factor scans scales as a ∝ k −1/2 during that epoch. Following the discussion of the previous section, the energy density of GWs that propagate inside the horizon at the time of big bang nucleosynthesis act as an additional radiation component increasing the background expansion rate, 3H 2 M 2 Pl = ρ rad + ρ GW . GWs produced prior to BBN do not alter BBN predictions and do not change the position and amplitude of the CMB acoustic peaks only if ρ GW /ρ γ ≤ (7/8)∆N ν,eff at BBN and recombination eras, respectively [89,90]. This is equivalent to a bound on the integrated energy density, For the lower limit of the integral the scale k m = k CMB ∼ 10 −2 Mpc −1 is chosen, that extends the integral to wider range of frequencies, compared to k m = k BBN , giving a stringent constraint. An inflationary era is a quasi-de Sitter phase and it is roughly H inf ∼ constant, slightly decreasing towards the end of inflation. This means that the tensor spectral tilt in Equation (42) is n t ∼ 0. We can thus assume a roughly constant value Ω GW (k, η 0 ) ≈ Ω GW ≈ 8×10 −16 r * , in the period between the moment k −1 CMB exits the horizon and the end of inflation. In terms of inflationary quantities it is k end = k * e N * (H end /H * ). k end can also be written in terms of reheating temperature and for a general equation of state according to the formula, For a kination domination regime that transits into radiation domination it is k end = k(M end , T rh , 1). The fractional energy density of GWs that impacts BBN is ≈ 8 × 10 −16 r * 5 × 10 12 T rh 10 6 GeV

+35 ln
T rh 10 6 GeV where α end ≡ H 2 end /H 2 * < 1 and we assumed a kination postinflationary phase until t rh . The CMB bound (45) implies that the reheating temperature for the kination domination models has to be larger than, The reheating temperature for a kination phase terminated by PBH evaporation is given by (27) and a lower bound on the formation rate is obtained, having assumed that roughly H end ∼ 0.1H * . Smaller values for β mean smaller amounts of Hawking radiation and a delayed reheating of the universe that experiences an extended kination domination phase. Comparing with β thresh (26) we see that it is β < β thresh only for small PBH masses or for small r * values. The bound (49) is relaxed if we add an initial amount of radiation in the postinflationary universe, due to, e.g., gravitational reheating, see Appendix A.
The aferomentioned results are modified if we depart from the simplistic approximation of w = 1. In that case an effective average equation of state during the entire reheating period should be considered,

BBN/CMB Constraints on Induced GWs from Kination Era
The induced GWs are sourced gravitational waves. The leading source is the scalar perturbations that evolve according to the equation, Here we have implicitly assumed an negligible anisotropic stress and introduced the variable x ≡ kη. Φ is the scalar transfer function defined Φ k (η) ≡ Φ(x) φ k where φ k is the Fourier mode of the primordial scalar perturbations and Φ k of the Bardeen potential. For the kination domination scenario, w = 1, it is Φ(x) = 2 x J 1 (x), and if we consider a transition into RD the scalar transfer function is given by [91], where x rh corresponds to the moment of reheating. The coefficients C 1 and C 2 are determined by the continuity of the potential and its derivative at the point of the transition, after expressing J 3/2 and Y 3/2 in terms of spherical Bessel functions. Φ(x) remains constant as long as x 1. According to Equation (52) the gravitational potential starts decaying at the horizon crossing roughly as x −3/2 experiencing maximal pressure during KD. After the transition to the RD phase the decay is even faster, x −2 . At some time t c , that might be during the kination or radiation era, the source Φ(x) has decayed and the production of IGWs ceases.
The power spectrum of the induced gravitational waves is expressed as a double integral of the curvature power spectrum and a tensor transfer function, where the over-line denotes the oscillation average. The tensor transfer function is given by where I(u, v, η, k) is the oscillation average of a kernel function composed of the Green's function . For kination domination it is ν = 0 and the Bessel functions J 0 and J 1 do not have a closed form representation and can only be written as infinite series. Therefore the transfer function T in kination domination, contrary to radiation domination, cannot be written in closed analytical expressions and the computation has to be conducted numerically.
For the scenario of the kination domination that transits into radiation and for a monochromatic power spectrum of scalar perturbations modeled by a delta-distribution peaked at k p , δ(k) = A 0 δ(ln(k/k p )), the energy density parameter at the moment t c is [91,92] where the unit step function Θ is for the conservation of momentum cutting off tensor modes with k > 2k p . For η c η rh no sharp peak in the IGWs appears even for a monochromatic P R (k). After the moment t c the GWs propagate freely. The present-day energy density parameter of the induced GWs is given by Equation (41). The energy density of induced GWs for a broad spectrum of curvature perturbations can be found from the expression (55) by the where is the width of a Gaussian distribution [91]. Similarly here, during kination domination the growth of induced GWs relative to the background is proportional to a 2 , given by the expression (44). In terms of conformal time it is This growth stops after the transition to the radiation domination phase. If it is η entry η rh we can neglect the effect from the RD era on the shaping of induced GWs, because by that time the gravitational potential sourcing the tensor modes is negligibly small. In the case of induced GWs, it is much different from the GWs produced by the de-Sitter stage of inflation. The spectrum of induced GWs has a prominent peak at wavenumbers relevant to PBH production. If there was not this amplification of the scalar power spectrum, the induced GWs would be irrelevant for our discussion. Notwithstanding, the existence of an early kination phase implies that the energy density of the produced induced GWs, associated with PBH production, becomes enhanced and might backreact on the geometry.
The induced GWs do not spoil BBN/CMB if they satisfy the bound (45). In particular, the increase of the expansion rate does not alter the position and amplitudes of the acoustic peaks of CMB if A numerical integration of a spectrum of induced GWs with peak at k ∼ 10 22 Mpc −1 and reheating temperature T rh ∼ 10 9 GeV associated with PBHs with mass M PBH ∼ 10 5 g, generated by a sharp peak at the scalar spectrum, violates the above bound five orders of magnitude. In order to perform the integration and demonstrate the conflict of induced GWs with BBN predictions, we will assume a rough, top-hat approximation about the peak for the spectrum of induced GWs. Let us assume that the energy of GWs is stored mainly in a narrow wave band (k 1 , k 2 ) with central wavenumber k p that enters the horizon at η entry . The top-hat approximation for the GW spectrum with amplitude Ω IGW = A IGW in the interval (k 1 , k 2 ) around k p gives, A IGW is related to the scalar power spectrum amplitude and we can consider that Ω GW (k p ) ∼ 10 −2 A 2 R , justified by our analytic and numerical study, where A R ≡ P R (k p ). Assuming a width k 2 /k 1 ∼ 10, we can express the above bound in terms of the reheating temperature, where we used that k rh ≈ 2 × 10 7 (T rh /GeV) Mpc −1 and Equation (46) to express the wavenumbers. This bound is much stringent than the bound coming from the inflationary GWs, Equation (48). For reheating following a kination domination phase β has to be larger than and at the same time less than β thres . This happens only for too small M PBH masses, less than 0.1 grams, which contradicts the bound (35). In other words, β has to be large enough so that it is always β > β thesh and a PBH domination phase is unavoidably realized.

Reheating after a PBH Domination Era
For a PBH domination phase the above bounds alter significantly. PBH domination is realized for β > β thresh and such values are obtained for a slight increase of the curvature power spectrum, see Equation (3). A PBH domination phase has been discussed in different contexts and has profound implications on GW signals, see, e.g., [8,[93][94][95][96]. In the following we will mainly discuss the evident impact of a matter domination phase on the GW energy density parameter, which is a a −1 decrease, considering inflationary and induced GWs.

BBN/CMB Constraints on GWs from Inflation and a PBH Domination Phase
A PBH domination era suppresses the energy density of GWs with respect to the background. The reheating temperature of a PBH domination phase is larger than that of a kination phase for the same PBH mass. Moreover, after inflation, the pre-radiation phase is partitioned between kination and PBH domination, see Figure 3. The ratio of the scale factor at BBN and at the end of inflation is analyzed as where kination domination lasts N KIN ≡ ln(a end /a eMD ) efolds, PBH domination (early matter domination) lasts N eMD ≡ ln(a eMD /a rh ) efolds, and a pre-BBN radiation domination N eRD ≡ ln(a rh /a BBN ) efolds. Accordingly, the GW energy density during BBN is decomposed as For inflationary GWs, that are roughly scale invariant, Ω GW (k, η 0 ) ≈ Ω GW , the above expression is written in terms of efolds, PBHs form during kination at t form about N form efolds after inflation, For M PBH ∼ 10 5 g it is N form ∼ 4 and for M PBH ∼ 10 9 g it is N form ∼ 7. PBHs dominate the energy budget of the early universe after ∆N dom efolds, and it is N KIN = N form + ∆N dom . PBH domination lasts N eMD folds until reheating that coincides with PBH evaporation. It is a eMD /a form ∼ (t eMD /t form ) 1/3 = (γβ) −1/3 , a rh /a eMD ∼ (t evap /t eMD ) 2/3 and Summing up the expressions for the amount of efolds, the bound on kination era (64) is recast into a bound on β value, This bound selects a wide part of the PBH domination (M PBH , β) parameter space where the tilted GWs from inflation do not change BBN predictions. Note that the lower limit of integration interval (62) [k BBN , k rh ] slightly modifies the result, thus the above bound applies also if it is extended to the CMB scale, [k CMB , k rh ].

BBN/CMB Constraints on Induced GWs from Kination Era and a PBH Domination Phase
The early matter domination era realized by PBHs suppresses the relative enhancement of induced GWs relative to the background energy density, see Figure 5. For the simplistic tophat approximation for the energy density spectrum of induced GWs the suppression is given by the exponential in front of the brackets, where we considered modes that propagate inside the horizon at the time of BBN/CMB. Induced GWs, associated with the PBH formation, do not change BBN predictions for This is the most stringent constraint posed on the scenario of runaway inflation reheated via PBH evaporation. We illustrate this constraint in 6. Note that the constraint (70) encompasses the constraint (68) from inflationary GWs.    Figure 5: In the left panel the early spectrum of induced GWs is depicted, produced during a kination domination phase by a monochromatic P R (k) and associated with PBHs with mass M PBH = 10 5 g. In the right panel the present day frequency spectrum is depicted for M PBH = 10 5 g (blue curve) and M PBH = 10 8 g (orange curve) and for β ∼ 10 −10 . Note the change of the shape of induced GW spectrum due to kination and early matter domination phase. The dotted curves depict the corresponding isocurvature-induced GWs, following Ref. [96]. The colored areas are enclosed by the sensitivity curves of LIGO/Virgo, Einstein Telescope, LISA and DECIGO.

Induced GWs and Detection by LIGO/Virgo, Einstein Telescope
The induced GWs associated with the production of mini PBHs have large frequencies that can be probed only from LIGO/Virgo experiment and the designed Einstein Telescope. Focusing on the frequency band around 100 Hz, that corresponds to wavenumbers about k EX ∼ 10 16 Mpc −1 , the induced GWs must have an amplitude, Ω(t 0 , k EX ), below the sensitivity curve (SC) of the operating experiments, Ω −1 SC . It is exciting that current or future GW experiments can rule out or support this cosmological scenario, probing part of its parameter space. Approximately, a GW wave experiment can probe the production of mini PBHs during a kination phase followed by PBH domination and reheating if Advanced LIGO/Virgo will reach a sensitivity if h 2 Ω −1 SC ∼ 10 −10 and the forecast for the Einstein Telescope is h 2 Ω −1 SC ∼ 10 −13 . Although the β values required for detection (71) are too small and already ruled out by BBN/CMB arguments it is an exciting perspective that testing experimentally exotic cosmological scenarios is in principle possible, see Figure 5. , that gives viable cosmology, lies in the PBH domination region. The left part of the parameter space (gray area) is excluded by the energy scale of inflation that determines the minimum PBH mass, shown for two different r * values. The right part of the parameter space is excluded by the requirement of radiation domination during BBN (red area). The induced GW constraints due to BBN/CMB exclude a great part (green area) of the parameter space as described in the text. The red dotted lines indicate the reheating temperature.

Isocurvature-Induced GWs
Ref. [95] noted that the isocurvature perturbation of PBHs (and of dark matter in general [97]) induces additional tensor modes. As discussed in Section 2, PBH formation is a rare event and takes place at those regions of space where the density perturbation is above the threshold δ c . Their space distribution is random and follows a Poisson type distribution [95,98,99]. For distances much larger than their mean separation,r = ((3/4π)n −1 PBH ) 1/3 , the PBH gas can be described as a continuous fluid that acts as an isocurvature perturbation S = 3(ζ PBH − ζ ϕ ) where ζ = −Φ + Hδρ/ρ is the curvature perturbation on uniform-energy density hypersurfaces and Φ is the Bardeen potential [100]. For the kination and PBH fluid they, respectively, ζ ϕ = −Φ + (1/6)δ ρ ϕ/ρ ϕ , ζ PBH = −Φ + (1/3)δρ PBH /ρ PBH , thus, at the formation time and with β 1 it is where the second equality follows from the isocurvature property δρ ϕ + δρ PBH = 0. S and separately ζ ϕ , ζ PBH are conserved on superhorizon scales. Initially the curvature perturbation receives its leading contribution from the scalar field ζ ≈ ζ ϕ , it evolves with time and in the PBH domination phase it is ζ = ζ PBH ≈ ζ ϕ + S/3. The comoving curvature perturbation R is related to ζ by a gauge transformation and coincide in the limit k → 0. During PBH domination it is R = −ζ ≈ 5Φ/3. The evolution of the Bardeen potential is given by the equation where the sound speed is The early isocurvature sources tensor modes that have size given by the formula (53) after replacing P R → P S , where [95,96] is the initial isocurvature perturbation power spectrum. k UV = a/r is an ultra-violet cutoff below which the grainy structure of the PBH fluid becomes important.
According to [101,102], and applying the notation to our scenario that discusses a kination era, the Kernel I(u, v, η, k) (54) is split into three contributions, I = I KIN + I eMD + I RD , from kination era, PBH domination and radiation era, respectively. I RD gives the dominant contribution to isocurvature induced GWs due to non-zero potential Φ field at small subhorizon scales. We recall that the perturbation Φ decays during kination domination but it remains constant during PBH domination and decays completely only after reheating of the universe. Ref. [96] analyzed a similar to our case problem with the difference that an early radiation era precedes PBH domination. They find a peak of the isocurvature induced GWs around k ∼ k UV with amplitude that depends on β 16/3 and M As a rule of thumb we can apply the above bounds to kination domination scenario, see Figure 6, because isocurvature-induced GW receive their main contribution at the late stages of PBH domination and, particularly, at the beginning of radiation domination. We remark that the sudden transition [101] from PBH domination to radiation, that gives the maximal GW and detectable signal, Figure 5, is realized only for a universal sudden evaporation that is possible only for a monochromatic PBH mass spectrum. Moreover, as noted in [95,96], the density contrast of PBH fluctuations exceeds unity at the time of evaporation and the linear analysis should break down at some point. A non-perturbative methodology has been developed in Ref. [103] where the deviations from ashpericity due to the absence o pressure has been explicitly considered. An overdensity of size k −1 evolves non-spherically and eventually collapses at a pancake configuration at the time t col . The spectral energy density parameter is given in terms of the quadrupole tensor Q ij and a probability density function F D (α, β, γ), that describes the size of asphericity of the overdensity, with the precise expression found in Ref. [103]. We leave a detailed analysis of this particular problem for a future work.

PBHs from Runaway Inflation Models
The inflationary potential that we consider is a runaway one without a minimum. Its precise form plays a critical role for the determination of both the inflationary and postinflationary cosmic evolution. The number of efolds N * has to be large, compatible with a kination domination phase that follows inflation. Furthermore, the position and characteristics of the feature that amplifies the power spectrum specify the mass and the abundance of the PBHs. A runaway potential is acceptable only if it can realize a sufficient reheating of the universe and, here, this is realized by the production and evaporation of mini PBHs. The construction of runaway inflation models that induce PBH production and act as quintessence models today is very challenging. This unified description works only for particular choices of the parameters that affect horizontally the observables in the cosmic time range from 10 −35 s to 13.8 billion years. To be specific, we mention that the number of efolds N * , related to the CMB scale, depends on the reheating temperature which is not a free parameter but depends on β and M PBH . This is in sharp contrast with conventional inflation models where the inflaton field decays perturbatively about the vacuum. A particular PBH mass M PBH is related to a specific k, where P R (k) is amplified, only after N KIN , N eMD are known. Hence, the characteristics of the peak in the power spectrum determine the mass of the evaporating PBHs and the reheating temperature of the universe and the remnant abundance! Additionally, the tail of the potential might lead to the observed late time acceleration of the universe.
In single field models, the special feature that amplifies the power spectrum can be either an inflection point [104][105][106] or a step-like transition in the inflationary plateau [107,108] 4 . The presence of step-like, or generally sharp, features in the inflation potential results in amplification and, additionally, in oscillatory patterns in the curvature power spectrum [108,111]. The potential might feature more than one steps, so that the relevant enhancement in the power spectrum becomes strong enough.
In the following subsections we will discuss in some detail the inflection point scenario.
horizon and its value freezes out is After the numerical computation of the Mukhanov-Sasaki equation the P R at all scales is obtained, see Figure 8.  (86) and potential depicted in Figure 7. The power spectrum has a peak at large wavenumbers triggering the production of mini PBHs during kination regime.
The position of the peak determines the PBH mass, but the knowledge of the duration of kination and PBH domination phase is also required. According to the discussion in Section 5.2 the peak of the power spectrum is at the wavenumber, where k rh ≈ 2 × 10 7 T rh GeV −1 Mpc −1 . Inflation ends at the wavenumber k end = k * e N * H end /H * .
A power spectrum with a peak at large wavenumbers, k k * , is welcome because it is not spoiling the n s and α s values measured at k * . Notwithstanding, shifting the peak at large wavenumbers does not render the P R (k) free from constraints. Hawking radiation might affect BBN and CMB observables and upper bounds on the P R (k) at large k-bands. Therefore, one has to be careful with the width and the tail of the power spectrum peak since PBHs with a distribution of masses might be generated. The population of (non-monochromatic) PBHs with mass M PBH 10 9 g, if not zero, has to be small enough. The stringent constraint is at the mass M PBH = 2.5 × 10 13 g, where the variance of the density perturbations has to satisfy [112], This bound is easily satisfied in our scenario that we demand a curvature power spectrum peak that generates ultra-light PBHs with mass M PBH 10 9 g, Figure 8, but this is not always true for wide power spectra that generate a distribution of PBH masses, in particular during an eMD era. For a wide curvature power spectrum with a peak close to the maximum k values the growth of P R (k) [113] together with the reheating temperature [112] have to be checked.

Inflection Point
Inflaton potentials with an inflection point can produce a power spectrum of curvature perturbations P R (k) with large hierarchies. At the region of the inflection point the amplitude can be amplified by many orders of magnitude due to the acceleration and deceleration of the inflaton [104][105][106]. Such a model may arise from the α-attractors models [32] as described in Ref. [31,35,55]. Successful CMB observables and a significant PBH population can be generated if the power spectrum peak is positioned at large k.
Furthermore, we note that peaks in the power spectrum can also be related to the existence of one or more sharp drops (steps) in the inflaton potential [107,108]. One could replace the inflection point with one or more sharp steps sufficiently close to each other and an analogous results can be attained.
The effective Lagrangian for the inflaton field ϕ in α-attractors reads where ReΦ = φ = √ 3 tanh(ϕ/ √ 6α) is a chiral superfield and we took M Pl = 1. Polynomial, trigonometric and exponential forms for the function f (φ) can feature an inflection point plateau sufficient to generate a significant dark matter abundance in accordance with the observational constraints [31]. An explicit working example is a combination of exponentials of the form, that generates the potential V (ϕ) = |f (φ/ √ 3)| 2 . In the above expressions we have taken α = 1, but other choices also work. The parameter f 0 is redundant in the sense that it can be absorbed by c 0 , c 1 and c 2 , but it is kept for numeric convenience.
For φ → √ 3 the potential drives the cosmic inflation and the CMB normalization gives the first constraint for the parameters. Furthermore, the n s and dn s /d ln k values have to be in accordance with the CMB data, given by Planck 2018 collaboration [33] as well as the recent bounds by BICEP [34]. Requiring PBH production with the right abundance chooses specific values for the parameters that shape the inflection point plateau. In the following we also demand the value of the runaway potential at the field value today to be slightly above zero in order to identify it as the dark energy.
As commented before, the precise determination of the parameter values is a subtle numerical process. The M PBH and the efold number values N * , ∆N dom , N KIN select the k-position of the P R (k) peak. The required β value determines the amplitude of the P R (k) peak which is found after a careful and precise selection of the potential parameters due to the exponential sensitivity of β to the amplitude. In our scenario particular attention has to be paid at the number of efolds N * that has to be in agreement with the postinflationary cosmic evolution that involves a kination and a PBH domination phase.

From Inflation to Dark Energy
A runaway inflaton has interesting implications not only for the early but, also, for the late universe cosmology. Regarding the early universe, the inflaton does not decay oscillating about a vacuum and a period of kinetic energy domination of the scalar field can be realized. It has distinct early universe phenomenology due to the stiff equation of state w ∼ 1 and the expansion rate is reduced. The most notable effect is that the tensor perturbations attain a spectrum with more power at small scales. Furthermore, the fact that a tiny trace of the inflaton field potential energy density might remain non-zero until the present-day universe draws the attention to these inflationary models that can play the role of quintessence. In the framework of α-attractors kination models have been constructed in [26][27][28].
After inflation the scalar field rolls down the runaway potential until it freezes due to Hubble friction at some value ϕ F . The residual potential energy at the frozen value acts temporarily like a cosmological constant and must be compatible with early and late time cosmological observations, V (ϕ F ) 10 −120 M 4 Pl . When the upper bound is saturated, the runaway inflaton field is identified as the quintessence field that drives the accelerated expansion today. The inflation runaway potential originating from the theory (86) can drive late time inflation after a proper choice of parameters. We ask for V (φ) → 0 as ϕ → −∞, or equivalently φ → − √ 3 in order that the potential is positive-definite. At the value ϕ F the potential is flat enough and with small enough slow-roll parameters to implement a wCDM quintessence model. The hierarchy of energy scales between α-attractors inflation and the present dark energy implies a tuning of the potential energy value at ϕ F , that is λ 1 ∼ 108 ln(10)/4 ∼ 62. This condition is a third constraint to the parameters of the potential, together with the CMB normalization and the positivity of the potential energy density. A working example is given by the set of parameter values c 0 = −8.7 × 10 −27 , c 1 = 0.1045, c 2 = −4 × 10 25 , λ 1 = 62.2 and λ 2 = −4430.97, f 2 0 = 3.115 × 10 −62 M 4 Pl , φ P = 0.995M Pl and initial field value ϕ CMB = 11.8M Pl . Details can be found in Ref. [35].
In summary, at the field values ϕ CMB and ϕ P inflation produces the seeds of CMB anisotropies and PBHs, respectively, and at φ end inflation ends and a kination stage commends. Later PBH form, dominate the energy density of the universe and evaporate reheating the universe. At the field value ϕ F the field freezes and its potential energy plays the role of the dark energy in the universe. It is remarkable that the simple theory (86) provides a complete example that can explain the main cosmological observations, for appropriate though highly tuned values of the parameters. Next, we will give an approximate quantitative description of the aforementioned stages of this cosmological model, that we also illustrate in Figures 9 and 10.

The Evolution of the Runaway Inflaton
Let us work out some simple analytic expressions and describe approximately but quantitatively the post-inflationary evolution of the field ϕ. After inflation ϕ fast-rolls the potential and a stage of kination domination takes over whereφ 2 /2 V (ϕ). For negligible potential energy the evolution of the scalar field is given by the system of equations From the Friedmann equation we obtain HM Pl ≈ ±φ/ √ 6, where the minus is for negative field velocity. After an integration the evolution of the field value is obtained, where t init can be defined as the moment thatφ 2 dominates the energy density. The initial velocity is found from the Friedmann equation and the Hubble parameter, . At the end of inflation t end it is V =φ 2 , but soon afterwards kination regime takes over, see Figure 9. We can approximate t init with the moment inflation ends, where we also considered negative initial field velocity. At the moment t form PBHs form and later at t evap they evaporate.
If there is no PBH domination phase, the universe becomes radiation dominated at the moment t rh and kination regime ends, where t end and t rh define the duration of the kination phase. Until reheating it is approximately a ∝ t 1/3 the reheating moment is found to be t rh = (Ω rad (t evap )) −3/2 t evap , where Ω rad (t evap ) = (3/2)γ 2 βM 2 PBH /m 2 Pl , given by Equation (25). Utilizing the above results, the field value at the reheating moment, ϕ(t rh ) ≡ ϕ rh , is found. After reheating, we can follow the same steps and assume that the scalar field keeps running away its potential with the Hubble parameter given by H ≈ 1/(2t) (radiation domination). Thus, in the reheated universe the scalar field evolves as where we approximated H 2 (t rh ) ≈φ 2 rh /3M 2 Pl at the reheating time. Notice that the field is not rolling fast the runaway potential but slows down and at late times t t rh the field freezes having covered a 2M Pl / √ 3 distance in field space from the value ϕ rh . Next we add a PBH domination phase, that is necessary for viable cosmology, according to the discussion in Section 5.

PBH Domination Phase
In the viable case that a PBH domination phase takes place, the kination phase ceases earlier, at the time t eMD , where eMD stands for early matter domination. There is the following hierarchy of moments, PBHs form with a nearly instantaneous collapse at the time M PBH /γ m 2

Conclusions and Discussion
In this work, we investigated the reheating of runaway inflation models via the evaporation of primordial black holes. A runaway inflaton potential has no minimum and the inflaton field does not decay but gradually loses potential energy. Since the standard mechanism of reheating cannot be realized in this type of models, another source of entropy production has to be introduced. Black holes are objects that exist in the present-day universe having a quite significant population and it is plausible to presume that they might have existed in the early universe as well. Mini black holes are motivated candidates because the early universe, although homogeneous at large scales, might be characterized by strong small scale fluctuations. If these fluctuations have size less than about k −1 ∼ 10 −19 Mpc and density contrast that exceeds a threshold value, gravitational collapse takes place and primordial black holes with an ultra light mass are formed. These mini black holes are ephemeral and promptly transform their entire mass into thermal degrees of freedom. Black hole evaporation is therefore an alternative reheating mechanism that can render, the otherwise problematic, runaway inflation models cosmologically viable. Primordial black holes with mass M PBH < 10 9 g have to be generated with a significant population in the early universe in order to implement a successful reheating. The abundance of the PBHs is determined by the parameter β that measures the fraction of the universe energy density that collapsed. At the same time with the PBH generation a substantial amount of gravitational radiation is produced, called induced or secondary GWs. Contrary to the ephemeral mini PBHs, the associated induced GWs propagate in the early and late universe and have significant observational implications. First, their energy density contributes to the expansion rate and impacts BBN and CMB. Second, the secondary GWs are in principle detectable by current and near future gravitational wave detectors.
The size of the β value plays a critical role and determines the details of the cosmological evolution. Small β values imply a small population of PBHs and an extended kination phase. Kination domination ends due to the radiation produced by PBH evaporation at a temperature T rh ∝ β 3/4 . Although reheating temperatures safely larger than O(10) MeV can be realized, the extended kination scenario is ruled out. The energy density of GWs, of both primary and secondary, becomes enhanced and alters BBN and CMB. Noteworthy, the stringent constraint on β comes from the secondary GWs. The β value has to be large enough, β 10 −18 − 10 −10 , so that PBHs dominate the energy density of the early universe, Equation (70). A PBH domination compensates the relative increase of the GW energy density that occurs during kination domination. In the present-day universe the induced GWs have a non-negligible amplitude but too large frequencies, which generally lay outside the sensitivity band of the current or near future detectors. A PBH domination phase comes with an isocurvature perturbation which, in turn, induces an extra component of secondary GWs [95] that constrains the maximum β value [96]. The final allowed parameter space (M PBH , β) is summarized in Figure 6. It is exciting that the allowed parameter space of our scenario can shrink further from data coming from GW experiments as well as from cosmological precision measurements of the effective relativistic degrees of freedom ∆N ν,eff during BBN or photon decoupling.
A very interesting implication of reheating runaway inflation models via PBH evaporation is that, on the one hand, the inflaton field has a non-zero potential energy in the presentday universe and, on the other, PBH evaporation might leave a remnant mass behind. The former form of energy can play the role of dark energy. The latter can constitute the dark matter, or part of it, in the galaxies. The mass of the remnant is expected to depend on the unknown physics that operates at the Planck energy scale or somewhere close to that scale. PBH remnants have a significant cosmological abundance for a particular remnant mass M rem which is proportional to M 5/2 PBH , see Equation (34), and lies in the range between 10 −15 m Pl and 10 7 m Pl . Due to the near one-to-one cosmological correspondence between M rem and M PBH a potential discovery of those exotic remnants will connect us, among others, with the early universe and PBHs. The reheating temperature of the universe, which also depends on the PBH mass as M −3/2 PBH , can range from few MeV up to 10 10 GeV. Let us note that the reheating temperature has important implications for other essential processes such as baryogenesis and dark matter production, topics which are outside the scope of this work.
Last, we engaged in an explicit inflationary model built in the framework of α-attractors [35] to implement this cosmological scenario. Mini PBHs are totally compatible with the inflationary observables, in particular the spectral index value, since the power spectrum has to be enhanced at very large wavenumbers far away from the CMB scale. It is remarkable that the minimal theory (86) of a single field with a non-canonical kinetic term and without the need of extra interactions, can explain the basic cosmological observations. It gives cosmic inflation, reheats the universe via the production and evaporation of mini black holes, predicts the presence of dark matter if the black holes leave a remnant behind, and drives the late acceleration of the universe via the residual vacuum energy of the very same scalar field. Let us comment that, although extremely economic, we advocate the presence of new physics that could give rise to such a type of quintessential inflationary field with inflection point or strong features in its inflationary trajectory.

A Gravitational and Instant (P)reheating
Let us briefly describe two other non-conventional reheating mechanisms: gravitational reheating and instant (p)reheating. The first mechanism is inefficient regarding the current observational constraints, while the second one depends on the coupling of the inflaton to another scalar field. This is why we regard them beyond the main line of the paper.

A.1 Gravitational Reheating
Gravitational particle production [12,26,65] is an inevitable but feeble reheating mechanism. Particle production of all light fields (with masses less than the Hubble scale) takes place due to the change of the spacetime metric at the end of inflation. This is essentially Hawking radiation in de Sitter space, which generates a radiation bath of temperature T H = H/2π. We assume that the energy density at the end of inflation is dominated by the kinetic energy of the scalar field (kination domination). The radiation density produced by gravitational reheating at the end of inflation is where g end * is the effective number of relativistic degrees of freedom at the energy scale of inflation and q ∼ 1 is an efficiency factor.
Using Equation (98), we have for the radiation density parameter at the end of inflation: where we have used that ρ tot = 3H 2 M 2 Pl . The universe is reheated when radiation takes over and dominates the kinetic density of the scalar field. This definitely happens, because the energy density of kination is red-shifted away much faster than that of radiation as ρ KIN ∝ a −6 and ρ rad ∝ a −4 . Using that a ∝ t 1/3 during kination, it is easy to show that the time when reheating occurs is t rh = (Ω end rad ) −3/2 t end .
As Ω rad = ρ rad /ρ KIN during kination, we can easily find that ρ rh KIN = (Ω end rad ) 3 ρ end φ . Thus, considering that radiation is thermalized by the time it comes to dominate the Universe, we find that the reheating temperature is which, using Equation (99), can be written as However, we can easily see that this reheating temperature does not obey the BBN constraint (48). For example, if we set H end ∼ 10 12 GeV, g end * = O(100) and g rh * = 10.75, Equation (102) gives T rh ≈ 10 4 GeV, while Equation (48) demands that T rh 10 7 GeV. Consequently, gravitational reheating is not efficient and we should rely on other mechanisms.

A.2 Instant (P)reheating
Instant preheating can work both in usual inflationary models where V (φ) has a minimum and in quintessential models such as those that we examine in this paper [13,27]. The basic assumption is that the inflaton field φ is coupled to another scalar field χ which is also coupled to a fermion field ψ. The interaction Lagrangian density is where g and h are perturbative coupling constants. In order for particle production to occur, we must have [27] |φ| which gives the range for φ: The density of the produced χ particles can be shown to be [13] ρ InP χ = g 5/2 |φ InP | 3/2 φ InP 8π 3 .
where the superscript "InP" denotes instant preheating values. In fact, we expect only the particles produced near the end of the particle production window in (105) to contribute significantly to ρ InP χ [27]. Thus, taking φ 0 = 0, we can set φ InP = |φ InP | g . Consequently Equation (106) can be written in the simpler form where we have considered that χ-particles decay to radiation instantaneously. For the reheating temperature in this model we have ρ rh rad = ρ rh φ = ρ InP φ (Ω InP rad ) 3 . Thus, where Ω InP rad can is found to be Ω InP rad = g 2 /(4π 3 ). Consequently in instant preheating T rh depends strongly on the value of the coupling constant g and takes larger values as we increase g. However, the region of g so that all cosmological constraints are satisfied is given approximately by [27] 10 −4 g 10 −2 .