Primordial black hole formation during slow-reheating: a review

1 Introduction

While the existence of supermassive and stellar-mass black holes today is thoroughly demonstrated, the possibility of a third species of black holes has been hypothesized in recent decades: the formation of primordial black holes (PBHs) during the early stages of evolution of the Universe [see, e.g., Refs. (Yu, 2010; García-Bellido, 2017; Sasaki et al., 2018; Carr and Kuhnel, 2020; Anne, 2021; Carr et al., 2021; Albert et al., 2022; Özsoy and Tasinato, 2023) for recent reviews]. These objects may be key to understand fundamental aspects of cosmology and particle physics.

One epoch in which primordial black holes might have emerged is the reheating phase that immediately followed the period of cosmic inflation [see, e.g., Refs. (Allahverdi et al., 2010a; Amin et al., 2014)]. The inflationary paradigm postulates that the early Universe experienced a period of rapid and exponential expansion in its earliest moments [see, e.g., Refs. (Linde, 1984; Keith, 1990; Lyth and Riotto, 1999; Liddle and Lyth, 2000; Baumann, 2011; Martin, 2020; Odintsov et al., 2023)]. This inflationary period played a crucial role in explaining the observed flatness of the Universe and the uniformity of the cosmic microwave background radiation, also offering an elegant explanation for the large-scale structure of the Universe [see, e.g., Ref. (Langlois, 2010; Alberto Vázquez et al., 2020)]. After inflation, the Universe entered a reheating phase characterized by the decay of the inflaton field, resulting in the transfer of its energy to matter and radiation. This process eventually led to the emergence of a hot and dense environment, providing the necessary conditions for the subsequent stages of cosmic evolution.

The process of reheating represents a crucial yet relatively understudied chapter in cosmology. It not only determines the thermal properties of the early Universe but also plays a fundamental role in particle production and the subsequent formation of cosmic structures [see, e.g., Refs. (Allahverdi et al., 2010a; Allahverdi et al., 2010b; Amin et al., 2014; El Bourakadi, 2021)]. Among the various possibilities that arise during reheating, the formation of PBHs has garnered significant attention due to their particular characteristics and potential cosmological implications [see, e.g., (Harada et al., 2016; Carr et al., 2017; Harada et al., 2017; Juan Carlos Hidalgo et al., 2017; Carr et al., 2018; Martin et al., 2020; Carrion et al., 2021; De Luca et al., 2022; Padilla et al., 2022; Bhattacharya, 2023; Harada et al., 2023; Hidalgo et al., 2023; Padilla et al., 2023)].

A key aspect of PBH formation during reheating lies in the collapse threshold for the density contrast compared to the formation process in a radiation-dominated background. During the inflationary phase the rapid expansion of space stretches small-scale quantum fluctuations to macroscopic and cosmological scales. Once inflation ceases and scales slowly re-enter the horizon, these fluctuations undergo substantial growth due to the dynamics of the reheating process. Consequently, localized regions with remarkably high density emerge. If the density within these regions surpasses a critical threshold, they may collapse and form PBHs.

Understanding the formation and characteristics of PBHs during the reheating phase poses a challenge spanning the fields of cosmology, particle physics, and astrophysics. By investigating the mechanisms of PBH formation, their mass spectrum, and cosmological abundance, we may gain insights into the fundamental physics that governed the early Universe. Moreover, such PBHs hold potential in elucidating intriguing astrophysical phenomena, including dark matter [see, e.g., (Ivanov et al., 1994; Frampton et al., 2010; Carr and Kuhnel, 2020; Anne, 2021; Villanueva-Domingo et al., 2021; Carr and Kuhnel, 2022)], gravitational waves [see, e.g., (Eroshenko, 2018; Franciolini, 2021; Papanikolaou et al., 2021; Ballesteros et al., 2022; Bavera et al., 2022; Papanikolaou et al., 2023)], and the origins of supermassive black holes [see, e.g., (Duechting, 2004; Kawasaki et al., 2012; Luis Bernal et al., 2018; Dolgov et al., 2019)]. All these possibilities further emphasize the significance of PBHs in our quest for a comprehensive understanding of the Universe.

In this paper we intend to present a comprehensive review of the formation of PBHs during the reheating epoch. Starting with the theoretical foundations, we discuss the various formation mechanisms. By studying the details of PBH formation during reheating, we aim to contribute to the understanding of the early Universe and shed light on the nature of black holes, providing potential avenues for future research and novel insights into the cosmic dark sector.

2 Inflation, preheating and the reheating epochs

2.1 Inflation setting the initial conditions for reheating

As mentioned above, cosmological inflation [see, e.g., (Alan, 1981; Langlois, 2010; Alberto Vázquez et al., 2020)] refers to a period of accelerated expansion of space. In the framework of general relativity, inflation usually stipulates the existence of a scalar field as the dominant energy content of the Universe during this period. In its simplest form, the inflationary scenario is described by the action

For the inflaton φ to drive the inflationary epoch, its energy must be dominated by a nearly constant potential energy V(φ). In this case, the inflaton field behaves effectively like a cosmological constant, causing the Universe to expand exponentially. When the kinetic part of the inflaton field is subdominant compared to the potential part V(φ), inflation kicks in, whereas when both quantities become comparable, the inflationary period ends. This requirement is typically expressed in the slow-roll conditions ϵ ≡ (1/2)[V′(φ)/V(φ)]² ≪ 1 and η ≡|V″(φ)/V(φ)|≪ 1 for the inflaton to produce inflation and ϵ ≃ 1 to end the inflationary epoch (see also Figure 1 for a sketch of the inflationary potential).

FIGURE 1. Schematic inflationary and reheating phases. Initially, the inflaton slowly rolls along its potential until it reaches a critical point where ϵ ≃1 at φ ≃ φ_end. Subsequently, the inflaton transits rapidly towards the bottom of the potential, where it oscillates rapidly at around the minimum. At such stage the process of reheating takes place.

The scalar perturbations during the inflationary epoch are of special importance since they are attributed the generation of the initial inhomogeneities that gave rise to the large scale structure in the Universe. Inflaton fluctuations may also get to form PBHs in the early Universe. That is, if an initial perturbation is dense enough when it reenters the cosmological horizon, it can collapse under its own gravity to form a black hole.¹ Hence, in order to assess the probability of formation of PBHs in the early Universe, it is necessary to determine precisely the average amplitude of scalar perturbations generated during inflation, that is, the primordial power spectrum.

The evolution of the scalar field is governed by the Klein-Gordon equation. In order to calculate the amplitude distribution or power spectrum of the field fluctuations δφ, the perturbed Klein-Gordon equation must be solved. Adopting the flat gauge, δφ is turned into the Sasaki-Mukhanov variable [see, e.g., Ref. (Malik and Wands, 2009)] and it proves convenient to define a new variable (in Fourier space)

where k is a comoving wavenumber scale. This allows us to rewrite the perturbed Klein-Gordon equation as the so-called Sasaki-Mukhanov equation (Sasaki, 1986; Viatcheslav, 1988):

where a prime (′) denotes a derivative with respect to the conformal time dτ = dt/a, $z\equiv a{\dot{\phi }}_b/H$, and the subscript b is used to refer to background quantities. Note that the comoving curvature perturbation ${\mathcal{R}}_k$ is given in terms of the perturbed scalar field (in flat gauge) δφ_k as

It is in terms of this quantity that the power spectrum of scalar perturbations is defined. Explicitly, the dimensionless primordial power spectrum of curvature perturbations is

The power spectrum is normalized to the amplitudes derived from the CMB at large scales, setting the normalization scale as the pivot scale at k* = 0.05Mpc⁻¹. The usual parametrization for the spectrum follows the evidence that at large scales the spectrum is almost scale-invariant. Thus

with the CMB normalization dictating $\ln (10^{10}{\mathcal{A}}_s)=3.044\pm 0.014$ and the spectral index n_s = 0.9649 ± 0.0042 (Akrami et al., 2020).

Such prescription for the power spectrum would produce a tiny amplitude of fluctuations and produce an extremely small number of PBHs [see, e.g., (Carr et al., 1994; Emami and George, 2018)]. The spectrum, however, may suffer modifications at small scales, where features in the potential may arise and impact the amplitude of fluctuations significantly (see, e.g., the setting in Section 7). It is precisely such possibilities what we aim to explore, and possibly constrain, through the probability of PBH formation in a fertile scenario; the reheating epoch.

2.2 Inflaton evolution after inflation

2.2.1 Reheating

