Interest in the primordial black holes (PBHs) is dramatically increasing since the gravitational waves discovery from the black holes mergers (Abbott, 2016). However, PBHs origin and possible formation mechanisms are still a topical issue of modern astrophysics and cosmology. The first ideas of such mechanisms had been proposed in (Zel’dovich and Novikov, 1967; Hawking, 1971; Carr and Hawking, 1974) and lately developed in many other works (see reviews and references within (Khlopov, 2010; Carr et al., 2020; Carr and Kühnel, 2020)). The different PBH spectra are used in papers (Carr, 1975; Dolgov and Silk, 1993; Sendouda et al., 2006; Clesse and García-Bellido, 2015; Garriga et al., 2016; Carr and Kühnel, 2019; Liu et al., 2020) depending on specific needs.

The phase transitions of the first (Hawking et al., 1982; Kodama et al., 1982; Konoplich et al., 1998; Jedamzik and Niemeyer, 1999; Konoplich et al., 1999) and the second order (Rubin et al., 2000; Rubin et al., 2001) might also underlay a mechanism of the PBH formation. In this paper, we continue elaboration of the model based on the second type phase transitions during the inflationary epoch (Rubin et al., 2000; Rubin et al., 2001; Khlopov et al., 2002). However, the described model has a flaw. It inevitably leads to a multiple production of small mass PBHs. That problem could not be avoided within the framework of the discussed scenario, and typical mass spectra have the falling form d N/d M ∝ M ^−α , α > 0 (see review (Belotsky et al., 2019)). Such a form of spectra could be unfavorable for explaining the observable effects. In addition, the overproduction of low-mass PBHs could contradict future experiments.

The solution to the above problems is to use the classical motion of massive scalar fields together with their quantum fluctuations. The idea was firstly studied in (Dokuchaev et al., 2010) to suppress the production of intermediate-mass black holes. In this research, we have elaborated this idea and obtain the analytical formula for the field distribution probability.

At present, there are a lot of models that contain a complicated form of scalar field potential. The latter is used in a variety of inflationary models predicting the potential landscape. In addition, inflation can be driven by the dynamics of several fields (Wands et al., 2007). For instance, the supergravity often produces more than one physical scalar field (Ketov and Starobinsky, 2012) and predicts nontrivial forms of inflaton potentials (Ketov, 2021). The string theory also predicts the landscape with a large number of vacua, local peaks, and saddle points (Susskind, 2003; Cline, 2005). Such a complex potential can have both random and quasi-periodic forms (Li et al., 2009) and leads to the multi-field inflation such as multi-stream one (Li and Wang, 2009; Duplessis et al., 2012), assisted inflation (Battefeld and Battefeld, 2009) or multi-field inflation with a random potential (Tye et al., 2009; Frazer and Liddle, 2012). Therefore, other non-inflaton scalar fields might have complicated potential as well.

In this paper, we adopt the mechanism of black hole formation (Rubin et al., 2000; Rubin et al., 2001; Khlopov et al., 2002; Belotsky et al., 2019) to the modern trends in the inflaton potential complication. A complex potential shape influences the classical motion of the fields so that the analytical form of the probability contains the classical trajectories as well as the quantum field contributions. The resulting black holes mass spectrum appears to be related to a shape of the scalar field’s potential. The tool elaborated here allows fixing the potential form by knowledge of black hole mass distribution.

This paper is organized as follows. In Section 2, we elaborate the way to involve the classical part of scalar fields into the expression for its fluctuations probability. The PBH spectrum depends on an initial position of the scalar field that allows us to adjust the model predictions to future observational data without inserting small parameters. The numerical results are represented in Section 3. Finally, Section 4 concludes the paper.

The discussed mechanism of PBHs production requires closed domain walls formation due to the quantum fluctuations of scalar fields at the inflation epoch (Rubin et al., 2000; Rubin et al., 2001). Let us take into account both the quantum and classical motion of fields. Consider the scalar field Φ of mass m and the standard action S = M Pl 2 2 ∫ d 4 x g 4 R + ∂ Φ 2 − m 2 Φ 2 . (1)

Here, M _Pl is the Planck mass. The scalar field could be the inflaton field as well as a spectator one. The field equation in the de Sitter space is represented as (Khlopov and Rubin, 2004) ∂ Φ ∂ t − 1 3 H e − 2 H t Δ Φ − ∂ V Φ ∂ Φ = y x , t ; y x , t ≡ − 1 3 H ∂ 2 ∂ t 2 − ∂ ∂ t + 1 3 H e − 2 H t Δ Q x , t . (2)

