Primordial black hole formation during slow-reheating: A review

In this paper we review the possible mechanisms for the production of primordial black holes (PBHs) during a slow-reheating period {in which the energy transfer of the inflaton field to standard model particles becomes effective at slow temperatures}, offering a comprehensive examination of the theoretical foundations and conditions required for each of formation channel. In particular, we focus on post-inflationary scenarios where there are no self-resonances and the reheating epoch can be described {by the inflaton evolving in} a quadratic-like potential. In the hydrodynamical interpretation of this field during the slow-reheating epoch, the gravitational collapse of primordial fluctuations is subject to conditions on their sphericity, limits on their spin, as well as a maximum velocity dispersion. We show how to account for all conditions and show that PBHs form with different masses depending on the collapse mechanism. Finally we show, through an example, how PBH production serves to probe both the physics after primordial inflation, as well as the primordial powerspectrum at the smallest scales.


I. 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.[1][2][3][4][5][6][7][8] 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.[9,10]).The inflationary paradigm postulates that the early universe experienced a period of rapid and exponential expansion in its earliest moments (see e.g.Refs.[11][12][13][14][15][16][17]).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. [18,19]).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.
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.
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.

II. INFLATION, PREHEATING AND THE REHEATING EPOCHS
A. Inflation setting the initial conditions for reheating As mentioned above, cosmological inflation (see e.g.[18,19,49]) 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 (1) 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 (φ)] 2 ≪ 1 and η ≡ |V ′′ (φ)/V (φ)| ≪ 1 for the inflaton to produce inflation and ϵ ≃ 1 to end the inflationary epoch (see also Fig. 1 for a sketch of the inflationary potential).
FIG. 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. 1 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. [53]) 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 [54,55]: where a prime ( ′ ) denotes a derivative with respect to the conformal time dτ = dt/a, z ≡ a φb /H, and the subscript b is used to refer to background quantities.Note that the comoving curvature perturbation 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 −1 .The usual parametrization for the spectrum follows the evidence that at large scales the spectrum is almost scaleinvariant.Thus with the CMB normalization dictating ln(10 10 A s ) = 3.044 ± 0.014 and the spectral index n s = 0.9649 ± 0.0042 [56].Such prescription for the power spectrum would produce a tiny amplitude of fluctuations and produce an extremely small number of PBHs (see e.g.[57,58]).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 Sec.VII).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.

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 Fig. 1) and is believed to have created the necessary conditions for the formation of primordial nuclei and structures within the universe. 2Historically, reheating was first treated perturbatively [60,61] and in this section we review a simple version of this process.
For transfering energy, it is usual to consider a nonminimal 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 [62] 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 [9] φb A typical form of the potential used for the reheating epoch is the quadratic potential, V (φ) = m 2 φ 2 /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 [56], 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 Fig. 2.This class of potentials are considered in this work as our main study cases.
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 ρ b = φb /2 + V (φ 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 ∼ √ Γm 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 14 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 slowreheating.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 [70] k 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 3 .
The preheating instability is a non-perturbative mechanism that arises in theories with non-minimal couplings between the inflaton and other fields, say χ [20].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 [73][74][75].Furthermore, its solutions exhibit exponential instabilities, that is, χ k ∝ exp(µ (n) k mt) within a series of resonance bands, located at specific frequency ranges ∆k (n) (here labeled by the integer index n).Such instabilities lead to an exponential growth in the occupation numbers of quantum fluctuation modes, denoted as nk(t) ∝ exp(2µ k mt), which are interpreted as the production of χ particles [60,71,73,76].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.

III. PBH FORMATION DURING PREHEATING
The production of PBHs during preheating was first studied in Ref. [77].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 nonminimal interaction between fields, the Preheating process has been modeled through lattice field theory simulations that evolve the scalar field equations on a homogeneous background [74,[78][79][80][81][82][83][84], though neglecting the backreaction of inhomogeneities on the local spacetime metric 4 .
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 [88][89][90][91][92].These factors encompassed the equation of state, the nature of the inflationary potential, and the presence of additional fields or interactions 5 .Importantly, even small inflaton self-interactions can accumulate over multiple oscillations, triggering the resonant growth of non-zero momentum inflaton modes, which constantly interact tive analysis of the system, where the resonance is obtained from a Mathieu-type equation [71,72]. 4Backreaction effects become important when overdensities grow large enough for gravity to be of the same order as self interactions in the field [85][86][87]. 5In fact, in Refs.[93,94], this dependence was investigated, illustrating within a spherically symmetric analysis the process of PBHs formation and their accretion for different equation state values.
with the homogeneous component.This process of selffragmentation 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.[95][96][97][98][99][100][101][102][103][104]).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.[105,106]).Specifically, Ref. [105] 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. [107] 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 [108].
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. [77] 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. [109] 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. [110] 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.

