Mechanisms of Producing Primordial Black Holes and Their Evolution

: Primordial black holes have become a highly intriguing and captivating ﬁeld of study in cosmology due to their potential theoretical and observational signiﬁcance. This review delves into a variety of mechanisms that could give rise to PBHs and explores various methods for examining their evolution through mass accretion


Introduction
The idea of primordial black hole (PBH) formation in the early universe was suggested in 1967 by Zeldovich and Novikov [1] and also by Hawking and Carr [2][3][4]. PBHs have been a subject of interest for fifty years. In particular, only PBHs, by construction, can be light enough for the so-called Hawking radiation to be essential to change the mass of a black hole. Today, the primordial origin of some discovered black holes (quasars at small z [5,6], BHs of intermediate masses detected by gravitational-wave observatories [7]) is hotly discussed [8,9]. It is also worth noting that black holes of primordial origin can be dark matter candidates [10]. At the moment, constraints are imposed on a wide range of masses of a PBH as a dark matter candidate [11]. It has also been shown [12][13][14] that primordial black hole formation mechanisms are able to form clusters of PBHs and then, some constraints imposed on black holes as dark matter candidates should be reconsidered [13,15,16].
The next aspect to be discussed is the accretion process. Accretion is the process in which a black hole can capture particles from a nearby source of fluid, and there is an increase in mass as well as angular momentum of the accreting object, see, e.g., [32]. There are different approaches to study this process and some of them are reviewed here. The most common assumption to calculate accretion rate is a spherical symmetry. This assumption is used in papers such as [33][34][35][36]. There have also been numerical estimations with relativistic hydrodynamics [37] and their results indicate that mass growth through radial accretion might be significant. It is also necessary to mention accretion disk models. The first realistic model of accretion disks around black holes was formulated in 1973 by Shakura and Sunyaev [38]. This approach was based on the consideration of matter rotating in circular Keplerian orbits around the compact object, losing angular momentum because of the friction between adjacent layers and spiralling inwards. In this process, gravitational energy was released, the kinetic energy of the plasma increased and the disk heated up, emitting thermal energy into space. For a relativistic analysis, see, e.g., [39]. Various models of accretion disks are considered in [40].
The work plan is organised as follows: In Section 2, we consider the mechanisms of black hole generation, including the traditional Zel'dovich and Novikov mechanism and various types of phase transitions. In Section 3, we discuss the mechanisms of mass accretion. We start with the conventional Bondi solution and one of its applications. Next, we consider a common marginal estimation of mass accretion, i.e., the Eddington limit. The section ends with a discussion of approaches based on the general theory of relativity. In Section 4, we give a brief description of the models considered in this review.

Initial Density Inhomogeneities
In this section, a brief review of [3] is presented. This mechanism is based on the collapse of primordial inhomogeneities in a hot plasma and arises within the framework of the standard Big Bang cosmology. Consider a region of the universe of radius R. It has a potential energy of self-gravitation ans also a kinetic energy of expansion where ρ is the energy density. At the radiation-dominated universe, the pressure and energy density are proportional to R −4 since p = 1 3 ρ. If the density is high enough in a region, gravitational forces may overcome the kinetic energy of expansion and pressure forces. As a result, the cosmological expansion is halted in that region. To overcome pressure forces, the gravitational energy must be greater than the internal energy, which, for p = 1 3 ρ, is U ∼ ρR 3 , so the necessary condition for collapse is ρR 2 >∼ 1.
A drawback of this model is that the mass spectrum in this model is close to monochromatic, so it is impossible to explain the existence of black holes of various masses. It is also impossible to produce clusters of PBHs within this model. This model does not contain a mechanism for the production of compact clusters, and consequently, it is not applicable to explain the merger rates observed by LIGO/VIRGO [9]. It must be stated that this is the first and "conventional" mechanism for black hole production. It does not require any assumption beyond standard Big Bang physics.

