PBH formation from spherically symmetric hydrodynamical perturbations: a review

Primordial black holes, which could have been formed in the very early Universe due to the collapse of large curvature fluctuations, are nowadays one of the most attractive and fascinating research areas in cosmology for their possible theoretical and observational implications. This review article presents the current results and developments on the conditions for primordial black hole formation from the collapse of curvature fluctuations in spherical symmetry on a Friedman-Lemaitre-Robertson-Walker background and its numerical simulation. We review the appropriate formalism for the conditions of primordial black hole formation, and we detail a numerical implementation. We then focus on different results regarding the threshold and the black hole mass using different sets of curvature fluctuations. Finally, we present the current state of analytical estimations for the primordial black hole formation threshold, contrasted with numerical simulations.


Introduction
One of the great mysteries in science is the composition of dark matter, which accounts for 27% of the present Universe [1]. Although there are different theories and candidates that try to explain it [2], still, the answer remains elusive. One of the most promising possibilities is Primordial Black Holes (PBHs), i.e., Black Holes (BH) generated at earlier than star formation times and therefore not of stellar origin. For the current observational status and constraints of PBHs in the form of dark matter, we refer the reader to [3][4][5][6][7][8].
PBHs were first considered in [9][10][11]. They could have formed in the very early Universe due to the gravitational collapse of cosmological perturbations. Under this scenario, PBHs could have been generated as a consequence of high non-linear peaks in the primordial distribution of density perturbations, leaving open the possibility that PBHs could constitute a significant fraction of Dark Matter (DM) .
PBHs with a size smaller than M PBH < 10 −20 M would have already evaporated due to Hawking radiation [39]. Therefore, those with higher masses can account for dark matter. In addition, they may be responsible for seeding supermassive black holes at the centre of galaxies [40,41], generating large-scale structures through Poisson statistics [42], or changing the thermal history of the Universe [43].
Although PBHs were theorised a long time ago, it has not been until now that their popularity has increased exponentially. The current big interest in PBHs (and their golden age) is due to the first detection of gravitational waves from a black hole merger by the LIGO/Virgo collaborations [44]. Soon after this remarkable discovery, it was suggested that the constituents of the black hole merger could have had a primordial origin [19,45]. Later on, several analyses of the data signals suggested that the population of black holes detected could have a primordial origin [27,[46][47][48][49][50] (see also the references therein). On the other hand, recent results from the NANOGrav experiment [51] have also been connected with PBHs [26,[52][53][54][55][56][57][58][59].
If formed by the collapse of inflationary perturbations, the abundance of PBHs is exponentially sensitive to the threshold of the gravitational collapse δ c [60] (where δ c is the minimum amplitude of the gravitational potential peak related to the perturbation undergoing gravitational collapse and leading to a BH). In order to obtain the necessary precision on δ c , a numerical analysis of PBH formation is an obvious way to go.
Numerical simulations of PBH formation started some time ago, the first one being performed in [61], where the non-linear behaviour of the gravitational collapse was studied. Later on, in [62,63], δ c was computed and a scaling behaviour for the PBH mass was found whenever the amplitude of the perturbation was close to δ c [62][63][64][65]. The value of the scaling exponent found was consistent with the one given in the literature and computed from a perturbative treatment [66,67] or from numerical simulations [68] in asymptotically flat spacetime. While the constant of proportionality appearing in the scaling law depends on the specific shape of the perturbation considered, the scaling exponent is a universal quantity only dependent on the type of fluid. Earlier simulations of PBHs were mainly based on the application of a Lagrangian hydrodynamic code based on finite differences and developed from an earlier work [69], solving Misner-Sharp equations [70], which describe the motion of a relativistic fluid in a curved spacetime. An inconvenience of this procedure is the appearance of a coordinate singularity after forming the apparent horizon, which breaks down the simulation. To resolve that issue, Hernandez-Misner equations [71] (which are the Misner-Sharp equations in null coordinates) are used to evade the formation of an apparent horizon and follow the subsequent evolution to determine the value of the black hole mass. The method is based on [72]. As was shown in [73], the formation of PBHs from the collapse of initial fluctuations at super-horizon scales does not form shocks (a discontinuous solution) in comparison with the earlier result of [64], where shock waves were found, but using initial conditions on sub-horizon scales.
On the other hand, even if a numerical code allows computing the threshold, the realisation of a numerical simulation can be quite expensive to be actually useful for statistical applications such as the calculations of PBH abundances. Therefore, an analytical expression would drastically simplify the problem. There are some analytical estimations of δ c (e.g., [74,75]) that were based on analytical models under certain assumptions. However, numerical simulations have shown that, even for a fixed equation of state, δ c is not universal [63,64,73,[76][77][78][79][80]. The main argument is that δ c is dependent on the specific detail of the initial curvature fluctuation [79], i.e., on the scale dependence or "shape" of the perturbation, something that these analytical estimations do not take into account.
Some works have tried to characterise the threshold for black hole formation in terms of the initial density perturbation profile. In [81], the amplitude of the curvature fluctuation and its second radial derivate at the centre were used to characterise the initial configuration profiles, showing that the probability of PBH formation is sensitive to them. Later on, in [82], numerical simulations were used to argue the formation of a PBH from an over-density peak in an FLRW Universe assuming the perturbation to be initially at super-cosmological scales, and it only depends on two master parameters: the integral of the initial curvature perturbation and the edge of the over-density distribution. The paper [79] recently refined the arguments of [82] by showing that these parameters may be more conveniently given in terms of the amplitude of the "gravitational potential" at its maximum (r = r m ), as already noticed in [77], and that it mainly depends on the functional form (shape) of the gravitational potential. More precisely, in [79], the threshold δ c was identified with the peak of the "compaction function" [77] at super-horizon scales. The compaction function closely resembles the Schwarzschild gravitational potential and is defined as twice the local excess-mass over the co-moving areal radius.
Following the aim to have a more accurate analytical estimation for the threshold of PBH formation, recently, in [83], it was argued that the threshold for primordial black hole formation should be quite insensitive to the physics beyond r m : the threshold is the amplitude for which a "black hole" with zero mass would be formed, taking an infinite time. This is an asymptotically and ideal limit. Therefore, all the over-density beyond r m will be diffused away while that just in the proximity of r m will block the collapse.
With those arguments, it was shown that during a radiation-dominated epoch (equation of state p = wρ with w = 1/3), with a very good approximation, there exists a universal threshold value for the volume-averaged compaction function, which is actually shape-independent. Taking into account that the volume average is dominated by scales near the maximum of the compaction function, in [83], it was shown that it is enough to parametrise the profile dependence of δ c by the curvature shape around the compaction function at its maximum. This result was used to make an analytical approximation for δ c taking into account the shape dependence on the curvature fluctuation, which matches that found in simulations to within a few percent [83].
From the results in [83], there was the question of whether or not this universality on the averaged critical compaction function is generic and not only an accident of radiation. Thus, in [84], the universality question beyond the case of a radiation fluid was addressed. In [84], the arguments of [83] were generalised to provide a new analytical estimation for the threshold in terms of the equation of state and the shape profile. This analytical formula was tested in front of numerical simulations and was found with an accuracy better than a 6%.
The other half-part of the problem is the determination of the final PBH mass. When a perturbation collapses, it will firstly form an apparent horizon with an associated mass M PBH,i . After that, due to the surrounding fluid of the FLRW background, there is an accretion process up to a final mass M PBH, f when the PBH is formed. As we have indicated, for fluctuations whose amplitude is near the threshold (the critical regime), a scaling law was found for the PBH mass in terms of the amplitude of the perturbation [62][63][64][65]. Although the importance of this result was pointed out in the context of critical collapse [85], currently, it is substantially relevant because, as indicated in the current literature [60], the PBHs with a higher probability of formation are those that are formed in the critical regime. The probability distribution function is exponentially suppressed beyond the threshold, and therefore, the scaling law can be directly used for estimations of PBH abundances.
On the other hand, some studies have also considered the maximum size of PBHs at the formation time and the accretion effect from the surrounding background in the final mass of the PBH M PBH, f . It was already noticed in [74,86] that the accretion effect for small PBHs should be small. Specifically, in [86], some analytical expressions were derived for the upper bound of a PBH formed by the collapse of hydrodynamical perturbations. In [87], the M PBH, f was computed numerically for different curvature profiles and the accretion effect was quantified, which could be substantially large for PBHs formed beyond the critical regime O(10%).
In conclusion, considering the recent developments in the field and future perspectives, we present this review paper to give the current state-of-the-art regarding the analytical, numerical, and theoretical results about PBH formation from the collapse of hydrodynamical perturbations on an FLRW background under spherical symmetry, currently the most common and popular scenario in the literature [39]. We hope that this review paper allows the scientific community to be introduced to the topic of PBH formation and understand the main insights.
We want to tell the reader that there are fascinating and valuable review papers on the PBH topic with a different perspective than the one considered here. For an interesting historical review of the dark matter candidates, we refer the reader to [2]. Regarding reviews that focus on gravitational waves applied to PBHs and the different constraints to account for dark matter, we have [3,14,15,39,[88][89][90][91]. Another interesting review with potential applications to PBHs is [92], which is focused on scenarios with deviations from radiation domination in the early Universe.
This review paper is organised as follows: In Section 2, we present the basics of PBH formation. In Section 3, we give some basic details about the statistical estimation of PBH abundances. In Section 4, we give the necessary ingredients to set up the formalism of PBH formation, in particular: the differential equations that describe the collapse of the perturbations, the initial conditions, the definition of the compaction function, and the threshold criteria. In Section 5, we discuss and comment on the different numerical techniques employed in the literature, and we focus on a specific method that is publicly available. In Section 6, we give numerical results regarding the threshold and PBH mass in terms of the specific profiles using numerical simulations. In Section 7, we discuss in detail the different developments for the analytical estimation of the threshold. Finally, in Section 8, we enumerate other scenarios of PBH formation not explicitly considered in this review.

Basics of PBH Formation
As we pointed out in the previous Section 1, PBHs could have been formed in the very early Universe during radiation domination due to the gravitational collapse of large curvature perturbations generated during inflation [9,10].
Under this scenario, the collapse or dispersion of those hydrodynamic perturbations depends on the perturbation's strength: If its greater than a given threshold, the perturbation will collapse and form a BH after the perturbation re-enters the cosmological horizon. Otherwise, if it is lower, it will disperse because of pressure gradients preventing the collapse (a schematic picture can be found in Figure 1). Both things could happen, i.e., the perturbation could undergo gravitational collapse and subsequently bounce. In this case, the fluid is dispersed and collapsed continuously, making rarefaction waves (we will see this in more detail in Section 6). This behaviour is particularly evident when the initial strength of the perturbation is very close to its threshold value. Therefore, pressure gradients play a crucial role in the collapse of the perturbation. If pressure gradients are strong enough, they will prevent gravitational collapse, which implies having a large threshold compared to the situation when the pressure gradients are small.  A big effort has been made during the past few decades to find the correct PBH formation criteria. The first estimation of a threshold was given by [9], using a Jeans length argument and Newtonian gravity. The criterion for the formation of a BH used was that the size of an over-density at the maximum expansion should be larger than the Jeans length, but also smaller than the particle horizon. This translates to the requirement that the peak value of the density contrast at scales smaller than the cosmological horizon must be at least w (where w is the equation of state of the perfect fluid that fills the FLRW Universe) in order to collapse.
A simple, intuitive picture of PBH formation and its threshold connected with the Jeans length argument can be found in [15], and we summarise it here. Consider a cosmological perturbation at super-horizon scales. The spacetime metric can be written as a closed FLRW Universe (we will see this in more detail in Section 4) with a radial dependence on K(r), where a(t) is the scale factor. If we ignore the radial derivative of K(r), the time-time component of the Einstein equations is given by, which basically corresponds to the FLRW equation, but with a small non-homogeneity term given by K(r).
The previous equation allows defining a density contrast as, where ρ b is the energy density of the background. A region will collapse instead of expanding when K > 0 if we ignore the radial dependence on K(r). This precisely happens when 3K/a 2 = 8πρ, which corresponds roughly to when the Hubble scale becomes equal to the length scale of the curved region. In this case, the assumption made with Equation (1) is no longer valid, but still can be assumed to obtain some threshold estimation.
When δρ(t c )/ρ b = 1 is the time when a homogeneous and isotropic Universe would stop expanding, we therefore assume that this epoch is the time of black hole formation t = t c , in particular, The criterion is that perturbations cannot collapse if they have smaller scales than the Jeans length, which immediately implies that k 2 = a 2 H 2 /c 2 s , where c s = √ w is the sound speed of the fluid (perfect fluid in our case) and k is the wavelength mode associated with the perturbation. Therefore, the condition that should satisfy the density contrast at the time when the perturbation re-enters the horizon ( which is basically the threshold value estimation obtained in [9]. The second analytical estimation came years later in [75]. The authors considered a compensated "three-zone" model, i.e., a central over-dense region followed by an under-dense layer (which compensates the over-density), and finally, the FLRW background, to estimate the threshold for BH formation. Following an argument about the sound waves at the maximum expansion of the perturbation, it was found that the threshold of the over-density (at r = 0) when it renters the horizon must be (3(1 + w)/(5 + 3w)) · sin 2 (π √ w/(1 + 3w)). Later on, however, it became clear that the threshold depends on the shape of the curvature perturbation [63,65,82]. Moreover, the threshold also depends on the specific equation of state of the fluid since pressure gradients play a role in its determination: the larger is w, the stronger are the pressure gradients, and therefore, the higher is the threshold.
In [77], a new criterion for PBH formation was introduced. Importantly, it was found that the peak of the compaction function (the average mass excess on a given volume) was a good criterion for PBH formation, and it currently has become the standard definition of the threshold, which we call δ c . Prior to this criterion, some other measures for the threshold "definition" were not empty of some issues. See [79] for a historical perspective about that.
In [79], simulations were performed for a radiation fluid with a set of curvatures profiles leading to a range of threshold values given by 0.41 δ c ≤ 2/3 using the criterion of the compaction function's peak as the threshold definition. Later on, in [83], the minimum threshold value was refined, obtaining 0.4 rather than 0.41 in the case of a radiation fluid. This difference was due to the choice of unphysical profiles in some extreme limit in [79].
Although δ c is a profile-dependent quantity, in [83,84], it was shown that there exists an (approximately) universal value given by the averaged critical compaction function, which allows building a sufficiently accurate analytical expression for δ c , which only depends on the type of fluid and the curvature around the peak of the compaction function.
On the other hand, if a perturbation is sufficiently strong, it will collapse and form an apparent horizon. Later on, it will follow a process of accretion of energy density from the surrounded FLRW background until reaching a stationary state, when the final mass of the PBH is achieved. As was shown in [62,63], the final PBH mass follows a critical scaling law when the amplitude of the perturbation is close to its threshold value (the critical regime), with values of γ consistent with the previous numerical computation [85,93] γ ≈ 0.356. The scaling exponent is universal and only depends on the type of fluid, not on the initial condition. M H is the mass of the cosmological horizon at the time of horizon crossing (when the length scale of the perturbations equals the cosmological horizon). The constant K depends on each initial condition used. As in the case of the threshold, the values of K are different depending on what the initial conditions are, but with values O(1) [94]. The simulations in [63] were performed for δ − δ c 10 −5 , which was not sufficient to test the critical regime up to very small values, i.e., up to machine precision δ − δ c ≈ 10 −15 . In [73], the scaling law up to machine precision was finally verified, and this was tested with an explicit example of a Gaussian perturbation.

The Abundance of PBHs
This section gives a brief review of the two approaches commonly used in the literature to perform PBH statistics (formed in a radiation-dominated Universe) and estimate their abundances. In this review, we do not focus on the statistical estimation of PBH abundances, for which there are extensive works in this direction, with different approaches and methods. If the reader is interested, we suggest checking these references and the ones therein [60,[95][96][97][98][99][100][101][101][102][103][104][105][106][107][108][109][110][111][112][113][114]. However, to show this topic at the basic level, we basically describe the Press-Schechter [115] and the peak theory [116] procedures. This section aims to clearly point out the relation between the threshold and PBH mass on the estimation of PBH abundances, especially why the threshold should be known with enough accuracy to be useful for statistics.
In both cases, we consider the amplitude of the peak of density contrast in the comoving slicing at super-horizon scales (usually denoted in the literature as ∆ ≡ δρ (r = 0, t)/ρ b (t) [60,117]), as a statistically distributed variable. The power spectrum associated with the density contrast is defined as, where δ D is the Dirac delta and the moments j of P ∆ are given by, In the Press-Schechter formalism, we make two assumptions: (i) the density contrast field ∆ is a Gaussian variable; (ii) perturbations with ∆ > ∆ c will collapse and form a PBH. Therefore, basically, we integrate the probability distribution P(∆) over the range ∆ c ≤ ∆ < ∆ max , where ∆ max is the maximum allowed value. In practice, we integrate up to ∆ max → ∞ since the probability distribution is a rapidly decreasing function above ∆ c , and therefore does not change the result and allows simplifying the computation. It is important to notice that ∆ c = δ c , since we are comparing the critical density contrast with the peak of the critical compaction function (as we will see later, it is related to the averaged density contrast).
Finally, consider a Gaussian distribution probability distribution, the abundance of PBHs can be computed as, where Γ U (a, b, c) is the confluent hypergeometric function, 1 and the factor 2 at the beginning of the integral put by hand is introduced to avoid the under-counting that happens in Press-Schechter theory. It has been assumed that the scaling law is accurate always even when ∆ ∆ c , something not true, as we will show in Section 6. However, since the probability distribution is exponentially smaller for large ∆, Equation (10) is accurate. In the case of considering a monochromatic mass spectrum, i.e., all the PBHs formed to have the same mass (in particular, the horizon mass M H ), the abundance is given by 2 , This can be obtained by simply taking the limit γ → 0 and K =K in Equation (10). The procedure in the peak theory approach is a bit different. It makes statistics on counting the numbers of over-threshold peaks on the over-density; another approach used the curvature fluctuation ζ peaks instead [102,118].
We summarise the procedure performed in [60] to estimate the abundances through peak theory. 1 The confluent hypergeometric function is defined as Γ U (a, b, z) = 1

Γ(a)
∞ 0 e −zt t a−1 (1 + t) b−a−1 dt. 2 The erf function is defined as erf(x) = 2 First of all, we define the variable ν = ∆/σ 0 , for which in the rare peak assumption, we have that ν 1. The number of rare peaks is given by [116], where t f is the time when the PBHs are formed. It is clear that we consider only that peaks higher than the threshold value ν c will contribute to the formation of PBHs. Therefore, the abundance of PBHs measured at formation with respect to the background energy-density can be computed as, The scaling law mass, in terms of ν, is given by, whereσ 0 = σ 0 a 2 H 2 and t m is the time of horizon crossing. Taking into account that numerical simulations have shown a f ≈ a m ≈ 3 [79], a f is weakly dependent on ν, as pointed out in [60], and using the saddle point approximation [119] ν s ≈ ν c + γ/ν c , finally, we obtain, where k * = σ 1 /( √ 3σ 0 ). From Equations (10), (11), and (15), we clearly see the exponential dependence on the threshold for PBH formation and the linear dependence on the constant K associated with the scaling law. This is the reason why an accurate determination of the threshold is essential, since the abundances depend exponentially on it. These two methods are however only an approximation of the true statistics, which consider the fact that each statistical realisation of the profiles has a different threshold. This was developed in [95], so we will not discuss it here as it is beyond the scope of this review paper. However, as can be seen in [95], one finds again the generic behaviour that β is exponentially sensitive to the threshold. We also mention that, in principle, a window function has to be used to estimate PBH abundances correctly. It is necessary to relate the power spectrum in Fourier space to the probability distribution function in real space by coarse-graining the cosmological perturbations. The necessity of considering a window function is remarkably important in the case of having a broad power spectrum (where several wavelength scales k are involved). Despite that, the choice of the window function is not unique, and the result depends on the choice used. We refer the reader to [105,111,120] to check the details.

Cosmological Setup for PBH Formation
This section reviews the needed ingredients to simulate the formation of PBHs in an FLRW Universe under spherical symmetry and filled by a perfect fluid. We will see what differential equations describe the gravitational collapse, the consistent initial conditions for PBH formation we need to use, and the boundary conditions to supply. Later on, we will study the formation of the apparent horizon and characterise it in terms of the expansion of the congruences. Finally, we will describe the process of the accretion of the PBH mass from the FLRW background and the estimation of the final mass.

Misner-Sharp Equations
The Misner-Sharp equations [70] describe the motion of a relativistic fluid with spherical symmetry. This corresponds to the Einstein field equations written in the comoving gauge (the gauge we use in this review paper). To obtain them, first of all, we need to consider a perfect fluid with the energy-momentum tensor, and with the following spacetime metric in spherical symmetry: where R(r, t) is the areal radius, A(r, t) is the lapse function, and dΩ 2 = dθ 2 + sin 2 (θ)dφ 2 is the line element of a two-sphere. The components of the four-velocity u µ are given by u t = 1/A and u i = 0 for i = r, θ, φ, since we are considering comoving coordinates (comoving gauge). Throughout the review paper, we use units G N = c = 1. Solving the Einstein field equations for this problem, the following quantities appear: where D t and D r are the proper time and distance derivatives, respectively. U is the radial component of the four-velocity associated with an Eulerian frame (so not comoving), which measures the radial velocity of the fluid with respect to the origin of the coordinates. The Misner-Sharp mass M(r, t) is introduced as: which is related with Γ, U and R though the constraint: where Γ is called the generalised Lorentz factor, which includes the gravitational potential energy and kinetic energy per unit mass. This also can be seen in the Newtonian limit, Finally, the differential equations governing the evolution of a spherically symmetric collapse of a perfect fluid in general relativity are: We refer the reader to [70] to see the details of the derivation. The boundary conditions in the Misner-Sharp equation are R(r = 0, t) = 0, leading to U(r = 0, t) = 0 and M(r = 0, t) = 0. Then, by spherical symmetry, we have D r p(r = 0, t) = 0, and also, we consider a perfect fluid p = wρ. At r → ∞, we should recover the solution of the FLRW background. However, in a numerical simulation, we have to handle the discretisation of time and space. Therefore, the condition D r p(r = r f , t) = 0 is used (where r f is the outer point of the grid) to match with the FLRW solution and avoid reflections from density waves. On the other hand, Equation (26) is called the Hamiltonian constraint, which is commonly used to check the correct solution of the numerical equations. For the case p = wρ, the lapse equation in Equation (27) can be solved considering A(r f , t) = 1 to match with the FLRW spacetime, where ρ b (t) = ρ 0 (t 0 /t) 2 is the energy density of the FLRW background and ρ 0 = 3H 2 0 /8π. Actually, the Misner-Sharp equations can be written in a more advantageous way to perform numerical simulations using the definitions of Equation (18), where the dotU represents the time derivate ∂U/∂t and ρ the radial derivative ∂ρ/∂r.

The Gradient Expansion Approximation, Initial Conditions, and Compaction Function
The gradient expansion method [121] (also called long-wavelength approximation) was used in [122] to set up consistent initial conditions for PBH formation on an FLRW background at super-horizon scales. First of all, let us consider a cosmological perturbation at super-horizon scales, i.e., with a length scale R m much larger than the Hubble horizon. We can define a parameter to relate the two scales: the Hubble horizon R H and the length scale of the perturbation R m , where R H = 1/H. It is clear that at super-horizon scales, we will have 1. The gradient expansion method allows expanding these inhomogeneities in the spatial gradient in terms of . In the limit → 0, the spacetime locally corresponds to the FLRW metric, when the perturbation is smoothed out at sufficiently large scales R m .
Consider a general spacetime metric in the 3 + 1 Arnowitt-Deser-Misner (ADM) formalism [123,124], in general, the spatial part of the metric can be decomposed into: where γ ij , β i , and α are the spatial metric, shift-vector, and lapse function, respectively. ζ is usually called the curvature fluctuation. Under the gradient expansion method, it was shown in [77,[125][126][127] that . Therefore, in spherical symmetry and in the limit → 0, the metric of Equation (35) can be written as, Notice the time independence of ζ(r) since at super-horizon scales, it is shown thatζ = O( 2 ). It is important to point out that although, we consider the comoving gauge, other gauges are possible. In particular, in [77], the constant mean curvature gauge was used instead. However, as was shown in [76,126], the differences between the two gauges in ζ are of O( 2 ). Finally, the metric Equation (37) can be also written in other coordinates as an FLRW metric with a non-constant curvature K(r), The change of coordinates r,r between the two metrics and the conversion between K(r), ζ(r) are given by the following relations, which were obtained in several works previously (not necessarily in the context of PBH formation for all of them) [79,81,109,128]: Therefore, the cosmological perturbation will be characterised by the curvature perturbation K(r), ζ(r). The Misner-Sharp equations can be solved at leading order in 1 using the long wavelength approximation as was performed in [122], where for → 0, we recover the (FLRW) solution. The perturbations of the tilde variables at first order in gradient expansion were computed in [122], which we summarise here: The perturbations in ther coordinate with ζ(r) are [79], From Equations (41) and (42), an expression for the density contrast at super-horizon scales in terms of the curvature fluctuations can be obtained [79], given by: The areal radius of the length scale R m is given by R m (t) = a(t)r m or also R m (t) = a(t)r m e ζ(r m ) with ther coordinate taking into account Equation (39). As was shown in [129], where the gradient expansion method beyond first order was applied using an iterative scheme (see [129] for more terms in the expansion in comparison with Equation (40)), the long-wavelength approximation at first order is accurate enough if a sufficiently small parameter is taken.
When we take = 0, we recover the background solution equations: We also consider as the initial condition that r m = ΛR H (t 0 ), where Λ is some number (usually taken as O(10) in simulations). A time scale given when (t m ) = 1 is commonly defined in the literature, which precisely corresponds to the time when the Hubble horizon equals the length scale of the perturbation (the time of horizon crossing) t m = t 0 (a 0 r m /R H (t 0 )) 1/(1−α) . Notice that t 0 is the initial time, and it is assumed to be a time where the perturbation is at super-horizon scales.
On the other hand, the amplitude of a cosmological perturbation can be measured as the averaged density contrast within a spherical region: where δρ = ρ − ρ b and V = 4πR 3 /3. Related to that, in [77], the compaction function C(r) was defined, which gives a measure of the mass excess on a given volume. In particular, where M b (r, t) is the mass of the FLRW background in the volume V. The compaction function can also be written in terms ofδ, as was shown in [76], At leading order in ,δ where: From the above definitions, we finally define the compaction function at super-horizon scales at leading order of gradient expansion as: The previous equation can be obtained introducing the initial condition of Equation (40) into Equation (45) and making an expansion in . Notice in Equation (50) the time-independence of the compaction function at super-horizon scales at leading order in gradient expansion, since the curvature fluctuations remain frozen at super-horizon scales, as we have mentioned before. The compaction function at super-horizon scales is an essential magnitude, since it allows defining the threshold. In particular, we now define the location of the maximum of C(r) as r m , and its value C max = C(r m ) is used as a criterion for PBH formation, as was proposed in [77], while r m (r m in the case of using ther coordinate) is considered the length scale of the perturbation. The threshold for PBH formation will correspond to the peak value of the critical compaction function, i.e., C c (r m ) = δ c . Therefore, we define the "amplitude" of the perturbation as δ m = C(r m ). Perturbations with δ m > δ c will collapse and form a PBH. In the opposite case, perturbations with δ m < δ c will disperse on the FLRW background with no black hole formation. Notice that from Equation (38), we have a bound on the amplitude value δ m at super-horizon scales given by δ m,max = f (w), since K max (r m )r 2 m = 1.
Because of the previous definitions, the value of r m is given by the solution of C (r) = 0: The definition of the compaction function at super-horizon scales using the ζ(r) perturbation instead of K(r) leads to a slightly different definition, which yields the following condition for the peakr m : and therefore, . Interestingly, in [79],δ(r m , t) = 3δρ(r m , t)/ρ b was found at super-horizon scales, which relates the averaged mass excess with the local density contrast at r m . In this review paper, we focus on type I PBHs, which are those that fulfil that the areal radius R is a monotonic function, i.e., R > 0. Instead, type II PBHs satisfy R < 0 [130]. Type II PBHs are still unexplored in the literature from a numerical point of view or even with a general analytical treatment.

Horizon Formation
If an initial perturbation at super-horizon scales has an amplitude δ m bigger than its threshold δ c , the perturbation will continue growing and, at some point, a trapped surface will be formed. This indicates the onset of gravitational collapse.
To identify when trapped surfaces are formed, we have to consider the expansion Θ ± of null geodesics' congruences k ± , orthogonal to a spherical surface Σ. The expansion Θ ± is defined as Θ ± ≡ h µν ∇ µ k ± ν , where h µν is the spacetime metric induced on Σ. There are two congruences: we call them inwards k + µ and outwards k − µ , whose components are k ± µ = (A, ±B, 0, 0) with k + · k − = −2. In the case of flat spacetime, Θ − < 0 and Θ + > 0, and these surfaces Σ are called normal surfaces. On the other hand, if Θ ± < 0, the surface is called trapped, while if both are positive Θ ± > 0, the surface is anti-trapped. In our case, we have that, In spherical symmetry, we can consider that any point (r, t) is a closed surface Σ with a proper radius R. These points can be classified as normal, trapped, and anti-trapped. Specifically, the transition from a normal to a trapped surface should satisfy Θ − < 0 and Θ + = 0, which corresponds to a marginally trapped surface, usually called the "apparent horizon". Taking into account that , the condition for the apparent horizon is given by In Figure 2, 2M/R and the congruences Θ ± for a supercritical (with an amplitude δ m > δ c ) Gaussian curvature fluctuation can be seen, after the first horizon has been formed. The same qualitative behaviour can be applied for any perturbation of type I PBHs (with w = 0).
Looking at Figure 2, we can see three points where 2M = R, and therefore three different surfaces that we can characterise: • Surface (A-B): In this region, we have a transition between Θ + > 0 and Θ − < 0 (normal region) to Θ + < 0 and Θ − < 0. This horizon, moves inwards of the computational domain, which eventually will encounter a singularity at r = 0; • Surface (B-C): In this region, we have a transition between Θ + < 0 and Θ − < 0 to Θ + > 0 and Θ − < 0 (normal region). This horizon, typically called "the apparent horizon", moves outwards of the computational domain. Actually, this horizon and the previous one emerge from a single marginally trapped surface. During the evolution, the motion of the horizon can be followed to obtain the PBH mass M(r * , t), where r * (t) is the location of the outer horizon in time. At the final stage of the gravitational collapse for very late times, this horizon will become static, i.e., r * = const; • Surface (C-D): In the last region, we have a transition between Θ + > 0 and Θ − < 0 (normal region) to Θ + > 0 and Θ − > 0. This case corresponds to an anti-trapped surface, which corresponds to the cosmological horizon, moving outwards. We should mention that this horizon does not correspond to the cosmological horizon R H from the FLRW background, since the perturbed medium (deviating from the FLRW solution) affects the evolution of the cosmological horizon.

PBH Mass and Accretion
Once the apparent horizon has formed, the initial PBH mass, i.e., the mass of the PBH M PBH,i at the moment of formation of the first apparent horizon t AH , will start to grow to arrive at a stationary value M PBH,f . This situation is different from the case of dust collapse w = 0, where the mass would increase forever since no pressure gradients would avoid the accretion of energy-matter.
The process of accretion of energy-matter from the FLRW background has been studied in the past [9,11,86,141,142] in different scenarios. To estimate the final mass of the PBH, one could follow the motion of the apparent horizon until very late times when the horizon remains static. However, since it could be computationally very expensive, we can use the Zeldovich-Novikov formula Equation (55) [11], which considers Bondi accretion [11,143,144], and in our case, assume that the energy density right outside the apparent horizon decreases as in an FLRW Universe. It is essential to indicate that this is not valid at the time t AH , since it neglects the cosmological expansion of the spacetime [145], but we can apply from sufficiently late time for t t AH considering an effective constant accretion efficiency rate F [143,144], which depends on the non-linear process of the accretion flow and the equation of state. This approximation was already used numerically in the context of PBH formation from domain walls in [146] and from curvature fluctuations in [80] with great success.
In particular, for very late times t t AH , the mass accretion follows, where F is commonly numerically found to be of order O(1). In the previous equation, the cosmological expansion is neglected, and also, it assumes a quasi-stationary flow onto a PBH. The inaccuracies of these assumptions are absorbed in F. We should also mention that this formula assumes that PBHs are at rest relative to the FLRW background or have negligible peculiar velocities, which is a good approximation in our case. We refer the reader to [142] for a discussion about the scenario with relative velocities. Taking into account the condition of the apparent horizon in spherical symmetry R PBH = 2M PBH , the previous equation is solved analytically as: where M a is the initial mass when the asymptotic approximation is used at the time t a . The value of F has been found to weakly depend on the specific shape of the perturbation and its amplitude, but strongly on the equation of state w. Specifically, F depends on the high non-linear process of accretion and the counterpart of pressure gradients. F can be obtained numerically by fitting the evolution of M PBH (t) (at sufficiently late times after the formation of the apparent horizon, when it is valid). Using that, the PBH mass is estimated as the asymptotic mass value at t → ∞, i.e., The final PBH mass will be dependent on the specific profile of the fluctuation, its initial amplitude δ m , and w. As was already shown in the past, in the critical regime where δ m is very close to the critical value δ c , the black hole mass follows a scaling law [63][64][65],  Table 1 in [67].
In Equation (58), the constant K depends on the specific shape of the profile considered and M H ≡ 1/2H(t m ) is the Hubble mass, which is calculated at the time t m . Specifically, M H in terms of the initial values of the perturbation can be written as: which will depend on the initial length scale r m and the equation of state w (we can consider gauge values a 0 = 1 and t 0 = 1). In the case of using ther coordinate, we just need to replace Λ → Λe ζ(ΛR H (t 0 )) in Equation (59). The scaling law starts to deviate at (δ m − δ c ) 2 × 10 −2 (for w = 1/3), as was indicated in [94]. It was shown explicitly in [80] that for large δ m beyond the critical case, the scaling law seems to deviate by order O(10%) (at least for w = 1/3), but it can be smaller or larger depending on the profile considered. The deviation from the critical regime is expected to be also found for other w's, but the quantification still is missing in the literature. As pointed out in [147], it is important to consider the scaling law behaviour Equation (58) for the predictions of the PBH mass spectra.

Numerical Techniques
Most of the times, Partial Differential Equations (PDEs) cannot be solved analytically, and numerical methods are needed. This is especially common in general relativity and is the case we have here. Currently, several methods and tools can be employed in numerical relativity [148][149][150][151][152][153][154][155][156]. In the case of PBH formation, there is the extra difficulty that the gravitational collapse happens on an FLRW background, so not on asymptotically flat spacetime.
Different authors have employed different numerical procedures to simulate PBH formation from the collapse of curvature fluctuations under spherical symmetry, as we mentioned in Section 1. From a historical perspective, several works have used the standard and very well-known finite differences approach to solve Hernandez-Misner-Sharp equations in the comoving gauge [62,63,65,78,79,82,94,122,157]. Alternatively, simulations were also performed at earlier times using a different gauge, in particular the mean curvature gauge [77]. The numerical results of both procedures with different gauges were verified to be consistent with each other in [76] at leading order in gradient expansion. A new numerical method came up recently using pseudospectral methods [80] to solve Misner-Sharp equations together with the use of an excision technique to obtain the final mass of the black hole.
In this review, we describe the pseudospectral collocation technique used in [80]. The basics of the numerical code are available here [158] (currently, the only one in the literature publicly available), and therefore, the reader can reproduce or extend some of the results presented in this review.
In what follows, we outline how to use the basics of the pseudospectral technique applied to our problem, which allows solving partially algebraically the differential equations describing the gravitational collapse. See also [177,178] for much more details and mathematical foundations. In particular, we use the Chebyshev collocation method (but there are also other formulations).
Let us consider a function f (x) and fit it with N cheb Chebyshev polynomials (any other orthonormal function could be used as well). The function f (x) can be composed by a sum of Chebyshev polynomials as, where T k (x) are the Chebyshev polynomials of order k. The coefficients c k with k = 0, 1, . . . , N cheb are then obtained by solving The points x k , which satisfy T k (x k ) = 0, are called Chebyshev collocation points. Therefore, wherec k = 2 if k = 0, N andc k = 1 in the other cases. The functions L k are the Lagrange interpolation polynomials.
The derivate of p-order of the function will be given by, We define the Chebyshev differentiation matrix as The components of the matrix are given by: To improve round-off errors in the numerical computations [178], the following identity of the diagonal components of the matrix D is used: In spectral methods, the error decays exponentially with N cheb . This represents a significant improvement in comparison with finite differences, where the error decays algebraically as 1/N v , where v > 0 and N is the number of points of the grid. The exponential convergence in spectral methods is obtained thanks to the fact that the derivate at a given point is computed taking into account all the other points of the grid, instead of finite differences, which only consider the neighbours.
The radial domain Ω in our case corresponds to Ω = [r min , r max ] where r max = N H R H (t 0 ) and r min = 0. N H is the number of initial cosmological horizons where we put the final point of the grid, for which it usually is enough to take N H ∼ 10 2 . The domain of the Chebyshev polynomials is [−1, 1]; therefore, we need to make a mapping between the physical domain and the spectral one. Different options are possible, but a linear mapping is commonly used: x k are the new Chebyshev points rescaled to our physical domain Ω. Furthermore, the Chebyshev matrix can be rescaled using the chain rule:D The imposition of boundary conditions in spectral methods is straightforward. In the case of the Dirichlet boundary condition at x k , such that f (x = x k ) = u D,bc , we just need to satisfy f N cheb (x = x k ) = u D,bc .
In the case of the Neumann boundary condition such that f (1) The numerical stability will depend on the values chosen for dt and N cheb . We do not have complete "freedom" to choose those values. If the number of points of the grid is increased, it will require a reduction of the time step dt, in order to avoid instabilities during the evolution. The Courant-Friedrichs-Lax (CFL) condition for hyperbolic PDEs already indicates that the maximum allowed time step should fulfil dt ∝ 1/N 2 .
There are situations where pressure gradients are substantially large during the numerical evolution. This, for instance, can happen when the curvature fluctuation has very sharp profiles or when w becomes large. In this situation, it is necessary to increase the numerical accuracy of the simulation. One possibility is to increase the number of points of the grid N cheb , but of course, this would imply reducing the time step dt to ensure stability. The other one (which is actually better and more clever) is to use a composite Chebyshev grid: split the full domain into several Chebyshev grids depending on the density of points needed in specific regions of the domain. In particular, the full domain Ω is split into M subdomains as Ω l = [r l , r l+1 ] with l = 0, 1 . . . , M. A mapping between the spectral and physical domain for each is required. We choose a linear mapping again as, wherex k,l are the new Chebyshev points re-scaled to the subdomain Ω l . The Chebyshev differentiation matrix is re-scaled again using the chain rule: The numerical evolution is performed independently in each subdomain, i.e., the spatial derivative is computed using the Chebyshev differentiation matrixD l associated with each subdomain, and the time integration is performed applying the Runge-Kutta 4 method on each subdomain.
To evolve across the different subdomains Ω l , it is necessary to impose boundary conditions. The approach followed is the one of [179], which is based on performing an analysis of the characteristics of the field variables. The prescription is that the incoming fields' derivative is replaced by the time derivative of the outgoing fields from the neighbouring domain at the boundaries. It was checked in [84] that all the fields of the Misner-Sharp equations are incoming except for the density field. Therefore, the boundary conditions that must be applied between each subdomain are: U(t, r l+1,i ) =U(t, r l, f ), where r l, f and r l+1,i mean the last grid point in the subdomain Ω l and the first grid point in the subdomain Ω l+1 , respectively. A diagram of the imposition of the boundary conditions between the subdomains can be found in Figure 4. The use of multigrids in spectral methods has a similar idea to the use of an Adaptive Mesh Refinement (AMR) in finite differences, to increase the resolution of the grid in some regions. Actually, the use of AMR was essential in [94] to verify the scaling law behaviour for the PBH mass in the critical regime. The procedure of the AMR is based on splitting the cell grids on a particular desired region and making subsequent subdivisions of the cells. A summarised and basic procedure is the following: (1) During the evolution, individual grid cells are tagged for refinement using some criteria (high-density regions have to be highly resolved, for instance). (2) All tagged cells are then refined, which means that a finer grid (inside each initial cell) is overlaid on the big one. (3) Finally, after refinement has been performed, individual grid patches on a single fixed level of refinement are evolved in time using the time-discretisation scheme, implemented from the differential equations of the problem. If the level of refinement implemented in a cell is greater than what is needed, the refinement of the cell can be reduced. The disadvantage as in spectral methods is that the stability of the method is modified, and a smaller time step will be required. Therefore, the AMR (and in general, any non-uniform grid procedure) allows resolving problems that are intractable due to the high amount of resolution needed for some regions of the domain.
A known weak point of spectral methods is when they have to handle discontinuous solutions, in particular with the presence of shocks, since the derivative at each point is computed globally, and therefore, the discontinuity is quickly propagated to the other parts of the grids. In this situation, it is better to use a finite differences approach, although usually, a modification of the numerical implementation is still needed to be able to capture shocks, such as the introduction of artificial viscosity.
However, as mentioned in Section 1, in our case, shocks are not formed in the simulations, since we should consider initial conditions at super-horizon scales [94].

Numerical Procedure
In all the simulations, we fixed the value of w, t 0 = 1 and a 0 = 1. From that, we automatically obtained H 0 = α, R H (t 0 ) = 1/H 0 , and recall that α = 2/3(1 + w). For the initial length scale of the perpetuation, r m = 10R H (t 0 ) (Λ = 10) is usually considered, which is sufficient to ensure the long wavelength approximation, as was indicated in [79].
First, we focus on describing the procedure to obtain the threshold δ c , and we leave the method to obtain the final mass M PBH,f for the following subsection. To obtain the threshold, a bisection method was used, which compares different regimes of δ m until finding the range in which the formation of the apparent horizon appears. The threshold δ c corresponds to the midpoint of this range.
An example of a successful procedure is the following: • Create the different grid pointsx k,l for the different subdomains Ω l with a given number of points N cheb,l , and compute the different Chebyshev differentiation matrices D l . Introduce also the initial time step dt 0 , which usually dt 0 ≤ 10 −3 is enough to ensure stability; • Specify a lower and an upper limit in δ m of the bisection. The maximum value δ m,max = f (w) can be chosen for the upper limit, and for the lower bound, δ m,min ≈ w (which corresponds to the threshold estimation using a Jeans length criteria). To reduce the computational time, a closer domain to δ c can be chosen if such a domain is known; • A value δ m is taken once the bisection starts. From that, the curvature profiles K(r) or ζ(r) can be built, and the location of the compaction function peak r m can be found,r m , using Equations (51) and (53) respectively. Subsequently, compute the perturbations of the hydrodynamical variables following Equations (40) and (41); • Integrate Equations (29) and (32) at each time step dt using the fourth-order Runge-Kutta method, imposing at each internal time step the boundary conditions at r = 0 and r = r f . Include the internal boundary conditions between the different Chebyshev grids Equation (72); • To distinguish numerically the formation of an apparent horizon or not, we used the peak value of the compaction function in time. When the peak approaches C max ≈ 1, this implies the formation of the apparent horizon 3 . In this case, the value of δ m chosen would imply the formation of a black hole, i.e., δ m,yes , and therefore, a lower value δ m closer to δ c must be found modifying the bound of the bisection in such way that δ c ∈ [δ c,min , δ m,yes ]. In the opposite case, when C max decreases continuously for some time, this will imply the dispersion of the perturbation, and the bisection bound should be modified as δ c ∈ [δ m,no , δ c,max ]; • The previous steps are iterated until the difference between δ m,yes and δ m,no is lower than the desired resolution, i.e., δ m,yes − δ m,no δ(δ c ). From that, the threshold value is defined as δ c = (δ m,yes + δ m,no )/2 ± δ(δ c ). When δ m approaches the critical value, in some cases, the grid resolution initially set up is not enough to follow the evolution; if so, then δ m is shifted according to δ(δ c ) to avoid being close to δ c . Another alternative would be the dynamical inclusion of more subgrids in the desired region where gradients start to be large.
During the time evolution, a conformal time step dt = dt 0 (t/t 0 ) α is used, which substantially improves the performance of the simulation. The two-norm of the Hamiltonian constraint Equation (26) is computed at each time step. It is expected that the constraint should remain constant from t = t 0 if the Einstein equations are correctly solved. In particular, where the label l refers to the sum of the different Chebyshev grids and k to the different points of the grid.

Numerical Estimation of the Final PBH Mass
To obtain the final mass of the PBH M PBH,f , we used an excision technique, and we describe it with detail in this section. Saying that, other alternatives in the literature have used null coordinates with the Hernandez-Misner equations [63,71,73] to obtain the final mass. We refer the reader to these references, but the primary approach followed can be summarised as follows: (i) The simulation was set up with some initial data using the Misner-Sharp equations. (ii) Starting the evolution, an outgoing null radial ray is sent from the centre to the outer point of the grid, and the values of the field variables are stored once the ray crosses each grid point. (iii) The Misner-Hernandez equations are evolved using the initial conditions obtained in the previous step, which are basically the values stored.
Let us focus now on excision. The main idea of an excision technique is based on the fact that nothing inside the event horizon can affect the physics outside. Therefore, the aim is to dynamically remove the segment of the computational domain inside the apparent horizon, to avoid large gradients appearing and breaking down the simulation.
There are different procedures to implement excision, but the one we used here is the following: We define two parameters, ∆r and dr (we always consider ∆r > dr), where ∆r is the distance separation between the apparent horizon and the excision surface that we establish after each redefinition of the excision surface, and dr is the maximum allowed displacement of the apparent horizon before we create a new excision boundary. At each time step, we locate the position of the apparent horizon, taking into account that 2M(r * , t)/R(r * , t) = 1, where r * (t) is the location of the apparent horizon in time. Since we were handling a discretised grid, we used a cubic spline interpolation to find the correct value of r * (t) (the difference of taking a quadratic spline interpolation is O(0.01%) for the PBH mass values).
In particular, the method used is the following: For a supercritical perturbation, when the peak of the compaction function C max ≈ 1.2 (the final result is not affected by the exact value that we take, as long as we consider C max ≈ O(1)), part of the computation domain is removed, creating an excision surface close to the location of the apparent horizon, whose separation from the excision boundary is basically given by ∆r. Once having done this, we created a dense, but small Chebyshev grid starting from the excised boundary. Another single Chebyshev grid can cover the other domain since large gradients are only developed near the excision surface.
We evolved in time the Chebyshev grids in the usual form independently, and we tracked the apparent horizon motion at each dt. Once the apparent horizon had been moved a distance greater than dr, we created a new excision surface separated a distance ∆r from the apparent horizon. During the evolution, we needed to continuously reduce the values ∆r,dr, since the velocity of the apparent horizon decreases in time. This is especially important for very small δ m − δ c . We decided to reduce ∆r,dr when the Hamiltonian constraint starts to be continually violated. In particular, it was found that considering ∆r ≈ 2dr ≈ O(10 −2 ) worked nicely for our purposes.
There is no physical boundary condition to apply at the excision surface. Still, it was found that freezing (leave it constant) the gradient of ρ (i.e., ρ ) at the excision surface, after each redefinition of the excision boundary, allowed increasing the stability of the method without changing the numerical results.

Numerical Results
In this section, we review some numerical results using the numerical procedure pointed out in the previous section. First of all, we introduce a set of curvature profiles that allows us to span all different possible thresholds for PBH formation. Later on, we will see the dynamical evolution of the collapse for three different cases in terms of the initial amplitude δ m and with a specific profile chosen. Finally, we show the results regarding the final PBH mass in terms of the different profiles, as well as the effect of the accretion from the FLRW background.
The numerical results were obtained using the procedure of Section 5 with a computer, i-7 core, 16 GB of RAM, and without parallelisation procedures. As an example of the efficiency, in the case of w = 1/3 and with a Gaussian profile, to obtain the threshold δ c with a resolution δ m − δ c ≈ 10 −3 , the running time was ∼3-5 min using the bisection method in Section 5.2. To estimate the final mass M PBH,f , we needed to run the excision method in Section 5.3 for ∼3 h for the largest PBHs. These values can differ when changing the profiles or the equation of state w. The typical size of the grid chosen was two or three grids with ∼100 points in each grid. However, more points could be needed for sharp initial profiles. With these configurations, the maximum resolution obtained for the threshold in the case of a Gaussian profile was δ m − δ c ≈ 10 −5 (in the case of w = 1/3). With sharper profiles or larger w, this resolution slightly reduced.

Curvature Fluctuations
The seed of the generation of the curvature fluctuation is the power spectrum, namely P(k), where k is the wavelength. We briefly review the procedure assuming the Gaussian statistics of [116] applied to our purposes. We considered that the amplitudes of the curvature perturbations are statistically Gaussian distributed accordingly to the specific inflationary model considered, and therefore, under this assumption, the rare over-densities leading to PBH formation are, with a very good approximation, spherically symmetric [116].
From the power spectrum P ζ (k), we can relate it to the curvature fluctuation ζ Equation (37) as, where ζ k is the Fourier mode k of the field ζ, Notice that in linear theory and at super-horizon scales, the density contrast δρ/ρ b (r) is related to ζ(r) as, and therefore: δρ From that, we can build the normalised two-point correlation function by, where W(k) is the window function, and the variance of the field ζ is given by, Therefore, the curvature fluctuation ζ(r) assuming peak theory can be built as, where µ is the peak value of ζ.
In conclusion, using Equation (81), we can relate the power spectrum to our curvature fluctuation ζ(r) and, from that, set up the initial conditions following Equations (41) and (42), or alternatively also build K(r) using the transformations of Equation (39). Although it goes beyond the scope of the review, we should mention that one could also consider the effect of the inclusion of non-Gaussianities in Equation (81), which can have important consequences for the PBH abundance estimation. There are vast and extensive works that consider its impact on the thresholds and PBH abundances [103,107,108,[180][181][182][183][184][185][186][187][188][189][190][191][192][193][194][195], including some of them that have used numerical simulations following the basics shown in this review [196,197]. Saying that, it is much more practical to already span a set of curvature profiles instead of considering different power spectrums when we are interested in studying the physical process of PBH formation. The reason is that setting up directly the curvature fluctuation allows us to modulate different shapes and therefore obtain all possible thresholds and estimate their effects on the PBH mass. In any case, some typical power spectrum templates can be found in [60,180].
In the literature, exponential families of profiles have been commonly used to modulate K(r) [77,79,81,82]. Another possibility is a polynomial family [84]. Both can be written in a convenient way as: The profiles of Equations (82) and (83) correspond to polynomial and exponential (for λ = 0) non-centrally peaked (for λ = 0) profiles, respectively. Notice that the profiles depend on the ratio r/r m . The profiles are written in a convenient way with the amplitude δ m , in such away that when computing C(r) using Equation (50), this automatically gives C(r m ) = δ m when evaluating r = r m . On the other hand, the parameter q is a dimensionless parameter, which relates the shape around the peak of the compaction function, which is defined as,

of 60
In terms of ther coordinate, the transformation in Equation (39) can be used, as was done in [198], to obtain, Interestingly, in [83], it was found numerically that profiles with the same q have the same threshold δ c upon a deviation of O(3-4%) in the case 1/3 ≤ w ≤ 1 (and O(1-2%) in the particular case of a radiation-dominated Universe), but we will discuss this in more detail in Section 7, since this result will be used to obtain analytical estimations on δ c . The benefit of Equation (82) in comparison with Equation (83) is that at r → 0, this fulfils the regularity conditions, whereas Equation (83) does not for some small values of q. Actually, the profile Equation (82) can be considered as a basis profile, since it allows spanning all possible threshold values, as we will see in more detail in Section 6.3. The profiles of C(r) for the polynomial profile can be found in Figure 5. When q 1, the peak of the compaction function is sharp, while when q 1, the peak is broad (which would be a homogeneous sphere). The opposite behaviour can be found for the density contrast (bottom panel of Figure 5). In this case, for small q, the density contrast profile is sharper, whereas for large q, the shape becomes broad up to r m (a homogeneous sphere). A profile with two modulating peaks in C was also considered in [87] with the aim to compare the other profiles. It is similar to the one considered in [78], where double-PBH formation was studied. This kind of profile would represent a perturbed region superposed on another one at much larger scales, as discussed in [78]. The compaction function profile is given by, where C b is equal to: the value of the second peak of C can be modulated as, C tt(peak,2) = (1 + q 1 )(r j + r m2 ) 2 δ 1 q 1 r 2 m1 + (r j + r m2 ) 2 (r j + r m2 /r m1 ) 2q 1 Different profiles of the three families are plotted in Figure 6. In particular, we plot the profiles of Equations (82), (83), and (86), and from them, we also plot the corresponding compaction function from Equation (50), as well as the density contrast Equation (41). The conversion from K(r) to ζ(r) can be made using Equation (39). The values taken for C tt (r) are δ 1 = δ c (q 1 ), q 2 = 3, r m1 = r m2 = 1, r j = 2r m1 , C tt(peak,2) = 0.3, with the corresponding δ 2 obtained from Equation (88) using the previous parameters, q 1 = 5 (violet) and q 1 = 1 (orange).

Non-Linear Behaviour of the Gravitational Collapse
In this subsection, we show the behaviour of the gravitational collapse of the curvature fluctuations for three different regimes. As an example, a Gaussian profile given by Equation (83) with q = 1 (λ = 0) and for different values of δ m and w = 1/3 is considered. However, the same qualitative behaviour can be found for other profiles and other w = 0. In particular, three different regimes can be considered: Figure 7 (supercritical), Figure 8 (subcritical), and Figure 9 (critical), we see the evolution of the variables ρ, Γ, U, and C.
In the case of a supercritical evolution, which is the case of Figure 7, we can observe that the compaction function grows during the evolution. The gravitational collapse increases the energy density on the central regions, and therefore, the mass excess grows. The formation of the two horizons is also obvious, as discussed in Section 4.3. The outer horizon moves outwards, and the inner one moves inwards, faster than the outwards one. The simulation breaks down when the inner horizon approaches r = 0. This is when the excision technique has to be applied to remove the singularity.
In Figure 8 (the subcritical case), the compaction function decreases continuously as the perturbation is diluted away due to the superiority of the pressure gradients against gravity. In general, for late times, the peak C max is pushed outwards due to the pressure gradients. This behaviour can also be seen with the peak value of ρ/ρ b . In the supercritical case, the peak value is continuously increasing. In contrast, in the subcritical case, the peak value starts to decrease at some moment (which means the perturbation starts to disperse). In addition, from Figures 7-9, it can be observed that Γ is not constant during the evolution. As expected, the gradient expansion approximation fails for sufficiently late times t t 0 .  In Figure 7 (supercritical case), we observe that U/Γ decreases quickly in time, and a negative velocity mainly dominates the profile before r m for late times, which indicates the collapse of the fluid. Instead, in Figure 8 (subcritical case), only a small negative value U/Γ is reached for early times, and beyond that, the values are positive, which indicates that the perturbation is dispersing, avoiding the gravitational collapse.
The most relevant behaviour is found in the critical regime, in Figure 9. Here, the fluid is divided into two parts, one going outwards (positive U) and one inwards (negative U), which generates an under-dense region. This under-dense region re-attracts the fluid with the net effect of a compression and rarefaction process, which increases its velocity in time. In this situation, the gradients become very large, and eventually, the simulation breaks down if a sufficient resolution has not been set up for the simulation.
It is interesting to mention that in the presence of secondary peaks in C(r) (not only one isolated one), with an amplitude roughly equal to or less than the first one, it was verified numerically in [197] that those secondary peaks decrease and disperse on the FLRW background, even in the case of supercritical perturbations. A different situation could happen when a secondary peak is sufficiently high. In such a case, the secondary peak could also collapse even if the first peak is not greater than the threshold value to collapse itself. In Figure 10 is shown the correctness of the numerical evolution using the Hamiltonian constraint of Equation (74) at each step. It corresponds to different evolutions of the gravitational collapse, both subcritical and supercritical. In the figure, it can be seen that the constraint starts to be violated for late times for (δ m − δ c ) ≈ O(10 −5 ), which gives a bound on the maximal resolution for δ c (in this particular example).
The convergence with the use of spectral methods can be found in Figure 11. It is clear that the convergence is exponential taking into account the fits performed (see the caption of the figure). In particular, the use of a multidomain grid improves the performance and accuracy of the simulations substantially, as is shown in the right panel of Figure 11, where the spectral accuracy is obtained with fewer points (in comparison with the left panel). Making the first grid Ω 1 from r = 0 up to r = 2r m allows capturing the region where large gradients are developed. Then, for the second grid, Ω 2 goes from r = 2r m up to r = r f , and only a few points are needed.
An interesting behaviour is found for t AH when w is changed, as noticed in [84]. Although the pressure gradients work against the gravitational collapse, since they are also a form of gravitational energy, it will mostly be favoured when the collapse is active [84]. Consequently, this implies a shorter time of formation for larger w, as can be observed in Figure 12.

Thresholds for PBH Formation
As we have mentioned, the threshold for PBH formation depends on the shape of the curvature fluctuation. Following the numerical procedure in Section 5.2, all the possible thresholds for PBH formation can be obtained using the basis polynomial profile of Equation (82). The results are shown in Figure 13 for different profiles. In the particular case of a radiation fluid, it was found in [79,83] that the range of all possible thresholds is δ c ∈ [0.4, 2/3]. The minimum threshold δ c,min = 0.4 corresponds to the limit case of q → 0 (broad profile in C). Instead, the maximum threshold δ c,max = 2/3 corresponds to q → ∞ (peaked profile in C). Increasing q, the pressure gradients are larger, and therefore, the threshold value is higher.
Remarkably, the relative deviation between the thresholds for the different profiles in the left panel of Figure 13 with the same q is within 2%. This was found for the first time in [83], and it was shown that different profiles with the same q (shape around the peak of the compaction function) have the same thresholds within a small deviation. In [84], simulations of PBH formation in the case of a perfect fluid with w ∈ (0, 1] were performed. The functional behaviour of δ c in terms of q for different w's is roughly similar, as is shown in the left panel of Figure 13 (the case with fixed w = 1/3), but quantitatively, the values are different. It is clear that increasing w, the gradients will be larger, and therefore, the threshold increases. The limits q → 0 and q → ∞ for the general case in terms of w are both amenable to further analysis: • The case q → 0 corresponds to a C that becomes approximately constant over a wide range of scales.
Simulations, in this case, become difficult due to the presence of conical singularities [84]. To solve this issue, in [84], a modified profile of Equation (82) was used to subtract the mass excess beyond r r m and fulfil that at the boundaries, the FLRW background is recovered. In Figure 14 are shown the numerical values of δ c (q → 0) in terms of w, which show a strong dependence on w. As expected, for w → 0, the threshold δ c → 0; • The case of large q → ∞ corresponds to the case of a sharply peaked profile for the compaction function. For this kind of profile, the pressure gradients acting against the gravitational collapse are maximal, and hence, the compaction function should be maximal as well. Analytically, this is shown by Equation (50) that C(r m ) < f (w) (at super-horizon scales). As was proven numerically in [79], δ c,max → f (w) in the case of radiation w = 1/3. The intuition is that the saturation should persist to larger values of w because larger values of w also imply larger pressures that fight against the collapse. Therefore, for w ≥ 1/3, the compaction function of a peaked profile must also saturate the bound f (w). This was verified numerically in [80]. This behaviour is no longer true for w < 1/3, something obvious in the case of dust w = 0, since f (0) = 0.6, but δ c (w = 0) = 0 always for any profile. This is shown in the left panel of Figure 14. It should be taken into account that beyond q > 40, the simulations do not obtain the threshold with enough resolution. This is the reason why the numerical value of the threshold for w = 1/3 does not exactly coincide with f (1/3). One expects that in the asymptotic limit, the threshold for q → ∞ should match with f (1/3), as was pointed out in [79] with the same argument. It is interesting also to notice the difference between considering K(r) and ζ(r) for the threshold and PBH mass. As we have seen in Section 4, there is a non-linear relation between ζ and C, which leads to some differences between the thresholds obtained for both profiles. Consider for instance a Gaussian profile in both cases, i.e., K(r) = A e −(r/r m ) 2 and ζ(r) = µ e −(r/r m ) 2 . Due to the non-linear relation between C and ζ at super-horizon scales, the δ m,ζ (for this example, let us call it δ m,K and δ m,ζ for both profiles) depends non-linearly on µ as, Due to this, the q ζ value in the case of the fluctuation ζ is dependent also on its amplitude µ and w, when for the K(r) in this example, it is just q = 1. A plot of the thresholds δ c for both profiles can be found in Figure 15 for different w's. Although both profiles are Gaussian, they have a different δ c due to the non-linear relation with ζ(r). The threshold for ζ(r) is higher than K(r) since the shape around the compaction function is sharper, which means a larger q, as can be observed in the bottom panel of the subplot in Figure 15. It also has an implication for the PBH mass, as we will see in the next subsection. A Gaussian family of profiles with ζ(r) = µe −(r/r m ) 2p was considered in [106] for the estimation of PBH abundances taking into account the non-linearity of ζ(r) for w = 1/3. There, a range of thresholds was found, 0.442 < δ c < 0.656 for 0.34 p 2.

PBH Masses
Once the apparent horizon is formed, an excision technique following Section 5.3 can be applied to remove the singularity. An example of the evolution of the PBH mass M PBH (t) for the profile of Equation (83) with q = 1, λ = 0, and different amplitudes δ m can be found in Figure 16. The initial time of the curves corresponds to the initial mass of the PBH M PBH,i at the time of formation of the apparent horizon t AH . As can be seen in Figure 16, when the amplitudes δ m increase, the PBH mass will be higher. The amount of time needed until reaching the stationary regime where the PBH mass remains constant could be substantial and therefore computationally expensive. As we have mentioned, to avoid that, Equation (56) can be used to estimate the final mass of the PBH after the accretion process accurately.
In order to find the domain in time where the approximation of Equation (56) is valid, in [80], the ratio of the increment in time of the black hole mass with respect to the Hubble scale Ψ =Ṁ/HM was computed, which is predicted to be Ψ < 1 when the evolution of the PBH mass satisfies this regime. During the application of the excision technique, the Hamiltonian constraint is computed using Equation (74) to check the correctness of the numerical evolution. In the simulations in [80,87], the constraint was fulfilled until very late times, as shown in Figure 17. Nevertheless, it seems that the evolution of the mass is not affected substantially by the violation of the constraint. Interestingly, in the bottom panel of Figure 17, a crossing for different evolutions of Ψ at a given time t * is observed. In Figure 18 is shown the evolution of M PBH (t) in time for different values of q and two different profiles. It is clear that the larger is q, the smaller is the final PBH mass M PBH,f , since the pressure gradients are stronger and prevent accretion from the FLRW background. Notice the differences of M PBH (t) between the two profiles considered in the figure. Even with the same q, the values for the mass are different since the full shape of the profile will contribute to the final mass.
To obtain the parameters t a , M a , and F, a non-linear fit using the function of Equation (56) with the data M PBH (t) coming from the simulations can be performed. The range of values to make the fit are those that fulfil Ψ 0.1. The values of F that were obtained in the literature [80] making the non-linear fit were in the range F ∈ [3.2, 3.8] (for the case of w = 1/3) in terms of the different profiles and amplitudes δ m . Some indicative results regarding the accuracy of the fits performed are: σ max ≈ 10 −2 (variance) s d (M a ) = 10 −4.2 , s d (t a ) = 10 −4.1 , and s d (F) = 10 −3.6 , where s d is the standard deviation of the parameters. This is consistent with what was already stated in [11,143,144]: the numerical coefficient F must be O(1). In particular, in the case of the collapse of the domain walls, F ≈ 3.8 was found for large black holes [146].  Interestingly, in [87], the effect of the accretion (from the FLRW background) to the growth of the PBH until reaching the stationary regime was quantified, namely M PBH,f /M PBH,i , for the case w = 1/3. As was pointed out, the accretion should not be important for small PBHs, i.e., with δ m near δ c . However, the situation differs for relatively large PBHs. This can be precisely seen in Figure 19. The accretion can be significantly important for large δ m .
Sharp profiles, corresponding to large q, have larger pressure gradients, and therefore, the ratio M PBH,f /M PBH,i is smaller, even for large δ m , since the gradients prevent the accretion. For low M PBH,f , the ratio M PBH,f /M PBH,i should be small, as expected [9,86]. When M PBH,f M H , i.e., for PBHs with a higher probability to form, M PBH,f 3M H was obtained. In comparison, a similar result was also found in the case of the collapse of a massless scalar field [141] and in [146,199] from BHs formed from the collapse of domain walls or vacuum bubbles, with a factor M PBH,f /M PBH,i ≤ 2.  Figure 20. In this case, the values of F differ substantially from w = 1/3. In particular, for the set of values w ∈ [0.2, 1/3, 0.5, 0.7], this corresponds to F ≈ [7.2, 3.8, 2.0, 1.37], respectively, which already clearly indicates the substantial accretion for small w. Moreover, the time needed to reach the asymptotic regime of Equation (56) is much higher the smaller w is. For example, in the case of δ m − δ c = 10 −2 for w = 0.2 and w = 0.5, a time of order t/t 0 10 4 and t/t 0 10 3 was needed, respectively, to reach the asymptotic regime. On the other hand, in [86], some analytical approximations were made regarding the maximum size of the PBHs formed from hydrodynamical perturbations at the time of the formation of the apparent horizon t AH ; in particular, it was found that: This analytical result (derived for 0 < w ≤ 1) was verified numerically in the case of a radiation fluid in [87], where an extensive study on the values M PBH,i , R PBH,i , and t AH was also performed for different profiles and amplitudes δ m .
To compute the value of the constant K in front, the scaling law Equation (58) is needed to estimate the final mass of the black hole M PBH, f using perturbations with an amplitude near the critical regime, i.e., δ m → δ c . A range of 10 −3 δ m − δ c 10 −2 is sufficient to estimate it, since the gravitational collapse of the fluctuations happens in the critical regime. Beyond that, the scaling law deviates. In particular, for a radiation fluid, the value γ = 0.357 [66,68] can be taken directly to obtain K from a given δ m − δ c (for other values of w, check Figure 3 and the caption).
In Figure 21 are shown the values of K for the profiles of Equations (82), (83), and (86) as an example of K in terms of the profile. The value of the PBH mass M PBH,f (and therefore, also K) is affected by the shape of the profile beyond the peak of C since the accretion process depends on the full shape of the fluctuation, in comparison with the threshold δ c , as was pointed out in Section 6.2, which mainly depends on the shape around the compaction function peak.  (86) (green). The parameters used for the profile C tt (r) are q 2 = 3, r j = 2r m1 , C tt(peak,2) = 0.3, and δ 2 , obtained from Equation (88) using the previous values with and q 1 = q.
In the case of the profile Equation (82), K tends to ≈ 3.5 for large values of q. The values of K tend to increase as q decreases, as shown Figure 21. Numerically, in [87], it was not possible to obtain the final mass M PBH,f for profiles q 0.5, due to conic singularities, as was already found in [84].
On the other hand, in Figure 22, an explicit example of M PBH, f can be seen for a Gaussian profile (in both K(r) and ζ(r)) in terms of M H for large values of δ m beyond the critical regime δ m 10 −2 with w = 1/3. As can be observed in the subplot of Figure 22, the scaling law deviates for the highest allowed values of δ m to O(15-20%) for both cases. For these particular cases, the maximum allowed mass of the PBH formed is M PBH, f ≈ 3.7M H and M PBH, f ≈ 2.3M H , respectively. It is expected in principle that this deviation from the scaling regime will not significantly modify the estimation of the PBH abundance, due to the unusualness of such perturbations since those PBHs with a higher probability to form have M PBH, f ≈ M H . It is important to also notice that the masses generated with the profile ζ(r) are lower than those with K(r), since although they are both Gaussian, they have q > 1 and q = 1, respectively, and sharper profiles have less mass due to stronger pressure gradients, as we see in Figure 18. Values of M PBH f /M H in terms of (δ m − δ c ) for two Gaussian profiles in K(r) and ζ(r). Dark points correspond to the numerical values for the profile in K(r). The solid red line corresponds to the scaling law behaviour with γ = 0.357, δ c = 0.49774 for K = 6.03 (pointed out in [94]) and the blue solid line with the numerical value obtained in the review with K = 5.91. The absolute value of the relative deviation d with respect to the numerical values and the ones coming from the scaling law are shown in the subplot. The orange vertical line is the value δ max = 2/3. On the other hand, the green dots correspond to the case of the Gaussian profile in ζ. The solid green line is the scaling law with γ = 0.357, δ c = 0.55257, and K = 3.98 (in [106], for this profile, K 4.0 was given). The vertical magenta line corresponds to the case δ max = 2/3 for this profile. In both cases, w = 1/3.

Analytical Estimation for the Threshold of PBH Formation
In this section, we review in detail the current analytical estimations used in the literature for the threshold of PBH formation from curvature fluctuations with an equation of state w. The first analytical estimation came from [74], where using a Jeans length approximation, it was shown that, Although this approximation is simple, it fulfils some properties: (i) at w = 0 (dust), any curvature fluctuation should collapse, i.e., δ c (w = 0) = 0; (ii) when w increases, the pressure gradients are larger, so the threshold should be higher. This approximation was obtained assuming the argument that the size of an over-density (that leads to PBH formation), at the maximum expansion, should be smaller than the particle horizon and larger than the Jean radius (Jeans criterion). In particular, where R PH 3/8πρ max , R J c s R PH . See also Section 2 for a simple derivation. Later on, in [75], the previous estimation of [74] was improved by using a sophisticated and interesting model (a three-zone model) of the collapse of a homogeneous over-dense sphere surrounded by a thin under-dense shell. We suggest the reader check [75] for more details on the model and the assumptions made. The threshold estimation that they obtained was, It was also provided with an upper bound δ + HYK and a lower bound δ − HYK on δ c , to take into account the different ways of applying the relativistic Jeans criterion for this model.
These previous estimations (in particular, Equation (94)) have been extensively used in the literature, but as has already been noticed [74,75], these estimations are not profile dependent and therefore do not take into account the specific shape of the curvature fluctuation.

Shape-Dependent Analytical Estimation of the Threshold for Radiation Fluid
In [83], it was proven that, to a very good approximation, the threshold for the w = 1/3 case only depends on the curvature of the compaction function at its maximum, under the assumption of a central over-dense peak in the density distribution. It was used to build an analytical formula for the threshold with a dependence on the profile considered. The formula was accurate enough and only deviated 2% in comparison with the numerical simulations.
The two main ingredients of the procedure used in [83] were: (i) parametrise the shape around the compaction function peak with a dimensionless parameter q; (ii) use the average of the critical compaction function. As was shown in [83] and as we have already mentioned, to a very good approximation, the threshold only depends on, which is a dimensionless measure of the curvature of C(r) at its maximum. This means that different profiles, but with the same q value, will have the same threshold upon a small deviation of <2%. The next point is to define a "basis" (or fiducial set of curvature profiles) such that, by varying q, this set covers the whole range of interesting thresholds and shapes with q ∈ (0, ∞) while also being regular at r = 0 and having ρ (r = 0, t) = 0. In [83], this basis was given in terms of the exponential functions Equation (83) (with λ = 0). Although it strictly does not satisfy regular conditions for q → 0 at r = 0, this limit was tested numerically with the polynomial basis of Equation (82).
Taking this into account, the averaged critical compaction function integrated from r = 0 up to the peak r = r m is defined as,C Introducing the profile considered into Equation (98), we can solve the integral analytically to obtain: where Γ(x) 4 is the gamma function, Γ(x, y) 5 is the incomplete gamma function, and we have used C c (r m ) ≡ δ c . In the case of a radiation fluid, when q → ∞, it was confirmed numerically that δ c (q → ∞) = f (w = 1/3) = 2/3 [79]. Taking into account this analytical limit into Equation (99), we obtain: The assumption made in [83] was to consider thatC c is a profile-independent quantity, and therefore, C c = 2/5 for any value of q. The assumption was tested extensively and numerically for different profiles and was found to be consistent with numerical simulations up to a small deviation of O(2%) in δ c . Inverting Equation (99) in δ c withC c = 2/5, we finally obtain the analytical threshold formula δ c in terms of the specific profile q, in the case of a radiation fluid, This remarkable result in [83] was the first analytical estimation shown in the literature for the threshold of PBH formation, taking into account the shape of the fluctuation, with enough accuracy to be used for estimations of PBH abundances [95]. A plot of Equation (101) in terms of q can be found in Figure 23.
Interestingly, a study of the thresholds for different power spectrums (leading to perturbations in ζ with the coordinater) was performed in [198] using Equation (101), with the corresponding q value obtained from Equation (85) and connecting with the density contrast of Equation (43).

Generalisation of the Analytical Threshold beyond the Radiation Fluid
Later on, the procedure used in [83] was generalised considering PBH formation within a perfect fluid with w ∈ (0, 1] on an FLRW background, and it was found that the threshold averaged compaction function within a concentric sphere of radius r m (1 − α(w)) ≤ r ≤ r m is, to a very good approximation, 4 The gamma function Γ(x) is defined as Γ(x) = ∞ 0 t z−1 e −t dt. 5 The incomplete gamma function is defined as Γ(x, y) = ∞ y t x−1 e −t dt.
again universal in the case 1/3 ≤ w ≤ 1, where α(w) is the location where we should start to integrate the compaction function, and it should depend on w. This was used to provide a new and improved analytical formula for the threshold that only depends on the normalised second derivative of the compaction function at its maximum and the equation of state w. This formula was proven to be accurate at order 6% compared to numerical simulations and reproducing the same result for the case of radiation w = 1/3. It was also again found that the δ c for different profiles with the same q are equal within a deviation less than 4%.
In [83], the basis profiles used were given in terms of the exponential functions; despite that, since the boundary conditions at r = 0 are violated for q < 1/2, the basis profile of Equation (82) was used instead. This set of profiles satisfies the appropriate boundary and regularity conditions.
The critical compaction function, averaged within a spherical shell extending from radius [1 − α(w)] r m to r m (the peak), is defined to be: where Inserting Equation (82) into with: and: where 2 F 1 is the hypergeometric function 6 .
The key point is that as was shown numerically in [84], the dependence of the averaged critical compaction function on the profile shape is weak enough to be ignored, and therefore universal. Then, Therefore, the analytic expression for the critical threshold value reads as: The functions α(w) andC c (w) were found using two different procedures: (I) a numerical double fit minimisation using the numerical results of the thresholds for different profiles and (II) making a single numerical fit minimisation already using the analytical result that δ c (q → ∞) for w ≥ 1/3. Both procedures gave equivalent results. Therefore, the assumption ofC c (w) being universal was found good enough 6 The hypergeometric function is defined as 2 and verified numerically for w 1/3. This result verified the previous one obtained for a radiation fluid w = 1/3 in [83], where α = 1 andC c = 4. In particular, using the procedure (I), it was found,

A Comparison between Analytical Estimates
In [84], the analytical estimation obtained with Equation (108) was compared with the previous ones in the literature, i.e., Equations (92) and (94). In Figure 24 is shown the comparison. The solid lines in Figure 24 show these approximations; the symbols with error bars show δ c (w) from the numerical simulations of the profiles of Equation (82) with q = 0.015, q = 0.1, q = 1, and q = 30. The dotted and dashed curves, which significantly better describe the result of the numerical simulations, show the result of inserting Equations (109) and (110) or Equations (111) and (112), into Equation (108). In both cases, Equation (108) (as well as the simulations) exceeds even the upper bound claimed by [75] at even lower w as q increases, which are Equations (95) and (96).
As can be observed, the discrepancy of the estimation Equations (92) and (94) for the case of small w is as large as 50% for the case q = 0.1. This is the limit that was considered successful and optimal for the approximations on which Equation (94) is based. This deviation is even bigger than the one found earlier because Equation (94) was compared with a Gaussian profile (Equation (83) with q = 1 and λ = 0). Actually, in the case w < 0.15, the solid black curve in Figure 24 give a satisfactory description of the simulations with q = 1. However, the top-hat profile with q → 0 (top-hat for the profile in the compaction function), which was the one used in the analytic estimations of [75], is quite better approximated by q 1. For q = 0.1, the formula Equation (94) does not describe the simulations particularly well, and the discrepancy at w < 1/3 is even worse when q = 0.015 is taken. It was claimed in [84] that this disagreement seems to suggest that the supposed agreement shown in Figure 3 of [75] is a result of numerical coincidences, with no physics implications. Figure 24. Dependence of threshold δ c on w when the initial profile is given by Equation (83) with q = 30, q = 1, q = 0.1, q = 0.015, and λ = 0 in all the cases (q = 0 would be a homogeneous sphere in the compaction function). The solid lines with dots and error bars show the results of the numerical simulations. The green and blue curves show the minimal and maximal bounds of δ c from [75], in particular Equations (95) and (96), respectively. The magenta line shows the approximation for δ Carr of Equation (92). The black curve shows δ HYK from Equation (94). The other curves show the approximation (Equation (108)) in which δ c depends both on w and q. Finally, the dotted curves use Equations (111) and (112) in Equation (108), whereas the dashed curves use Equations (109) and (110) in Equation (108). Finally, the polynomial basis profile of Equation (82) fulfils the regularity conditions for q → 0, as we have pointed out before, something not fulfilled with the exponential basis Equation (83) (with λ = 0). Precisely, this led in previous literature [79] to inferring that the minimum threshold for PBH formation (in the case of a radiation fluid) must coincide with the estimation of [75], i.e., δ c,min (w = 1/3) = δ HYK (w = 1/3) ≈ 0.41. In [79], the profile leading to the minimum numerical threshold of PBH formation was given by Equation (83) with q = 0.1 and λ = 0, since as was indicated there, the profiles were too sharp to take even smaller values of q and produce a numerical simulation. This actually looks related to the fact that the exponential profile Equation (83) does not satisfy the regularity conditions at r → 0 for q < 1/2. Instead, in the case of using the polynomial profile Equation (82), it allows performing simulations even with smaller values in q, up to q = 0.01, as done in [83]. Therefore, this led to the conclusion in [83] that, for a radiation fluid, the correct range of values for the threshold is 0.4 ≤ δ c ≤ 2/3, and therefore, the minimum value was lower than δ HYK . Indeed, the generalisation of δ c in w already indicates that the minimum threshold should not correspond necessarily to the one of δ HYK .

Other Scenarios of PBH Formation
In this review paper, we mainly focused on the study of PBH formation from the spherical collapse of primordial density perturbations, which is currently the case most studied in the literature and with more results and applications [39]. Despite that, many other mechanisms and scenarios could have led to black holes formed in the very early Universe [12]. In this last section, we try to briefly enumerate other scenarios beyond the ones considered in this review.
The case w = 0 is very special as the threshold of PBH formation, for an infinitely long matter-dominated era; it is zero if non-spherical effects are not considered. However, as was pointed out in [200], the dust phase could induce an angular momentum on the collapsing perturbation, therefore effectively increasing the threshold value for the black hole formation, making it non-vanishing. Thus, the study of abundances, where the Universe is matter-dominated for a finite time, differs greatly from the case of radiation. Some works have considered this scenario with scalar field types [201][202][203][204][205][206], specifically taking into account the oscillatory behaviour of the scalar field during the reheating. Recently, in [207], simulations were performed in full GR using scalar fields to modulate the perturbation and the expanding dust background of the Universe, showing that the PBH mass grows beyond the self-similar limit M PBH ∝ H −1 , at least initially. The dust collapse has also been studied in a more oriented perfect fluid scenario [208]. The PBH formation has also been considered in the context of the Lemaître-Tolman-Bondi (LTB) model [209,210]. In [211], the effect of inhomogeneity on PBH formation in the matter-dominated era and its consequences on the PBH production were studied. On the other hand, in [200,212], the non-spherical effects and spins for PBH formation in the case of a dust-dominated Universe were studied, showing this could have a significant effect.
A non-spherical effect, in principle, could dramatically change the threshold value and PBH mass [213]. Some analytical estimations taking into account an ellipsoidal collapse were performed in the past, showing that the effect could be important [209,214]. However, in [215], for the first time in the literature, the first numerical simulation in 3 + 1 spacetime (beyond spherical symmetry) for PBH formation from the collapse of curvature fluctuations was performed, with a small (perturbative) non-sphericity initial condition (with a radiation fluid). The main conclusion was that the effect of ellipticity on the threshold for PBH formation is negligible in the standard scenario of PBH formation in the radiation-dominated Universe. This situation could dramatically change if an equation of state different from radiation is considered, as indicated in [215]. In particular, for a matter-dominated Universe, the effect of any non-sphericity could be important [200,208,212]. Another possibility is the inclusion of angular momentum to analyse the spin generation to PBHs; some analytical estimations have been performed [216], but still, numerical simulations to test those arguments are necessary. Recently, PBH formation with an anisotropic perfect fluid has also been considered analytically [217].
The usual scenario of the collapse of cosmological density perturbations considers a constant equation of state w (mainly w = 1/3 in the literature). However, the Universe has a thermal history where the equation of state is not necessarily constant, i.e., w(t). A particular (and very popular) example applied to the PBH scenario is the QCD epoch [218][219][220][221]. During the QCD phase transition, the equation of state w could have been softened. Looking at the result of the simulations in Figure 24, it is already clear that the threshold δ c will be reduced, and therefore, this makes the formation probability of PBHs exponentially more likely. Several works have addressed this scenario with progressive refinements [218,[222][223][224][225][226][227], specifically with the new calculations regarding lattice QCD [228,229]. However, numerical simulations are still needed to fully verify the effect on how much the threshold, PBH mass, and mass function modifying the equation of state in time could change.
A bit more different mechanism of PBH formation comes from the collapse of Q-balls and oscillons, which are features of supersymmetric extensions of the standard model [230]. Specific realisations of those ideas have been made through theory-motivated scalar field potentials, e.g., the axion-monodromy potential [231]. The mechanism for PBH formation of all these cases is naively similar: small number densities of defects lead to large fluctuations relative to the background density. These fluctuations become gravitationally bound and collapse to form black holes once the relic density has come to dominate. Finally, the relics decay due to some instabilities. Precisely, it has been argued that some solitonic type solutions such as Q-balls and oscillons could produce a significant fraction of dark matter in the form of PBHs [232][233][234][235][236], without relying on any spectrum of density perturbations. Some recent numerical works on oscillons and Q-balls have obtained results in this direction [237][238][239].
On the other hand, scalar force instability can lead to the growth of structures and the formation of halos of interacting particles (even in a radiation-dominated Universe) [240][241][242][243][244]. Actually, in [243], it was theorised that these structures could have collapsed and formed PBHs, but later on, it was realised that they should remain as virialised dark matter [244]. Precisely, in [245], it was shown that the radiative cooling due to the scalar radiation allows those halos to collapse into black holes. This was shown with a simple model with fermions interacting via the Yukawa interaction. The general relativistic formulation of the interaction between a Fermi gas and a scalar field in cosmology was carefully studied recently in detail in [246].
All these mechanisms could induce small or large spins to the PBHs in terms of their formation evolution and the process of radiative cooling, as was pointed out in [213]. Another proposal is the one in [247], where the dilution of heavy quarks produced during inflation, which become confined by QCD flux tubes after horizon re-entry, could have led to PBH production.
Another mechanism is through the collapse of domain walls. Domain walls are topological defects that may form when a discrete symmetry is spontaneously broken in the very early Universe [248,249]. Several works have addressed this scenario and have shown the possibility of forming PBHs, with astrophysical consequences [146,[250][251][252][253]. Another possibility comes from the collapse of cosmic strings [254][255][256], which are 1 + 1 topological defects that are predicted beyond the Standard Model and could have led to the production of PBHs [257][258][259][260][261][262][263].
Other results have also been pointed out in the context of modified gravity, in particular in [275], where the PBH production in G-inflation was studied, as well as in [276], where this was studied with the inclusion of a Gauss-Bonnet term. In [277], the effects on the threshold in Eddington-inspired Born-Infeld gravity were considered.
Finally, a new mechanism appeared in [278], where PBHs could have been formed from the collapse of primordial CDM isocurvature fluctuations. Related to that, a new work has numerically shown the formation of PBHs [279] from isocurvature fluctuations generated by a massless scalar field.

Conclusions
In this review, we focused on explaining in detail the process of PBH formation from the collapse of curvature fluctuations under spherical symmetry on an FLRW background and its numerical implementation through a particular numerical technique. Specifically, we saw that numerical simulations are essential to capture the correct and full-non-linear behaviour of the collapse and formation of PBHs. This leads to a better and accurate determination of the threshold of PBH formation and their masses compared to analytical estimations, something essential for the precise analysis of PBH abundances. In this review, we saw a simple and efficient numerical procedure based on pseudospectral methods that allows us to fully solve the non-linear gravitational collapse of those perturbations and obtain the needed values of the threshold and PBH mass in the range of practical interest.
Two aspects in the literature have been remarkably important for the development of the field: (I) the results found from numerical simulations about the critical scaling regime for the PBH mass and (II) the analytical estimations for the threshold of PBH formation. Both allow us to estimate the abundances of PBHs and constrain them. In particular, we saw that the threshold and PBH mass are both dependent on the specific shape of the curvature fluctuation considered, and it is crucial to take this into account for the precise analysis of PBH abundances instead of using a constant threshold or PBH mass value. Moreover, we saw that the gravitational collapse is highly dependent on the equation of state of the fluid w. The larger w is, the smaller is the mass of the PBHs and the accretion from the FLRW background, but the higher the threshold is.
It is interesting to mention that the compaction function, firstly used in [77] years ago, has become currently essential, not only for the definition of the threshold of PBH formation, but also to obtain an analytical estimation of it, using its average.
There is still much effort to be made in other directions and different scenarios of PBH formation to consider, as the ones already pointed out in Section 8. Future theoretical developments and the implementation of numerical simulations to test and check those scenarios could lead to new possibilities for PBHs to account for dark matter and update the abundance estimations. In particular, it would be interesting to compare the results already pointed out in this review in other scenarios.
We hope that in the next few years, several scientific results in this and other directions will help to fully understand and constrain the PBH scenario and its consequences in our Universe, especially taking into account future gravitational wave observations. Funding: AE is currently supported by "bourse de post-doctorat" from Université Libre de Bruxelles (ULB), and was supported by contract PID2019-106515GB-I00/AEI/10.13039/501100011033 (Spanish Ministry of Science and Innovation).
Acknowledgments: I want to thank Vicente Atal, Guillem Domenech, Cristiano Germani, Claudia Gonzalez-Boquera and Marc Oncins for extremely useful comments about the manuscript. I would also like to thank discussions on the PBH topic to Vicente Atal, Sebastien Clesse, Guillem Domenech, Jaume Garriga, Cristiano Germani, Carsten Gundlach, Ilia Musco, Marc Oncins, Javier G. Subils, Ravi K. Sheth, Antonio Enea Romano, Yuichiro Tada, Shuichiro Yokoyama and Chulmoon-Yoo. Finally, I also want to thank the anonymous referees for the constructive and valuable comments that have helped to improve this work.