Primordial black holes from slow phase transitions: a model-building perspective

We investigate the formation of primordial black holes (PBHs) through delayed vacuum decay during slow cosmic first-order phase transitions. Two specific models, the polynomial potential and the real singlet extension of the Standard Model, are used as illustrative examples. Our findings reveal that models with zero-temperature scalar potential barriers are conducive to the realization of this mechanism, as the phase transition duration is extended by the U-shaped Euclidean action. We find that the resulting PBH density is highly sensitive to the barrier height, with abundant PBH formation observed for sufficiently high barriers. Notably, the phase transition needs not to be ultra-supercooled (i.e. the parameter $\alpha\gg1$), and the commonly used exponential nucleation approximation $\Gamma(t)\sim e^{\beta t}$ fails to capture the PBH formation dynamics in such models.


Introduction
Due to their profound cosmological and astrophysical implications, primordial black holes (PBHs) have attracted significant research attention [1][2][3][4].Those hypothetical black holes originate in the early Universe shortly after the Big Bang, long before the existence of stars and galaxies.Cosmic first-order phase transitions (FOPTs) provide a favorable environment for PBH formation: the Universe transitions from a higher energy false vacuum to a lower energy true vacuum via tunneling, resulting in the true vacuum bubbles nucleate and expand in the false vacuum background.This may lead to the formation of PBHs, if the bubble collisions compress a large amount of energy into a small volume [5][6][7], or if some false vacuum remnants form significant overdensities via delayed vacuum decay  or trapping fermions [29][30][31][32][33][34][35][36][37][38][39][40][41][42][43][44][45][46][47].While FOPTs might be generic during the thermal history of the Universe [48], the formation of PBHs during such transitions is highly nontrivial.It typically imposes additional constraints on the temperature, strength, duration, and bubble expansion velocity of the FOPT, as well as the particle content of the underlying model.Depending on scenarios, the PBHs from FOPTs have various mass, spin, and abundance distributions.
In this work, we focus on the mechanism based on the delayed vacuum decay [11] (also see a similar scenario in Ref. [8]).During the FOPT, the energy stored in the false vacuum regions is released and converted to radiation energy.Because of the randomness of the vacuum tunneling process, the timing of bubble nucleation varies across Hubble patches.Some patches experience delayed nucleation compared to the average time, and become overdense over time as they have larger vacuum energy fractions, whose density maintains a constant, while the overall Universe is dominated by radiation that dilutes its energy density via expansion.Those patches eventually decay and release vacuum energy to radiation energy, and when the overdensity reaches a certain threshold, the delayed-decay patches collapse into PBHs.
This mechanism has been extensively studied in the literature on its detailed calculations, applications, and phenomenology [12][13][14][15][16][17][18][19][20][21][22][23][24][25].Most studies are agnostic about the underlying particle models, employing the FOPT parameters such as α and β/H * as input parameters.In particular, the exponential nucleation rate approximation Γ(t) ∼ e βt is commonly adopted.In this paper, we focus on the implementation of the mechanism in particle physics models beyond the Standard Model (SM), and for the first time analyze its general relationship with the scalar potential structure of the model.By investigating models that incorporate zero-temperature potential barriers, we find that such models typically exhibit a U-shaped Euclidean action S 3 (T )/T , which results in slow FOPTs, conducive to delayed vacuum decay and hence PBH formation.The PBH abundance is very sensitive to the height of the zero-temperature barrier.We also confirm that the Γ(t) ∼ e βt approximation is not suitable due to its reliance on the linear time-expansion of the Euclidean action, which fails to capture the characteristics of a U-shaped action in such models [49].
To illustrate our findings, we consider two concrete models.The first model utilizes a polynomial potential, serving as a prototype or toy model.The second model is the Z 2 -preserving real singlet extension of the SM (dubbed as the Z 2 -xSM), representing a realistic scenario concerning the first-order electroweak phase transition.In Section 2, we describe our calculation framework in detail, and also compare different approximation methods for the nucleation rate.Subsequently, we apply this framework to analyze the polynomial potential in Section 3 and the Z 2 -xSM model in Section 4, presenting our main results.The conclusion will be given in Section 5.