Here, Q (x, t) is the “quick” part of the Fourier field decomposition. This equation was simplified: we have omitted the second time derivative due to a slow roll approximation and have neglected higher powers of the function y (x, t). The latter is supposed to be small so that we may find a solution to the equation in the form Φ = Φ cl + ϕ . (3)

The deterministic part of the classical field Φ_cl is governed by the equation ∂ Φ cl ∂ t − 1 3 H e − 2 H t Δ Φ cl − ∂ V Φ cl ∂ Φ cl = 0 , (4) while its random part ϕ depends strictly on quantum fluctuations according to the linear equation ∂ ϕ ∂ t − 1 3 H e − 2 H t Δ ϕ − V ″ Φ cl ϕ = y x , t . (5)

Here, we consider the limit Φ_cl ≫ ϕ which is valid if the random “force” y (x, t) is small. Let us denote m 2 t ≡ V ″ Φ cl . (6)

The parameter m is positive if we are near the bottom of potential and is imaginary if we are near the potential maximum. It is supposed that m(t) varies slowly during inflation.

We are interested in the super horizon scales where the fluctuations do not depend on the space coordinates. The uniform distribution Φ = Φ(t) is governed by the more simple equation ∂ Φ cl ∂ t + 1 3 H ∂ V Φ cl ∂ Φ cl = 0 , (7) ∂ ϕ ∂ t + μ t ϕ = y t ; μ t ≡ m 2 t 3 H , (8) provided that H(t) = const. The correlator of the random function y(t) may be approximated as follows (Rey, 1987) y t 1 y t 2 = D x , t 1 ; x , t 2 = H 4 π 2 δ t 1 − t 2 . (9)

The delta function in the rhs of this expression indicates that the random function y(t) is distributed according to the Gauss law with the density W y = const exp − 1 2 σ 2 ∫ y 2 t d t , σ = H 3 / 2 2 π . (10)

The probability distribution of the function ϕ is proportional to that of the function y(t) due to their linear relationship (8). It means that the probability to find the specific value ϕ(t) inside some small interval is equal to (Feynman et al., 2010) d P ϕ = const D ϕ exp − 1 2 σ 2 ∫ ∂ ϕ ∂ t + μ t ϕ 2 d t . (11)

Let’s obtain the probability to find a quantum part of the field ϕ ₂ at an instant t ₂ provided that a value ϕ ₁ at an instant t ₁ is known. Evidently, we have to integrate over all values of the field inside the interval (t ₁, t ₂) and come to the expression d P ϕ 2 , t 2 ; ϕ 1 , t 1 = const d ϕ 2 ∫ ϕ 1 ϕ 2 D ϕ exp − 1 2 σ 2 ∫ t 1 t 2 ∂ ϕ ∂ t + μ t ϕ 2 d t . (12)

The constant factor in this equation is determined by normalization condition ∫ − ∞ ∞ d P ϕ 2 , t 2 ; ϕ 1 , t 1 d ϕ 2 d ϕ 2 = 1 . (13)

Functional integral (12) can be calculated in the standard manner by finding an extreme trajectory of the integral in the exponent ϕ ̈ − μ 2 t ϕ = 0 , (14) where the term μ ̇ is neglected due to slow variation of μ(t). The boundary conditions for (14) are as follows ϕ t 1 = ϕ 1 ; ϕ t 2 = ϕ 2 . (15)

Exact solution to this equation is ϕ t = A ⁡ exp M t + B ⁡ exp − M t ; M t ≡ ∫ t 1 t μ t ′ d t ′ ; (16) A = ϕ 2 − ϕ 1 e − M t 2 2 ⁡ sinh M t 2 , B = − ϕ 2 − ϕ 1 e M t 2 2 ⁡ sinh M t 2 . (17)

Notice that M(t ₁) = 0 by definition.

Substituting this solution into the integral in the exponent of the expression (12) one obtains the desired probability in the saddle point approximation d P ϕ 2 , t 2 ; ϕ 1 , t 1 = const d ϕ 2 exp − 2 A 2 σ 2 ∫ t 1 t 2 μ 2 t e 2 M t d t = const d ϕ 2 exp − q 2 ϕ 2 − ϕ 1 e − M t 2 2 , (18) q 2 ≡ 1 2 σ 2 ⁡ sinh 2 M t 2 ∫ t 1 t 2 μ 2 t e 2 M t d t . (19)

It describes the probability to find specific value of “quantum” part of the field (3). The “classical” part of the field Φ_cl is incorporated into the function M(t). The probability for the field value F (the distribution function f) is easily obtained by substitution ϕ(t) = Φ(t) − Φ_cl(t) into the formula above. f Φ 2 , t = d P Φ 2 , t ; Φ 1 , t 1 d Φ 2 = q 2 π exp − q 2 Φ 2 − Φ cl t 2 − Φ 1 − Φ cl t 1 e − M t 2 2 . (20)