The reheating of the Universe refers to the process through which the energy stored in the inflaton field, responsible for driving the inflationary expansion, is transferred to other particles present in the Universe. This energy transfer takes place at the end of the inflationary period (see Figure 1) and is believed to have created the necessary conditions for the formation of primordial nuclei and structures within the Universe.² Historically, reheating was first treated perturbatively (Abbott et al., 1982; Albrecht et al., 1982) and in this section we review a simple version of this process.

For transferring energy, it is usual to consider a non-minimal coupling of the inflaton with, say, a second scalar field χ through an interaction in the lagrangian. That is,

where g is a dimensionless coupling constant and Σ is a mass term. The decay rate of the inflaton field into χ particles is thus given by (Greene and Kofman, 1999)

where m is the “effective” inflaton mass. The energy loss of the inflaton through its conversion to χ particles can be approximated by the following Klein-Gordon equation (Allahverdi et al., 2010a)

A typical form of the potential used for the reheating epoch is the quadratic potential, V(φ) = m²φ²/2, since, in order to have efficient reheating, it is necessary for the inflaton to oscillate around its global minimum. This precise form of the potential is not necessarily followed during the inflationary epoch and, in fact, for inflation to yield the correct observations given by Planck data (Akrami et al., 2020), a power law of the form V(φ) ∼ φ^α, with α < 1 is required. However, there are many realizations that meet this condition at large values of the field and that converge to the simple quadratic form at small field values. Some examples of such potential are shown in Figure 2. This class of potentials are considered in this work as our main study cases.

FIGURE 2. Different normalized inflationary potentials as a function of φ. The dotted red line corresponds to the typical chaotic potential V(φ)∼ φ². The brown, blue, and black lines correspond to a polynomial (Kallosh and Linde, 2022; Iacconi and David, 2023), T-models (Ellis et al., 2013; Kallosh et al., 2013; Kallosh and Linde, 2013), and generalized α-attractor T-models (Germán, 2021; Francisco, 2023), respectively. In the plot, we used M =1.

To proceed with the analysis, let us note that in the case of an small coupling constant (Γ ≪ H), the interaction term Γ can be neglected and the equation of motion of the inflaton reduced to

This form suggests that in the limit in which m ≫ H (a limit that is fulfilled during the reheating epoch) the field φ_b experiences damped oscillatory motion about φ_b = 0,

Here all quantities are displayed in units of the Planck mass M_Pl. In terms of the scale-factor and averaging over several oscillations we obtain

where the background energy density is ${\rho }_b={\dot{\phi }}_b/2+V({\phi }_b)$, and we have used the subscript _end to refer to quantities evaluated at the end of inflation. The above result shows a pressure-less matter behavior, which in the background can be adopted while the inflaton dominates the overall energy density.

Once the Hubble expansion rate decreases to values comparable to Γ, the χ-particle production becomes efficient, and the energy associated with the inflaton is transferred to the field χ. The temperature at the time at which Γ = H is known as the reheating temperature and is given by $T_R\sim \sqrt{\mathrm{\Gamma }{m}_{\text{Pl}}}$.

As shown in Eq. 8, Γ is proportional to the square of g and, since typically g ≪ 1, we should expect that the reheating temperature occurs at low energy scales (low compared to the energy scale of inflation, at around 10¹⁴ GeV). This temperature can be as low as the scale of Big Bang Nucleosynthesis (BBN, at around 10 MeV), since this is the maximum energy scale at which we have evidence of a radiation-dominated Universe. This allows us to consider a scenario where reheating could have lasted a few e-foldings. This so called slow-reheating process is the scenario we explore in this review. We mainly look at the implications of a slow-reheating brings on the formation of PBHs, the conditions for collapse in this scenario, and the kind of inflationary models that can be constrained with this observable.

To start, it is important to define the stages of the evolution of overdensities of the inflaton field during slow-reheating. In the canonical mechanism, the k-modes associated with the quantum fluctuations of the inflaton field have a fixed amplitude when they stretch beyond the Hubble horizon during the accelerated expansion phase. Once inflation ends, these modes reenter the cosmological horizon after

e-folds of expansion, where the subscript _HC refers to quantities evaluated at the horizon crossing time. Inside the cosmological horizon, two regimes are distinguished, which are separated by the scale k_Q, usually referred to as the quantum Jeans scale or simply the Jeans scale, given by (Suárez and Chavanis, 2015)

Inhomogeneities can be characterised by the density contrast δ = δρ/ρ_b, for which the time evolution (neglecting the decaying mode) is

Here k_H = aH is the scale associated to the size of the cosmological horizon. Thus, density fluctuations with a characteristic scale k > k_Q undergo damping via oscillations, while at scales k_H < k < k_Q fluctuations experience a growth in amplitude proportional to the scale factor³.

2.2.2 Preheating

The preheating instability is a non-perturbative mechanism that arises in theories with non-minimal couplings between the inflaton and other fields, say χ (El Bourakadi, 2021). The dynamics of the inflaton during the preheating process can be described by the Mathieu equation, which is related to periodic or quasi-periodic oscillating systems (Norman, 1964; Kofman et al., 1994; Bassett et al., 2006). Furthermore, its solutions exhibit exponential instabilities, that is, ${\chi }_k\propto \exp ({\mu }_k^{(n)}mt)$ within a series of resonance bands, located at specific frequency ranges Δk⁽ⁿ⁾ (here labeled by the integer index n). Such instabilities lead to an exponential growth in the occupation numbers of quantum fluctuation modes, denoted as $n_{\bar{k}}(t)\propto \exp (2{\mu }_k^{(n)}mt)$, which are interpreted as the production of χ particles (Abbott et al., 1982; Kofman et al., 1997; Bassett et al., 2006; Lyth and Liddle, 2009). In short, preheating describes the process through which the energy density of the created particles, calculated within the above formalism, is extracted from the energy density of the oscillating inflaton field.

3 PBH formation during preheating

The production of PBHs during preheating was first studied in Ref. Anne and Malik (2001). In such work, the authors studied a two-field chaotic inflation model and found that for a wide range of parameters the resonant amplification of modes during preheating leads to an overproduction of PBHs, before backreaction terminates the resonance.

In order to handle the non-perturbative and non-minimal interaction between fields, the Preheating process has been modeled through lattice field theory simulations that evolve the scalar field equations on a homogeneous background (Kofman et al., 1994; Khlebnikov and Tkachev, 1996; Felder and Kofman, 2001; Micha and Igor, 2004; Felder and LATTICEEASY, 2008; Khlebnikov et al., 2012; Jo and Rubio, 2016; Dux et al., 2022), though neglecting the backreaction of inhomogeneities on the local spacetime metric⁴.

With the use of lattice simulations came the development of computational techniques that captured the non-linear aspects of the problem and showed the importance of various factors on PBH formation (Niemeyer and Jedamzik, 1999; Shibata and Sasaki, 1999; Harada et al., 2013; Nakama et al., 2014a; Nakama et al., 2014b). These factors encompassed the equation of state, the nature of the inflationary potential, and the presence of additional fields or interactions⁵. Importantly, even small inflaton self-interactions can accumulate over multiple oscillations, triggering the resonant growth of non-zero momentum inflaton modes, which constantly interact with the homogeneous component. This process of self-fragmentation in the inflaton field not only redistributes the initial energy density but also leads to the formation of localized soliton-like structures referred to as oscillons [see, e.g., (Bogolyubsky and Makhankov, 1976; Linde, 1990; Gleiser, 1994; Kolb and Tkachev, 1994; Copeland et al., 1995; Honda and Choptuik, 2002; Fodor et al., 2006; Amin et al., 2012a; Antusch et al., 2018; Hong et al., 2018)]. The presence of large numbers of oscillons has motivated the search for PBH production and the implications for reheating, placing constraints on various single-field inflation models and other models accommodating oscillon solutions [e.g., (Cotner et al., 2018; Cotner et al., 2019)]. Specifically, Ref. Cotner et al. (2018) shows that the fragmentation of the inflaton into oscillons can give rise to the formation of PBHs in single-field inflation models or other models permitting oscillon solutions. Subsequently, Ref. Suyama et al. (2005) showed that PBH production in preheating does not exceed astrophysical bounds because the mass of PBHs is small enough to evaporate before BBN. As a result, these PBHs are not constrained by observation, even if they are overproduced, unless they leave behind Planck mass relics (Barrow et al., 1992).