First-Order Phase Transitions
Let us consider the mechanism for producing PBHs by means of scalar field dynamics. First-order phase transition as a mechanism to produce black holes was proposed by Khlopov et al. [41]. For a realization of this mechanism, it is necessary that the field potential has at least two minima, and one of them must be false. Let the field be in a false vacuum φ 1 at the initial moment of time, and let us denote the true vacuum by φ 0 (Figure 1). As a result of quantum tunnelling, in one region of space, the field will have a value of φ 1 and in another region of space, φ 0 . These regions are called bubbles. In this formulation, the free energy of a bubble consists of two parts, the volume and the surface energy. We denote the surface energy density µ and the difference of potential values at minima ∆V = E(φ 0 ) − E(φ 1 ). In this case, the free energy of a bubble with radius R and the surface energy density µ can be written as Obviously, the dependence (5) has a maximum at the point R cr = 2µ/∆V, after which it becomes energetically advantageous for the bubble to expand infinitely. Then, the expansion of bubbles of true vacuum in the region of false vacuum becomes possible, while the potential energy of the false vacuum is converted into a kinetic energy of the walls, which leads to an ultrarelativistic speed of expansion in short time. When a pair of bubbles of true vacuum collide, a new bubble of false vacuum can arise in the region of false vacuum. If its size is smaller than its gravitational radius r g = 2GM (G is a gravitational constant), it becomes a black hole for a distant observer. However, if the bubble shrinks to a size comparable to the wall thickness d ∼ δ > r g , the collapse in the PBH will not occur, since the bubble will oscillate, lose energy, and finally decay by tunnelling.
For some detailed analysis how PBHs could be formed within this mechanism see, e.g., [42,43] and the references therein.
The advantage of this mechanism is a broad mass spectrum [44]. However, one must pay attention to its drawbacks, e.g., it leads to strong inhomogeneities [45], and there is no possibility to form clusters of PBHs. This model also does not contain a mechanism for compact cluster production.