The limit m → 0 restores the textbook formula.

The next section aims to demonstrate how the obtained formulas can be applied to a particular scalar field potential. It is assumed that the potential may possess many extremes of a different kind. In our consideration, we choose a part of phase space containing two maxima and at least one saddle point.

In this section, we show that the classical motion of fields together with their quantum fluctuations influence the PBH mass spectra. To this end, we have to find the classical trajectory and use the probability (20) derived above.

The fields move between potential local maxima that lead to complicated spectra of fluctuations. The latter are discussed in papers (Wands et al., 2007; Ketov and Starobinsky, 2012). At the same time, the presence of saddle points is the reason for the closed domain walls formation, see details in (Gani et al., 2018; Murygin et al., 2021). In the following, they could collapse to black holes (Rubin et al., 2001; Belotsky et al., 2019).

Let us consider the model of two real scalar fields with the Lagrangian L = 1 2 ∂ μ Φ ∂ μ Φ + ∂ μ X ∂ μ X − V Φ , X . (21)

We choose the potential possessing n peaks and saddle points V Φ , X = m ϕ 2 2 Φ 2 + m χ 2 2 X 2 + ∑ i = 1 n δ V i Φ , X , δ V i Φ , X = Λ i 4 ⁡ exp − Φ − ϕ i Δ i 2 − X − χ i Δ i 2 . (22)

Here, δV _i describes the ith local maximum. The global minimum of the potential is located at the point (ϕ _min, χ _min) = (0, 0) with exponentially small errors. Hereinafter, all variables are taken in the Hubble units H where H ≈ 10¹³ GeV at the inflationary epoch.

For our estimates, we choose the fields masses m _ϕ = 0.4 and m _χ = 0.5. For simplicity, we consider the potential with two peaks (n = 2) with the coordinates ϕ ₁ = − 9.0, χ ₁ = 3.0 and ϕ ₂ = − 1.7, χ ₂ = 4.5. The parameters corresponding to the peaks heights are Λ₁ = 3.0 and Λ₂ = 1.5, and the peaks widths are set with Δ₁ = 0.5 and Δ₂ = 1.5. The initial fields values are ϕ _in = − 8.0 and χ _in = 45.0. Note, all chosen parameters have the values ∼ O ( 1 ) .

Following Section 2, the first step consists of finding the classical trajectory Φ_cl(t), X_cl(t) of the fields Φ, X. Starting from the initial values (ϕ _in, χ _in) at the inflation epoch, the scalar fields tend to the potential minimum. The process is described by the classical motion equations Φ cl , t t + 3 H Φ cl , t + ∂ V Φ cl , X cl ∂ Φ cl = 0 , X cl , t t + 3 H X cl , t + ∂ V Φ cl , X cl ∂ X cl = 0 . (23)

The classical evolution of the fields and the form of the specific potential are shown in Figure 1. At the same time, quantum fluctuations lead to fields “diffusion” during inflation. The probability density f to find the fields F or X in any point of the physical space is given by (20).

Both F and X distributions depend on a classical position of the fields at the instant t. In our estimates, we suppose that the probability function for the quantum parts ϕ or χ of the fields is separated into two independent fluctuation processes f t , ϕ , χ = f ϕ t , ϕ f χ t , χ = d P d ϕ d P d χ , (24) where distribution functions of each field ϕ, χ can be defined as in (18).

As noted above, the fields should reach a saddle point of potential for domain wall formation. Suppose that some quantum fluctuation crosses a saddle point and ends up at a point of the area O. As was shown in (Gani et al., 2018; Murygin et al., 2021), it causes nontrivial field solutions of the system (23) characterized by a nonzero winding number. Such configurations might lead to the domain walls formation after the inflation is finished. Detailed explanation might be found in (Gani et al., 2018). The calculation of an exact shape of the area O is a separate, quite complicated task, so that we limit ourselves with the following approximation. Let us assume that the area O is bordered by two lines χ _cl-sp(Φ) and χ _{min  −sp}(Φ) in the phase space. The first line connects the classical value at the instant t and the saddle point (ϕ _sp, χ _sp) χ cl - sp Φ = χ sp + Φ − ϕ sp ϕ cl t − ϕ sp χ cl t − ϕ sp . (25)

The second one connects the vacuum value and the saddle point χ min − sp Φ = χ sp + Φ − ϕ sp ϕ min − ϕ sp χ min − ϕ sp . (26)