It has also been argued that the assumption of a Gaussian probability distribution for density perturbations at horizon crossing is crucial. Since the density perturbations that lead to PBH formation are very rare and sensitive to the tail of the distribution, on average, PBHs are not overproduced during the violent non-equilibrium phase of preheating that follows the end of many inflationary models. In Ref. Anne and Malik (2001) a linear approximation estimated the time when backreaction becomes significant and when the amplitude of density perturbations surpasses a certain threshold, treating them separately. This leads to a criterion for PBH formation. It is important to note that even a small error in determining the backreaction time can lead to incorrect conclusions due to the exponential growth of perturbation amplitude. This discrepancy highlights the sensitivity of the results to the precise timing of backreaction and its potential impact on the predictions for PBH production or overproduction. On the other hand, Ref. Torres-Lomas and Arturo Ureña-LAlpez (2013) modified a version of HLattice to numerically solve the relevant equations of motion and analyze the mass variance as a means to explore the formation of structures during the preheating phase. The study revealed that preheating has the potential to generically produce PBHs. However, the results highlighted the influence of the smoothing scale values and emphasized the need for backreaction, to confirm the obtained results. Subsequently, Ref. Angelo et al. (2022) found strong backreaction effects in the system, invalidating the standard perturbation theory approach. They also observed that the non-Gaussianity of the comoving curvature perturbation is large in the linear regime but gets suppressed as the dynamics become nonlinear. This suppression of non-Gaussianity relaxes the bounds on PBH overproduction, allowing instead for an observable gravitational wave signal at interferometer scales.

4 PBH formation during slow-reheating: direct gravitational collapse

The transition of the Universe to the standard Big Bang cosmology prior to the BBN epoch can be achieved through a variety of mechanisms. One possibility is the fragmentation of the inflaton condensate into its own quanta, triggered by self-resonance (Amin, 2010; Amin et al., 2010; Lozanov and Amin, 2018; Fukunaga et al., 2019). In the same context, as mentioned above, particles coupled to the inflaton can be resonantly produced (Kofman et al., 1994; Kofman et al., 1997) leading to prompt thermalization (Lozanov and Amin, 2017) or a potential oscillon dominated epoch (Amin, 2010; Amin et al., 2010; Amin et al., 2012b; Lozanov and Amin, 2018). During this epoch, the formation of PBHs is possible due to the gravitational collapse of the perturbations that were resonantly amplified.

In contrast, if parametric resonance is not present, the particle production occurs through a more gradual, perturbative processes [as described in, e.g., (Abbott et al., 1982; Albrecht et al., 1982)]. As discussed in Section 2.2, these processes can take place during a relatively long period of expansion when the Universe is governed by a nearly homogeneous condensate, in a (nearly) φ² potential, and when Γ ≪ H. Eventually, the condensate fragments due to the gravitational growth of perturbations (Jedamzik et al., 2010b; Easther et al., 2011), in what is dubbed a primordial structure formation process, where inhomogeneities virialize. When the matter concentration is sufficiently high, the hoop conjecture prescribes that some of these structures may instead collapse and form PBHs. The mass of such black holes should be close to the mass of the cosmological horizon evaluated at the time of horizon crossing, conveniently expressed as:

Here, γ is a constant that encrypts the efficiency of the collapse. The precise value of γ ought to be determined through numerical calculations currently in progress for the reheating scenario [partial progress has been reported in (Padilla et al., 2021a; Eloy de Jong et al., 2022)]. For the sake of the argument, we will adopt the value γ = 1. In this and the following two sections, we review the conditions for which PBH formation may occur in a slow-reheating scenario.

Since the inflaton during reheating behaves like pressureless matter, one may be tempted to extrapolate the perfect fluid criterion lim_w→0δ_th → w (Harada et al., 2013) and deduce a copious formation of black holes, following the Press-Schechter formalism usually employed in the standard caculation, and as exemplified for this case in Appendix A. Such approach is an over-simplification of the problem since, as we have mentioned in Section 5.3–and will discuss in more depth in the following, the inflaton presents an effective pressure due to its quantum nature, which prevents the formation of PBHs from overdensities of infimum amplitude. Moreover, the collapse criteria used for dust-like overdensities can be adopted and complement those derived for a cosmological scalar field. In the following we review the diverse criteria with the aim of assessing PBH formation in a slow-reheating scenario.

4.1 The sphericity criterion

One of the most widely used criteria to describe the formation of PBHs in a slow-reheating scenario was discussed early in the develpment of the theory (Khlopov and Polnarev, 1980; Polnarev et al., 1985). Physically, this criterion limits the configurations of initial pressureless overdensities to be sufficiently close to spherical symmetry, so as to collapse onto a black hole. Inspired by this, Ref. Harada et al. (2016) presents a more detailed analysis of the so-called sphericity criterion for the collapse of overdensities. The latter article investigates the formation of PBHs in the matter-dominated phase of the Universe, where nonspherical effects in gravitational collapse play a crucial role. The authors apply the Zel’dovich approximation (Zel’dovich, 1970), Thorne’s hoop conjecture (Wheeler and Klauder, 1972), and Doroshkevich’s probability distribution (Doroshkevich, 1970) to derive the production probability β₀ of PBHs. In summary, in the limit of a small variance of δ evaluated at the horizon crossing time, the relation

approximates the initial abundance of PBHs as a function of the variance σ² in a dust dominated Universe. Note that this is not directly linked to a threshold amplitude, but instead this prescription alone may overproduce PBHs even for a scale-invariant power spectrum, if reheating lasts long enough (Juan Carlos Hidalgo et al., 2017).

4.2 Conservation of angular momentum criterion

The initial angular momentum of overdensities plays an important role in the formation of PBHs. In Ref. Harada et al. (2017) it was shown that this can lead to a significant suppression of the production rate. In particular, it was found that the limit on the ammount of angular momentum allowed for collapse provides a threashold value for PBH formation which complements the sphericity criterion. Specifically, the production probability β₀ is restricted to

where $\mathcal{I}$ is a parameter of order O (1) and f_q (q_c) is the fraction of mass with a level of quadrupolar asphericity q smaller than a threashold q_c. Comparing this result with (Eq. 17) and assuming $\mathcal{I}=f_q(q_c)=1$ [as assumed in Ref. Harada et al. (2017)] it is evident that the angular momentum criterion of Eq. 18 is more stringent than the sphericity criterion of Eq. 17 for a standard deviation σ ≲ 0.005 (the relevance and matching of these constraints is illustrated in Figure 6).

We conclude this section by mentioning that according to the study conducted in Ref. Eloy de Jong et al. (2023). on the formation of spinning PBHs during an early matter-dominated era, the efficiency of mass transfer was found to be approximately 10%, while the efficiency of angular momentum transfer was estimated to be around 5%. That reference further suggests that unless the matter era is short, the final dimensionless spin of PBHs is expected to be negligible.

4.3 Reheating time criterion

Another employed criterion for characterizing the generation of PBHs arising from an slow-reheating epoch can be found in Ref. Martin et al. (2020) (see also, e.g., (Goncalves, 2000; Carr et al., 2018), for earlier work that considers this criterion of PBH formation). Unlike the previous criteria that consider the morphology and angular momentum of the initial inhomogeneity, this criterion is mainly focused on studying the time required for the collapse of configurations. In this way, one can impose a condition on the contrast density evaluated at the horizon crossing time that a perturbation must meet in order to form a PBH. In this section, we will review the most important results of such work, considering that the reader interested in the details of the calculations presented here will be able to review the original papers.

As previously argued in Section 2.2, modes within the resonance band k_H < k < k_Q are expected to behave as standard pressureless matter fluctuations, which may therefore collapse to form a primordial structure such as a PBH. The time t_c required for such collapse in the spherical collapse model is given by Goncalves (2000):

Here the subindex _bc is related to quantities evaluated at the time when the mode k transitions across the instability band. Note that there are two ways in which a k − mode may enter the instability band. First are those scales which exited the cosmological horizon during inflation and reenter (entering from above) after the inflationary epoch. The second possibility is for those scales that enter the instability band from below, due to the fact that the quantum Jeans scale decreases in size as ∼ a^−1/4. Following the work of Ref. Martin et al. (2020), in this section we shall only concentrate on the scales which enter the instability band from above and thus we equate t_bc = t_HC. Additionally, we can also re-express Eq. 19 in terms of the number of e-foldings. Considering that during a matter-dominated Universe we have H = 2/(3t) and assuming $3\pi /(2{\delta }_{\text{HC}}^{3/2}(k))\gg 1$, which is a good approximation in the perturbative regime, we obtain

We can identify the last mode to enter into the instability band from above as k_Γ = a_ΓH_Γ, where subindex _Γ indicates quantities evaluated at the time the inflaton field decays. We can then relate the scale k_Γ with the scale at the end of inflation, k_end by $k_{\mathrm{\Gamma }}/k_{\text{end}}={({\rho }_{\mathrm{\Gamma }}/{\rho }_{\text{end}})}^{1/6}$. Consequently, the scales that can collapse gravitationally during the slow-reheating epoch to form a PBH are those within the interval

