# Revisiting formation of primordial black holes in a supercooled first-order phase transition

Marcos M Flores, Alexander Kusenko, Misao Sasaki

► To cite this version:

Marcos M Flores, Alexander Kusenko, Misao Sasaki. Revisiting formation of primordial black holes in a supercooled first-order phase transition. Physical Review D, 2024, 110 (1), pp.015005. 10.1103/PhysRevD.110.015005. [hal-04491710](https://hal.science/hal-04491710v1)

HAL Id: hal-04491710

https://hal.science/hal-04491710v1

Submitted on 29 Aug 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Revisiting formation of primordial black holes in a supercooled first-order phase transition

Marcos M. Flores<sup>1,2</sup> Alexander Kusenko<sup>1,3</sup> and Misao Sasaki<sup>3,4,5</sup>

<sup>1</sup>Department of Physics and Astronomy, University of California, Los Angeles, California 90095-1547, USA

<sup>2</sup>Laboratoire de Physique de l'École Normale Supérieure, ENS, Université PSL, CNRS, Sorbonne Université, Université Paris Cité, F-75005 Paris, France

<sup>3</sup>Kavli Institute for the Physics and Mathematics of the Universe (WPI), UTIAS The University of Tokyo, Kashiwa, Chiba 277-8583, Japan

<sup>4</sup>Center for Gravitational Physics and Quantum Information, Yukawa Institute for Theoretical Physics, Kyoto University, Kyoto 606-8502, Japan

<sup>5</sup>Leung Center for Cosmology and Particle Astrophysics, National Taiwan University, Taipei 10617, Taiwan

<figure><img src="image_2.png" /></figure>

(Received 4 March 2024; accepted 4 June 2024; published 3 July 2024)

We reexamine production of primordial black holes in a supercooled phase transition. While a mere overdensity associated with a surviving false-vacuum patch does not imply formation of a black hole, it is possible for such a patch to evolve and create a black hole, thanks to the gradient energy stored in the bubble wall.

DOI: 10.1103/PhysRevD.110.015005

# I. INTRODUCTION

Primordial black holes (PBHs) are an attractive candidate for dark matter [1–52] and also could be responsible for various astrophysical phenomena, including seeding supermassive black holes [14,53–56], synthesis of heavy elements [57–59], explaining nonrepeating fast radio bursts [57,60,61], and G objects discovered in the Galactic Center [62].

A relatively slow, supercooled first-order phase transition, in which the false-vacuum “islands” can linger, possibly creating inhomogeneities, could provide the conditions for the formation of PBHs [63–70]. We will examine the viability of this scenario, while avoiding unreliable and incorrect assumptions that have appeared in the literature. We will ultimately show that there is a way to form PBHs from shrinking domains of the false vacuum, thanks to the gradient energy stored in the bubble wall.

# II. INAPPLICABILITY OF A SIMPLE OVERDENSITY CONDITION

The physics of supercooled phase transitions is well understood. False-vacuum decay, being a stochastic

process, allows for the possibility that large patches of space-time remain in false vacuum, even after the percolation of bubbles in the majority of the Universe. Regions that transitioned to true vacuum subsequently reheated and have become radiation dominated. Eventually, the delayed patch may also transition into radiation. Given the difference in transition times, the radiation density in the majority of the Universe has redshifted more relative to the delayed patch. The relative difference in radiation densities generates a local overdensity.

In determining whether an overdensity of this kind might lead to formation of a PBH, it may be tempting to apply a simple overdensity criterion for black hole formation adopted from the inflationary PBH formation discussion. However, in most cases this approach cannot succeed. Indeed, let us define a density contrast

$$ \delta(x) := \frac{\rho(x) - \bar{\rho}}{\bar{\rho}}. \qquad (1) $$

Let us suppose that some region of size L with a density contrast δ remained in the false-vacuum phase at a relatively late time in the transition, and it converted to the true vacuum at a time when $\bar{\rho}$ has decreased to a lower value than $\rho(x)$ averaged over the region L.

Is there a critical value $\delta_c$ such that a patch with $\delta > \delta_c$ collapses to a black hole? In a flat universe, the answer is no. Indeed, if a patch of radiation dominated Universe has a higher density, it expands in accordance with the Friedmann equations, at the rate that is higher for larger values of $\rho$.

<footnote>Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI. Funded by SCOAP³.</footnote>
The difference with the widely discussed scenarios for PBH formation due to collapse of curvature perturbations is rooted in the fact that a curvature perturbation on some superhorizon scales makes the Universe locally nonflat. Let us briefly review the perturbation scenario and juxtapose it with the phase transition.

In inflation-based formation mechanisms of PBHs, cosmological inflation produces a locally perturbed region with comoving size larger than the horizon size. The density contrast on this comoving slice can be related to the curvature K associated with the comoving slice [28]. In particular, K is a function of the conserved comoving curvature perturbation, typically denoted $R_c$ [71].

Approximating $K(r) = K = \text{const}$, one obtains the Friedmann equation

$$H^2 + \frac{K}{a^2} = \frac{8\pi G}{3} \rho, \qquad (2)$$

where $H = \dot{a}/a$, and $a$ is the scale factor. One can now relate the density contrast to the curvature K,

$$\Delta := \frac{\rho - \bar{\rho}}{\bar{\rho}} = \frac{3K}{8\pi G \bar{\rho} a^2} = \frac{K}{H^2 a^2}. \qquad (3)$$

During radiation domination, $\bar{\rho} \propto a^{-4}$, which signals that it is the curvature perturbation that induces the density perturbation. Now the perturbed patch of the Universe evolves as a nonflat universe with a curvature.

Once the perturbation grows to $\Delta = 1$, the expansion stops, and the patch can collapse to a black hole. Specifically, if the size of the perturbation is larger than the Jeans length, then collapse will occur. Using [28]

$$\frac{c_s^2 k^2}{a^2} = H^2, \qquad (4)$$

one obtains

$$1 = \Delta(t_c) = \frac{K}{c_s^2 k^2}, \qquad (5)$$

i.e., $K = c_s^2 k^2$, where $t_c$ is the time when the density perturbation $\Delta$ is one. The criterion for PBH formation is that the comoving slice density contrast reenters the Hubble horizon with a value greater than $\delta_c = c_s^2$. This physical understanding of the critical density criterion will be sufficient for our purposes, but much work as been done to refine the value of $\delta_c$ needed for collapse [72–76].

As discussed in the context of the inflationary scenario, all PBH formation scenarios require a clear criterion that demarcates the formation of a black hole. The simplest requirement, which we will call the Schwarzschild criterion, requires that the mass enclosed in a given volume fit within that same mass’s Schwarzschild radius. This can be generalized to the so-called hoop conditions [77], but spiritually resembles the Schwarzschild criterion, so we

will not differentiate between the two. On the other hand, as we have seen in the inflationary context, there are other conditions that signal the formation of a black hole.

While the critical density criterion established in inflationary scenarios appears to be applicable generically, there are numerous subtleties that forbid its application to overdense regions that did not originate from inflationary perturbations. The generation of large subhorizon overdensities, including those larger than the critical density criterion, do not signal the formation of a black hole. The existence of structure in our present-day Universe makes this point obvious. Naturally, subhorizon scale black holes can form but the criterion for formation should instead be replaced by the Schwarzschild criterion.

It may also seem reasonable to apply the critical density criterion to horizon-sized density perturbations. This too proves problematic. In the inflationary picture, the generation of an overdensity is linked to the existence of a curvature perturbation. It is this curvature perturbation that provides a contracting solution for the Friedmann equations and the subsequent formation of a black hole. Without this curvature component, the equations of motion describing any Hubble-sized patch containing radiation provide only expansionary solutions. Via this argument, the generation of a horizon-sized radiation overdensity is not sufficient to signal the formation of a black hole.

In this paper, we will demonstrate that PBHs can form from supercooled phase transitions. Instead of relying on the critical collapse criteria, we will instead use the Schwarzschild criterion. A false-vacuum region surrounded by true vacuum and radiation naturally implies the existence of a domain wall separating the two regions. The pressure differential between these two regions pushes the domain wall boundary inward. During this contraction, the energy density of the volume remains constant. However, the domain wall energy per unit volume scales as $1/R$. Eventually, the domain walls will fit within the Schwarzschild radius associated with the false-vacuum region leading to the formation of a PBH.

In Sec. III, we review false-vacuum decay and specify the scalar potential used in our study. In Sec. IV, we will discuss the background evolution of the Universe and the dynamics of bubble formation and evolution. In Sec. V, the domain wall equation of motion will be discussed as well as the formation of PBHs in this scenario. Section VI will discuss the probability of generating large supercooled regions and present a formula for the PBH abundance. Finally, in Sec. VII we will provide benchmark points of parameter space that allow for the generation of PBHs relevant for explaining all of dark matter (DM) or which may be consistent with microlensing observations.

### III. FALSE-VACUUM DECAY RATE

False-vacuum decay is a well-studied phenomena with a plethora of applications to early Universe physics [78–81].
The tunneling rate from the false-vacuum state to the true vacuum will be the primary focus of this section. The tunneling rate determines the nucleation rate of true-vacuum bubbles. Tunneling can occur either thermally or via quantum tunneling. In the case of our interest, the thermal tunneling rate is appropriate,

$$ \Gamma = T^4 \left( \frac{S_3}{2\pi T} \right)^{3/2} e^{-S_3/T}, \qquad (6) $$

where $\Gamma$ denotes the tunneling rate per unit volume, and $S_3$ is the three-dimensional Euclidean action.

In this work, we will assume the nucleated bubbles are spherically symmetric and produced by a single scalar field. Under these assumptions, the three-dimensional, O(3)-symmetric Euclidian action is given by

$$ S_3 = 4\pi \int dr r^2 \left[ \frac{1}{2} \left( \frac{d\phi}{dr} \right)^2 + V(\phi) \right]. \qquad (7) $$

Extremizing the action given by Eq. (7) requires solving the equations of motion,

$$ \frac{d^2\phi}{dr^2} + \left(\frac{2}{r}\right) \frac{d\phi}{dr} = \frac{\partial V}{\partial\phi} \qquad (8) $$

subject to the boundary condition that

$$ \phi(r \to \infty) = \phi_-, \qquad \left. \frac{d\phi}{dr} \right|_{r=0} = 0, \qquad (9) $$

where $\phi_-$ corresponds to the value of $\phi$ at the true minimum of the potential. We have also assumed that the metastable minimum lies at the origin. The solution of Eq. (8), called a thermal instanton, may be substituted into the action to find the tunneling rate.

For our purposes, we will examine a simple scalar theory whose properties allow for strong phase transitions. In particular, we will consider a finite-temperature effective potential of the form

$$ V(\phi, T) = \frac{1}{2} m^2(T) \phi^2 - \frac{\delta(T)}{3} \phi^3 + \frac{\lambda(T)}{4} \phi^4. \qquad (10) $$

It is useful to define the dimensionless parameter,

$$ \kappa(T) \equiv \frac{\lambda(T) m^2(T)}{\delta^2(T)}, \qquad (11) $$

and introduce the dimensionless field and the coordinates such that $\phi \to m^2\phi/\delta$ and $r \to \rho/m$. The potential is then reexpressed as

$$ \begin{aligned} V(\phi, T) &= \frac{m^6(T)}{\delta^2(T)} \hat{V}(\phi, T), \\ \hat{V}(\phi, T) &= \frac{1}{2}\phi^2 - \frac{1}{3}\phi^3 + \frac{\kappa(T)}{4}\phi^4. \end{aligned} \qquad (12) $$

Here the parameter $\kappa(T)$ lies between $-\infty < \kappa(T) < \kappa(T_c) = \kappa_c = 2/9$, where $T_c$ corresponds to the value of $T$ where the two minima of the potential are degenerate. The action $S_3$ is written as

$$ \begin{aligned} S_3(T) &= \frac{m^3(T)}{\delta^2(T)} \hat{S}_3(T), \\ \hat{S}_3(T) &= 4\pi \int d\rho \rho^2 \left[ \frac{1}{2} \left( \frac{d\phi}{d\rho} \right) + \hat{V}(\phi, T) \right]. \end{aligned} \qquad (13) $$

Determining the action for any $T$, or equivalently any $\kappa$, amounts to first numerically evaluating Eq. (8), subject to its boundary conditions, Eq. (9), then evaluating the action for this solution. This is well-trodden ground, and numerous software packages have been developed that are capable of tackling precisely this problem [82–85]. For simplicity, we will use fits constructed in Ref. [86] for the one-parameter potential given in Eq. (12). The resulting action may then be approximated by

$$ S_3(T) \approx \frac{m^3(T)}{\delta^2(T)} \begin{cases} \frac{2\pi}{3(\kappa-\kappa_c)^2} \bar{B}_3(\kappa) & \kappa > 0, \\ \frac{2\pi}{2} \left( \frac{1+\exp(-|\kappa|^{-1/2})}{1+|\kappa|/\kappa_c} \right) & \kappa < 0, \end{cases} \qquad (14) $$

where $\bar{B}_3(\kappa)$ is a fitting function defined in Appendix A and which is normalized so that $\bar{B}_3(\kappa_c) = 1$. This fitting function reproduces the analytic solutions obtained in the $\kappa \to 0$ limit and also accommodates the correct $\kappa \to \kappa_c$ behavior.

To move forward, we must specify the functional form of $m(T)$, $\delta(T)$, and $\lambda(T)$. Following Ref. [87], we define

$$ \begin{aligned} m^2(T) &= 2D(T^2 - T_0^2), \\ \delta(T) &= 3(ET + A), \\ \lambda(T) &= \lambda. \end{aligned} \qquad (15) $$

It is convenient to define the zero-temperature minimum $\phi_-(T=0) \equiv v$. This enables us to express the temperature $T_0$ as

$$ T_0^2 = v^2 \left( \frac{\lambda - 3A/v}{2D} \right). \qquad (16) $$

We can also express zero-temperature mass $m_\phi$ as

$$ m_{\phi}^2 = \frac{\partial^2 V}{\partial \phi^2}(v, 0) = 2\lambda v^2 - 3Av. \qquad (17) $$
False-vacuum decay can only occur once the minimum $\phi_-(T) < 0$. This occurs at the critical temperature $T_c$,

$$T_c = \frac{1}{\lambda D - E^2} \left[ AE + \sqrt{\lambda D(A^2 + (\lambda D - E^2)T_0^2)} \right]. \quad (18)$$

Additionally, the value of the field $\phi_-(T_c) \equiv \phi_c$ is given by

$$\phi_c = T_c \left[ \frac{2(E + A/T_c)}{\lambda} \right]. \quad (19)$$

For the specific set of functions $\{m(T), \delta(T), \lambda(T)\}$ in (15), the parameter $\kappa$ is given by

$$\kappa(T) = \frac{2\lambda D(T^2 - T_0^2)}{9(ET + A)^2}. \quad (20)$$

Given that $\kappa$ alone determines the nondimensional action $\hat{S}_3$, it is a simple exercise to demonstrate that the action depends only on the parameters $D$, $E$, and $\lambda$ and the ratios $A/v$ and $T/v$. The ratio $A/v$ plays a significant role in the strength of the phase transition allowed by the potential Eq. (12). Therefore, for the remainder of the paper, we will fix [87]

$$E = 1/16, \quad D = 2.725, \quad \lambda = \frac{1}{8} + \frac{3A}{2v}, \quad (21)$$

where $\lambda$ is chosen so that the zero-temperature mass $m_\phi$ is given by $m_\phi = v/2$. The requirement that the vacuum energy vanishes at $T = 0$ gives

$$\rho_V = \frac{v^4}{4} \left( \lambda - \frac{2A}{v} \right) = \frac{v^4}{4} \left( \frac{1}{8} - \frac{A}{2v} \right). \quad (22)$$

Note that our choice of the parameter $D$ differs from Ref. [87] in order to ensure that a period of supercooling occurs.

Another important aspect of our scenario is the physical properties of the domain walls that will separate regions of true and false vacua. For a temperature-dependent potential, this calculation is nontrivial, especially considering the fact that temperature changes across the domain wall.

We will approximate the size of the domain wall $\ell$ and its surface energy density $\sigma$ using the values at $T = T_c$, namely, when the minima of the potential are degenerate. This is also known as the “thin-wall regime.” The solution for the potential equation (10) in the degenerate vacuum limit is

$$x = \int_0^\phi \frac{d\varphi}{[2V(\varphi, T_c)]^{1/2}}, \quad (23)$$

which gives

$$\phi(x) = -\frac{1}{2}\phi_c \tanh\left(\frac{\lambda\phi_c}{\sqrt{8}}x\right). \quad (24)$$

Therefore, the characteristic scale associated with the domain wall is

$$\ell \sim \frac{8^{1/2}}{\lambda\phi_c}. \quad (25)$$

The surface energy density is given by

$$\begin{aligned} \sigma &= \int dx \left[ \frac{1}{2} \left( \frac{d\phi}{dx} \right)^2 + V(\phi, T_c) \right] \\ &= \int_0^{\phi_c} [2V(\phi, T_c)]^{1/2} d\phi = \frac{\phi_c^3}{6} \left(\frac{\lambda}{2}\right)^{1/2}. \end{aligned} \quad (26)$$

These quantities will play an important role in our discussion of the domain wall equation of motion and the conditions for the PBH formation.

## IV. BACKGROUND EVOLUTION AND BUBBLE DYNAMICS

The formation of PBHs is determined by the evolution of energy densities within different causal regions of the Universe. In the supercooled framework, there are two relevant energy components, namely, radiation and the false-vacuum energy

$$\rho_{\text{tot}} = \rho_R + \rho_V, \quad (27)$$

where generally the radiation component will include the energy density stored in relativistic bubble-wall propagation, scalar waves, and plasma excitation [69]. For simplicity, we will group these effects into a singular energy component $\rho_R$ and will leave a holistic study for future work.

The evolution of the scale factor, radiation energy density, and false-vacuum energy density are controlled by the differential equations

$$H^2 = \left(\frac{1}{a}\frac{da}{dt}\right)^2 = \frac{1}{3M_{\text{Pl}}^2}[\rho_R(t) + \rho_V(t)], \quad (28)$$

$$\dot{\rho}_R + 4H \frac{d\rho_R}{dt} = -\frac{d\rho_V}{dt}. \quad (29)$$

For the sake of evaluation, it is useful to define the dimensionless time and energy densities

$$\tau \equiv \frac{H_{\text{eq}}}{\sqrt{2}} t, \quad \hat{\rho}_i(\tau) = \frac{2\rho_i(\tau)}{3M_{\text{Pl}}^2 H_{\text{eq}}^2}, \quad (30)$$

where $H_{\text{eq}}$ is the Hubble parameter at equality, i.e., when $\rho_R(\tau_{\text{eq}}) = \rho_V(\tau_{\text{eq}})$. We will leave a detailed derivation of the nondimensionalized differential equations in Appendix B.
It is also useful to define a temperature $T_{eq}$, associated with the equality time, namely,

$$H_{\text{eq}}^2 = \frac{2\pi^2}{90} g_{\star}(T_{\text{eq}}) \frac{T_{\text{eq}}^4}{M_{\text{Pl}}^2} = \frac{2\rho_{\text{V}}}{3M_{\text{Pl}}^2}. \quad (31)$$

From here, we see that the equality temperature is determined by $v$ and $A/v$, assuming $\lambda$ is fixed by Eq. (21).

While the energy density $\rho_V$ is initially constant, as soon as the temperature falls below $T_c$ thermal tunneling and bubble nucleation begin to occur. Once a bubble forms at $t'$, its radius at time $t$ is given by [88]

$$R(t; t') = \int_{t'}^{t} d\bar{t} \frac{v_w(\bar{t})}{a(\bar{t})}. \quad (32)$$

The expected volume of true-vacuum bubbles per unit volume of space at time $t$, $I(t)$, dictates the probability that a given point of space remains in the false-vacuum state at time $t$. Given a bubble of true vacuum with coordinate radius $R(t; t')$, then [88]

$$I(t) = \frac{4\pi}{3} \int_{t_c}^{t} \Gamma(t') a^3(t') R^3(t; t') dt', \quad (33)$$

where $t_c$ corresponds to the time in which $T = T_c$ or, in other words, when phase transitions become energetically feasible. The time dependence of vacuum energy density $\rho_V(t)$ is then given by the product

$$\rho_V(t) = \rho_{V,eq} \exp[-I(t)], \quad (34)$$

where $\rho_{V,eq}$ is given by the constant value presented in Eq. (22).

Equations (28)–(34) describe a two-dimensional coupled, integrodifferential equation. In Appendix B we outline how to transform this system into a set of first-order differential equations that may be evaluated with any off-the-shelf numerical method software.

## V. FORMATION OF PBHS

The formation of PBHs will follow the collapse of a false-vacuum patch due to the pressure differential between the true- and false-vacuum regions. This scenario was first explored in Refs. [63,64] and later refined by Ref. [65]. The formalism established by Refs. [63–65] has also been utilized to describe the possible formation of PBHs from bubble nucleation during inflation [89,90].

References [63–65] consider a circular patch of false vacuum surrounded by a region of true vacuum. This physical scenario presents two apparent paradoxes.

First, for a sufficiently large false-vacuum patch, an observer deep within the false-vacuum region would unambiguously observe an inflating space-time. However, an observer outside the patch would observe an unstable,

collapsing region as the pressure differential between the two regions induces an inward force.

Second, if the two observers sat arbitrarily close to each other, but separated by a domain wall between the two regions, each would again observe different phenomena. The inner, false-vacuum observer would see the domain wall radius of curvature increasing. However, the outer true-vacuum observer would not see an increase in the radius of curvature for the same reasons as discussed in the previous paragraph. Given that general relativity will ensure continuity between the two regions, the observers must agree on the observed domain-wave radius of curvature.

The resolution to these apparent paradoxes was first identified in Ref. [63] and further explored in Ref. [65]. General relativity and the non-Euclidean geometry of space-time allow for a consistent solution across the de Sitter-Schwarzschild barrier. References [63–65] highlighted the possibility that bubble universes might detach from the original space-time. Observers inside this detached region would observe an inflating universe, while those outside would simply see a black hole.

On the other hand, small regions of false vacuum will never reach a state of eternal inflation. Alternatively, these smaller false-vacuum regions may simply collapse into “conventional” black holes after their domain walls, and their associated energy, enters into their respective Schwarzschild radii. It is this scenario that is applicable to the formation of PBHs from supercooled phase transitions.

The equation of motion for the domain wall dividing the true- and false-vacuum regions can be derived using the Israel junction conditions [65,91,92]. In particular, the junction conditions lead to an equation of motion for the wall radius [65]

$$M = \frac{4\pi}{3} \rho_V r^3 - 8\pi^2 G \sigma^2 r^3 + 4\pi\sigma r^2 \left[ 1 - \left( \frac{8\pi G}{3} \rho_V \right) r^2 + r^2 \right]^{1/2}, \quad (35)$$

where $\cdot \equiv d/d\tau$ and $\tau$ is the proper time derivative as measured along the domain wall trajectory. It is convenient to define the Hubble scales $H_V$ and $H_\sigma$ as

$$H_V^2 = \frac{\rho_V}{3M_{\text{Pl}}^2} = \frac{H_{\text{eq}}^2}{2}, \quad (36)$$

$$H_{\sigma} = \frac{\sigma}{2M_{\text{Pl}}^2}, \quad (37)$$

as well as the ratio of the two,

$$\eta \equiv \frac{H_{\sigma}}{H_{V}}. \quad (38)$$
For the specific potential considered in Sec. III,

$$ \eta \sim \frac{v^3/2M_{\text{Pl}}^2}{\sqrt{v^4/3M_{\text{Pl}}^2}} \sim \frac{v}{M_{\text{Pl}}} \ll 1. \quad (39) $$

We will also introduce

$$ \chi^2 \equiv H_V^2 + H_\sigma^2 = H_V^2(1 + \eta^2). \quad (40) $$

This allows us to define the nondimensional quantities

$$ z^3 \equiv \frac{\chi^2}{2GM} r^3 \quad (41) $$

and

$$ \tau' \equiv \frac{\chi^2 \tau}{2H_\sigma}. \quad (42) $$

Following Ref. [65] this allows us to write Eq. (35) as

$$ \left(\frac{dz}{d\tau'}\right)^2 + U(z) = E, \quad (43) $$

where

$$ \begin{aligned} E &= \frac{4H_{\sigma}^{2}}{(2GM)^{2/3}\chi^{8/3}} \\ &= \frac{4\eta^{2}}{(2GMH_{V})^{2/3}(1 + \eta^{2})^{4/3}} \end{aligned} \quad (44) $$

and

$$ U(z) = -\left(\frac{1-z^3}{z^2}\right)^2 - \frac{\gamma^2}{z}, \quad \gamma \equiv \frac{2\eta}{(1+\eta^2)^{1/2}}. \quad (45) $$

Note that $0 \le |\gamma| \le 2$ and that only $E$ depends on the mass $M$. Therefore, the wall equation of motion, Eq. (35), has been reduced to the motion of a particle moving in one dimension under the influence of the potential given above.

The implicit solution to Eq. (43) is

$$ \tau' = \int_z \frac{dz'}{\sqrt{E - U(z')}}, \quad (46) $$

The domain of integration is determined by the initial conditions of the bubble wall and can lead to dramatically different physical outcomes.

Fortunately, we are only interested in the collapse timescale and so the solution Eq. (46) is sufficient. However, to establish an intuition for the process of collapse we will first examine properties of the potential, Eq. (45). The potential $U(z)$ is both negative and concave down for all $z$. Furthermore, the potential has one

maximum $U_m$, located at

$$ z_m^3 = \frac{1}{2} \left\{ \left[ 8 + \left(1 - \frac{1}{2}\gamma^2\right)^2 \right]^{1/2} - \left(1 - \frac{1}{2}\gamma^2\right) \right\}, \quad (47) $$

which gives

$$ U_m \equiv U(z_m) = - \frac{3(z_m^6 - 1)}{z_m^4}. \quad (48) $$

The maximum of the potential also allows us to define the critical mass $M_{cr}$ as $E(M_{cr}) = U_m$ given by

$$ M_{cr} = \bar{M} \frac{\gamma^3 z_m^6 (1 - \gamma^2/4)^{1/2}}{3\sqrt{3} (z_m^6 - 1)^{3/2}}, \quad (49) $$

where

$$ \bar{M} \equiv \frac{4\pi M_{\text{Pl}}^2}{H_V}. \quad (50) $$

Observe that, in the limit that $\gamma \to 0$, $M_{cr} \to \bar{M}$. In connection with our original parameters, the $\gamma \to 0$ limit is equivalent to $v \to 0$. Further, $\bar{M}$ acts as the minimum mass needed for a region to be able to collapse into a black hole.

Finally, it is convenient to express the Schwarzschild radius $r = r_s = 2GM$ and the de Sitter horizon $r = H_V^{-1}$ in terms of $z$. It is a simple exercise to show that

$$ z_s = \frac{\gamma^2}{|E|}, \quad (51) $$

$$ z_H = \frac{\sqrt{|E|}}{\gamma(1 - \gamma^2/4)^{1/2}}. \quad (52) $$

The key properties of the potential, Eq. (45), can be observed in Fig. 1. With Eq. (39) in mind, we note that

$$ z_s(\gamma \to 0) = (M/\bar{M})^{2/3}, \quad (53) $$

$$ z_H(\gamma \to 0) = (M/\bar{M})^{-1/3}. \quad (54) $$

This fact and the maximum mass defined in Eq. (49) limit the values of $z_s$ and $z_H$ in the $\gamma \to 0$ or, equivalently, the $\eta \to 0$ limit.

For the purposes of forming PBHs from supercooled first-order phase transitions, regions of false vacuum will originally be generated with $z > 0$ ($r > 0$) and with an inward velocity, i.e., moving toward $z = 0$. For a region to collapse, its radius must lie on the left side of the potential. If the initial radius lies instead on the right side of the potential, solutions to Eq. (43) indicate that the radius of the region would simply expand to infinity after reflecting off of the rightmost side of the potential. The region below the curve defined by Eq. (45) cannot be physically realized.
<figure><img src="image_3.png" /><figcaption>FIG. 1. The nondimensional potential from Eq. (45). The critical mass $M_{cr}$ is labeled, along with the lines that specify the Schwarzschild or Hubble radius $H_V^{-1}$ for a fixed energy $E$. The purple line and arrow illustrate the motion of a wall with initial size $r_{TP}$ which falls into its Schwarzschild radius to form a black hole. For illustrative purposes, we set $\gamma = 1.3$ in this figure, though $\gamma \ll 1$ in realistic models.</figcaption></figure>

Therefore, when searching for a viable set of parameters, only initial radii left of the potential $U(z)$ will be included.

We will define the radius $z_{TP}$, and correspondingly $r_{TP}$, as the turning point radius in which a system with fixed mass $M$, and therefore fixed “energy” $E$, intersects the leftmost wall of the potential. Based on our previous discussion, the regions that can collapse into PBHs must have initial radii between $z_s$ and $z_{TP}$. As an example, the purple line in Fig. 1 illustrates the path taken by a wall with $E = -6$, beginning $z_{TP}$, and falling inward toward $z_s$.

Finally, in order for a PBH to form in this scenario, a number conditions must hold. The first condition is related to the collapse time $\tau'$ or, alternatively, its dimensionful counterpart $\tau$. We require that $\tau$ be much smaller than the Hubble time, i.e.,

$$ \frac{\tau}{H_V^{-1}} = \frac{2\eta}{(1 + \eta^2)} \tau' \ll 1. \qquad (55) $$

To further avoid complications during collapse we will also require that bubble nucleation rate is slow compared to the rate of collapse. Given that $\Gamma$ is the rate per volume, the corresponding nucleation timescale within a patch of radius size $z$ at time $t$ is

$$ \tau_{\text{nuc}}(z; t) = \left[ \frac{4\pi}{3} r_s^3 \left(\frac{z}{z_s}\right)^3 \Gamma(t) \right]^{-1}. \qquad (56) $$

Therefore, we will also require that

$$ \frac{\tau}{\tau_{\text{nuc}}(z; t)} \ll 1, \qquad (57) $$

where the values $z$ and $t$ will be specified in Sec. VI.

The final PBH formation condition we will specify involves the physical size of the surrounding domain wall. In particular, we require that the size of the domain wall $\ell$ fits within the Schwarzschild radius of a black hole with mass specified by Eq. (35). This leads to the condition

$$ 1 \ge \frac{\ell}{2GM} \approx 6.7 \times 10^{-7} \left( \frac{10^{-13} M_{\odot}}{M} \right) \left( \frac{1 \text{ GeV}}{v} \right), \qquad (58) $$

which is quite easy to satisfy. To reiterate, $\ell$ is the thin-wall approximated value of the domain wall size. However, it is difficult to believe that any correction to the size of the domain wall will exceed the 6 orders of magnitude needed to overcome the above inequality.

So long as the above conditions are met, a PBH will form following the collapse of the false-vacuum region as described by Eq. (43).

## VI. THE POSSIBILITY OF AN ISLAND

Crucial to the supercooled phase transition paradigm is the existence of large islands of the false vacuum in which the phase transition is delayed. As is expected, the majority of bubble nucleation will occur when

$$ \Gamma(t_n) \approx H^4(t_n) \qquad (59) $$

where $t_n$, and the corresponding temperature $T_n$, satisfies the above equality. In our realization of PBH formation from supercooled phase transitions, the probability of collapse corresponds to the probability that a patch with comoving radius $\mathcal{R}$ can exist until some late time. Specifically,

$$ \mathcal{P}_{\text{coll}} \equiv \mathcal{P}_{\text{surv}}(t_n; \mathcal{R}), \qquad (60) $$

where $\mathcal{P}_{\text{surv}}(t_n; \mathcal{R})$ specifies the probability that a patch with comoving radius $\mathcal{R}$ has survived up to time $t_n$. As one expects, $t_n < t_{n_i}$.

The collapse probability can be written explicitly as [63,64]

$$ \ln \mathcal{P}_{\text{surv}}(t_{n_i}; \mathcal{R}) = - \int_{t_c}^{t_{n_i}} dt' \Gamma(t') a^3(t') \mathcal{V}(t'; t_{n_i}), \qquad (61) $$

where

$$ \mathcal{V}(t'; t_{n_i}) = \frac{4\pi}{3} [\mathcal{R}(t_{n_i}) + \mathcal{R}(t_{n_i}; t')]^3 \qquad (62) $$

is the volume corresponding to the spacelike slice at $t'$ of a patch with comoving radius $\mathcal{R}$ at time $t_{n_i}$. It is convenient to define the ratio

$$ \xi = \frac{\mathcal{R}(t_{n_i})}{[a(t_{n_i})H(t_{n_i})]^{-1}}, \qquad (63) $$
which allows us to directly compare the size of the comoving radius $\mathcal{R}$ and the comoving Hubble radius at time $t_{n_i}$. In particular, for $\xi < 1$ we are dealing with subhorizon patches, while for $\xi > 1$ the delayed patch is larger than the horizon size. For the parameters relevant for our discussion, we will always be well within the sub-horizon case.

As discussed in Sec. V, the initial radius of the false-vacuum region must lie within $z_s$ and $z_{TP}$ for a conventional black hole to form. This limits the acceptable values of $\xi$ in this circumstance. For simplicity, we will perform our calculations for regions size $z_{TP}$ at $t_{n_i}$. It is these values of $z$ and $t$ that will be used in the nucleation rate condition specified by Eq. (57).

Calculating the PBH abundance from supercooled phase transitions is nontrivial. In general, one has to consider the collapse of false-vacuum patches with varying sizes. This will lead to a nonmonochromatic distribution of PBHs. As is standard, we express the present-day PBH abundance as

$$f_{\text{PBH}} \equiv \frac{\rho_{\text{PBH},0}}{\rho_{\text{DM},0}}. \qquad (64)$$

In this particular context,

$$f_{\text{PBH}} = \frac{1}{\rho_{\text{DM},0}} \left(\frac{4\pi}{3} H_0^{-3}\right)^{-1} \times \int M N_{\text{patches}}(M) P_{\text{surv}}(M) g(M) dM. \quad (65)$$

Here, $N_{\text{patches}}$ denotes the number of patches within the present-day horizon that have the associated mass $M$. Equation (35) provides the relation between the mass $M$ and the initial size of the patch. The function $g(M)$ describes the distribution of patches. In the absence of an analytic expression for $g(M)$, we will utilize the monochromatic limit

$$g(M) = \delta(M - M_{\text{PBH}}), \qquad (66)$$

which will approximate the peak of the true distribution. This leads to

$$f_{\text{PBH}} = P_{\text{surv}} \frac{M_{\text{PBH}} N_{\text{patches}}}{\rho_{\text{DM},0} \frac{4\pi}{3} H_0^{-3}}. \qquad (67)$$

Here,

$$N_{\text{patches}} \approx \frac{1}{\xi^3} \cdot \left( \frac{a_{\text{eq}} H_{\text{eq}}}{a_0 H_0} \right)^3, \qquad (68)$$

where we used the fact that

$$a(t_{n_i}) H(t_{n_i}) \approx a(T_{\text{eq}}) H(T_{\text{eq}}) \qquad (69)$$

within the delayed, inflating patch. The expression for $N_{\text{patches}}$ can be further simplified by utilizing conservation of entropy, namely,

$$g_{*S}(T_{\text{eq}}) a_{\text{eq}}^3 T_{\text{eq}}^3 = g_{*S}(T_0) a_0^3 T_0^3. \qquad (70)$$

This allows us to express $N_{\text{patches}}$ as

$$N_{\text{patches}} \approx \frac{1}{\xi^3} \cdot \frac{g_{*S}(T_0)}{g_{*S}(T_{\text{eq}})} \left( \frac{T_0 H_{\text{eq}}}{T_{\text{eq}} H_0} \right)^3, \qquad (71)$$

where we used $g_{*S}(T_0) = 3.94$, $T_0 = 2.7$ K = $2.4 \times 10^{-13}$ GeV, and $H_0 = 1.4 \times 10^{-42}$ GeV.

This leads to a parametrization of $f_{\text{PBH}}$ of the form

$$f_{\text{PBH}} \approx \left( \frac{7.3 \times 10^{-5}}{\xi} \right)^3 \left( \frac{P_{\text{surv}}}{3 \times 10^{-8}} \right) \times \left( \frac{M_{\text{PBH}}}{10^{-13} M_{\odot}} \right) \left( \frac{g_*(T_{\text{eq}})}{100} \right)^{1/2} \left( \frac{T_{\text{eq}}}{150 \text{ MeV}} \right)^3. \quad (72)$$

## VII. RESULTS AND DISCUSSION

The above formalism enables us to identify numerous points in parameter space that could have important phenomenological implications. For the scope of this work, we will highlight a few points of interest and leave an exhaustive parameter search for later work. Table I presents three sets of parameters that offer PBH populations relevant for dark matter, microlensing observations, and gravitational wave detectors. A summary of our results can also be seen in Fig. 2 alongside the present-day landscape of observational constraints for PBHs.

Given that our scenario relies on the existence of a false-vacuum region surrounded primarily by true vacuum we require that the ratio $\rho_V^{\text{bkg}}/\rho_{\text{rad}}^{\text{bkg}} \sim \mathcal{O}(0.1)$ at $t = t_{n_i}$.

Figure 3 shows the nucleation rate associated with $M_{\text{PBH}} = 10^{-13} M_{\odot}$. A notable feature of the decay rate is that $(\Gamma/H^4)_{\text{max}} \lesssim 10^4$. In our parameter search, this feature has generally been required for the successful formation of PBHs. Adjusting the ratio $A/v$ can lead to extremely large values of $(\Gamma/H^4)_{\text{max}}$. After equaling 1, the ratio $\Gamma/H^4$ increases extremely rapidly. Since the nucleation rate is so large, the probability to survive any appreciable time after $t_n$ immediately becomes small [see Eq. (61)]. Given that the

<table><caption>TABLE I. Parameters that generate PBHs with masses relevant for dark matter or microlensing experiments.</caption><thead><tr><th>Model</th><th><i>v</i> (GeV)</th><th><i>A</i>/<i>v</i></th><th><i>M</i><sub>PBH</sub> [<i>M</i><sub>☉</sub>]</th><th><i>t</i><sub>n<sub>i</sub></sub>/<i>t</i><sub>eq</sub></th><th><i>f</i><sub>PBH</sub></th></tr></thead><tbody><tr><td>I</td><td>1</td><td>0.1052</td><td>10<sup>-13</sup></td><td>1.4106</td><td>1</td></tr><tr><td>II</td><td>10</td><td>0.1040</td><td>10<sup>-8</sup></td><td>1.3305</td><td>10<sup>-3</sup></td></tr><tr><td>III</td><td>50</td><td>0.1037</td><td>10<sup>-6</sup></td><td>1.7332</td><td>2 × 10<sup>-2</sup></td></tr></tbody></table>
<figure><img src="image_4.png" /><figcaption>FIG. 2. PBH masses generated from our realization of the supercooled phase transition scenario compared with observational constraints. A summary of the constraints can be found in Refs. [93,94]. Here, HE specifies Hawking evaporation constraints. The remaining constraints are due to microlensing observations by the Hyper Suprime-Cam (HSC) [95,96], Kepler (K) [97,98], EROS [99], OGLE [100], and Icarus (IC) [101].</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIG. 3. The nucleation rate per unit volume normalized according to the nondimensionalization specified in (B4). The maximum value of the rate is relatively small, which allows for the majority of the Universe to transition to radiation domination before the false-vacuum island collapses into a black hole. The parameters utilized correspond to model I in Table I.</figcaption></figure>

delay times needed for PBH formation occur significantly after $t_n$, we also note that the approximation

$$ \Gamma(t) \sim \Gamma_n e^{-\beta(t-t_n)}, \qquad (73) $$

where

$$ \beta \equiv HT \frac{dS_3}{dT} \bigg|_{T=T_n} \qquad (74) $$

does not hold in this scenario. This is demonstrated by Fig. 3, which shows that the time $t_{n_i}$ needed for PBH formation is relatively far from $t_n$. Therefore, a linear approximation, as expressed through Eq. (74), will fail to accurately approximate the rate $\Gamma(t)$ when $t = t_{n_i}$.

In summary, we have presented the possibility that supercooled phase transitions can lead to the formation of PBHs. Large islands of false vacuum have the possibility to persist after a short intermediate de Sitter phase. The domain waves separating the false- and true-vacuum regions are driven inward due to the pressure differential between the two phases. These domain waves fall within the Schwarzschild radius of the false-vacuum region, leading to the formation of a black hole. This framework is significantly different from existing literature, which instead relied on the critical density condition to signal the formation of PBHs [66–70,102,103]. We will also point to existing literature, such as Refs. [104–106], which also previously examined the formation of PBHs from supercooled phase transitions, but which have instead opted to use the Schwarzschild formation condition.

While not emphasized here, it is well known that first-order phase transitions can lead to sizable gravitational wave signals due to collision of bubble walls. The amplitude for such a signal is determined by the ratio [107–110]

$$ \Omega_{GW} \sim 10^{-6} \left(\frac{H}{\beta}\right)^2, \qquad (75) $$

while the peak frequency is given by

$$ f_{\text{peak}} \sim 1 \text{ mHz} \left(\frac{\beta}{H}\right) \left(\frac{T_{\text{RH}}}{100 \text{ TeV}}\right), \qquad (76) $$

where $T_{\text{RH}}$ is the temperature when the majority of the Universe is restored to a radiation dominated era.

For the parameters we have explored, $\beta/H \sim 100$ and $T_{\text{RH}} \lesssim 0.1$ GeV leading to $\Omega_{GW} \sim 10^{-10}$ with $f_{\text{peak}} \sim 10^{-7}$ Hz. Frequency is relevant for future pulsar timing array measurements. That being said, it would be difficult to associate any gravitational waves in this frequency range directly with the formation of PBHs. However, such a gravitational wave signal might help in distinguishing this PBH formation mechanism from the expected signals of others.

## ACKNOWLEDGMENTS

We thank K. Petraki for helpful discussions. This work was supported by the European Union's Horizon 2020 research and innovation program under Grant Agreement No. 101002846, ERC CoG CosmoChart. The work of A. K. was supported by the U.S. Department of Energy (DOE) Grant No. DE-SC0009937. The work of A. K. and M. S. was also supported by World Premier International Research Center Initiative (WPI), MEXT, Japan, and by
Japan Society for the Promotion of Science (JSPS) KAKENHI Grants No. JP20H05853 and JP24K00624.

## APPENDIX A: BOUNCE DETAILS

Calculating the single parameter bounce action, as defined by Eqs. (7)–(14), is a well-understood procedure. As discussed in the main text, a potential of the form of Eq. (10) can be redefined through field and coordinate definitions to obtain an equation of motion that depends only on one temperature-dependent parameter $κ$. In order to obtain the bounce solution, one must then solve Eq. (8) subject to the boundary conditions given in Eq. (9). Generally, this must be done numerically and many software packages have been written that can solve this simple bounce problem and its generalizations [82–85].

For simplicity, we will utilize a fitting function for the $κ$-dependent bounce action $S_3$ obtained in Ref. [86]. In particular, for $κ > 0$,

$$S_3(T) \approx \frac{m^3(T)}{\delta^2(T)} \frac{2\pi}{3(\kappa - \kappa_c)^2} \bar{B}_3(\kappa), \quad (A1)$$

where $\bar{B}_3(\kappa)$ is obtained from numerical fits and given by [86]

$$\bar{B}_3(\kappa) = \frac{16}{243} \left[ 1 - 38.23 \left(\kappa - \frac{2}{9}\right) + 115.26 \left(\kappa - \frac{2}{9}\right)^2 + 58.07 \sqrt{\kappa} \left(\kappa - \frac{2}{9}\right)^2 + 229.07 \kappa \left(\kappa - \frac{2}{9}\right)^2 \right]. \quad (A2)$$

We verified that this fitting function is in agreement with the results obtained in Ref. [87], as well as with our own numerical implementation of the one-parameter bounce solution.

## APPENDIX B: NUMERICAL APPROACH

In this section, we hope to illustrate one possible methodology for solving (28)–(34). Our first step is rewriting the background evolution equations (28) and (29) in terms of the dimensionless time variable $τ$ [see Eq. (30)]. This simple exercise yields

$$a'(\tau) = a(\tau)[\hat{\rho}_R(\tau) + \hat{\rho}_V(\tau)]^{1/2}, \quad (B1)$$

$$\hat{\rho}'_R(\tau) + 4 \frac{a'(\tau)}{a(\tau)} \hat{\rho}_R = -\hat{\rho}'_V(\tau), \quad (B2)$$

where $\tau \equiv d/d\tau$. We will normalize the scale factor such that $a(\tau_{eq}) = 1$. In addition, we must also nondimensionize Eqs. (32) and (33). It is a simple exercise to show that this leads to

$$\hat{R}(\tau; \tau') = \int_{\tau'}^{\tau} d\tilde{\tau} \frac{v_w(\tilde{\tau})}{a(\tilde{\tau})}, \quad (B3)$$

$$I(\tau) = \frac{4\pi}{3} \int_{\tau_c}^{\tau} \hat{\Gamma}(\tau') a^3(\tau') \hat{R}^3(\tau; \tau') d\tau'. \quad (B4)$$

The nondimensional radius $\hat{R}$ and rate $\hat{\Gamma}$ are related to their dimensionful quantities by

$$R(\tau; \tau') = \frac{\sqrt{2}}{H_{\text{eq}}} \hat{R}(\tau; \tau'), \quad (B5)$$

$$\Gamma(\tau) = \frac{H_{\text{eq}}^4}{4} \hat{\Gamma}(\tau). \quad (B6)$$

The next step in our method will be to define

$$r(\tau) \equiv \int_{\tau_c}^{\tau} d\tilde{\tau} \frac{v_w(\tilde{\tau})}{a(\tilde{\tau})} \quad (B7)$$

such that

$$\hat{R}(\tau; \tau') = r(\tau) - r(\tau') \quad (B8)$$

and

$$r'(\tau) = \frac{v_w(\tau)}{a(\tau)}, \quad r(\tau \le \tau_c) = 0. \quad (B9)$$

Furthermore, we can define

$$v_i(\tau) = \int_{\tau_c}^{\tau} \hat{\Gamma}(\tau') a^3(\tau') r^i(\tau') d\tau' \quad (B10)$$

such that

$$v_i'(\tau) = \hat{\Gamma}(\tau) a^3(\tau) r^i(\tau), \quad v_i(\tau \le \tau_c) = 0. \quad (B11)$$

This enables us to write $I(\tau)$ as

$$I(\tau) = \frac{4\pi}{3} [r^3(\tau)v_0(\tau) - 3r^2(\tau)v_1(\tau) + 3r(\tau)v_2(\tau) - v_3(\tau)]. \quad (B12)$$

Thus, we have extended our two-dimensional system of integrodifferential equation to a seven-dimensional system of first-order ordinary differential equations (ODEs),

$$\{a(\tau), \rho_R(\tau)\} \rightarrow \{a(\tau), \rho_R(\tau), r(\tau), v_i(\tau)\}. \quad (B13)$$

The initial conditions are

$$\{a(\tau), \rho_R(\tau), r(\tau), v_i(\tau)\}_{\tau=\tau_c} = \{1, 1, 0, 0\}. \quad (B14)$$

The cost of solving seven-dimensional first-order ODEs is minimal, and so this method provides quick, easy to implement strategy to examine background evolution of a supercooled phase transition.
[1] Y. B. Zel'dovich and I. D. Novikov, Sov. Astron. **10**, 602 (1967).

[2] S. Hawking, Mon. Not. R. Astron. Soc. **152**, 75 (1971).

[3] B. J. Carr and S. W. Hawking, Mon. Not. R. Astron. Soc. **168**, 399 (1974).

[4] M. Khlopov, B. Malomed, and I. Zeldovich, Mon. Not. R. Astron. Soc. **215**, 575 (1985).

[5] A. Dolgov and J. Silk, Phys. Rev. D **47**, 4244 (1993).

[6] J. Yokoyama, Astron. Astrophys. **318**, 673 (1997).

[7] E. L. Wright, Astrophys. J. **459**, 487 (1996).

[8] J. Garcia-Bellido, A. D. Linde, and D. Wands, Phys. Rev. D **54**, 6040 (1996).

[9] M. Kawasaki, N. Sugiyama, and T. Yanagida, Phys. Rev. D **57**, 6050 (1998).

[10] A. M. Green, A. R. Liddle, K. A. Malik, and M. Sasaki, Phys. Rev. D **70**, 041502 (2004).

[11] M. Y. Khlopov, Res. Astron. Astrophys. **10**, 495 (2010).

[12] B. J. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama, Phys. Rev. D **81**, 104019 (2010).

[13] P. H. Frampton, M. Kawasaki, F. Takahashi, and T. T. Yanagida, J. Cosmol. Astropart. Phys. 04 (2010) 023.

[14] S. Clesse and J. Garcia-Bellido, Phys. Rev. D **92**, 023524 (2015).

[15] M. Kawasaki, A. Kusenko, Y. Tada, and T. T. Yanagida, Phys. Rev. D **94**, 083523 (2016).

[16] B. Carr, F. Kuhnel, and M. Sandstad, Phys. Rev. D **94**, 083504 (2016).

[17] K. Inomata, M. Kawasaki, K. Mukaida, Y. Tada, and T. T. Yanagida, Phys. Rev. D **95**, 123510 (2017).

[18] S. Pi, Y.-l. Zhang, Q.-G. Huang, and M. Sasaki, J. Cosmol. Astropart. Phys. 05 (2018) 042.

[19] K. Inomata, M. Kawasaki, K. Mukaida, Y. Tada, and T. T. Yanagida, Phys. Rev. D **96**, 043504 (2017).

[20] J. Garcia-Bellido, M. Peloso, and C. Unal, J. Cosmol. Astropart. Phys. 09 (2017) 013.

[21] J. Georg and S. Watson, J. High Energy Phys. 09 (2017) 138.

[22] K. Inomata, M. Kawasaki, K. Mukaida, and T. T. Yanagida, Phys. Rev. D **97**, 043514 (2018).

[23] B. Kocsis, T. Suyama, T. Tanaka, and S. Yokoyama, Astrophys. J. **854**, 41 (2018).

[24] K. Ando, K. Inomata, M. Kawasaki, K. Mukaida, and T. T. Yanagida, Phys. Rev. D **97**, 123512 (2018).

[25] E. Cotner and A. Kusenko, Phys. Rev. Lett. **119**, 031103 (2017).

[26] E. Cotner and A. Kusenko, Phys. Rev. D **96**, 103002 (2017).

[27] E. Cotner, A. Kusenko, and V. Takhistov, Phys. Rev. D **98**, 083513 (2018).

[28] M. Sasaki, T. Suyama, T. Tanaka, and S. Yokoyama, Classical Quantum Gravity **35**, 063001 (2018).

[29] B. Carr and J. Silk, Mon. Not. R. Astron. Soc. **478**, 3756 (2018).

[30] C. Germani and I. Musco, Phys. Rev. Lett. **122**, 141302 (2019).

[31] U. Banik, F. C. van den Bosch, M. Tremmel, A. More, G. Despali, S. More, S. Vegetti, and J. P. McKean, Mon. Not. R. Astron. Soc. **483**, 1558 (2019).

[32] A. Escrivà, C. Germani, and R. K. Sheth, Phys. Rev. D **101**, 044022 (2020).

[33] C. Germani and R. K. Sheth, Phys. Rev. D **101**, 063520 (2020).

[34] F. Hoyle and R. A. Lyttleton, Proc. Cambridge Philos. Soc. **35**, 405 (1939).

[35] E. Cotner, A. Kusenko, M. Sasaki, and V. Takhistov, J. Cosmol. Astropart. Phys. 10 (2019) 077.

[36] A. Kusenko, M. Sasaki, S. Sugiyama, M. Takada, V. Takhistov, and E. Vitagliano, Phys. Rev. Lett. **125**, 181304 (2020).

[37] M. M. Flores and A. Kusenko, Phys. Rev. Lett. **126**, 041101 (2021).

[38] J. de Freitas Pacheco and J. Silk, Phys. Rev. D **101**, 083022 (2020).

[39] V. Takhistov, G. M. Fuller, and A. Kusenko, Phys. Rev. Lett. **126**, 071101 (2021).

[40] V. De Luca, G. Franciolini, and A. Riotto, Phys. Rev. Lett. **126**, 041303 (2021).

[41] G. Dvali, F. Kühnel, and M. Zantedeschi, Phys. Rev. D **104**, 123507 (2021).

[42] M. Biagetti, V. De Luca, G. Franciolini, A. Kehagias, and A. Riotto, Phys. Lett. B **820**, 136602 (2021).

[43] Y.-F. Cai, X.-H. Ma, M. Sasaki, D.-G. Wang, and Z. Zhou, Phys. Lett. B **834**, 137461 (2022).

[44] S. Pi and M. Sasaki, Phys. Rev. D **108**, L101301 (2023).

[45] B. Carr and F. Kuhnel, SciPost Phys. Lect. Notes **48**, 1 (2022).

[46] P. Lu, K. Kawana, and A. Kusenko, Phys. Rev. D **107**, 103037 (2023).

[47] T. Harada, K. Kohri, M. Sasaki, T. Terada, and C.-M. Yoo, J. Cosmol. Astropart. Phys. 02 (2023) 038.

[48] G. B. Gelmini, J. Hyman, A. Simpson, and E. Vitagliano, J. Cosmol. Astropart. Phys. 06 (2023) 055.

[49] M. Chen, G. B. Gelmini, P. Lu, and V. Takhistov, arXiv: 2312.12136.

[50] M. Chen, G. B. Gelmini, P. Lu, and V. Takhistov, Phys. Lett. B **852**, 138609 (2024).

[51] B. Carr, S. Clesse, J. Garcia-Bellido, M. Hawkins, and F. Kuhnel, Phys. Rep. **1054**, 1 (2024).

[52] G. Domènech and M. Sasaki, arXiv:2401.07615.

[53] R. Bean and J. Magueijo, Phys. Rev. D **66**, 063505 (2002).

[54] M. Kawasaki, A. Kusenko, and T. T. Yanagida, Phys. Lett. B **711**, 1 (2012).

[55] M. M. Flores, Y. Lu, and A. Kusenko, Phys. Rev. D **108**, 123511 (2023).

[56] Y. Lu, Z. S. C. Picker, and A. Kusenko, arXiv:2312.15062.

[57] G. M. Fuller, A. Kusenko, and V. Takhistov, Phys. Rev. Lett. **119**, 061101 (2017).

[58] V. Takhistov, Phys. Lett. B **789**, 538 (2019).

[59] V. Takhistov, Phys. Lett. B **782**, 77 (2018).

[60] M. A. Abramowicz, M. Bejger, and M. Wielgus, Astrophys. J. **868**, 17 (2018).

[61] K. Kainulainen, S. Nurmi, E. D. Schiappacasse, and T. T. Yanagida, Phys. Rev. D **104**, 123033 (2021).

[62] M. M. Flores, A. Kusenko, A. M. Ghez, and S. Naoz, Phys. Rev. D **108**, L061301 (2023).

[63] K. Sato, H. Kodama, M. Sasaki, and K.-i. Maeda, Phys. Lett. **108B**, 103 (1982).
[64] H. Kodama, M. Sasaki, and K. Sato, Prog. Theor. Phys. **68**, 1979 (1982).

[65] S. K. Blau, E. I. Guendelman, and A. H. Guth, Phys. Rev. D **35**, 1747 (1987).

[66] J. Liu, L. Bian, R.-G. Cai, Z.-K. Guo, and S.-J. Wang, Phys. Rev. D **105**, L021303 (2022).

[67] M. J. Baker, M. Breitbach, J. Kopp, and L. Mittnacht, arXiv:2110.00005.

[68] K. Kawana, T. Kim, and P. Lu, Phys. Rev. D **108**, 103531 (2023).

[69] Y. Gouttenoire and T. Volansky, arXiv:2305.04942.

[70] I. Baldes and M. O. Olea-Romacho, J. High Energy Phys. 01 (2024) 133.

[71] D. H. Lyth, K. A. Malik, and M. Sasaki, J. Cosmol. Astropart. Phys. 05 (2005) 004.

[72] M. Shibata and M. Sasaki, Phys. Rev. D **60**, 084002 (1999).

[73] I. Musco, J. C. Miller, and L. Rezzolla, Classical Quantum Gravity **22**, 1405 (2005).

[74] I. Musco, J. C. Miller, and A. G. Polnarev, Classical Quantum Gravity **26**, 235001 (2009).

[75] T. Harada, C.-M. Yoo, and K. Kohri, Phys. Rev. D **88**, 084051 (2013); **89**, 029903(E) (2014).

[76] T. Harada, C.-M. Yoo, and Y. Koga, Phys. Rev. D **108**, 043515 (2023).

[77] K. S. Thorne, in Magic Without Magic: John Archibald Wheeler, edited by J. R. Klauder (Freeman, San Francisco, 1972), p. 231.

[78] S. R. Coleman, Phys. Rev. D **15**, 2929 (1977); **16**, 1248(E) (1977).

[79] C. G. Callan, Jr. and S. R. Coleman, Phys. Rev. D **16**, 1762 (1977).

[80] A. D. Linde, Phys. Lett. **70B**, 306 (1977).

[81] A. D. Linde, Nucl. Phys. **B216**, 421 (1983); **B223**, 544(E) (1983).

[82] A. Masoumi, K. D. Olum, and B. Shlaer, J. Cosmol. Astropart. Phys. 01 (2017) 051.

[83] R. Sato, Comput. Phys. Commun. **258**, 107566 (2021).

[84] P. Athron, C. Balázs, M. Bardsley, A. Fowlie, D. Harries, and G. White, Comput. Phys. Commun. **244**, 448 (2019).

[85] V. Guada, M. Nemevšek, and M. Pintar, Comput. Phys. Commun. **256**, 107480 (2020).

[86] N. Levi, T. Opferkuch, and D. Redigolo, J. High Energy Phys. 02 (2023) 125.

[87] A. Megevand and S. Ramirez, Nucl. Phys. **B919**, 74 (2017).

[88] M. S. Turner, E. J. Weinberg, and L. M. Widrow, Phys. Rev. D **46**, 2384 (1992).

[89] J. Garriga, A. Vilenkin, and J. Zhang, J. Cosmol. Astropart. Phys. 02 (2016) 064.

[90] H. Deng and A. Vilenkin, J. Cosmol. Astropart. Phys. 12 (2017) 044.

[91] W. Israel, Nuovo Cimento B **44S10**, 1 (1966); **48**, 463(E) (1967).

[92] W. Israel, Phys. Rev. **153**, 1388 (1967).

[93] B. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama, Rep. Prog. Phys. **84**, 116902 (2021).

[94] B. Carr and F. Kuhnel, Annu. Rev. Nucl. Part. Sci. **70**, 355 (2020).

[95] H. Niikura et al., Nat. Astron. **3**, 524 (2019).

[96] N. Smyth, S. Profumo, S. English, T. Jeltema, K. McKinnon, and P. Guhathakurta, Phys. Rev. D **101**, 063005 (2020).

[97] K. Griest, A. M. Cieplak, and M. J. Lehner, Phys. Rev. Lett. **111**, 181302 (2013).

[98] K. Griest, A. M. Cieplak, and M. J. Lehner, Astrophys. J. **786**, 158 (2014).

[99] P. Tisserand et al. (EROS-2 Collaboration), Astron. Astrophys. **469**, 387 (2007).

[100] H. Niikura, M. Takada, S. Yokoyama, T. Sumi, and S. Masaki, Phys. Rev. D **99**, 083503 (2019).

[101] M. Oguri, J. M. Diego, N. Kaiser, P. L. Kelly, and T. Broadhurst, Phys. Rev. D **97**, 023518 (2018).

[102] Y. Gouttenoire, arXiv:2311.13640.

[103] A. Conaci, L. Delle Rose, P. S. B. Dev, and A. Ghoshal, arXiv:2401.09411.

[104] R. Jinno, J. Kume, and M. Yamada, Phys. Lett. B **849**, 138465 (2024).

[105] M. Lewicki, P. Toczek, and V. Vaskonen, J. High Energy Phys. 09 (2023) 092.

[106] M. Lewicki, P. Toczek, and V. Vaskonen, arXiv:2402.04158.

[107] T. Konstandin, J. Cosmol. Astropart. Phys. 03 (2018) 047.

[108] D. Cutting, M. Hindmarsh, and D. J. Weir, Phys. Rev. D **97**, 123513 (2018).

[109] D. Cutting, E. G. Escartin, M. Hindmarsh, and D. J. Weir, Phys. Rev. D **103**, 023531 (2021).

[110] M. Lewicki and V. Vaskonen, Eur. Phys. J. C **80**, 1003 (2020).