Deriving the PBH profiles
In the early Universe, the scalar potential's shape changes with temperature.Initially, the Universe occupies the minimum of the potential, known as the vacuum.However, below the critical temperature, the potential develops another deeper minimum, which becomes the true vacuum.The Universe remains stuck in the initial false vacuum due to a barrier that prevents a smooth transition to the global minimum.Thermal tunneling allows the Universe to transition to the true vacuum with a probability per unit volume and per unit time [50] Γ(T ) ∼ T 4 S(T ) 2π where S(T ) = S 3 (T )/T , and S 3 (T ) is the O(3)-symmetric Euclidean action evaluated from the scalar potential.As the temperature T depends on cosmic time t, we can express the action and decay rate as S(t) = S (T (t)) and Γ(t) = Γ (T (t)), respectively.Denoting t c as the time corresponding to the critical temperature T c where the two vacua are degenerate, we have S(t c ) = ∞ and Γ(t c ) = 0. We now describe the evolution of the Universe in a FOPT.For a given Hubble patch, according to the first Friedmann equation, the Hubble constant is given by where M Pl = 1.22 × 10 19 GeV is the Planck scale, a(t) is the scale factor, and ρ(t) = ρ r (t) + ρ v (t) is the energy density that can be decomposed into the radiation and vacuum components where the factor g * is the number of effective degrees of freedom, and ∆V (t) is the energy density difference between the false and true vacua, which is 0 for t ⩽ t c and positive for t > t c .The false vacuum volume fraction is defined as where t d is the time that the Hubble patch starts to nucleate, and v w is the bubble expansion velocity.
As the FOPT proceeds, F (t) decreases from 1 to 0, and the vacuum energy is converted to radiation.This can be clearly seen from the variation of the second Friedmann equation, Given t d , with the initial condition ρ r (t c ) = π 2 g * T 4 c /30 and ρ v (t c ) = 0, one is able to resolve Eqs.(2.2), (2.4), and (2.5) consistently by iteration, and obtain the evolutions of the physical observables, such as T (t), ρ r,v (t), H(t), and F (t), etc.As mentioned before, t d varies in different patches since vacuum tunneling is probabilistic.For simplicity, we categorize the Hubble patches into two types: normal ones and delayed-decay ones.Normal patches initiate nucleation after the Universe cools to the critical temperature, i.e. t d = t c , which is the most prevalent case.Delayed-decay patches, on the other hand, have a late nucleation time, i.e. t d > t c .We use the subscripts "out" and "in" to label the two types of patches and the corresponding physical observables, respectively.The normal patch can be resolved once the underlying particle model is set, while the delayed-decay patch is determined by one extra parameter t d .
The delayed-decay patches form overdensities with respect to to normal patches because the latter receive more radiation contributions, and ρ r ∝ a −4 redshifts rapidly while ρ v ∝ a 0 remains a constant.The overdensity is defined as As time t increases from t c , δ(t) increases from 0 to a maximum value, δ max , and then decreases back to zero as t → ∞.If δ max exceeds a threshold value δ c , the delayeddecay patch collapses into a PBH at t pbh satisfying δ(t pbh ) = δ c .In this paper, we take δ c = 0.45 as the PBH formation criterion [51,52].This criterion can be obtained under the assumption that the Universe is radiation dominant and that the overdensity region is spherically symmetric.If the overdensity region occurs due to vacuum energy, an inflationary expanding baby Universe is realized there [53][54][55][56], and the region appears to the outside observer as a PBH.As it has not been numerically shown what value of δ c should be taken in such situations, we take the well-known value δ c = 0.45 as a benchmark in our study.
We can establish a mapping between t d and δ max .Generally, δ max increases with t d since a patch with later vacuum decay forms a larger energy contrast.However, the probability of "a patch remains in the false vacuum until t d " decreases rapidly as t d increases.Thus, we infer that the smallest value of t d allowing for collapse dominates the abundance of PBHs.This corresponds to the value of t d where δ max = δ c , and the moment when this maximum value is reached defines the PBH formation time t pbh , specifically δ(t pbh ) = δ max = δ c .By this procedure, we can determine t d and t pbh simultaneously (note t pbh > t d ).The resultant PBH mass is where V H (t pbh ) = (4π/3)H −3 in (t pbh ) is the Hubble volume.The factor γ represents the ratio between the PBH mass and the Hubble mass.We here take γ ≈ 0.2 [57].
To determine the abundance of PBHs, it is necessary to evaluate the probability of a patch nucleating at time t d .That is to find a region with a Hubble volume (4π/3)H −3  in (t pbh ) that remains bubble-free until t d .Within the time interval t c < t < t pbh , this region has a volume of (4π/3)[a in (t)/a in (t pbh )] 3 H −3 in (t pbh ), and the probability of no bubble formation in this region during the time interval [t, t + dt] is given by dP = 1 − 4π 3 Consequently, the probability of no bubble nucleation during the time interval t ∈ [t c , t d ] is given by This expression demonstrates that P (t d ) is exponentially dependent on t d , emphasizing the significance of setting the PBH formation time and mass through δ(t pbh ) = δ max = δ c .This lower limit on t d governs the PBH density.With P (t d ) in hand, we can calculate ρ pbh .Suppose the numbers of the delayed-decay and normal Hubble patches are N in and N out , respectively, and then the energy density of PBH is where V H is again the Hubble volume.The delayed-decay probability P (t d ) implies and substituting this relation back to Eq. (2.10) yields where the approximate equality holds because P (t d ) ≪ 1.The present-day energy density of PBHs is given by ρ 0 pbh = ρ pbh s 0 /s, where s ≈ 2π 2 g * s T 3 out (t pbh )/45 and s 0 ≈ 2891.2 cm −3 [58] the entropy at PBH formation and today, respectively.The relic abundance of PBHs can be expressed as Ω pbh = ρ 0 pbh /ρ 0 , where ρ 0 = (3M 2 Pl /8π)H 2 0 with H 0 ≈ 67.4 km/(s • Mpc) the current Hubble constant [58].We further define the PBH fraction of dark matter as where Ω dm h 2 ≈ 0.12 is the dark matter relic abundance.
We have described our framework to determine the formation time, mass, and relic abundance of PBHs.The formalism applies to general particle models with various shapes of potentials.Our calculation disregards the FOPT reheating effect, assuming a constant product of temperature (T ) and scale factor (a).This simplification may not hold for ultra-supercooled FOPTs with the parameter α ≫ 1, where the Universe's temperature can be increased by a factor of approximately (1 + α) 1/4 after the transition.However, our analysis focuses only on the mild-supercooled case when α is moderate, making the reheating effect negligibly small.More discussions on calculations involving very strong FOPTs can be found in Refs.[16,17].