The condition adopted in (Martin et al., 2020) to determine if a scale k can collapse to form a PBH was to require that the time for the collapse of the perturbation be smaller than the total period that the slow-reheating epoch spans. Noting that such a time can be expressed as

and using Eq. 19, we obtain that the condition for a perturbation to collapse is given by

where

and the maximum value in this context is determined by the requirement for the horizon-crossing fluctuation to be within the perturbative regime (Martin et al., 2020), imposing this upper bound should not significantly impact the abundance of PBHs, considering that the amplitudes of density contrasts are exponentially suppressed, preventing an overly abundant production (Niemeyer and Jedamzik, 1998; Harada et al., 2017).

To calculate the abundance of PBHs β₀ formed in this context, we can adopt the Press-Schechter formalism (see Appendix A for details):

where $\mathbb{P}[\delta >{\delta }_{\text{th}}]$ is the probability that a smoothed density field exceeds the threshold value δ_th and is given by

In the above expression R = 1/k, erfc(x) = 1 − erf(x) is the complementary error function, and σ(R) is, as previously defined, the variance of δ evaluated at the horizon crossing time.

It is important to note that this formation criterion assumes that all perturbations with a collapse time shorter than the remaining time for reheating will gravitationally collapse to form PBHs. Consequently, this criterion is likely to overestimate the actual abundance of PBHs formed, β₀. This is due to the omission of significant physical effects discussed throughout this review. In the following we focus on the physics of fluctuations from an oscillating scalar field, while in Section 6 we address the PBH formation criteria in such environment.

5 Dynamics and structure formation in a slow-reheating epoch

In the preceding section, we described the process of PBH formation during a period of slow-reheating, assuming a dust-like behavior for the governing inflaton field. However, such approximation oversimplifies the dynamics of a cosmological scalar field, as it neglects its inherent quantum nature. In order to describe more accurately the evolution of the post-inflationary epoch one must solve the Einstein-Klein-Gordon (EKG) system of equations in a cosmological background. However, such task presents the practical difficulties described below.

The post-inflationary regime for an (almost) free field is subject to the condition m ≳ H_end. As indicated by Eq. 11, the oscillation frequency of the homogeneous inflaton field is precisely m. On the other hand, the Universe evolution is characterized by the Hubble time 1/H ∼ a^3/2. Thus, a few e-folds after the end of inflation, the condition turns into 1/H ≫ 1/m. This thus stipulates two dissimilar characteristic time scales in the numerical evolution of the complete EKG equations, which turn the evolution over several Hubble times computationally unfeasible. Efforts in this direction are found in Refs. Alcubierre et al. (2015); Rekier et al. (2016); Eloy de Jong et al. (2022).

Here we follow a different approach, acknowledging the similarities between a slow-reheating epoch and the phase of structure formation in the scalar field dark matter (SFDM) model. Such parallelism was initially proposed in (Nathan et al., 2020) and has extended in subsequent studies (Juan Carlos Hidalgo et al., 2017; Jens, 2020; Carrion et al., 2021; Eggemeier et al., 2021; De Luca et al., 2022; Eggemeier et al., 2022; Padilla et al., 2022; Chavanis, 2023; Eggemeier et al., 2023; Hidalgo et al., 2023; Padilla et al., 2023). Taking advantage of the extensive literature on SFDM models, we aim for a better understanding of the Universe right after inflation. In the upcoming sections, we will review various SFDM results that can be adapted to the slow-reheating scenario. With such results at hand, in the following section, we formulate the conditions under which the primordial structures of this period may undergo gravitational collapse, resulting in the formation of PBHs.

5.1 The Schrodinger-Poisson picture

Let us look at a few approximations to simplify the treatment of the post-inflationary Universe. For instance, when a particular scale becomes non-linear, it is expected to be well within the Hubbet horizon and with a bulk motion exhibiting a nonrelativistic behavior, with large occupation numbers. These specific characteristics in the system are the precise requirements to apply the Newtonian approximation, where the EKG system can be reduced to the Schrödinger-Poisson (SP) system of equations. In such approximation, matter is represented by the non-relativistic wavefunction ψ, and the (Newtonian) gravitational potential Ψ is determined by solving the Poisson equation. Below we quickly review how the Newtonian description of the system is reached.

In the slow-reheating era, it is possible to describe gravity through the weak-field approximation. Well within the cosmological horizon, we adopt a spatially flat background metric and deal with scalar perturbations in the Newtonian gauge (see, e.g., Ref. Malik and Wands (2009)),

Considering that the anisotropic stress of a minimally coupled scalar field vanishes, we can identify the Newtonian potential as Ψ = −Φ. This allow us to write the Einstein-Hilbert action for subhorizon scales (k ≫ aH) and by considering quantities at first order in the potential and at second order in spatial derivatives as (Niemeyer, 2019)

Further simplification is achieved through the following considerations; while φ oscillates at a frequency m, the density field changes slowly within the nonrelativistic regime. To account for the rapid oscillations, we can introduce the complex field ψ, as follows:

With such definition, disregarding oscillatory terms that involve powers of exp (±imt/ℏ), and subsequently incorporating the simplifying assumptions above justified, namely, $\dot{\psi }\ll m\psi$ and m ≫ H, Eq. 28 is simplified to

where we write explicitly the ℏ factors to restore units and emphasize the quantum nature of the system here described.

After varying the above action S with respect to the gravitational potential Ψ and the field ψ, we finally arrive at the SP system of equations:

Here ⟨ρ⟩ is the smooth background value of the density of the scalar field, ⟨ρ⟩ = m⟨ψψ*⟩.

5.2 Quantum hydrodynamics equations

One of the great advances in the study of the SP system was the realisation that this pair of equations can be reformulated like classical hydrodynamics. The hydrodynamic version of the SP system introduces an additional quantity Q, known as the “quantum potential” or, in some instances, the “Bohm potential” (Bohm, 1952a; Bohm, 1952b). As a result, these equations are commonly named as the “quantum hydrodynamics” (QHD) equations (we recommend (Wyatt, 2005) for a textbook presentation of these equations). It is important to emphasize that both the hydrodynamic and the field expressions of the SP system, are equivalent and offer robust methods for addressing the nonlinear dynamics of a cosmological scalar field. Let us outline in this section two methods for deriving these QHD equations.

5.2.1 The Madelung–Bohm formulation of quantum hydrodynamics

We can consider a Madelung transformation (Erwin, 1927) of the form

If we define the bulk flow velocity of the field v as v = ∇θ, we can reexpress the SP system of equations as

where Q is given by

This set of equations constitute the QHD system. The first of these equations, Eq. 33a, corresponds to a continuity equation, typically found in classical fluid dynamics. Such an equation is used to describe the conservation of mass within the system. Additionally, the second equation, Eq. 33b is an Euler-like equation stating momentum conservation. However, instead of the conventional terms associated with a fluid pressure gradient, in these QHD equations, we encounter a novel potential Q that encapsulates the quantum properties of the scalar field.

5.2.2 Phase space formulation

We can also obtain the QHD equations by taking momentum moments of the SP system of equations (Takabayasi, 1954). Following Ref. Taha et al. (2021), we sketch the steps in the following.

We define the Wigner function (Wigner, 1932) as

From this expression, the number density at a point in coordinate space is determined by the integral of W over momentum space,

Moreover, the local average of a quantity A over momentum space is computed through

As an example, if we use the bulk velocity v = p/m and we perform the above integration, we arrive at the same bulk velocity in terms of ∇θ as before. One can also calculate the velocity dispersion tensor as follows

To derive the equation of motion for the Wigner function, we compute the partial time derivative (∂W/∂t) and incorporate it into the SP system of equations. The outcome is known as the Wigner-Moyal equation, which bears resemblance to the collisionless Boltzmann equation (CBE), also known as the Vlasov equation. For this particular case, the Wigner-Moyal equation introduces additional terms that encode the quantum characteristics of the scalar field.

When computing momentum moments of the Wigner-Moyal equation, it is possible to derive, from the 0^th moment the continuity equation, Eq.33a and from the first moment an Euler-like equation from classical fluid mechanics. Once again, in this formalism a new pressure-like term emerges, referred to as the “quantum pressure” tensor, Π_ij, which is linked to the velocity dispersion tensor as

It is worth mentioning that, in general, there is not a real difference between the first moment equation and Eq. 33b, since both equations coincide once the quantum potential Q is defined as

which reduces to the explicit form of Eq. 33c.

Note that Eqs. 38, 39 stipulate that the force originating from the quantum potential term in the momentum equation corresponds to the effective “pressure” present in the momentum flux density associated with the internal distribution of momentum in the phase space derivation. In both situations, these “quantum” terms arise from the kinetic term of the SP equations and represent the wave-like behavior of the scalar field. In practice this represents a force opposing gravitational collapse.