Second-Order Phase Transitions
In contrast to first-order phase transitions, in this mechanism, we do not face collisions of growing domain walls. The closed domain walls shrink and eventually collapse to form black holes. For a realization of this mechanism, it is necessary that the field potential possesses at least two vacuums of the same energy.
The idea is to generate domain walls which could collapse into a PBH after crossing the cosmological horizon. There are two ways to implement this. The first case is based on spontaneous symmetry breaking. The second case is based on the inflationary stage quantum fluctuations.
Let us qualitatively review the first case. Its key feature is a spontaneous symmetry breaking with a consequent change of potential. After the temperature drops below a certain level, the potential acquires several equal-energy minima. This can lead to the generation of shrinking domain walls; see, e.g., [20], where the case of spontaneous symmetry breaking is considered.
It is necessary to mention a model elaborated by A. Dolgov in 1993 [46] with a lognormal mass spectrum. The proposed model was the first where the cosmological inflation and Affleck-Dine baryogenesis [47] were applied to the PBH formation. Mass distribution of PBHs was given by the expression: where γ is a dimensionless constant, according to [46]. It was predicted that M 0 ∼ 10M [48] and this result was in good agreement with observations [11,49]. This is the only mass spectrum of PBHs which was tested by observations and found in good agreement with the model. The second mechanism is based on multiple quantum fluctuations generated during the inflationary stage. There is no symmetry breaking in this case, but such mechanism also generates closed domain walls.
Let the field be near the maximum of the potential at the initial moment of time. Classical motion of the field would lead to the field rolling into one of the minima; however, as a result of multiple quantum fluctuations at the inflationary stage (in a field with frozen classical motion), the field can gradually flip over the maximum. As a result, it rolls into different minima in different regions of space. These distinct regions are connected by closed domain walls. The shrinking of such domain walls can lead to the formation of PBHs. One of the first studies dedicated to the elaboration of this mechanism of PBH was performed in [20].
In this case, there is a very broad mass spectrum of PBHs, see, e.g., [13,50]. Nonetheless, the mass spectrum in that case strongly depends on the initial field value, so fine-tuning is required, although fine-tuning is a common problem for inflation theories [51].
There are common features in these two approaches. The characteristic scale of nonvanishing fluctuation at the inflationary stage is H −1 inf . If it is formed at the moment t during inflation so by the end of the inflation-dominated phase, it would be e N inf −H inf t times bigger. Its further evolution depends on the relation between two timescales: t σ = 1/2πGσ, where σ is the surface energy density of the wall and G is, again, a gravitational constant, and the moment at which the wall crosses the cosmological horizon, t H .
If t σ t H , the wall is called subcritical and a black hole forms much earlier, before the wall can become dominant in the universe. However, in the case where t σ t H , the wall is called supercritical, and there is a wormhole formation as a way to "export" the problem of wall domination to a baby universe; see [13,[52][53][54] and references therein for a detailed analysis.
The advantage of this mechanism is the same as mentioned above, a broad mass spectrum, and it could also be coupled with external processes, e.g., baryogenesis [46]. In contrast with the case considered above, there is no explicit mechanism to form compact clusters of black holes.

PBH Production in f (R)-Gravity
The idea of the proposed mechanism is based on the known possibility of formation of domain walls during cosmological inflation followed by their collapse into primordial black holes [13,20]. The formation of such domain walls requires a scalar field with a nontrivial potential containing several vacuums. It is this effective scalar field that arises in multidimensional f (R)-models in the Einstein frame [25,55,56]. This field controls the size of compact additional space, and its different vacuums correspond to different universes.
The model considered contains quadratic and tensor corrections to the scalar curvature: The multidimensional space is represented as a direct product M = M 4 × M n , where M 4 is a four-dimensional space, M n is a compact extra space, which is chosen in the form of a multidimensional sphere with n dimensions: One can obtain the scalar curvature: where R, R 4 , R n are the scalar curvatures for M, M 4 , and M n , respectively. As shown in [25], at the limit of the effective field theory: one obtain an effective field theory where the scalar curvature of extra space is considered a scalar field.
The effective potential V(φ) happens to be appropriate for the domain wall formation. The kinetic term and the potential are given by the expressions [56]: The right minimum of potential (12) (Figure 2) corresponds to the observable universe [55,56], while the left minimum corresponds to the case of macroscopic extra dimensions, but due to the nontrivial kinetic term (11), it is impossible to reach the left minimum of the potential in finite time.
One can simplify the Lagrangian (10) by substituting and obtain a differential equation for the field evolution in a simple form. As mentioned above, due to nontrivial kinetic term (11), field φ and consequently ψ are not able to reach the left minimum, so the lower limit in (13) is set to be φ 0 φ min , where φ min corresponds to the right minimum of the potential. It is assumed that cosmological inflation is an external process, so there are inflationary constraints on the presented mechanism. As a result of repeated quantum fluctuations during cosmological inflation, the field ψ (φ) can be flipped from the rolling down to the right minimum to the region of rolling down to the left minimum in some area of the inflationary universe [20]. During inflation, the ψ field is frozen near the potential maximum. After the end of inflation, the field tends to one minimum inside the bubble and another minimum outside of it. Increasing the energy density gradually forms the domain wall around the bubble.  Figure 2. Graphs of potential (12) and kinetic terms (11) with the parameters chosen in [24]: n = 6, c 1 = −8000, c 2 = −5000, a 2 = −500.
It must be stated that f (R)-gravity adds a new possibility for the formation of black holes without any involvement of matter fields. It is also worth to mention that Planck collaboration data are well described by models of f (R)-gravity, e.g., [56,57], although the mass spectrum in this particular theory and possible observable properties are yet to be studied.
It also must be stated that since this model is based on the generation of domain walls via quantum fluctuations during the inflationary era, it has the same drawbacks and advantages as stated in the previous section. In that case, domain walls form due to quantum fluctuations during inflation, we also have a broad mass spectrum, see e.g., [13,50], and the mass spectrum in that case also strongly depends on the initial field value, so finetuning is required.
Further research should shed light on the mass spectrum in this theory. For example, it seems that the theory of modified gravity presented in this section could lead to very large and dense domain walls, so that they might be able to form baby universes.

Bondi Accretion
One of the first solutions to the accretion problem was that of Hoyle and Littleton [58], but the analytical formula was derived by Bondi [33]. Below is an overview of Bondi's solution.
The Bondi accretion is a spherically symmetric accretion on a compact object (Figure 3). The accretion rate is assumed to beṀ ≈ πR 2 ρv, where R is the capture radius or impact parameter, ρ is the density of the surrounding matter, and v is the the relative speed. The capture radius can be determined from the equality of the escape velocity and some characteristic velocity of matter. It is usually assumed to be equal to the speed of sound in the surrounding matter, then the accretion rate is obtained byṀ ≈ πρG 2 M 2 /c 3 s , where c s is the speed of sound in the matter surrounding the compact object. The obvious drawback of the model is that it does not take into account possible relativistic effects, it (the model) is Newtonian. However, it is possible in this approach to take into account the expansion of the universe [59]. The Bondi problem can be posed as follows:Ṁ = 4πr 2 ρv, where β(z) is the viscosity coefficient of the plasma around the accretor due to the interaction of electrons with photons, basically the Compton effect. The viscosity is given by the expression β(z) = 2.06 · 10 −23 x e (1 + z) 4 c −1 , where x e is the electron fraction in the plasma, and z is the redshift. Term M(< r) represents the mass contained within the radius r. The equation of state parameters are K and γ, where γ is a dimensionless constant. If we go to coordinates r = a(t)x, the Hubble expansion is added to the viscosity, and the effective viscosity is β e f f = β(z) + H. Thus, the rapid expansion of the universe reduces the accretion rate, and hence the accretor cannot significantly increase its mass. It is necessary to note that Bondi's solution is designed to describe a stationary nonrelativistic accretion. There are also approaches that take into account the luminosity of the accreting object if one radiates [60] and that consider supermassive black hole accretion within the Bondi solution [61].

Accretion Inside Neutron Star
In the paper [62], the accretion by a small black hole trapped inside a neutron star is considered. With this mechanism, it is possible to give constraints on the population of small-mass black holes, i.e., to use neutron stars as "dark matter detectors". The observed population of neutron stars imposes constraints on the mass range of PBHs: 10 −15 M M PBH 10 −9 M . As stated above, it is assumed that a black hole is trapped inside a neutron star and begins to absorb it. The equation of state of the matter inside the star is chosen in the form P = Kρ Γ , which gives the speed of sound near the centre of the neutron star a c ≈ ΓP c ρ c 1/2 . Then, the trapping radius of the black hole can generally be defined as: If we set m(r c ) ≈ 4π 3 r 3 c ρ c , then we can write the accretion rate as: Using (17) and (18), we can write the characteristic accretion time: Note that in the limit m(r c ) M BH , it is reduced to the Bondi accretion. This approach is mostly dedicated to impose a constraint on a specific mass range of PBHs, although it has a drawback: it is model-dependent since in order to calculate the speed of sound inside the star, one has to make an assumption about the equation of state of the plasma inside the star. That means in order to impose constraints on PBHs, one has to make an assumption about the equation of state and its parameters.
Neutron stars should be studied in order to figure out the equation of state inside them.

Eddington Limit
The Eddington luminosity or Eddington limit is the maximum luminosity that an object can achieve given with an equilibrium (see Figure 4) between the gravitational force and the radiative pressure force. This state is called a hydrostatic equilibrium. The gravitational force acting on the test particle with mass m is F GRAV = GMm R 2 , and the radiation pressure is given by the formula , where Φ is the flux and the denominator is where the integrand function characterises the distribution of matter within the object. Denote the plasma transparency by κ, then the radiation pressure force on the test object is F RAD = P RAD κm. Equating the radiation pressure force to the gravitational force and assuming a spherical symmetry, that the main component of the plasma are protons, and the dominant process is the Thomson electron scattering, we obtain The luminosity created by the accretion of matter can be represented in the following form L acc = εṀc 2 , where ε is the radiative efficiency. Equating the Eddington limit to the luminosity from the accretion of matter, we obtain: It is easy to get the dependence of mass on time: Formula (23) is usually used as a marginal estimate of the mass change due to accretion, but in reality, this limit can be circumvented by the lack of spherical symmetry [63,64]; however, in the process of observing early quasars at redshifts z∼6 [65], it was found that they were accreting at nearly the Eddington limit. It is also worth noting that the radiative efficiency ε is generally referring to a function of the angular momentum of the accreting object. This dependence can significantly affect the black hole mass growth; see, e.g., [66], where the influence of the angular momentum on the accretion rate was considered in some detail.
Although originally it was derived for stars to limit their maximum mass, it is also applicable as an estimation of the accretion rate for accretion disks [67,68]. It also must be stated that realistic accretion disks models are too complicated and are utilized only via numerical simulations.

Accretion in Schwarzschild Spacetime
As an approach related to the general theory of relativity equations, the Schwarzschild metric is given by the expression: where dΩ 2 = dθ 2 + sin 2 θdφ 2 . In [34,69], the (24) metric and the associated equality to zero of the covariant derivative of the energy-momentum tensor (EMT) were considered. The EMT was chosen as an ideal fluid as follows: Omitting the details of the derivation given in [34,69], let us write down the expression for the accretion rate: where A = const, which, generally speaking, depends on the parameter of the equation of state ω (p = ωρ), and ρ ∞ and p ∞ are the energy density and pressure at infinity, respectively. The Equation (26) does not change its form for arbitrary ω (only the constant A changes). For example, in [70], the RD (radiation-dominated epoch of the universe evolution with the equation of state parameter ω = 1/3) stage was considered and the authors obtained the same equation, but with a different constant A. Equation (26) is easily integrated if we put ρ(t) = ρ 0 (t 0 /t) 2 , and we obtain: The obvious drawback of this model is that the Schwarzschild metric is not asymptotically Friedmann-Robertson-Walker (FRW), which casts doubt on the effectiveness of the resulting formula in the early stages of the evolution of the universe, when the characteristic accretion time is comparable to or greater than the Hubble time, i.e., the characteristic expansion time of the universe. There is an approach to match the cosmological Friedmann-Robertson-Walker metric with the local Schwarzschild metric in the context of evaluating the radiation accretion in an expanding universe [70].

McVittie Solution
In 1933, McVittie studied the influence of cosmic expansion on local physics [36] and derived a specific metric given by (28). The mechanism of generation of PBH proposed in [24] within the modified gravity framework possibly allows one to produce them practically immediately after the cosmological inflation. The question about the growth of mass of a PBH with the evolution of the universe naturally arises. Since there are no observational data on this distant period of evolution of the universe, it is necessary to estimate the influence on the final result of the free parameters of the problem, namely, the parameter ω (p = ωρ) and the moment of the end of reheating t reh . Of course, the initial mass of a black hole is also a free parameter; however, taking into account the fact that the early universe is still very rapidly expanding, it is necessary to consider an accretion model of matter taking into account the expansion of the universe. In [35,71], exact solutions to the accretion problem were considered. The solution with the McVittie metric and the energy-momentum tensor of a nonideal fluid with a radial flow is of most interest for the evaluation of the accretion in that case. In that model, it is possible to calculate how many times the mass of the PBH will change, as the mass of the PBH grows with the universe, depending on the state parameter ω.
In the paper [35], the following metric was considered: where . This metric is asymptotically FRW, and when the universe "stops" expanding, it passes into the Schwarzschild metric in isotropic coordinates via the radial coordinate substitutionr = ar. Obviously, it describes a strongly gravitating object. In this spacetime, the physically relevant [72][73][74] mass will be the quasi-local mass m H (t) = m(t)a(t), and m(t) is only the coefficient of the metric. The EMT of matter around a black hole is in which where q c describes the radial energy flux and u a is the four-velocity of the surrounding fluid. Then, the field equations are as follows: where A = dθdϕ √ g Σ = 4πa 2 A 4 r 2 and C =ȧ a +ṁ Ar . From the last two equations of the system, we obtain: Thus, we obtain the accretion rate: Consider Formula (36). It is useful to find the dependence of the accretion rate on cosmological parameters; to do so, following [71], let us take the r → ∞ limit for expression (36): It is interesting to compare this formula with others, such as the Babichev-Eroshenko-Dokuchaev formula [34,69,70]: Formula (38) can be derived using the stationary Schwarzschild metric or the nonstationary Schwarzschild metric, in which metric coefficient M is now a function of the time M(t) [75]: Since the metric (28) changes into (39) at the "stop" expansion of the universe (in other words, we insertȧ = 0 into the field equations) and replaces the radial coordinate, one would expect that the accretion rate (36) would change into Formula (38), but in general, this is not the case. The reason is that in deriving this formula, the assumption lim r→∞ (ur 2 ) = −2AG 2 M 2 was made. We continue to consider the limit r → ∞ for Equation (36). For a large r, we obviously have: u(r; t) = u ∞ (t)/r 2 + O(1/r 3 ).
Then, by substituting these approximations into the field equations and by combining the conservation laws, one can eventually obtain the dependence of m H (t) on the scale factor. See [71] for details of the derivation. For the quasi-local mass, we obtain a second-order differential equation: The last expression contains a constant of integration m 0 , which has not yet been assigned a definite meaning. The origin of this term is as follows: the differential equation for the first coefficient of the main part of the Laurent series for the density (40) can be derived from conservation laws. The equation for ρ 1 is as follows: from which we get Thus, an additional assumption about the value of m 0 at some initial point in time is required. Hence, it follows that (44) requires not two but three initial conditions to find a partial solution.
Thus, the dependence of the quasi-local mass on the scale factor is as follows: in which C 1 and C 2 are determined from the initial mass m H (t 0 ) and the initial accretion rateṁ H (t 0 ). An example of using (46) is given in Figure 5. This model takes into account the cosmic expansion through the metric of the spacetime, although this solution cannot be considered as realistic, since the fluid velocity on the sphere r = Gm(t)/2 is superluminal (it is easy to see from Equation (36)).

Discussion
In the present review, we considered some mechanisms of the formation of PBHs such as the conventional mechanism based on the collapse of density fluctuations, phase transitions, e.g., first-order phase transitions which are strongly constrained by now. Moreover, the inflationary mechanism for domain wall production with their consequent collapse was considered, as well as the case of symmetry breaking with a consequent wall collapse. A recently elaborated mechanism based on the dynamic of extra space in a model of f (R)-gravity was also considered. They are mentioned in the following with a brief list of advantages and disadvantages: • Initial density inhomogeneities: Advantages: a broad mass spectrum. Disadvantages: it leads to strong inhomogeneities and there is no possibility to form clusters of PBHs. • First-order phase transitions: Advantages: it does not require any assumption beyond standard Big Bang physics. Disadvantages: The mass spectrum in this model is close to monochromatic, so it is impossible to explain the existence of black holes of various masses. It is also impossible to produce clusters of PBHs within this model. • Second-order phase transitions: Advantages: it has a broad mass spectrum and it also could be coupled with external processes, e.g., baryogenesis. Disadvantages: a fine-tuning of the initial conditions is required.
The best-known black hole accretion models were also considered, such as the conventional Bondi solution (Newtonian approach), as well as a special case of the latter-accretion inside neutron stars. The Eddington limit as an approach of the marginal estimation of the accretion efficiency was also discussed. Next, a generalization of the Bondi accretion to the case of a curved Schwarzschild spacetime was considered. We also considered the accretion on a McVitty black hole, which was cosmologically coupled to the FRW universe. Here are the models and a brief list of their advantages and disadvantages: The topic of primordial black holes at the moment is a subject of great interest. They are potentially capable of explaining various cosmological problems, e.g., the origin of supermassive black holes, early structures' formation, and so on; thus, it is essential to elaborate and study the possibilities of their formation and evolution.