IV. 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 [111][112][113][114].In the same context, as mentioned above, particles coupled to the inflaton can be resonantly produced [71,74] leading to prompt thermalization [115] or a potential oscillon dominated epoch [111][112][113]116].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.[60,61]).As discussed in Section II B, 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) φ 2 potential, and when Γ ≪ H.Eventually, the condensate fragments due to the gravitational growth of perturbations [117,118], 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 [119,120]).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 [90] 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 Sec.V C-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.

A. 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 [121,122].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. [31] 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 [123], Thorne's hoop conjecture [124], and Doroshkevich's probability distribution [125] to derive the production probability β 0 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 σ 2 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 [26].

B. Conservation of angular momentum criterion
The initial angular momentum of overdensities plays an important role in the formation of PBHs.In Ref. [32] 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 β 0 is restricted to where 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 (17) and assuming I = f q (q c ) = 1 (as assumed in Ref. [32]) 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 Fig. 6).
We conclude this section by mentioning that according to the study conducted in Ref. [126] 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.

C. Reheating time criterion
Another employed criterion for characterizing the generation of PBHs arising from an slow-reheating epoch can be found in Ref. [29] (see also, e.g.[30,127], 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 Sec.II B, 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 [127]: 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. [29], 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 matterdominated universe we have H = 2/(3t) and assuming 3π/(2δ 3/2 HC (k)) ≫ 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 Γ /k end = (ρ Γ /ρ 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 [29] 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 [29], 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 [32,128].
To calculate the abundance of PBHs β 0 formed in this context, we can adopt the Press-Schechter formalism (see Appendix A for details): where P[δ > δ 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, β 0 .This is due to the omission of significant physical effects discussed throughout this review.In the follwing we focus on the physics of fluctuations from an oscillating scalar field, while in Sec.VI we address the PBH formation criteria in such environment.

V. 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.[119,129,130].
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 [131] and has extended in subsequent studies [22-26, 34, 132-136].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.

A. 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. [53]), 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 [137] Further simplifcation 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 ψ ≪ mψ 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⟨ψψ * ⟩.

B. 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" [138,139].As a result, these equations are commonly named as the "quantum hydrodynamics" (QHD) equations (we recommend [140] 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.

The Madelung-Bohm formulation of quantum hydrodynamics
We can consider a Madelung transformation [141] 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.

Phase space formulation
We can also obtain the QHD equations by taking momentum moments of the SP system of equations [142].Following Ref. [143], we sketch the steps in the following.
(34) 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 1 st moment an Euler-like equation from classical fluid mechanics.Once again, in this formalism a new pressurelike 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 1 st 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) and (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.

C. 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 [137]).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 6 .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⟩ = √ 2v, 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 N .Such enhancement is due to the phenomenon of Bose-Einstein stimulation (see [137]) with Correspondingly, the scattering time is given by [145] τ where the momentum-transfer cross section σ g for Rutherford scattering is given by 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.

D. Soliton solutions
If we set the scale factor a = 1 and assuming ρ ≫ ⟨ρ⟩ 7 , the SP system of equations ( 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 equations (31) and the above ansatz has numerous solutions satisfying appropriate initial and boundary conditions [146,147].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 [148][149][150][151][152][153][154][155]).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 [156]).In addition, it has been also shown in [145] 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 halfmass radius R 1/2 and the virial velocity v vir are given by (see for example appendix B in Ref. [157]): 7 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.
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, ∂ t v = 0 = v.If we set for simplicity ∇ ∼ 1/R sol , where R sol is a characteristic radius of the soliton profile (typically λ dB ), we obtain, And the condition of equilibrium configuration yields, and ∂ρ ∂t = 0.
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 mantains the same parameter dependence found in the exact numerical treatment, Eq. ( 44).
A general-relativistic regime.-Wecan anticipate from the relation M sol ∝ R −1 1/2 (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)8 (see Fig. 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 [152][153][154][155] 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.
FIG. 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 [159][160][161][162][163].However, as previously discussed, these excited states have been found to be unstable and undergo gravitational cooling, ultimately transitioning to the ground state solution.

E. 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 [143,[164][165][166][167] we choose to utilize the Husimi representation of ψ, which is a smoothed phase space representation.This representation, introduced in [168], 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 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 F in a manner akin to the Wigner-Moyal equation.This involves computing ∂F/∂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 F with the phase space distribution function of its collisionless Nbody 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.

The fluid approximation
Similarly than in Sec.V B 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 [143], which the reader should consult for a more comprehensive discussion.However, in particular the 0 th and 1 st 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 ≡ ρσ 2 ij , where σ 2 ij 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 V B 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. [22] (see also [27]).In a subsequent section, we describe the necessary conditions under which the collapse of inflaton perturbations can indeed take place and form PBHs.

F. 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.[169,170].In that work, the authors solved the SP system of equations, Eqs.(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 Sec.V D) surrounded by an NFW-like envelope from a radius determined by the incoherent fluctuations of the scalar field 9 .A schematic plot of such configurations is shown in Fig. 4.This solitonenvelope structure has been confirmed to exist by several studies that considered simpler scenarios [171][172][173], e.g., in the nearly simultaneous merger of several soliton configurations [174].In the context of reheating, such core-envelope structures have also been reported [134].

FIG. 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 Sec.V D.
Numerical simulations prescribe a soliton-halo 10 mass relation given by M sol ∝ M 1/3 halo /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 [175], while others suggest that their specific kinetic energy is what is comparable [176].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 [174,177].Regardless of the interpretation, the above studies agree in the following relation: In the above expression M halo = (4π/3)ρ 200 (a NL )R 3 halo 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, ρ 200 (a NL ) = 200ρ b (a NL ), and subindex NL refers to quantities evaluated at the time the k-mode becomes non-linear 11 .We can re-express Eq. ( 54) more conveniently as follows: 10 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. 11The number of e-folds after horizon crossing necessary for a perturbation to become nonlinear can be simply expressed as [22].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.
Or equivalently, with some more algebra, In the above expressions m 5 ≡ m/(10 −5 m Pl ) and ρ 11 (a) ≡ ρ 200 (a)/(10 11 GeV) 4 .This is an important result since, for a given m 5 , 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 5 = 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 −4 g.

VI. PBH FORMATION DURING SLOW-REHEATING: COLLAPSE FROM PRIMORDIAL STRUCTURES
A. 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.[22].

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 ∆N NL (k) = ln[1.39δ−1 HC (k))] (see footnote 11), we can finally obtain the condition for PBH formation:

Soliton collapse
In accordance with our previous discussion in Sec.V D, 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 δ (soliton) th beyond which the central soliton becomes unstable and may collapse to form a PBH.The condition for PBH formation in this scenario is given by12 : B. 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.Fig. 5 summarizes this general timeline, which we will discuss in more detail below.
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 NL (k) = [2/(3H end )][e NHC(k) 1.39/δ 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.Once the halo undergoes virialization, the condensation of the inflaton at its core generates a solitonlike 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.
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 IV).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, β 0 , in terms of the variance by assuming a Gaussian probability density.We thus bring all these effects together in Fig. 6, which shows limits to the abundance at the time of formation, β 0 , 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.
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: 64) 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 efolds 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 Fig. 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 acheived 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. [25].Illustrating this mechanism in a plot equivalent to Fig. 6 is a task to be tackled in future work.17) and ( 18), respectively.The solid black and dashed cian lines were plotted by using the standard Press-Schechter formalism, Eq. ( 25), and the threashold values ( 59) and ( 60), respectively.For comparison, we also included the PBH density fraction β 0 calculated in a standard radiation-dominated universe, where we took δ (rad) th = 0.41 (see for example [90]).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.