5.3 Some important scales

As shown above, the hydrodynamic formulation of the SP equations is a system that closely resembles the description of the nonlinear dynamics of a pressure-less fluid, with the exception of an additional term that accounts for the quantum properties of the scalar field. In this formulation it is easy to see that the quantum potential term becomes important only when its contribution is comparable to the kinetic and gravitational potential. This happens for scales that fulfill R ∼ λ_dB, where R (λ_dB ≡ ℏ/(mv)) is the characteristic scale (de Broglie wavelength) of the configuration [see, for example, (Niemeyer, 2019)]. At scales R ≫ λ_dB the expected behavior of the scalar field configurations should be very similar to that of a dust-like component, while for scales with R ∼ λ_dB the quantum contribution must be considered.

In order to estimate the timescales at which the quantum characteristics of the scalar field become significant, one may look at the gravitational scattering time for wave scattering within a condensate⁶. In the absence of external influences, the scattering rate Γ_s, inversely proportional to the time interval τ, is dependent on the scattering cross section σ_g, the average relative velocity $⟨v⟩=\sqrt{2}v$, and the number density n = ρ/m. Namely, Γ_s ∝ σ_g⟨v⟩n. However, when the final state experiences macroscopic occupation, the rate is further enhanced by the scalar field phase space density, often referred to as the occupation number $\mathcal{N}$. Such enhancement is due to the phenomenon of Bose-Einstein stimulation [see (Niemeyer, 2019)] with

Correspondingly, the scattering time is given by Levkov et al. (2018)

where the momentum-transfer cross section σ_g for Rutherford scattering is given by σ_g ≃ πG²m²v⁻⁴ log  Λ, with Λ ≃ R/λ_dB.

This implies that over period of order O(τ) we may expect effect due to the quantum nature of the scalar field. One of such effect, which we shall discuss in more detail later, is the formation of solitonic structures through the Bose-Einstein condensation.

5.4 Soliton solutions

If we set the scale factor a = 1 and assuming ρ ≫⟨ρ⟩⁷, the SP system of Eq. 31 admits solutions of the form

where r is the radial coordinate and E is the energy associated to the configuration. The system described by the SP Eq. 31 and the above ansatz has numerous solutions satisfying appropriate initial and boundary conditions (Siddhartha Guzman and Arturo Ureña-Lopez, 2004; Siddhartha Guzman and Arturo Ureña-Lopez, 2006). These solutions, often referred to as Newtonian boson stars (NBS), are characterized by the number of nodes present in ψ before the solution asymptotically decays. The solution without nodes, the soliton solution, is the ground state of the SP system, presenting the lowest energy. Accordingly, solutions with nodes are called the excited NBSs.

The soliton solution is the most widely studied in the literature [see, for example, (Ruffini and Bonazzola, 1969; Seidel and Suen, 1991; Seidel and Suen, 1994a; Alcubierre et al., 2002; Siddhartha Guzmán and Arturo Ureña López, 2004; Siddhartha Guzmán and Arturo Ureña-López, 2006; Chavanis, 2011; Alvarez-Ríos and Francisco, 2022)]. This is because the soliton is the attractor solution of the SP system: scalar field configurations with arbitrary initial conditions tend to migrate through a “gravitational cooling” mechanism to the ground state solution of the SP system [see (Seidel and Suen, 1994b)]. In addition, it has been also shown in Levkov et al. (2018) that initially homogeneus scalar fields with Gaussian-distributed initial conditions in momenta evolve to form localized soliton profiles by Bose-Einstein condensation in a timescale t ∼ τ. After its formation, a soliton accretes mass according to

Note that the condensation time τ dictates the timescale at which the soliton structures form and evolve.