Thus, the probability for the fields to attain the area O where domain walls might form is calculated by integrating (24) P t = ∫ ∫ Ω f ϕ , χ , t d X d Φ = ∫ ϕ sp + ∞ f ϕ ϕ , t d Φ ∫ χ cl - sp Φ χ min − sp Φ f χ χ , t d X . (27)

Here, both distribution functions f _ϕ and f _χ are defined in (20) and according to (3) ϕ = Φ − Φ_cl, χ = X − X_cl. The algorithm for calculating the probability (27) discussed above is more accurate then that used in the previous papers.

Now, let us find the mass spectra of primordial black holes. In the considered model, they are formed due to the collapse of domain walls. Here, we briefly reproduce the idea, while details may be found in the review (Belotsky et al., 2019). As was shown in (Gani et al., 2018; Murygin et al., 2021), domain walls might be formed due to the quantum fluctuations in field models with potential possessing at least one saddle point and a local maximum. The proto-soliton is formed if the fields achieve a saddle point (in our model, we have noted this area as Ω). These proto-soliton field configurations are quickly expanded during inflation. The final scale of such configuration depends on an e-fold number N. More definitely, the configuration scale is stretched in the factor ∼ e N inf − N to the end of inflation. The soliton is quickly formed after the end of inflation. The total mass of the field configuration is proportional to its area. It could collapse into a black hole after the end of inflation (Rubin et al., 2000; Rubin et al., 2001; Khlopov et al., 2002).

The regions number where the fields reach the critical values ϕ _cr and χ _cr belonging to O can be found as n t = P t e 3 H t . (28)

Here, the term e ^3Ht is the number of causally independent regions of the size H ⁻¹ at the instant t from the beginning of inflation. After the end of inflation at t = N _inf H ⁻¹ = 60H ⁻¹, the size of each region is expanded r 0 t = H − 1 ⁡ exp N inf − H t ≫ r h , 0 . (29)

Here, N _inf ≈ 60 is the total e-folds number, and r _h,0 is the horizon size at the end of inflation. At the radiation stage (RD), each region expands as ∝ τ while the horizon size is r _h = H ⁻¹(τ) = 2τ. Here, τ is time after the beginning of the RD epoch. After a domain wall goes under the horizon, its collapse begins (the details of the process taking into account detachment from the Hubble flow may be found in (Khlopov et al., 2005; Dokuchaev et al., 2008; Belotsky et al., 2019)). Thus, the maximal size of a domain wall can be written as the function of the instant t r t ≈ 1 2 e 2 N inf − H t H N inf . (30)

After eliminating of t from (28) and 30, one can finally get the distribution n(r) of closed walls sizes which can be rearranged into the mass spectrum.

Next, we have to find masses of PBHs. For simplicity, we assume total energy of domain wall converts to a black hole mass during its collapse and neglect nonsphericity of a domain wall and losses caused by gravitational waves. The energy density of a domain wall might be found by a common way. The energy momentum tensor for the Lagrangian (21) is given by T ν μ = ∑ i = 1,2 ∂ μ φ i ∂ ν φ i − 1 2 ∂ α φ i ∂ α φ i δ ν μ + V δ ν μ , (31) where φ ₁, φ ₂ correspond to the fields Φ_cl and X_cl, respectively. Then, the energy density of a domain wall is found to be ε x = T 0 0 = 1 2 ∑ i = 1,2 ∂ t φ i 2 + ∂ x φ i 2 + V . (32)

Upon integrating (32) over the all possible values of x (infinite interval), the surface energy density of a domain wall σ may be found. Finally, masses of black holes are M r ( t ) ≃ 4 π σ r 2 ( t ) . Taking into account (28) and (30), one can find the mass spectrum of primordial black holes. The PBH mass distribution for the parameters of the Lagrangian (21) is shown in Figure 2. Note, the mass spectrum has the non-power form due to taking into account both the quantum and classical motion of scalar fields and the compound form of the potential leading to nontrivial classical fields trajectory. The obtained spectrum is free from the overproduction of light PBHs and, therefore, is much more adaptable to new observational effects. We leave a detailed analysis of these possibilities for future research.

In this note, we have shown that the mechanism of PBHs formation in the second-order phase transitions of scalar fields might produce a wide variety of the PBH mass spectra. It is expected that observations (e.g. cosmic gamma-rays, gravitational waves spectrum (Sakharov et al., 2021), gravitational lensing (Toshchenko and Belotsky, 2019)) will help to select an appropriate one. The key point is the classical field motion which was taken into account together with the quantum fluctuations at the inflationary stage. We show here that the probability to find a particular value of the scalar field at a space point depends on its classical dynamics. We have derived the appropriate analytical formula and have applied it to obtain one of the PBH mass spectra. The elaborated method is the useful tool to fit an observable spectrum in the near future.