VII. 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 −1 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 Σ 2 p .For this example we shall consider the model parameters A s = 2.099 × 10 −9 , n s = 0.9634, B s = 0.084, Σ p = 0.03k end , k p = 0.6k end , and k end = 0.346 m −1 .The power spectrum generated for this set of parameters is sketched in Fig. 7.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's worth noting that this figure bears a striking resemblance to Fig. 2 in [23], 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 [23].
From the bound in the tensor-to-scalar ration (r ≤ 0.032) [179], 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, ρ end = (1.368× 10 16 GeV) 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). 13aving computed that, in Fig. 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.

FIG. 8:
The density fraction β 0 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).
For each one of the scales that manage to form halos or solitons, we compute the probability of PBH formation.For this purpose, we employ 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) and (60).The result obtained is shown in Fig. 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.
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. [23].
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 grams 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 Fig. 8, we have included the respective constraints applicable to each scenario.To derive these constraints, we used the publicly available software PBHBeta [180], 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 (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. [24]).

VIII. 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 nonspherical 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 pescribes a bound to the abundance of PBHs, as illustrated in Fig. 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 Fig. 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 Fig. 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 [181,182], curvaton [183,184] or multiple fields [185,186], 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) oscilating 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 Sec.III, 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.

FIG. 5 :
FIG. 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 solitonlike 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.

FIG. 6 :
FIG. 6:The PBH density fraction β 0 plotted as a function of the variance of the density contrast σ.The blue and green lines are plotted by following Eqs.(17) and (18), respectively.The solid black and dashed cian lines were plotted by using the standard Press-Schechter formalism, Eq. (25), and the threashold values (59) and (60), respectively.For comparison, we also included the PBH density fraction β 0 calculated in a standard radiation-dominated universe, where we took δ

FIG. 7 :
FIG. 7: Power spectrum as a function of the scale wavenumber k.