The exponential approximation and beyond
Once the underlying particle model is built, one is able to derive S(t) and also Γ(t), and proceed the calculation as described in Section 2.1.To simplify the calculation, most previous studies have chosen to expand the action at a specific moment t * , and have written down where T * is the temperature at t * , and Γ * = T 4 * (S(t * )/(2π)) 3/2 .Here the Taylor expansion coefficients are or they can be rewritten in the dimensionless forms and transferred into the derivatives on temperature, which are more explicitly derived from the particle model: where H * is the Hubble constant at t * , and S(T ) = S 3 (T )/T .Regarding the expansion Eq. (2.14), the commonly used approach is to consider only the linear β-term, leading to the exponential nucleation rate approximation Γ(t) ≈ Γ * e β(t−t * ) , where β −1 is interpreted as the time scale of the phase transition.Therefore β/H * can be understood as the ratio of Hubble time scale to FOPT duration.Besides β/H * , there is another important parameter α describing the strength of the FOPT, defined as i.e. the ratio of phase transition latent heat to the radiation energy density, evaluated in the normal patches.The parameters α and β/H * , together with the transition temperature T * and bubble wall velocity v w , can be used to calculate other observables of the FOPT dynamics, such as the stochastic gravitational wave (GW) spectrum [48,[59][60][61].
As for the expansion moment t * , there are typically two choices.One is the nucleation time t n , the moment that the bubble number in the normal patch reaches 1, i.e.
The other is the percolation time t p , which is the moment when the bubbles in the normal patch form an infinite connecting cluster, i.e., F (t p ) = 0.71 [62].In general t p > t n , and when the FOPT is prompt (β/H * ≫ 1) and moderate (α ≲ 1), t p ≈ t n .However, if the FOPT is very strong or lasts for a long time, t p differs significantly from t n .In that case, it is proposed to adopt t * = t p , and replace β/H * with (8π) 1/3 v w /(H * R) in the GW calculation, where R is the mean separation of bubbles [49,63,64], which can be adopted as [n b (t * )] −1/3 , with the bubble number density in the normal patches.Furthermore, one should also check the FOPT completion condition [49,65] such that the physical volume of the false vacuum is decreasing at percolation.The exponential approximation is widely used in FOPT research and is generally effective.It is also commonly employed in the study of PBH formation from delayed vacuum decay [11][12][13][14][15][17][18][19][20][21][22][23][24][25][26].However, PBH formation predominantly occurs in the small β/H * region, which corresponds to slow FOPTs with prolonged durations.For instance, previous studies indicate β/H * ≲ 10 for abundant PBH formation [11,16,17].When the first-order expansion coefficient β of S(t) is small, it becomes necessary to check if truncating the expansion Eq. (2.14) at linear order provides a satisfactory approximation.This can be determined by examining whether the ratio of the second-order expansion to the first-order one, In our study, we directly compare the results obtained using three different methods for Γ(t): (1) the full expression, (2) the first-order expansion with only the β-term (exponential approximation), and (3) the second-order expansion with both the β-and ζ-terms.We find that the first-order expansion is inadequate, and it results in significant differences in PBH mass and abundance compared to the results obtained from the full expression case.By including the quadratic term, the results are much improved and closer to those of the full expression case, although visible differences still persist, particularly for relatively mild FOPTs.Our results demonstrate the importance of utilizing the full expressions for Γ(t) in the PBH formation calculation in realistic models.3 The model with a polynomial potential

Model description and PBH formation
Let ϕ be the scalar that triggers the FOPT, and V (ϕ, T ) its potential in the thermal bath of the early Universe at temperature T , we parametrize the potential as and restrict the value of µ 3 to be within the range At T = 0, V (ϕ, 0) has two minima at ϕ = 0 and w, and the latter is the true vacuum.The ϕ boson mass square is m 2 ϕ = d 2 V (ϕ, 0)/dϕ 2 ϕ=w .Between the two minima, there is a local maximum at ϕ = w b , corresponding to the zero-temperature barrier (3.3)It exists when µ 3 > m 2 ϕ /w.If we increase µ 3 , the barrier height increases accordingly, and V (w, 0) is also shifted.When µ 3 reaches 3m 2 ϕ /w, the ϕ = w vacuum becomes degenerate with the ϕ = 0 one.Therefore, we consider µ 3 in the range in Eq. (3.2) to ensure there is a barrier at zero-temperature, and ϕ = w is the true vacuum.
The finite temperature correction to the potential is described by the c-term in Eq. (3.1), where we only include the leading T 2 -terms of the thermal integrals coming from light degrees of freedom.The subdominant −c ′ T ϕ 3 term from the bosonic degrees of freedom is omitted.The Universe stays in the false vacuum ϕ = 0 initially and it acquires a probability of decaying to the true vacuum when T < T c .The decay rate Γ(T ) ∼ e −S 3 (T )/T is defined in Eq. (2.1).At T = T c , the two vacua are degenerate and hence S 3 (T c ) = ∞; when the Universe cools down from T c , S 3 (T )/T also decreases.If there exists a zero-temperature barrier, then S 3 (T ) approaches a finite value for T → 0, resulting in lim T →0 S 3 (T )/T = ∞.Therefore, in such a case we have a U-shaped S(T ) = S 3 (T )/T curve, which efficiently suppresses the vacuum decay rate and hence results in a long FOPT.To illustrate this, we adopt m ϕ = 300 MeV, w = 900 MeV, and c = 0.11 as a benchmark, and test two different µ 3 values: m 2 ϕ /w = 100 MeV and 1.5m 2 ϕ /w = 150 MeV.The former corresponds to vanishing barrier at zero-temperature and hence S 3 (T )/T → 0 when T → 0, while the latter corresponds to a barrier that results in a U-shaped S 3 (T )/T curve, as clearly illustrated in Fig. 1.For the S 3 (T )/T calculation, we have used the semi-analytical expression [66] S 3 (T ) where (3.5) The approximation works well when u ∈ [0, 1].More discussions on the FOPT dynamics on the U-shaped S 3 (T )/T can be found in Refs.[67,68], while here we only focus on the impacts on PBH formation.For a given set of {m ϕ , w, µ 3 , c}, we apply the methodology outlined in Section 2.1 to obtain the FOPT profiles and calculate the PBH formation.We set v w = 1.0 as a fixed parameter.An illustrative example is presented in Fig. 2 with µ 3 = 154.1 MeV.In the left panel, we depict the evolution of the false vacuum fraction F (t). From the figure it can be clearly seen that delayed-decay patches initiate nucleation at later times compared to normal patches.Consequently, this results in distinct energy component evolutions, as shown in the middle panel.In both cases, the vacuum energy ρ v (t) is converted into radiation energy ρ r (t); however, delayed-decay patches exhibit higher values of ρ r (t) at later times.The right panel displays the overdensity δ(t), which progressively increases with time until it reaches the maximum value δ max = δ c .This defines the PBH formation time t pbh = 4.70×10 −5 s.The corresponding PBH mass is m pbh = 2.72×10 33 g = 1.37 M ⊙ , where M ⊙ is the solar mass.For reference, the nucleation time t n and percolation time t p of the normal patches are also indicated in the figures, and we can see that t pbh > t p > t n , which is a typical feature of this PBH mechanism.
In Fig. 3, we present the scan results in the c-µ 3 plane, with fixed values of m ϕ = 300 MeV and w = 900 MeV.By varying µ 3 and c around the center point at µ 3 = 154 MeV and c = 0.11, we explore the parameter space for FOPT and PBH formation.We observe that the PBH fraction f pbh is highly sensitive to these parameters, exhibiting a narrow band where it varies from 10 −100 (the white dotted line) to ∼ 10 7 (the upper edge of the color shaded region).This variation occurs within a small range, with µ 3 changing by only about ∼ 0.6% and c changing by approximately ∼ 9%.Therefore, f pbh is most sensitive to µ 3 , which describes the height of the zero-temperature potential barrier.This can be understood as the duration of FOPT is determined by the barrier, and f pbh is extremely sensitive to this duration, which is a main feature of the PBH formation mechanism through delayed vacuum decay.The resulting PBHs form at around 10 −5 s after the Big Bang, and have a mass of ∼ 10 33 g (∼ M ⊙ ).Within this mass region, f pbh is constrained to be ≲ 10 −2 �  due to microlensing observations [1], and the blue mesh region is already excluded.
We evaluate the FOPT parameters, namely α, β/H * , and (8π) 1/3 v w /(H * R), and observe that their contours closely align with those of f pbh .To maintain clarity in the figure, we present the corresponding FOPT parameter values for different f pbh in Table 1 instead of showing more contours in Fig. 3.We note that α at t p tends to be larger than at t n due to the greater dilution of radiation energy, as t p > t n .We do not have α ≫ 1, and hence the FOPT is not very strong.Regarding parameters related to the FOPT duration, we find that β/H * at t n and (8π) 1/3 v w /(H R) at t p decrease with longer FOPT durations.Conversely, (8π) 1/3 v w /(H * R) at t n remains insensitive to parameter variations, while β/H * at t p can even become negative.The reasons for these observations will be discussed in the following subsection, and we will demonstrate that (8π) 1/3 v w /(H * R) at t p can measure the ratio of Hubble time to FOPT duration.
The GWs generated by the FOPT are computed using numerical formulae [59], and their spectral peak falls within the 10 − 100 nHz range.This may provide an explanation for the recent excess observed in pulsar timing array (PTA) experiments, such as NANOGrav [69], CPTA [70], EPTA [71], and PPTA [72].Following the analysis conducted by the NANOGrav collaboration [73], we define "favored by the data" as a GW spectrum that matches the first 14 frequency bins of the NANOGrav-15yr dataset, highlighting this parameter space in red of Fig. 3.The parameter spaces accessible to future PTA experiments, such as IPTA [74] and SKA [75], are depicted in gray and lighter gray shades, respectively.Note that while much of the parameter space capable of explaining the NANOGrav excess has been ruled out by PBH detections, a narrow region remains unexplored, awaiting further investigation.We refrain from delving into the interplay between the PTA data, the FOPT, and PBHs, as the primary focus of this study is the relationship between the structure of particle models and PBH formation. 1

Comparison among different expansions of S(t)
In this subsection, we discuss different treatments of the Euclidean action S(t): linear and quadratic approximations, and the full expression.To begin, we need to determine an appropriate expansion point t * .Although the nucleation moment t n and percolation moment t p are potential choices, the previous subsection revealed that the region of parameter space where PBH formation is abundant has a negative value of β/H * at t p .This makes it unsuitable for expansion since it corresponds to a decreasing Γ(t) ∼ e βt at the linear level approximation.Therefore, we select t n as our expansion point.
Assuming parameters m ϕ = 300 MeV, w = 900 MeV, c = 0.11, and µ 3 = 154.1 MeV, Fig. 4 displays the evolutions of the Euclidean action and vacuum decay rate, with the argument switched to temperature for convenience.The left panel shows that the three treatments intersect at T n , but when T deviates from T n , they exhibit distinct shapes, particularly the linear approximation.We now also easily understand why β/H * at T p is 1 For more discussions on the MeV-scale FOPT explanation of the PTA data, see Refs.[19,[76][77][78][79][80][81][82][83][84][85][86][87][88][89][90][91][92] and the references therein.Due to the extensive literature on this topic, we have only included those directly relevant to the June 2023 PTA excess announcement.negative: the full S(T ) = S 3 (T )/T is U-shaped, and hence when the bubbles percolate at late time, dS/dT will become negative.The discrepancy in S(T ) results in differences in Γ(T ), as illustrated in the right panel.The exponential nucleation rate approximation fails completely in capturing the peak shape of the decay rate originating from the U-shaped action, whereas the quadratic approximation is much closer to the full expression, albeit with some visible differences.As a result, both the quadratic approximation and the full expression yield m pbh ∼ 2.7 × 10 33 g (∼ 1.36 M ⊙ ) and f pbh ∼ 2 × 10 6 , while the linear approximation has m pbh ∼ 1.6 × 10 34 g (∼ 8.1 M ⊙ ) and f pbh ∼ 10 −6 .The behavior of the linear and quadratic approximations can be understood by assuming H ∝ T n near T * , then and hence Eq. (2.14) can also treated as an expansion on temperature T , and we can have up to quadratic approximation level.Eq. (3.7) gives Since typically n ∈ [0, 2], the linear approximation always has d 2 S/dT 2 * < 0, which is qualitatively different from the full U-shaped expression that has d 2 S/dT 2 * > 0 as . The curves of α pbh (left), β/H * (middle), and (8π) 1/3 v w /(H * R) (right), scanning over µ 3 for different treatments of S(t).The same parameter setup as Fig. 5. illustrated in Fig. 4.Such a discrepancy is not pronounced when the FOPT is prompt, as it completes near T * .However, for PBH formation during a prolonged FOPT, the completion temperature deviates from T * .This deviation leads to a significant distinction between the full S(T ) and its linear time expansion.Therefore, the exponential approximation Γ(t) ∼ e βt cannot be applied to calculate PBH formation in scalar potentials with zerotemperature barriers when considering delayed vacuum decay.It's worth noting that this argument also applies to more general models with S(T ) ∝ T r and r > 1.
To further explore the differences among the three treatments, we maintain m ϕ = 300 MeV, w = 900 MeV, and c = 0.11, while varying µ 3 , and derive the curves of f pbh and m pbh in Fig. 5.The quadratic approximation exhibits a similar trend to the full expression and provides a close quantitative match in the large µ 3 region.However, for small µ 3 , the values of f pbh can differ by several orders of magnitude.In contrast, the linear approximation yields entirely different curves.This discrepancy arises because the decay rate from the linear approximation is generally higher than the other two treatments, as illustrated in the right panel of Fig. 4. Therefore, it results in faster FOPTs, earlier PBH formation, lower values of f pbh , t pbh , and m pbh .
We obtain the FOPT parameters α, β/H * , and (8π) 1/3 v w /(H * R) at t n and t p , and present their dependence on µ 3 in Fig. 6.The left panel demonstrates the FOPT is mild in three different treatments.In the middle panel, we observe that β/H * derived at t n or t p becomes negative for sufficiently large µ 3 , indicating it cannot be interpreted as a ratio of time scales.Conversely, the right panel demonstrates that (8π) 1/3 v w /(H * R) at t p is positive and decreases with increasing FOPT duration in both the full expression and quadratic approximation cases.Therefore, it serves as a measure of the ratio between the Hubble time and the FOPT time scales.
Before closing this subsection, we emphasize the importance of verifying Eq. (2.21) when studying the formation of PBHs resulting from delayed vacuum decay.The FOPTs being considered here have a long duration, necessitating the confirmation of their completion.For instance, in Figs. 5 and 6, the curves representing the full expression evaluation end at µ 3 ≈ 154.18 MeV, indicating that the FOPT cannot be fully accomplished beyond this point.If Eq. (2.21) is not satisfied, the physical volume of the false vacuum does not decrease at percolation, even if the volume fraction decreases; consequently, the FOPT remains incomplete [65].4 The Z 2 -symmetric singlet extension of the SM This section discusses a more complicated and realistic model with double-field FOPT dynamics: the Z 2 -symmetric xSM [93][94][95].The scalar sector contains the SM Higgs doublet H and a Z 2 -odd real gauge singlet s, and the joint potential at finite temperature can be written as where h and s are the neutral Higgs and singlet background field, respectively, and the coefficients are derived from the high-temperature expansion of the thermal collections, with g, g ′ , and y t the SU (2) L , U (1) Y gauge couplings, and the top Yukawa, respectively. 2e require λ h > 0, λ s > 0 and √ λ h λ s + λ hs > 0 to ensure the potential bounded from below, and such that at T = 0 there exists two local minima at (h, s) = (v, 0) and (0, w), with v = −µ 2 h /λ h and w = −µ 2 s /λ s , and the former is the true vacuum, preserving the Z 2 symmetry.The two minima are separated by a tree-level barrier induced by the λ hs term.Since s does not acquire a vacuum expectation value at zero-temperature, it does not mix with h, and the physical masses of the scalars can be easily obtained as, Given m h = 125 GeV and v = 246 GeV, we only have three free parameters, which can be adopted as {m s , λ hs , λ s }.Note that m s > m h /2 is required by phenomenology to evade the h → ss bound from collider experiments.It is well-known that when T > 0 the Z 2 -symmetric xSM can realize a two-step phase transition via the trajectory (0, 0) → (0, w * ) → (v * , 0), in which the second step is a firstorder electroweak phase transition.Under the thermal potential Eq. (4.1), this is the only possible FOPT pattern, and its necessary condition is [97]  derived from the existence condition of T c .Within this region, we resolve the FOPT dynamics of the model with the assistance of the Python package cosmoTransitions [98], which helps to evaluate S 3 (T )/T .In Fig. 7, the left panel displays the contours of the potential V (h, s, T ) at T = 0, with fixed values of m s = 218 GeV, λ s = 2, and λ hs = 1.108.Two distinct local minima, located at (0, w) (red point, false vacuum) and (v, 0) (blue point, true vacuum), are separated by a barrier resulting from the λ hs term.This barrier is clearly illustrated in the middle panel, where the field values are taken along the green dashed line in the left panel.The barrier gives rise to a U-shaped S 3 (T )/T , as depicted by the blue line in the right panel of Fig. 7. Additionally, we show the corresponding linear and quadratic approximations of S 3 (T )/T .We again see the linear expansion fails to capture the shape of the complete expression.
To provide a comprehensive overview of the parameter space, we conducted a 2dimensional scan for m s and λ hs while keeping λ s as different fixed values, as depicted in the top panel of Fig. 8.In the top-left panel, we specifically plotted the parameter space for f pbh ∈ [10 −100 , 1] and observed that it consists of narrow regions resembling "lines".These "lines" have endpoints determined by ensuring the triviality bound [99].It is well known that large scalar coupling constants are in general required to realize the FOPT.However, such large scalar coupling constants are strongly constrained by the triviality bound [100].We take into account the constraint from the triviality bound by employing mass dependent beta functions [101].We require that the Landau pole of the scalar couplings remains above 10 TeV to satisfy the triviality bound.
The top-right panel of Fig. 8 shows the details of the parameter space for λ s = 2, near the vicinity of m s = 218.4GeV and λ hs = 1.106.The PBH formation parameter space corresponds to strong GW production, and we have evaluated the corresponding GW spectra and verified that they are easily detected by the future LISA [102], TianQin [103], Taiji [104], and DECIGO [105] interferometers.Moreover, we illustrate the deviations in the hhh and hW W/hZZ couplings as purple and green dashed contours, respectively.These deviations can be precisely measured at future collider experiments like HL-LHC [106], CEPC [107], ILC [108], and FCC-ee [109].Those future collider experiments can offer complementary and correlated probes to this PBH formation mechanism.
In the bottom-left and -right panels of Fig. 8, we present the variations of f pbh and m pbh with respect to λ hs , considering different treatments for S(t) while keeping m s = 218 GeV and λ s = 2 fixed.We observe that the PBH fraction exhibits rapid changes as a function of λ hs , primarily due to the significant influence of the barrier height on the probability of delayed vacuum decay.Conversely, the PBH mass ∼ 10 28 g (∼ 10 −5 M ⊙ ) changes gradually as it is predominantly determined by the temperature of the FOPT.It can be seen in the figure that the linear expansion fails to accurately describe the dynamics of PBH formation, as expected.

Conclusion
We explore the PBH formation through the delayed vacuum decay in the slow FOPTs, also known as the "late-blooming" mechanism, from a model-building perspective.We for the first time investigate this mechanism in models with zero-temperature potential barriers, taking the polynomial potential and the Z 2 -symmetric xSM as two examples.Our findings reveal that a U-shaped Euclidean action significantly prolongs the phase transition, leading to abundant PBH formation with sufficiently high barriers.Unlike models with classically conformal invariance [18][19][20], the FOPT in these models is moderate with an O(1) α parameter.Therefore, our study demonstrates a general type of particle models that can realize the PBH formation without requiring ultra-supercooling.The PBH abundance is highly sensitive to model parameters, which is a distinctive feature of this mechanism and has been demonstrated in previous literature [16][17][18][19].Additionally, we explicitly show that the exponential nucleation approximation does not work in our models due to the U-shaped action, emphasizing the necessity for the full expression of the decay rate in particle-based studies.
Although our focus in this study is on particle models with zero-temperature barriers, the discussions presented in Section 2, especially in Section 2.2 regarding the validity of the exponential approximation, also have implications for more general models, including the classically conformal ones.In those models, where S 3 (T )/T ∼ T r with r > 1, satisfying Eq. (2.23) becomes challenging if one desires a small β/H * to facilitate PBH formation.Consequently, the exponential approximation may also fail in PBH studies of such models, resulting in significant errors of several orders of magnitude.It is also important to note that for FOPTs characterized by long-lasting or ultra-supercooling regimes, it is necessary to verify the fulfillment of Eq. (2.21) to ensure the completion of the phase transition, which is not realized in some relevant literature.
Our work can be extended in several ways.While we have chosen v w = 1 as a benchmark value, it is important to consider the friction force induced by plasma particles, which typically leads to bubble wall velocities smaller than 1 [110][111][112][113][114][115].By varying v w as different constants in our calculations, we find that lower values of v w result in later percolation and larger PBH abundance, with differences spanning several orders of magnitude.Therefore, accurately determining the velocity within a specific model is crucial for deriving the PBH abundance.Different definitions of overdensity δ(t) exist in the literature.In our research, we adopt the definition given by Eq. (2.6), which leads to a decrease in overdensity after reaching its maximum value δ max .This is consistent with Refs.[14,17].However, an alternative definition of δ(t), as used in Ref. [11], is given by δ(t) = ρ in (t) ρ out (t) a 4 in (t) a 4 out (t) − 1, which tends to be constant after reaching δ max .It would be interesting to further explore and compare these two definitions from a theoretical perspective.Furthermore, the finite temperature potential expression can be improved by including the full one-loop thermal corrections and beyond, and the calculation framework can be enhanced to incorporate more detailed considerations of PBH formation possibilities [15,16], or by conducting numerical simulations of FOPTs [116,117].Given the sensitivity of PBH formation to FOPT features in this mechanism, a more comprehensive treatment is expected to significantly impact the PBH abundance and the viable parameter space of a given model.However, the key qualitative conclusions drawn from our research will remain unchanged, including the relationship between zero-temperature barriers and PBH formation, as well as the breakdown of the exponential approximation.

Figure 2 .
Figure 2. The FOPT and PBH profiles for m ϕ = 300 MeV, w = 900 MeV, c = 0.11, and µ 3 = 154.1 MeV.t n , t p , and t pbh are the nucleation, percolation, and PBH formation time, respectively.H c and ρ c denote the Hubble constant and radiation energy density at t c , respectively.Left: the false vacuum fraction F (t). Middle: the energy components ρ r,v (t).Right: the overdensity δ(t).

Figure 3 .
Figure3.The scan result in the c-µ 3 plane for m ϕ = 300 MeV and w = 900 MeV.The white region cannot realize the FOPT.The white solid, dashed, and dotted lines represent f pbh = 1, 10 −10 , and 10 −100 , respectively.The red, gray, and lighter gray regions correspond to the parameter space that produces GWs favored by the NANOGrav excess, reached by future IPTA and SKA, respectively.The blue mesh region is excluded by microlensing experiments.

Figure 7 .
Figure 7. Left: the zero-temperature potential contours V (h, s, 0), with the blue and red points denoting the true and false vacua, respectively.Middle: the field value along the green dashed line of the left panel.Right: the S 3 (T )/T curves for different treatments.The parameters are adopted as m s = 218 GeV, λ s = 2, and λ hs = 1.108.
[16]ughly measures the duration of FOPT over the Hubble time scale, and according to the calculations in Appendix A, it can be replaced by O(10) × (β/H * ) −1 .In addition, we can generally assume H ∝ T n and S ∝ T r near T * , and then Eq. (2.22) reduces to PBH formation requires the right-hand side (r.h.s.) of Eq. (2.23) to be ∼ 0.2 × 10 ∼ 2. However, in ultra-supercooled and mild FOPTs with vacuum and radiation domination, the left-hand side (l.h.s.) typically takes values of n ≈ 0 and n ≈ 2 respectively, while a realistic particle model usually has r > 1.Consequently, it is challenging to satisfy the exponential approximation condition in the parameter space favored by the PBH mechanism.If the l.h.s. is not sufficiently small, including additional terms in the expansion, such as the quadratic ζ-term[16], becomes necessary.

Table 1 .
−10, and 10 −100 , respectively.The red, gray, and lighter gray regions correspond to the parameter space that produces GWs favored by the NANOGrav excess, reached by future IPTA and SKA, respectively.The blue mesh region is excluded by microlensing experiments.The FOPT parameters corresponding to the scan in Fig.3.The subscript "n" and "p" represent the values derived at nucleation (t n ) or percolation (t p ), respectively.