Once a soliton structure is virialized, its properties are mostly related to its mass M_sol. For example, the half-mass radius R_1/2 and the virial velocity v_vir are given by [see, for example, appendix B in Ref. Lam et al. (2017):

Additionally, the coherent length λ_dB for the virial velocity is:

which implies that the solitons formed from an scalar field present sizes of the order of the de Broglie wavelength associated to the configuration.

From the QHD system, Eq. 33, the physical nature of the soliton profile can be clearly understood. Since soliton solution is a static configuration, ∂_tv = 0 = v. If we set for simplicity $\nabla \sim 1/{\mathcal{R}}_{\text{sol}}$, where ${\mathcal{R}}_{\text{sol}}$ is a characteristic radius of the soliton profile (typically λ_dB), we obtain,

And the condition of equilibrium configuration yields,

In the above expression ν is a dimensionless constant.

From the above expression we can see then that the soliton profile can be understood as a result of the equilibrium between the forces generated by the quantum and the gravitational potentials. Using Eq. 47a it is also evident that the soliton profile must fulfill the condition

which maintains the same parameter dependence found in the exact numerical treatment, Eq. 44.

A general-relativistic regime.- We can anticipate from the relation $M_{\text{sol}}\propto R_{1/2}^{-1}$ (from Eq. 44) that for certain soliton masses we should expect that a general relativistic treatment must be necessary to describe the configurations adequately. In fact, when general relativistic effects are incorporated to the system a different mass-radius relation for the soliton profile is obtained for small radius (large masses)⁸ (see Figure 3). In particular, a limiting maximum mass is predicted to exist for the soliton profile. Such critical mass has been widely studied in the literature (Ruffini and Bonazzola, 1969; Seidel and Suen, 1991; Seidel and Suen, 1994a; Alcubierre et al., 2002) and is given by:

Above this critical mass, no stable soliton solutions are expected to exist. This is because, for soliton structures with larger masses, we anticipate that the force generated by the quantum potential is insufficient to balance the gravitational force resulting from the self-gravity of the system. This phenomenon is equivalent to the Chandrasekhar mass limit for white dwarf stars but associated to soliton structures.

FIGURE 3. Mass-radius relation for the soliton profile obtained from the Newtonian approximation, Eq. 44, and the General Relativistic treatment.

In the case of configurations that include excited states of a scalar field with a specific mass, it has been demonstrated that the resulting configurations can possess larger masses (Seidel and Suen, 1990; Hawley and Choptuik, 2003; Ureña-López, 2009; Bernal et al., 2010; Ureña-López and Bernal, 2010). However, as previously discussed, these excited states have been found to be unstable and undergo gravitational cooling, ultimately transitioning to the ground state solution.

5.5 The Schrödinger-Vlasov correspondence

As mentioned earlier, we can use either the SP or the QHD equations to explain the process of structure formation during the slow-reheating epoch. However, this approach necessitates addressing the characteristic length scale λ_dB, which is significantly much smaller than the cosmological horizon in the post-inflationary Universe. Consequently, modeling the gravitational collapse of a structure using either of the previously discussed formulations becomes a challenging task in general.

Nonetheless, there exists an alternative approach, although it is approximate in nature, which allows us to simplify the resolution of the λ_dB scale and capture the significant effects arising from the wave mechanics of the field on large scales (much larger than λ_dB). This method involves smoothing out the intricate details of the dynamics occurring at small scales (less than or approximately equal to λ_dB) governed by the SP equations. We shall elaborate more on the details of this description in this section.

Based on the studies conducted in (Skodje et al., 1989; Lawrence, 1993; Kopp et al., 2017; Mocz et al., 2018; Taha et al., 2021) we choose to utilize the Husimi representation of ψ, which is a smoothed phase space representation. This representation, introduced in (Kôdi, 1940), involves smoothing out ψ using a Gaussian window with a width parameter η and subsequently performing a Fourier transform of the form

The Husimi distribution function, which defines an smoothed mass density structure of phase space, is defined as

The methodology for handling this distribution closely mirrors the techniques used in conventional Wigner function analysis. However, in this context, we replace the Wigner function with the newly introduced distribution function $\mathcal{F}$. Consequently, we can apply the same methods as before to derive local number density, bulk velocity, and velocity dispersion.

Furthermore, we can calculate the equation of motion for $\mathcal{F}$ in a manner akin to the Wigner-Moyal equation. This involves computing $\partial \mathcal{F}/\partial t$ and then replacing it into the SP equation. The resulting equation, when smoothed over scales significantly larger than λ_dB (η ≫ λ_dB), simplifies to the CBE or Vlasov equation:

This last equation is the same equation that is typically used to describe the structure formation process for a dust-like component, as it is the case, for example, of the CDM model for dark matter. The only difference is that in the case of dust, we would replace $\mathcal{F}$ with the phase space distribution function of its collisionless N-body particles. Other than that, both formalisms are entirely equivalent, implying that at scales larger than λ_dB, the dynamics of dust and the dynamics of a scalar field must be entirely identical.

5.5.1 The fluid approximation

Similarly than in Section 5.2.2, it is also possible to find fluid equations from this description by computing momentum moments of the CBE. Such procedure is well explained in (Taha et al., 2021), which the reader should consult for a more comprehensive discussion. However, in particular the 0^th and first momentum of Eq. 52 are found to mimic the continuity and momentum equations of classical hydrodynamics

The quantity P_ij is defined as $P_{ij}\equiv \rho {\sigma }_{ij}^2$, where ${\sigma }_{ij}^2$ represents the phase space velocity dispersion. Remarkably, this P_ij serves as an effective “pressure” term, analogous to the quantum pressure tensor Π_ij found in the exact QHD equations.

Comparing the results from Section 5.2.2 with the system in Eq. 53, we note that the process of smoothing over scales much greater than λ_dB reduces the “quantum pressure” tensor Π_ij to an effective velocity dispersion tensor P_ij derived from the CBE. Thus, this tensor accounts for the effects of the quantum potential/pressure on large scales. This velocity dispersion plays a crucial role in the process of structure formation, as it is essential for the stability of galactic systems, opposing the gravitational collapse and maintaining dynamical equilibrium in the system. Clearly, a small velocity dispersion can result in collapse, while excessive dispersion may cause the system to disintegrate. In summary, the effective pressure resulting from the velocity dispersion, generated by the quantum potential of the scalar field, may act to prevent PBH formation at scales above the characteristic de Broglie wavelength. Such effects were investigated in detail in Ref. Padilla et al. (2022). [see also (Harada et al., 2023)]. In a subsequent section, we describe the necessary conditions under which the collapse of inflaton perturbations can indeed take place and form PBHs.

5.6 Soliton cores and its halo-like exterior

In the context of dark matter, the first realistic simulations of a scalar field considering cosmological initial conditions were conducted by Refs. Schive et al. (2014a); Schive et al. (2014b) In that work, the authors solved the SP system of equations, Eq. 31, considering a scalar field as the only constituent of the Universe. One of their key findings was that the final structures formed from this cosmological scalar field can be well-described by an inner soliton profile (described by the theory we reviewed in Section 5.4 surrounded by an NFW-like envelope from a radius determined by the incoherent fluctuations of the scalar field⁹. A schematic plot of such configurations is shown in Figure 4. This soliton-envelope structure has been confirmed to exist by several studies that considered simpler scenarios (Jan and Niemeyer, 2016; Schwabe et al., 2016; Alvarez-Ríos et al., 2023), e.g., in the nearly simultaneous merger of several soliton configurations (Mocz et al., 2017). In the context of reheating, such core-envelope structures have also been reported (Eggemeier et al., 2022).

FIGURE 4. In dotted blue is plotted a halo-like structure that is formed in cosmological simulations whereas in dashed black is plotted the soliton profile, predicted by the theory reviewed in Section 5.4.

Numerical simulations prescribe a soliton-halo¹⁰ mass relation given by $M_{\text{sol}}\propto M_{\text{halo}}^{1/3}/m$. Several hypotheses have been proposed to justify this relation. For instance, it has been suggested that the specific energy of the central soliton and the host halo are of the same magnitude (Bar et al., 2018), while others suggest that their specific kinetic energy is what is comparable (Bar et al., 2019). There is also the suggestion that the velocities of the soliton and the host halo, such as the circular, virial, or dispersion velocity, are equivalent (Mocz et al., 2017; Chavanis, 2019). Regardless of the interpretation, the above studies agree in the following relation:

In the above expression $M_{\text{halo}}=(4\pi /3){\rho }_{200}(a_{\text{NL}})R_{\text{halo}}^3$ is the total mass of the halo structure, which should coincide with the mass of the cosmological horizon evaluated at the horizon crossing time (see Eq. 16),

R_halo is the virial halo radius, ρ₂₀₀ (a_NL) = 200ρ_b (a_NL), and subindex _NL refers to quantities evaluated at the time the k-mode becomes non-linear¹¹. We can re-express Eq. 54 more conveniently as follows:

Or equivalently, with some more algebra,

In the above expressions m₅ ≡ m/(10^–5 m_Pl) and ${\rho }_{11}(a)\equiv {\rho }_{200}(a)/{(10^{11}\mathrm{G}\mathrm{e}\mathrm{V})}^4$. This is an important result since, for a given m₅, we anticipate that the mass of the solitons formed in the post-inflationary Universe (during reheating) must be closely related to the amplitude of the perturbations δ_HC(k) and independent of the mass of its host halo. For example, in the case where m₅ = 1 and δ_HC ∼ 10^–5 (typically predicted by CMB observations) we have that the solitons formed in this context should have a mass of M_sol(k) ≃ 6.54 × 10⁻⁴ g.

6 PBH formation during slow-reheating: Collapse from primordial structures

6.1 New criteria of PBH formation in a slow-reheating scenario

As we previously shown, if reheating lasts long enough, two types of structures could form during this phase. Firstly, the formation of halo-like structures that resemble a typical NFW profile when smoothing over scales larger than λ_dB. On the other hand, when considering scales R ∼ λ_dB, the formation of soliton-like structures is expected. In this section, we will investigate the conditions under which we can expect the gravitational collapse of both types of structures onto PBHs, closely following Ref. Padilla et al. (2022).

6.1.1 Halo collapse

If the halo-like structures that arise following inflation are massive enough, they could become gravitationally unstable and collapse, resulting in the formation of PBHs. Specifically, a reliable indicator of such collapse is when overdensities reach a state of virialization within a radius for the associated halo comparable to, or smaller than, the Schwarzschild radius; that is, R_Sch ≡ 2GM_halo. In other words, if R_halo ≤ R_Sch, the collapse to PBHs is likely to occur. By comparing the halo’s virial radius with the Schwarzschild radius, we can determine that PBH formation will always take place when the following condition is satisfied:

When this inequality is combined with Eq. 12 and the mass of the halo (which coincides with the mass of the cosmological horizon at the horizon crossing time), it reduces to the following condition:

By using the relationship $\mathrm{\Delta }{N}_{\text{NL}}(k)=\ln [1.39{\delta }_{\text{HC}}^{-1}(k))]$ (see footnote 11), we can finally obtain the condition for PBH formation:

6.1.2 Soliton collapse

In accordance with our previous discussion in Section 5.4, there is a maximum possible mass for soliton structures, as described by Eq. 49. This maximum mass is reached when the quantum pressure resulting from the Heisenberg uncertainty principle is insufficient to counterbalance the self-gravitational forces within the soliton structure. By substituting this maximum mass into Eq. 56, we can determine a critical threashold value ${\delta }_{\text{th}}^{(\mathrm{s}\mathrm{o}\mathrm{l}\mathrm{i}\mathrm{t}\mathrm{o}\mathrm{n})}$ beyond which the central soliton becomes unstable and may collapse to form a PBH. The condition for PBH formation in this scenario is given by¹²:

6.2 Overview of the mechanisms of PBH formation in slow-reheating

Let us describe the integral picture of the variety of mechanisms leading to the formation of PBHs during reheating. Figure 5 summarizes this general timeline, which we will discuss in more detail below.

FIGURE 5. Evolutionary stages of density fluctuations during the reheating period as a function of the number of e-folds after inflation concludes (N_end). After the onset of reheating, a mode with a specific wavenumber k reenters the horizon N_HC(k) e-folds later, reaching a nonlinear amplitude at N_NL(k). This nonlinear state gives rise to the formation of a halo-like structure (or a PBH) after a Hubble time, at N_halo(k). Once the halo undergoes virialization, the condensation of the inflaton at its core generates a soliton-like structure (or a PBH), emerging at N_soliton(k). Reheating is expected to reach completion and achieve thermalization at N_reh. However, if this final event occurs earlier, the sequential progression is disrupted, and only some or none of the k-dependent processes may occur.

As previously stated, during the slow-reheating phase we expect a few perturbation modes (the ones at the smallest scales in the spectrum) to reenter the horizon. The number of e-folds at horizon reentry N_HC(k), after inflation ends, can be calculated using the following equation:

Once perturbations reenter the horizon, they grow as δ ∼ a (see Eq. 15) and may reach a nonlinear stage. The number of e-folds necessary to reach this regime depends on the wave number k and the density contrast amplitude at horizon crossing δ_HC(k), and it can be expressed as (see footnote 11):

When inhomogeneities reach a nonlinear amplitude, halo-like structures or PBHs are expected to form within a Hubble time. Expressed in terms of e-folds, this happens at:

where $t_{\text{NL}}(k)=[2/(3H_{\text{end}})]{[e^{N_{\text{HC}}(k)}1.39/{\delta }_{\text{HC}}(k)]}^{3/2}$. The collapse time (or the number of e-folds up to collapse) approximately corresponds to the time derived in Eq. 20, indicating that this collapse criterion is necessary but not sufficient for PBH formation.

During the collapse, three important effects may hinder the gravitational collapse into PBHs. To wit, the sphericity of the configurations, the conservation of angular momentum, and the velocity dispersion. The first two criteria limit the abundance of PBHs as a function of the variance of fluctuations (as discussed in Section 4). The latter effect can be expressed in terms of a threshold amplitude for collapse, as shown in the previous section. With the aid of the Press-Schechter integral, we can express the abundance at the time of formation, β₀, in terms of the variance by assuming a Gaussian probability density. We thus bring all these effects together in Figure 6, which shows limits to the abundance at the time of formation, β₀, that should be carefully interpreted. In the consideration of direct collapse, the velocity dispersion criterion (the halo collapse/black curve), seems to impose the most stringent bound to the production of PBHs. However, if configurations do not virialize, one could follow the evolution of fluctuations as that of dust, which subjects the abundance of collapsed objects to the sphericity (blue line) criterion at larger variance values. Note that the distribution of spin in the initial fluctuations yields the limit imposed by the green line, which is never the most stringent bound to the production of PBHs.

FIGURE 6. The PBH density fraction β₀ plotted as a function of the variance of the density contrast σ. The blue and green lines are plotted by following Eqs. 17, 18, respectively. The solid black and dashed cian lines were plotted by using the standard Press-Schechter formalism, Eq. 25, and the threashold values (Eq. 59) and (Eq. 60), respectively. For comparison, we also included the PBH density fraction β₀ calculated in a standard radiation-dominated universe, where we took ${\delta }_{\text{th}}^{(\mathrm{r}\mathrm{a}\mathrm{d})}=0.41$ [see, for example, (Harada et al., 2013)]. In the plot we have also included the brown and purple circles, which denotes the points at which the sphericity criterion becomes more restrictive for PBH formation (for larger σ) than the criterion related to velocity dispersion or the case in which the PBH formation happens in a radiation-dominated universe, respectively.

The above is, however, not the full story. If reheating extends sufficiently to reach t_soliton(k) = t_NL(k) + τ(k), a soliton-like structure gets to form at the core of virialized haloes. Then, a new type of PBHs may emerge at the center of the halo-like structures via the collapse of the Bose-Einstein condensate. The condensation process initiates when the inhomogeneity becomes nonlinear, and it can be described by the following equation:

In the above expression, we use the condensation time τ from Eq. 41 and reexpress it in terms of halo quantities. The condensation time thus determines the number of e-folds necessary for soliton/PBH structures to form:

The criterion for discriminating PBHs from solitons is given by the threshold value of the overdensity at horizon crossing, presented in Eq. 60. The associated abundance as a function of the variance is presented in Figure 6 by the dashed line. This is an alternative route to the direct collapse and the abundance is thus not subject to the sphericity or spin criteria. Let us emphasize that if reheating is terminated early enough, the soliton collapse will not take place.

In concluding this section, we highlight a final possibility, not been previously mentioned in this paper. In the process of formation of structures, such as halos or solitons, with a mass below their critical collapse threshold, the process of accreting matter from their surroundings can be achieved if reheating lasts long enough. Consequently, the structures can grow in mass until they reach the critical value required for them to collapse and the formation of black holes. This particular possibility has been explored in Ref. De Luca et al. (2022). Illustrating this mechanism in a plot equivalent to Figure 6 is a task to be tackled in future work.

7 Testing the mechanisms in a simple setting

Let us apply the presented hypothesis to a relevant inflationary scenario. We consider a power spectrum parameterized by the following expression

where k* = 0.05 Mpc⁻¹ is a pivot scale. We thus approximate the power spectrum with the usual slow-roll approximation plus a Gaussian peak located at k_p and with a variance ${\mathrm{\Sigma }}_p^2$. For this example we shall consider the model parameters ${\mathcal{A}}_s=2.099\times 10^{-9}$, n_s = 0.9634, ${\mathcal{B}}_s=0.084$, Σ_p = 0.03k_end, k_p = 0.6k_end, and k_end = 0.346 m⁻¹. The power spectrum generated for this set of parameters is sketched in Figure 7.

FIGURE 7. Power spectrum as a function of the scale wavenumber k. The red, blue, and cyan lines represent the minimum k-modes relevant to the reheating epoch, the formation of halo-like structures, and the formation of solitonic structures. It is worth noting that this figure bears a striking resemblance to Figure 2 in (Hidalgo et al., 2023), with the sole distinction being that it considers a reheating period lasting for N_reh =30 e-folds of expansion, as opposed to the 40 e-folds used in (Hidalgo et al., 2023).

From the bound in the tensor-to-scalar ration (r ≤ 0.032) (Tristram et al., 2022), we can impose a constraint on the Hubble parameter evaluated at the horizon crossing time of the pivot scale k*:

Assuming that H* ≳ H_end and using the Friedmann equation, we limit the energy scale at witch the end of inflation took place:

For definiteness we adopt the maximum value allowed, ${\rho }_{\text{end}}={(1.368\times 10^{16}\mathrm{G}\mathrm{e}\mathrm{V})}^4$.

We turn our attention to the value of the density contrast δ_HC(k) evaluated at the horizon crossing time. If we consider the mean amplitude of density perturbations as that of the configurations to collapse, taken as,

we can compute the number of e-folds necessary for each particular scale to reenter the cosmological horizon, Eq. 61, then form halo-like structures, Eq. 63, and finally form soliton-like structures, Eq. 65.¹³

Having computed that, in Figure 7 we mark the minimum k’s (largest scales) that undergo each of these three processes. Let us stress that we have assumed a reheating period that lasted for 30 e-folds of expansion.

For each one of the scales that manage to form halos or solitons, we compute the probability of PBH formation. For this purpose, we use the Press-Schechter formalism [see Appendix A], where we calculate the variance of the contrast density following Eq. (A3) and use the threshold values given by Eqs. 59, 60. The result obtained is shown in Figure 8. For comparison, we also included the abundance of PBH predicted by the sphericity criterion, Eq. 17, and the conservation of angular momentum criterion, Eq. 18, that would modulate the PBH production in the case of pure dust. As expected, the velocity dispersion criterion reduces the abundance of PBHs with respect to the other criteria.

FIGURE 8. The density fraction β₀ as a function of M_PBH for the collapse of solitons (cian) and halo-like structures (black). We have also included the sphericity criterion (blue) and the conservation of angular momentum criterion (green). For comparison, we have included the abundance of PBHs formed in a radiation-dominated universe. We also showed the current constraints on the abundance of PBHs for each of the scenarios studied and the mass range of PBHs formed. In particular, the dashed cian, black, and red lines are constraints that apply for PBHs that form due to the collapse of solitons, halos, and perturbations during a radiation-dominated epoch, respectively (see main text for further discussion).

On the other hand, we can observe that due to the small value of the threshold for collapse of soliton structures, in this example we would have a much larger abundance of PBHs due to the collapse of the central soliton of the halos. It is worth mentioning that in the case of PBHs formed by solitons we would not have a one to one relationship of scale vs. mass of the PBH. This is due to how the soliton’s mass depends on the total mass of its host halo (see Eq. 56). For a more detailed discussion of this, we refer the reader to Ref. Hidalgo et al. (2023).

It is evident that the given power spectrum may result in the formation of PBHs with different mass ranges, depending on the duration of reheating and the dominant matter content. Specifically, our analysis shows that PBHs with masses around 0.74, 202, or 17 g are more prominently produced when they originate from the collapse of scalar field solitons, scalar field halos, or overdensities in a radiation-dominated Universe, respectively. Each of these PBH populations influence specific phenomena at astronomical and cosmological levels, thus meeting different constraints depending on their mass.

In Figure 8, we have included the respective constraints applicable to each scenario. To derive these constraints, we used the publicly available software PBHBeta (Gomez-Aguilar et al., 2023), which enabled us to determine constraints for PBHs at various mass ranges in extended reheating scenarios. These constraints were adjusted for two cases: one where the reheating period lasts for N_reh = 30 e-folds of expansion (to calculate the constraints in the case of collapsing halo- and soliton-like structures) and another with N_reh = 0 to obtain constraints for the collapse during the radiation-dominated scenario. Interestingly enough, our analysis shows that in the case of PBH formation during a radiation epoch, the abundance of PBHs formed would be well below the limits imposed by observations. On the other hand, for PBHs formed from the collapse of Inflaton halos, the abundance is much closer to the limit value. Finally, for the same primordial power spectrum, we would observe an overabundance of PBHs resulting from the collapse of solitonic structures, surpassing the limits imposed by observations.

To conclude this section, it is necessary to mention that the abundances obtained in this example would indicate that models with features similar to that in Eq. 66 would be strongly ruled out in a scenario of extended reheating where reheating lasts long enough. However, in a short reheating period, where PBHs from the peak in the primordial power spectrum do not have time to form, or where the peak in the primordial power spectrum becomes smaller, the simple model (Eq. 66) still satisfies the constraints imposed by the Planck mass relics and cannot be ruled out. Of course, to draw more realistic conclusions for different inflationary models existing in the literature, it would be necessary to compare such models with the criteria for PBH formation discussed here [as was done, for example, in Ref.Padilla et al. (2023)].

8 Summary and outlook

In this article, we have reviewed the criteria for the formation of Primordial Black Holes (PBHs) in the context of a slow-reheating scenario. We focus on an extended reheating stage dominated by an oscillating field in a quadratic potential. Specifically, we have examined how the gravitational collapse of primordial inhomogeneities re-entering the cosmological horizon can be influenced by three significant effects: the effects related to the morphology of the initial perturbation (how spherical or non-spherical it is), the possible angular momentum it might have presented, and the effects due to velocity dispersion. We have described in detail how the latter emerges from the quantum nature of the dominating scalar field once it is averaged over scales much larger than the associated de Broglie scale. Each of these effects prescribes a bound to the abundance of PBHs, as illustrated in Figure 6. In particular, we have found in the case in which σ ≤ 0.04 that the criterion for PBH formation due to velocity dispersion effects is more restrictive than the criteria considering the morphology of the perturbation and its angular momentum. As a result, this criterion is the most important to consider for those values of σ. On the other hand, in the case in which σ > 0.04, the most important effect preventing the collapse of perturbations into PBHs is the one related to the morphology of the system.

At scales comparable to the de Broglie wavelength of the scalar field, we have seen that the formation of a solitonic-like structure (due to the Bose-Einstein condensation process) is expected at the center of virialized configurations that form in the post-inflationary Universe. In this article, we have also reviewed the necessary condition for the formation of PBHs due to the gravitational collapse of these solitons and we have calculated the abundance of PBHs that should be expected by this mechanism (which can be also seen in Figure 6). In the figure it is easy to see that the collapse of solitons into PBHs is much more likely to occur than in the case of the collapse of the total perturbation.

We emphasize that in order to achieve the formation of PBHs due to either of these two formation mechanisms (collapse of the total perturbation or solitonic center), reheating should last sufficiently long to allow for each of these processes to take place. The timeline for this process is presented in Figure 5. Another important remark is that we have assumed the dominating scalar field during reheating to be the inflaton field itself, oscillating around the bottom of a quadratic potential. The extension to scalar fields of other nature, such as axion monodromy (Silverstein and Westphal, 2008; McAllister et al., 2010), curvaton (Enqvist and Martin, 2002; Lyth and Wands, 2002) or multiple fields (Iacconi et al., 2022; Iacconi and Mulryne, 2023), is still work in progress.

This work has been carried out as a review article aimed at guiding the reader through the different criteria for PBH formation in the context of a slow-reheating scenario. Our study considers exclusively the case of a scalar field (inflaton or another) oscillating around the minimum of a quadratic potential and dominating the Universe prior to the standard radiation domination. Extensions to (self-)interacting fields have been mentioned in Section 3, where the computation of PBH formation requires numerical analysis. This alternative reheating scenario shall thus be addressed elsewhere.

Our intention is to provide the proof of principle and the tools for the reader to adapt the mechanism to specific models of the early Universe, in order to test different inflationary models with this PBH formation and explore the detection window. To this end, we have also included a simple example that considers a primordial power spectrum with a Gaussian peak at small scales, showing the possibility of overproduction of PBHs. This is an exciting possibility that we shall study in more detail elsewhere.

Author contributions

LP: Conceptualization, Formal Analysis, Investigation, Methodology, Writing–original draft, Writing–review and editing. JH: Conceptualization, Supervision, Validation, Writing–review and editing. TG-A: Investigation, Writing–original draft. KM: Conceptualization, Supervision, Writing–review and editing. GG: Conceptualization, Supervision, Writing–review and editing.

Funding

The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. LP and JH acknowledge sponsorship from CONAHCyT Network Project 304001 “Estudio de campos escalares con aplicaciones en cosmología y astrofísica”, and through grant CB-2016-282569. The work of LP is also supported by the DGAPAUNAM postdoctoral grants program, by CONAHCyT México under grants A1-S-8742, 376127 and FORDECYT-PRONACES grant 490769. JCH acknowledges financial support from PAPIITUNAM programme Grant IG102123 “Laboratorio de Modelos y Datos (LAMOD) para proyectos de Investigación Científica: Censos Astrofísicos”. KM is supported in part by STFC grants ST/T000341/1 and ST/X000931/1.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

¹In this review we focus on PBHs formed at horizon entry. For PBHs formed inside the horizon see however Refs. Lyth et al. (2006); Zaballa et al. (2007); Torres-Lomas et al. (2014).

²This is the simplest and most widely accepted mechanism by which it is assumed that the reheating of the Universe took place. There are however, alternative scenarios such as a reheating generated by moduli fields, a massive metastable particle, or reheating from PBH evaporation [see, for example, (Allahverdi, 2020)]. In this paper we will focus on the simplest case of reheating due to the inflaton field, although many of the results shown here can be easily extended to the other mechanisms.

³This resonance band appears naturally when doing a perturbative analysis of the system, where the resonance is obtained from a Mathieu-type equation (Kofman et al., 1997; Jedamzik et al., 2010a).

⁴Backreaction effects become important when overdensities grow large enough for gravity to be of the same order as self interactions in the field (Frolov, 2008; Figueroa et al., 2021; Figueroa et al., 2023).

⁵In fact, in Refs. Nadezhin et al. (1978); Novikov and Polnarev (1980), this dependence was investigated, illustrating within a spherically symmetric analysis the process of PBHs formation and their accretion for different equation state values.

⁶A Bose-Einstein condensate is an exotic state of matter where particles clump together and behave as a single quantum entity. The idea in this studied context is that the inflaton must condense during the post-inflationary Universe, forming structures (such as solitons) as a result of this process.

⁷Normalizing a = 1 we are taking periods of evolution of the system that are not very large compared to the period of evolution of the Universe. On the other hand, the condition ρ ≫⟨ρ⟩ would be met for virialized structures that form in the post-inflationary Universe.

⁸See Ref (Guzmán, 2009) for a recent review in how to obtain the soliton solution for the complete EKG system of equations.

⁹Note that this configuration aligns completely with the Schrödinger-Vlasov correspondence discussed in Section 5.5; if we smooth the profile derived from simulations over scales larger than λ_dB, we should obtain a final distribution similar to that of an NFW profile, typically obtained in the conventional dust-like scenario.

¹⁰For the sake of clarity, we use here the term ‘halo’ to define the composite structure consisting of the central soliton together with the NFW-like envelope.

¹¹The number of e-folds after horizon crossing necessary for a perturbation to become nonlinear can be simply expressed as $\mathrm{\Delta }{N}_{\text{NL}}(k)\equiv N_{\text{NL}}(k)-N_{\text{HC}}(k)=\ln (1.39{\delta }_{\text{HC}}^{-1}(k))$ (Padilla et al., 2022). The discrepancy between this value and Eq. 20 is less than one e-fold of expansion. In this paper we will interchangeably use both quantities.

¹²In the context of the SFDM, this mechanism has been also proposed to explain the formation of supermassive black holes in the model (see, for example, (Padilla et al., 2021b)).

¹³In a more accurate approximation, δ_HC(k) follows a Gaussian distribution, with a specific number of e-folds required for the collapse of each amplitude. By taking the period required for collapse from the mean amplitude, we are underestimating the mass of primordial structures and overestimating the collapse time of the tail of the distribution.

References

Appendix A Press-Schechter formalism

In the Press-Schechter formalism (William, 1974), the likelihood of having collapsed objects with masses greater than M is analogous to the probability that a density field, after being smoothed, surpasses the threshold value δ_th:

Adopting that the probability distribution function of δ, P(δ), follows a Gaussian distribution, i.e.,

with σ(R) the variance of δ evaluated at the horizon crossing time,

$\tilde{W}(kR)=\exp (-k^2{R}^2/2)$ the Fourier transform of the window function used to smooth the density contrast over a scale R = 1/k, and ${\mathcal{P}}_{\delta }$ the power spectrum of density perturbations, we can rewrite Eq. (A1) as

In the above expression erfc(x) = 1 − erf(x) is the complementary error function. We can finally compute the abundance of PBHs of a given mass M at the time of formation, β₀(M), by using the expression

where the factor 2 is included to fit the cloud in cloud correction.