Sublunar-Mass Primordial Black Holes from Closed Axion Domain Walls

TL;DR

This paper proposes a model where sublunar-mass primordial black holes originate from closed axion domain walls near the QCD scale, linking axion parameters to PBH abundance and dark matter contributions.

Contribution

It introduces a novel formation mechanism for PBHs from closed axion domain walls and constrains axion parameters based on PBH observational limits.

Findings

PBHs in the mass range $10^{20}$-$10^{22}$ g can account for up to 1% of dark matter.

The model favors axion mass around the meV scale ($f_{a} extasciitilde 10^{9}$ GeV).

PBH abundance is sensitive to the formation efficiency of closed axion domain walls.

Abstract

We study the formation of primordial black holes (PBHs) from the collapse of closed domain walls (DWs) which naturally arise in QCD axion models near the QCD scale together with the main string-wall network. The size distribution of the closed DWs is determined by percolation theory, from which we further obtain PBH mass distribution and abundance. Various observational constraints on PBH abundance in turn also constrain axion parameters. Our model prefers axion mass around the meV scale ($f_{a}\sim 10^{9}$ GeV). The corresponding PBHs are in the sublunar-mass window $10^{20}$-$10^{22}$ g (i.e., $10^{-13}$-$10^{-11}M_{\odot}$), one of few mass windows still available for PBHs contributing significantly to dark matter (DM). In our model, PBH abundance could reach $\sim1\%$ or even more of DM, sensitive to the formation efficiency of closed axion DWs.

Equations

$$V_a=m_a^2(T)f_a^2[1-\cos (\varphi /f_a)]$$

$$m_a(T)=\{\begin{array}{ll}{f}_a^{-1}{\chi }_0^{1/2},& T\le T_c\\ f_a^{-1}{\chi }_0^{1/2}(T/T_c{)}^{-\beta },& T\ge T_c\end{array}$$

$$m_a(t_1)t_1\simeq 1.$$

$$\xi (T)\simeq m_a^{-1}(T)$$

$$n_s\propto s^{-\tau }\exp (-\lambda s^{2/3}).$$

$$f(r_1)=f_0\cdot r_1^{2-3\tau }\cdot {\mathrm{e}}^{\lambda (1-r_1^2)}$$

$${\int }_1^{\infty }d{r}_{1}4\pi (\xi r_1{)}^{2}f(r_1)\simeq \gamma \cdot \frac{1}{\xi }.$$

$$t_2\simeq r_1^2{t}_1.$$

$$\left[{\partial }_t^2+\frac{3{\partial }_t}{2t}-\frac{{\partial }_{\mathcal{R}}^2}{a^2(t)}-\frac{2{\partial }_{\mathcal{R}}}{a^2(t)\mathcal{R}}\right]\tilde{\varphi }+m_a^2(t)\sin \tilde{\varphi }=0$$

$$\tilde{\varphi }(t=t_2,\mathcal{R})=$$

$$+{\tan }^{-1}[{\mathrm{e}}^{m_a(t_2)(-\mathcal{R}-R_2)}]\}$$

$$\frac{R_s}{R}=\frac{2GE(t,R)}{R}\gtrsim 1\Rightarrow S(t,R)\gtrsim m_{\mathrm{P}}^2$$

$$S_{\max }=k_1[m_a(t_2)R_2{]}^{k_2}\cdot f_a^2$$

$$r_1^2\gtrsim \frac{m_a(t_1)}{m_a(t_2)}{\left(\frac{m_P^2}{k_1{f}_a^2}\right)}^{1/k_2}.$$

$$M_{\mathrm{PBH}}\simeq 4\pi R_2^2\sigma (t_2)\simeq 4\pi r_1^4\cdot m_a^{-2}(t_1)\cdot \sigma (r_1^2{t}_1)$$

$$\frac{d{\rho }_{\mathrm{PBH}}(t)}{d{M}_{\mathrm{PBH}}}=M_{\mathrm{PBH}}(r_1)\cdot f(r_1)\cdot {\left[\frac{T(t)}{T_1}\right]}^3\cdot \frac{d{r}_1}{d{M}_{\mathrm{PBH}}}$$

$$\frac{d{\Omega }_{\mathrm{PBH}}(t_{\mathrm{eq}})}{d{M}_{\mathrm{PBH}}}=\frac{M_{\mathrm{PBH}}(r_1)\cdot f(r_1)}{{\rho }_{\mathrm{cr}}(t_1)}\cdot \frac{T_1}{T_{\mathrm{eq}}}\cdot \frac{d{r}_1}{d{M}_{\mathrm{PBH}}}$$

$${\Omega }_{\mathrm{PBH}}={\int }_{r_{1,\min }}^{\infty }\left(\frac{M_{\mathrm{PBH}}(r_1)\cdot f(r_1)}{{\rho }_{\mathrm{cr}}(t_1)}\cdot \frac{T_1}{T_{\mathrm{eq}}}\right)d{r}_1.$$

$$⟨M_{\mathrm{PBH}}⟩=\frac{{\int }_{r_{1,\min }}^{\infty }d{r}_1{M}_{\mathrm{PBH}}(r_1)f(r_1)}{{\int }_{r_{1,\min }}^{\infty }d{r}_{1}f(r_1)},$$

$$T_1\simeq 1\mathrm{GeV}\cdot {\left(\frac{10^{12}\mathrm{GeV}}{f_a}\right)}^{1/6}.$$

$$t_{\mathrm{w}}\simeq \mu (t_{\mathrm{w}})/\sigma (t_{\mathrm{w}}),$$

$$T_{\mathrm{w}}\simeq 600\mathrm{MeV}\cdot {\left(\frac{10^{12}\mathrm{GeV}}{f_a}\right)}^{1/6}$$

$${\int }_{t_{\mathrm{w}}}^{t_{\mathrm{hit}}}\frac{vdt}{a(t)}\simeq \frac{t_{\mathrm{w}}}{a(t_{\mathrm{w}})}$$

$$t_{\mathrm{hit}}\simeq 5.1t_{\mathrm{w}},T_{\mathrm{hit}}\simeq 0.44T_{\mathrm{w}}\simeq 0.26T_1,$$

$$T_{\mathrm{hit}}\lesssim T_{\mathrm{decay}}$$

$$T_{\mathrm{hit}}\ll T_{\mathrm{decay}}.$$

$$\frac{{\partial }^2\tilde{\varphi }}{\partial {\tilde{t}}^2}+\frac{3}{2\tilde{t}}\frac{\partial \tilde{\varphi }}{\partial \tilde{t}}-\frac{1}{a^2(\tilde{t})}\left(\frac{{\partial }^2\tilde{\varphi }}{\partial {\tilde{r}}^2}+\frac{2}{\tilde{r}}\frac{\partial \tilde{\varphi }}{\partial \tilde{r}}\right)+\frac{m_a^2(\tilde{t})}{m_a^2({\tilde{t}}_2)}\sin \tilde{\varphi }=0,$$

$$\tilde{\varphi }(\widetilde{t_2},\tilde{r})=4\{{\tan }^{-1}[{\mathrm{e}}^{(\tilde{r}-{\tilde{r}}_2)}]+{\tan }^{-1}[{\mathrm{e}}^{(-\tilde{r}-{\tilde{r}}_2)}]\},$$

$${\frac{\partial \tilde{\varphi }(\tilde{t},\tilde{r})}{\partial \tilde{t}}|}_{\tilde{t}=\widetilde{t_2}}=0$$

$$a(\tilde{t})={\left(\frac{t}{t_2}\right)}^{1/2}={\left(\frac{\tilde{t}}{{\tilde{t}}_2}\right)}^{1/2}.$$

$$\frac{m_a(\tilde{t})}{m_a({\tilde{t}}_2)}={\left(\frac{t}{t_2}\right)}^{\beta /2}={\left(\frac{\tilde{t}}{{\tilde{t}}_2}\right)}^{\beta /2}.$$

Full text

††journal: Physics of the Dark Universe

1 Introduction

Primordial black holes (PBHs) have long been considered as viable dark matter (DM) candidates, see Refs. [1, 2, 3] for recent reviews. Despite various observational constraints, some mass windows remain valid in which PBHs could significantly contribute to DM: sublunar-mass range $\mathcal{O}(10^{20}{\rm g})$ and intermediate mass range $\mathcal{O}(10M_{\odot})$ [1, 2, 4]. In addition to the frequently studied mechanism of PBH formation from the collapse of overdense regions in the early universe [1, 2], PBHs could also be formed from the collapse of topological defects [5, 6, 7, 8, 9, 10, 11, 12, 13, 14].

QCD axion was originally proposed as a solution to strong CP problem [15, 16, 17, 18, 19, 20, 21]. As Peccei-Quinn (PQ) symmetry gets spontaneously broken at PQ scale $T_{\rm PQ}\sim f_{a}$ in the early universe, axion strings are formed. If PQ symmetry is broken after inflation ($f_{a}\lesssim H_{I}$, post-inflationary scenario), axion domain walls (DWs) will be formed later near QCD scale $T_{1}\sim$ GeV with the pre-existing strings as boundaries, which we call the string-wall network [22, 23]. Otherwise, in the pre-inflationary scenario, the pre-existing strings are ‘blown away’ and the axion field gets homogenized by inflation, so no DWs can be formed at $T_{1}$. Propagating axions generated from misalignment mechanism and topological decays are also DM candidates [24, 25].

Recently, Refs. [26, 27] have studied PBH formation from the collapse of closed axion DWs. The PBH mass obtained in Ref. [26] is $\sim 10^{-8}M_{\odot}$ ($10^{25}$ g), but much heavier in Ref. [27] $\sim 10^{4}$-$10^{7}M_{\odot}$ since an extra bias term is considered there lifting the energy enclosed by DWs. Closed DWs in Refs. [26, 27] are related to the network fragment which could occur much later than $T_{1}$, and PBH formation there is significantly affected by the fragment time which is however very hard to determine [28, 29, 30, 31, 32].

In this paper, however, we study the closed axion DWs initially formed at $T_{1}$ together with the main string-wall network. The closed DWs thus evolve independently of the network fragment. Also, we focus on $N_{\rm DW}=1$ case. The size distribution of $N_{\rm DW}=1$ closed DWs initially formed at $T_{1}$ is well predicted by percolation theory, from which we can further calculate the PBH mass distribution and abundance. Another advantage is that $N_{\rm DW}=1$ model naturally avoids the known DW problem that arises in $N_{\rm DW}>1$ models leading to a DW-dominated universe [24, 33]. The DW problem in $N_{\rm DW}>1$ cases can also be avoided with a bias term introduced, which is adopted in Ref. [27], although there is only little room in parameter space for this term [24].

In our model, for axion decay constant $f_{a}\sim 10^{9}$ GeV, PBHs formed from the collapse of closed axion DWs are in the sublunar-mass window $\sim 10^{20}$-$10^{22}$ g, one of few allowed windows constrained by observations. In addition to the propagating axions generated from misalignment mechanism and topological decays as conventional DM candidates, PBH abundance in our model could reach $\sim 1\%$ of DM, sensitive to the formation efficiency of closed DWs at $T_{1}$. Additionally, various observational constraints on PBH abundance in turn could constrain QCD axion parameter space.

The paper is organized as follows. In Section 2, we briefly review the formation of axion DWs and discuss the size distribution of $N_{\rm DW}=1$ closed axion DWs predicted by percolation theory. In Section 3, we study the criterion for a closed DW to collapse into a black hole. In Section 4, we present the PBH mass distribution and abundance obtained in our model, in comparison with the constraints from astrophysical observations. Also, the constraints on PBH abundance in turn are used to constrain QCD axion parameter space. We draw the conclusions in Section 5.

2 Size distribution of closed axion DWs

We start with a brief review of axion DWs formation. Non-perturbative QCD effects induce an effective potential for the axion field $\phi$ [24, 25]:

[TABLE]

with $0\leq\phi/f_{a}\leq 2\pi N_{\rm DW}$ where $N_{\rm DW}$ is the model-dependent chiral anomaly coefficient [34] that also represents the number of degenerate vacua locating at $\phi/f_{a}=2k\pi$. The axion mass is [35, 36]

[TABLE]

where $T_{c}\simeq 150$ MeV is the QCD transition temperature, $\chi_{0}=(75.6{\rm~{}MeV})^{4}$ is the zero-temperature topological susceptibility and $\beta\simeq 4$ [35, 37].

$V_{a}$ is unimportant until $m_{a}(T)$ increases to the scale of the inverse of Hubble radius $H\sim t^{-1}$ at $t_{1}$ [24]

[TABLE]

We say axion mass effectively turns on at $t_{1}$. The corresponding temperature is $T_{1}\sim 1$ GeV, much lower than PQ scale. In the post-inflationary scenario, axion DWs start to form due to Kibble-Zurek mechanism [38, 39] at $T_{1}$ when different regions of the universe fall into different vacua. The typical length of each region is the correlation length $\xi$ (see e.g. Refs. [40, 41]):

[TABLE]

Using Eq. (3), we further get $\xi(T_{1})\simeq t_{1}$, i.e. the correlation length at DW formation point $t_{1}$ is approximately the Hubble radius.

If $N_{\rm DW}=1$, the topology of vacuum manifold has two discrete values, $\phi/f_{a}=0,2\pi$, corresponding to the same physical vacuum. It is known that DWs can be formed in this case as $\phi$ interpolates between the two topological branches [math] and $2\pi$ [24, 42], and they could live long enough against tunnelling process to have important implications [42, 43]. If we ignore the pre-existing strings at $T_{1}$ (the effects of which will be discussed later), $N_{\rm DW}=1$ model can be treated as $Z_{2}$ model, for they have identical topology of vacuum manifold: both have two discrete values [44]. The formation of such walls in the early universe has been widely studied in the literature (see e.g. Refs. [45, 46]): different ‘cells’ (typical length $\xi$) fall into one of the two values randomly with equal probability. Two or more neighbouring cells falling into the same value form a finite cluster (closed DW). A mathematical theory known as percolation theory studies the size distribution of such clusters, which gives [45]:

[TABLE]

$n_{s}$ is the number density of finite clusters with size $s$ (number of cells within a cluster). $\tau=-1/9$ and $\lambda\approx 0.025$ are two coefficients from percolation theory 111$\lambda$ is obtained indirectly. In percolation theory, $\lambda^{-1}$ is the crossover size where $\lambda^{-1}\simeq\left|p-p_{c}\right|^{\rm-1/\sigma}$ valid for $\left|p-p_{c}\right|\ll 1$ (see e.g. Refs. [47, 48, 49]). $p$ is the probability of each cell choosing one of the two topological branches, so $p=0.5$ in our case; $p_{c}=0.31$ for cubic lattice and $\sigma=0.45$ in 3D [50], so $\lambda\approx 0.025$ for $\left|p-p_{c}\right|\ll 1$ well satisfied. The other coefficient $\tau=-1/9$ for $p>p_{c}$ is obtained in a field theoretical formulation of the percolation problem [51, 50].. Although Eq. (5) is originally obtained with the assumption $s\gg 1$, it can be extrapolated down to the smallest clusters $s=1$ with high accuracy [52].

Eq. (5) can be translated into DW language straightforwardly. Finite clusters are closed DWs with volume $R_{1}^{3}\simeq s\xi^{3}$, where $R_{1}$ is introduced as the radius of closed DWs. We can write $n_{s}$ in differential form as $n_{s}=d\mathfrak{n}/ds$ where $\mathfrak{n}$ denotes the number density of finite clusters with size smaller than $s$. Then, Eq. (5) becomes

[TABLE]

where $r_{1}\equiv R_{1}/\xi$, $f(r_{1})\equiv d\mathfrak{n}/dr$. $f_{0}\equiv f(r_{1}=1)$ is the distribution at the smallest size $R_{1}=\xi$.

Closed DWs are indeed observed in computer simulations. In $Z_{2}$-system, closed DWs account for $\gamma\sim 13\%$ of total wall area [45]. We expect the proportion is lower in $N_{\rm DW}=1$ models with strings present, because the presence of strings makes less space available to form closed DWs. This has also been seen in simulations [45, 53]. But it’s hard to determine the strings effects exactly. One difficulty is that simulations are sensitive to simulation size [45] and may not be properly applied to the universe at $T_{1}$. Another difficulty is that simulations only apply to DWs formed soon after strings formation [45] which contradicts the realistic case $T_{1}\ll T_{\rm PQ}$. Despite simulation difficulties, we can absorb the strings effects on closed DWs at $T_{1}$ into $\gamma$ (defined as the proportion of closed DWs area in total wall area [40]), implying $\gamma\lesssim 13\%$ with strings present. Additionally, in contrast with the traditional view, $N_{\rm DW}=1$ DWs could also be formed in the pre-inflationary scenario ($f_{a}\gtrsim H_{I}$) based on the argument that different topological branches cannot be separated by inflation [54, 40] 222$N_{\rm DW}=1$ closed axion DWs formed in the pre-inflationary scenario are crucial in Refs. [54, 40]. The closed walls there accumulate baryons or anti-baryons inside. They finally evolve into the axion quark nuggets (AQNs) which have many intriguing astrophysical and cosmological implications. See the original paper [54] and recent developments [40, 55, 56, 57, 58, 59, 60, 61, 62, 63, 44, 64] for details.. In that scenario, the pre-existing strings are blown away by inflation, so they cannot affect the formation of closed DWs at $T_{1}$, implying that $\gamma\sim 13\%$, the same as $Z_{2}$ case.

We can also interpret the correlation length $\xi$ as the average distance among DWs, to get

[TABLE]

The best information we have about $\gamma$ in the post-inflationary scenario is $\gamma\lesssim 13\%$ (but nonzero, since closed DWs are observed with strings present [45, 53]). One might worry that closed DWs could be destroyed by intercommuting with walls bounded by strings in the late time evolution after $T_{1}$, but our analysis shows that closed DWs will survive, see A for details.

3 Collapse into PBHs

Closed DWs with size $r_{1}>1$ (i.e. $R_{1}>\xi(T_{1})$) are super-Hubble structures since $\xi(T_{1})\simeq t_{1}$. They do not collapse until the size is surpassed by Hubble horizon. We emphasize that super-Hubble DWs are formed not because $\phi$ is physically correlated in super-Hubble scale, but a natural result of random combinations of self-correlated cells predicted by percolation theory.

Instead of contraction, super-Hubble closed DWs first expand due to the universe’s expansion with the scale factor $a(t)\propto T^{-1}\propto t^{1/2}$ (radiation-dominated era). However, the Hubble horizon $H^{-1}\sim t$ increases faster, implying that some time after $t_{1}$ (labeled as $t_{2}$), $H^{-1}$ will catch up with the closed DWs size, $R_{2}\simeq t_{2}$. $R_{1}$ and $R_{2}$ are connected by the universe’s expansion, $R_{2}/R_{1}\simeq(t_{2}/t_{1})^{1/2}$. Recalling that $r_{1}\equiv R_{1}/\xi(T_{1})\simeq R_{1}/t_{1}$, we have

[TABLE]

Closed DWs start to collapse at $t_{2}$ as the DW tension overcomes the universe’s expansion.

The collapse of closed DWs is dominated by the axion Lagrangian $\mathcal{L}=1/2(\partial_{\mu}\phi)^{2}-V_{a}$ with $V_{a}$ from Eq. (1). The equation of motion (EoM) is

[TABLE]

where we have incorporated the universe’s expansion. $\mathcal{R}=R/a(t)$ is the co-moving distance. Also, the axion field is redefined as $\tilde{\phi}=\phi/f_{a}$ (dimensionless). For simplicity, we treat closed DWs as nearly spherical, so the EoM is written in the spherically symmetric form. We can use the kink-antikink pair as the initial configuration of spherical DWs [26, 41]

[TABLE]

where the initial scale factor is set as $a(t_{2})=1$. We also assume walls initially at rest, $\dot{\tilde{\phi}}(t=t_{2},\mathcal{R})=0$.

Following the procedure of Ref. [26], we define $E(t,R)$ as the energy contained within a sphere of radius $R$ at time $t$ during collapse of a closed DW. If for some $t$ and $R$, we have $R$ smaller than the corresponding Schwarzschild radius $R_{s}=2GE(t,R)$, a black hole will be formed. The above criterion can be expressed as [26]

[TABLE]

where $S(t,R)\equiv 2E(t,R)/R$ and $m_{\rm P}$ is the Planck mass. By numerically solving the EoM (9) with the initial conditions above, we can obtain the evolution of $S(t,R)$. The detailed numerical calculations are shown in B. The key result is that the maximum $S(t,R)$ is related to the initial collapse size $R_{2}$ by

[TABLE]

where $k_{1}\approx 3.1\times 10^{3}$ and $k_{2}\approx 2.76$. This should be compared with a similar relation in Ref. [26] where $k_{1}\approx 21.9$ and $k_{2}\approx 2.7$. The crucial difference is that in our model closed DWs are originally formed at $T_{1}$ together with the main network and the collapse point $T_{2}$ could be earlier than the QCD transition $T_{c}$ (i.e. $T_{c}<T_{2}<T_{1}$), so the full expression of axion mass Eq. (2) where $m_{a}(T)$ increases rapidly with $T$ before $T_{c}$ must be included in solving the EoM (9). Additionally, our EoM includes the universe’s expansion. In comparison, Ref. [26] considered collapse of fragments from the string-wall network. The fragment process could occur later than $T_{c}$, so $m_{a}$ is treated as a constant there.

Also, fragments in Ref. [26] inherit angular momentum from strings motion, which could significantly suppress PBH formation. However, our model does not suffer from this suppression. Closed DWs have no initial angular momentum at $T_{1}$ since they are formed independently of the main network, and the simple assumption of spherical shape guarantees no angular motion later but only radial motion.

Substituting Eq. (12) into Eq. (11) and using Eq. (8), we can finally express the criterion of PBH formation in terms of $r_{1}$:

[TABLE]

The classical window of current axion mass is $10^{-6}~{}{\rm eV}\lesssim m_{a,0}\lesssim 10^{-2}~{}{\rm eV}$ [65], implying $10^{8}~{}{\rm GeV}\lesssim f_{a}\lesssim 10^{12}~{}{\rm GeV}$ [Eq. (2)]. $r_{\rm 1,min}$ is the minimum radius satisfying the criterion Eq. (13). With $f_{a}$ known, $t_{1}$ and $t_{2}$ are also known from Eqs. (2), (3) and (8), so $r_{\rm 1,min}$ is merely determined by $f_{a}$. In Fig. 1, we plot the relation $r_{\rm 1,min}$-$f_{a}$ (see also B for more numerical details).

4 PBHs as DM

Eq. (13) roughly determines whether a closed axion DW could collapse into a PBH. To exactly calculate the PBH mass, however, we need to answer many complicated questions, e.g. how the PBH as the core alters the wall dynamics and the fraction of the wall falling into the PBH, etc. For simplicity, we estimate the PBH mass as the energy initially stored in the closed wall at $t_{2}$ when it starts to collapse:

[TABLE]

where $\sigma=8f_{a}^{2}m_{a}$ is the DW tension [42].

The PBH mass distribution is related to the size distribution of closed axion DWs Eq. (6) via

[TABLE]

where $\rho_{\rm PBH}(t)$ is the mass density of PBHs. $[T(t)/T_{1}]^{3}$ is the matter density decrease with the universe expanding. We further define $\Omega_{\rm PBH}(t)=\rho_{\rm PBH}(t)/\rho_{\rm cr}(t)$ where $\rho_{\rm cr}(t)=3H^{2}(t)/8\pi G$ is the critical density. $\Omega_{\rm PBH}(t)$ remains constant after the epoch of matter-radiation equality $T_{\rm eq}\approx 0.8$ eV, so the present mass distribution of PBHs is

[TABLE]

By integrating Eq. (16), the present PBH abundance is

[TABLE]

The average mass of PBHs can be calculated as

[TABLE]

which does not change with the universe’s expansion. There is a one-to-one correspondence between $\left<M_{\rm PBH}\right>$ and $f_{a}$. In Fig. 2, we plot PBH mass distributions for different $f_{a}$. We see that PBHs are generally within the mass range $10^{19}$-$10^{29}$ g, but the distribution for each $f_{a}$ is quite narrow centering at $\sim\left<M_{\rm PBH}\right>$ and heavy PBHs are greatly suppressed due to Eq. $(\ref{eq:fr})$.

We emphasize that PBH mass reaching the scale $10^{19}$-$10^{29}$ g is due to the large size of closed DWs which is inversely proportional to the axion mass at $T_{1}\sim$ GeV, i.e. $\xi\simeq m_{a}^{-1}(T_{1})$, rather than the current axion mass $m_{a,0}$. There is a huge difference between $m_{a,0}$ and $m_{a}(T_{1})$. For example, for $m_{a,0}$ as large as $10^{-4}$ eV, we have $m_{a}(T_{1})\sim 10^{-8}$ eV [Eq. (2)]. Another factor contributing to closed DWs size is $r_{1}$ predicted by percolation theory. See also Eq. (14) where $m_{a}^{-1}(T_{1})$ and $r_{1}$ enter the PBH mass expression.

PBHs surviving today contribute to DM with the trivial constraint $\Omega_{\rm PBH}\leq\Omega_{\rm DW}$. Furthermore, various astrophysical observations constrain $\Omega_{\rm PBH}$ for a wide mass window [1, 2]. Most of the valid constraints assume the PBH mass function is monochromatic. Although PBHs in our model have a mass distribution, it is narrow as we see in Fig. 2. If we approximate our model as one which has the monochromatic mass function $M_{\rm PBH}=\left<M_{\rm PBH}\right>$ with the same abundance $\Omega_{\rm PBH}$, the astrophysical constraints on $\Omega_{\rm PBH}$ can be roughly applied to our model.

$\Omega_{\rm PBH}$ in Eq. (17) depends on $f_{a}$ which determines the DWs formation point $t_{1}$ and also the DW tension $\sigma$. Another parameter that also significantly affects $\Omega_{\rm PBH}$ is $\gamma$ [contained in $f(r_{1})$, via Eqs. (6), (7)], $\Omega_{\rm PBH}\propto\gamma$. In Fig. 3, we plot $\Omega_{\rm PBH}/\Omega_{\rm DM}$, the present fraction of PBHs in DM, as a function of $\left<M_{\rm PBH}\right>$ (or $f_{a}$ in the second x-axis, one-to-one corresponding to $\left<M_{\rm PBH}\right>$) for different $\gamma$, with various observational constraints. We see that for $f_{a}\sim 10^{9}$ GeV, PBHs are in the sublunar-mass window $\left<M_{\rm PBH}\right>\sim 10^{20}$-$10^{22}$ g, one of few allowed windows 333Like many other discussions (e.g. Refs. [70, 68]), Fig. 3 does not include the constraint from observations of neutron stars [71] which depends on the controversial assumption of PBHs as DM existing in globular clusters. Many observations disfavor DM existing in such regions, see e.g. Ref. [72].. For the typical value $\gamma=0.1$, PBHs could account for up to $\sim 1\%$ of DW in this mass window. If closed DWs are formed more efficiently, PBHs could contribute more to DM.

We can in turn constrain QCD axion parameter space using the constraints on $\Omega_{\rm PBH}$. Fig. 3 shows that $f_{a}\gtrsim 10^{10}$ GeV is almost excluded, although extremely small $\gamma\lesssim 10^{-3}$ is still plausible resulting in $\Omega_{\rm PBH}\lesssim 10^{-3}\Omega_{\rm DM}$. For $f_{a}\lesssim 10^{8}$ GeV, PBH abundance is very tiny ($f_{a}\lesssim 10^{8}$ GeV is actually excluded by independent observations of supernovae cooling [73]). Our model prefers $f_{a}\sim 10^{9}$ GeV corresponding to $m_{a,0}\sim$ meV (see a similar result in Ref. [27] but depending on a totally different mechanism). Additionally, PBH formation mechanism suggested in this work can also be applied to axion-like particles (ALPs) where $m_{a}$ and $f_{a}$ are not linked. In the ALP case, PBH formation could even be more efficient due to the larger DW sizes since the ALP mass could be lower than $10^{-12}$ eV [74].

5 Conclusions and discussions

We have studied PBH formation from the collapse of closed QCD axion DWs naturally arising when axion mass effectively turns on. PBH mass distribution can be obtained from the size distribution of closed DWs predicted by percolation theory. Our model prefers axion mass at the meV scale (several experiments can detect axion in this mass range, see Ref. [75] for a review). The resulting PBHs are in the sublunar-mass window $10^{20}$-$10^{22}$ g, one of few allowed windows constrained by observations. PBH abundance in our model could vary a lot and it could reach $\sim 1\%$ of DM, where the formation efficiency $\gamma$ of closed DWs plays a key role.

Sublunar-mass PBHs have other significant implications. Ref. [70] suggests that their interactions with neutron stars could solve the long-standing puzzle of r-process nucleosynthesis, which might get indirect supports from aLIGO, aVirgo and KAGRA experiments [76, 77, 78] in the near future. In Fig. 3, r-process is denoted as the dashed line, the region above/below which is the parameter space that fully/partially explains r-process observations [70]. Ref. [79] discussed the possibility of detecting gravitational waves generated by sublunar-mass PBH binaries. Ref. [80] proposed the sublunar-mass PBHs detection through the diffractive microlensing of quasars in long wavelengths with sublunar-mass PBHs as lenses, which could also detect the PBH mass distribution. These experiments might support or exclude our proposal of PBH formation.

Acknowledgments

The work was initiated in the conference IPA 2018 (Interplay between Particle and Astroparticle Physics) in Cincinnati, USA. I thank IPA organizers for this excellent conference. I also thank Ariel Zhitnitsky for useful comments on the work. This work was supported in part by the National Science and Engineering Research Council of Canada and the Four Year Doctoral Fellowship (4YF) of UBC.

Appendix A Survival of the closed axion DWs in the pre-collapse evolution

As we discussed in the main text, closed axion DWs are formed at $T_{1}$ and start to collapse at $T_{2}=T_{1}/r_{1}$ when their sizes are surpassed by the Hubble horizon. The minimum $r_{1}$ required to collapse into PBHs is about $4$ to $14$ for different $f_{a}$ as we see in Fig. 1 in the main text. The pre-collapse evolution refers to the evolution of closed axion DWs from $T_{1}$ to $T_{2}$. During this period, in addition to closed DWs, walls bounded by strings (which we call string-wall objects) are also copiously present in the system (post-inflationary scenario), whose intercommuting with closed DWs might destroy closed DWs [22]. In this section, we are going to study how string-wall objects affect closed DWs and demonstrate that closed DWs will survive against these effects.

The string-wall objects are formed at $T_{1}$ as strings become boundaries of walls. They are like pancakes or large walls with holes [53]. $T_{1}$ can be obtained from Eqs. (2) and (3):

[TABLE]

Another critical time is the time when the domain wall tension dominates over that of strings. We denote the time as $t_{\rm w}$, which is defined by [53, 28, 81]

[TABLE]

where $\sigma\simeq 8f_{a}^{2}m_{a}$ is the wall tension and $\mu\simeq\pi f_{a}^{2}{\rm ln}(f_{a}/m_{a})$ is the energy per unit length of strings [53]. Solving Eq. (20), we get [53]

[TABLE]

which is below $T_{1}$. After $T_{\rm w}$, the dynamics of string-wall objects is dominated by walls, whereas, before $T_{\rm w}$ it is dominated by strings [53]. Thus, the evolutions of the string-wall objects are totally different before and after $T_{\rm w}$, so we should should discuss their effects on closed walls separately.

Before $T_{\rm w}$. In this stage, we have $t<\mu(t)/\sigma(t)$ and strings dominate the dynamics of string-wall objects. The evolution of strings in this stage is no qualitatively different from that before $T_{1}$ when walls have not been formed yet [81]. The main source of strings is closed loops (or wiggles on long strings) with the typical size $t$ [53]. These strings move relativistically and are likely to hit closed walls, which will create holes on walls [22]. However, the holes that are formed in this stage (before $T_{\rm w}$) will shrink and disappear [81]. This is because the force of tension in a string $\sim\mu(t)/t$, is greater than the wall tension $\sigma(t)$, for $t<\mu(t)/\sigma(t)$ [81]. We thus conclude that although the relativistically moving strings may create holes on walls, these holes will disappear themselves as the tension in a string loop can easily overcome the wall tension in this stage.

On the other hand, at $T_{1}$, closed walls with string holes on them could also be formed initially with strings present. This is one of the reasons why $\gamma\lesssim 13\%$ compared to the case without strings. But as we discussed above, these holes tend to disappear themselves in the initial stage, and thus these holey walls initially formed at $T_{1}$ may become closed, which actually brings $\gamma$ closer to $13\%$. This is another thing we can learn from $t<\mu(t)/\sigma(t)$.

After $T_{\rm w}$. The wall tension becomes greater than that of strings. In this stage, if strings hit closed walls and create holes on them, these string holes will inevitably increase in size pulled by the walls, which may significantly decrease the rate of closed walls collapsing into PBHs. However, compared with the first stage, the crucial difference is that the motion of a string after $T_{\rm w}$ is greatly constrained by its own wall originally attached, for the walls dominating the dynamics of the string-wall objects. Also, the string-wall objects will quickly decay into axions [53]. As we will see below, string-wall objects cannot reach the nearest closed walls before these string-wall objects totally decay.

In the first stage (before $T_{\rm w}$), the strings move at relativistic speeds [53]. If a string and a wall collide, the intercommuting probability is very high (close to $1$) [22, 81, 82]. Thus, large closed walls will eat the incoming string-wall objects quickly and efficiently in the first stage (the holes created will disappear as discussed above). With the surrounding regions cleared up, the typical distance between a closed wall surface and the neighbouring string-wall object is the Hubble scale $\sim t$, saturating the requirement of causality444This is also commonly assumed in many related studies of topological defects where the interactions are efficient, see e.g. Refs. [24, 83]. This is also consistent with the numerical simulations of string-wall objects where the wall area parameter $\mathcal{A}\lesssim 1$ [28], implying on average there is one or less horizon-size string-wall object per horizon.. The equilibrium will be kept until $T_{\rm w}$ when the dynamics of string-wall objects is greatly altered. Now at $T_{\rm w}$, for string-wall objects, more energy is stored in walls rather than strings and thus the bulk motion of string-wall objects is determined by walls. We should check what will happen to the system. The simulation result of walls speed is $v\sim 0.4c$ [83]. At $T_{\rm w}$, the distance between a string-wall object and its nearest closed wall surface is $\sim t_{\rm w}$. Then, the time needed for the string-wall object to hit the closed wall can be estimated as

[TABLE]

from which we get

[TABLE]

To obtain $T_{\rm hit}\simeq 0.26T_{\rm 1}$, we also used $T_{\rm w}\simeq 0.6T_{1}$ [Eqs. (19) and (21)].

$T_{\rm hit}$ should be compared with the temperature at which the string-wall objects totally decay. Soon after $T_{\rm w}$, string-wall objects will decay into axions, as the strings pulled by the wall tension quickly unzip the attached walls [53]. Recent simulations show that string-wall objects totally decay at $T_{\rm decay}\simeq T_{1}/3$ [29]555 It is $T_{\rm decay}\simeq T_{1}/4$ obtained in Ref. [30]. However, the exact value of $T_{\rm decay}$ is not essential for us. As we will see below, in the realistic case that $m_{a}(t)$ increases rapidly with time, the wall speed is much lower, which finally leads to Eq. (25).. The crucial point for us is that

[TABLE]

which implies that string-wall objects cannot reach the nearest closed walls before these string-wall objects totally decay into free axions. In other words, closed domain walls will not be destroyed by the string-wall objects after $T_{\rm w}$.

One more comment is that the wall speed $v\sim 0.4c$ obtained in Refs. [83] is relatively high, because they did not consider that the axion mass $m_{a}(t)$ increases with time drastically. With the time-dependent $m_{a}(t)$ taken into consideration, the bulk speed is expected to be lower (even non-relativistic). This could be possibly explained as follows. The speed $v$ is related to the ratio of kinetic energy to rest energy $E_{\rm kin}/E_{\rm rest}$ [83, 84] where $E_{\rm kin}\sim\left<\frac{1}{2}\dot{\phi}^{2}\right>$ and $E_{\rm rest}\sim\left<\frac{1}{2}(\nabla\phi)^{2}+m_{a}^{2}(t)\right>$. With $m_{a}(t)\propto T^{-\beta}$ increasing rapidly, the ratio becomes much lower and so does the wall speed $v$. We could see this picture more intuitively in Fig.2 of Ref. [28], where the simulations show that the string-wall objects are constrained “locally” to decay with almost no bulk motion (close to zero)666The bulk motion should not be confused with the strings motion pulled by the walls. After $T_{\rm w}$, due to the wall tension, a string is accelerated to relativistic speed in the direction of the wall to which it is originally attached (“unzip”) [53]. So the strings motion is constrained locally by the position of walls in the string-wall objects (see e.g. Fig.2 of Ref. [28]). However, the bulk speed of the string-wall objects is low as we have discussed.. Thus, Eq. (24) is quite conservative, and actually we should have

[TABLE]

We conclude this section that closed walls will survive the pre-collapse evolution. Therefore, $\gamma$ formed at $T_{1}$ remains unaffected and becomes important in calculating the PBH abundance.

Appendix B Numerical details of the collapse of closed axion DWs

In this section, we are going to show the details of numerically solving the collapse of closed axion DWs, including how we get the expression of $S_{\rm max}$ as shown in Eq. (12) and also the relation between $r_{\rm 1,min}$ and $f_{a}$ as plotted in Fig. 1 in the main text.

For the convenience of numerical calculations, we define $\tilde{r}=\mathcal{R}/m_{a}^{-1}(t_{2})$ and $\tilde{t}=t/m_{a}^{-1}(t_{2})$ as dimensionless variables, then the EoM Eq. (9) and the initial conditions (Eq. (10) and $\dot{\tilde{\phi}}(t=t_{2},\mathcal{R})=0$) can be written as

[TABLE]

[TABLE]

[TABLE]

where $\tilde{r_{2}}=R_{2}/m_{a}^{-1}(t_{2})$ and $\tilde{t}_{2}=t_{2}/m_{a}^{-1}(t_{2})$ are respectively the rescaled initial radius and rescaled initial time at the starting point of the collapse of closed DWs, consistent with the definitions of $\tilde{r}$ and $\tilde{t}$. Note that $\tilde{r}_{2}=\tilde{t}_{2}$ since $R_{2}=t_{2}$. As we mentioned in the main text, the initial scale factor is set as $1$, $a(\tilde{t}_{2})=1$. In the radiation-dominated era, we have

[TABLE]

If PBHs are formed before the QCD transition $T_{c}$, according to Eq. (2) the axion mass that enters Eq. (26) is

[TABLE]

Later, we will discuss the effect of QCD transition on the collapse of closed axion DWs. As we mentioned in the main text, $\beta\simeq 4$. One of the most recent calculations on axion mass is given by Ref. [35] based on lattice QCD method which shows that the exact value is $\beta=3.925$ 777Ref. [35] does not give the value of $\beta$ directly, but the Supplementary Information of that paper provides the related data. By fitting the data provided, we get $\beta=3.925$..

$E(t,R)$ is defined as the energy contained within a sphere of radius $R$ at time $t$ during collapse of a closed DW, which can be calculated as

[TABLE]

We add the prefactor $1/f_{a}^{2}$ in LHS because $\phi$ is redefined as a dimensionless variable $\tilde{\phi}=\phi/f_{a}$ as we mentioned in the main text. Now, the term $S(t,R)$ related to the criterion of PBH formation can be expressed as

[TABLE]

The maximum value of $S(\tilde{t},\tilde{r})$ during the collapse is

[TABLE]

We see that $S_{\rm max}/f_{a}^{2}$ is a function of $\tilde{r}_{2}$.

We then study the collapse of closed axion DWs by numerically solving Eqs. (26)-(30), from which we obtain the evolution of $S(\tilde{t},\tilde{r})$ (based on Eq. (32)) and further $S_{\rm max}$. We do numerical calculations for different values of the initial radius $\tilde{r}_{2}$, and finally we obtain the relation between $S_{\rm max}/f_{a}^{2}$ and $\tilde{r}_{2}$ which is plotted in Fig. 4. We see that $S_{\rm max}/f_{a}^{2}$ linearly depends on $\tilde{r}_{2}$ in the log-log scale, consistent with Ref. [26] which however did the numerical calculations for a constant $m_{a}$. By fitting the numerical results in Fig. 4, we get

[TABLE]

where $k_{1}=3106.28$ and $k_{2}=2.7626$. In Fig. 5, we also plot the relation between $t_{\rm max}$ and $\tilde{r}_{2}$ where $t_{\rm max}$ is the time when $S(\tilde{t},\tilde{r})$ reaches its maximum value $S_{\rm max}$. The numerical results show that

[TABLE]

We see that the collapse is a very fast process, with the scale factor $a(t)$ only enlarged by $(t_{\rm max}/t_{2})^{1/2}\approx 1.76$ times from $t_{2}$ to $t_{\rm max}$. Similar to Ref. [26], we also observed that $S_{\rm max}$ is reached when the wall collapses to the radius close to zero. So the speed of collapse can be estimated as $(t_{\rm max}/t_{2})^{1/2}t_{2}/(t_{\rm max}-t_{2})\approx 0.84$, close to the speed of light.

Substituting Eq. (34) into the criterion Eq. (11), and using Eqs. (3) and (8), the criterion of PBH formation can be expressed in terms of $r_{1}$:

[TABLE]

Taking equal sign in Eq. (36), we obtain the lowest limit of the size of closed axion DWs at the formation point $t_{1}$ which could finally collapse into PBHs, denoted as $r_{\rm 1,min}$.

However, Eq. (34) is only applicable when the axion mass relation Eq. (30) works, which assumes that $S_{\rm max}$ is reached before QCD transition, i.e. $t_{\rm max}<t_{c}$. Using Eqs. (8) and (35), this condition ($t_{\rm max}<t_{c}$) becomes a constraint on the size of closed DWs at the formation point:

[TABLE]

The interpretation of this relation is straightforward. The larger a closed DW is at $t_{1}$, the later it will collapse according to Eq. (8), so a sufficiently large closed DW will collapse after $T_{c}\simeq 150$ MeV. If Eq. (37) is satisfied, we can substitute the axion mass relation Eq. (30) into Eq. (36) to get

[TABLE]

We see that $r_{\rm 1,min}$ is merely determined by $f_{a}$. The relation between $r_{\rm 1,min}$ and $f_{a}$ is plotted in Fig. 6, denoted as line 1.

For the case $t_{2}>t_{c}$, i.e. closed axion DWs start to collapse after QCD transition, the axion mass that enters the EoM is a constant according to Eq. (2). $t_{2}>t_{c}$ corresponds to the condition $r_{1}>T_{1}/T_{c}$. Ref. [26] numerically solves the collapse of closed axion DWs with $m_{a}$ constant, in which $S_{\rm max}$ has the same form as Eq. (34) but with $k_{1}\approx 21.9$ and $k_{2}\approx 2.7$ 888Although Ref. [26] does not incorporate the effect of the universe’s expansion into the EoM, the results of that paper can still be applied here for constant axion mass. This is because the universe’s expansion plays only a minor role as we see in Eq. (35) where the scale factor is only enlarged by $1.76$ times during the collapse which is a very fast process.. Then, from Eq. (36) we can derive $r_{\rm 1,min}$ in this case:

[TABLE]

We also plot $r_{\rm 1,min}$ in this case as a function of $f_{a}$ in Fig. 6, denoted as the dashed line.

In Fig. 6, we also plot $T_{1}/T_{c}$ and $0.57(T_{1}/T_{c})$ in comparison with Eqs. (38) and (39). Region I (between line 1 and line 2) is the parameter space where the condition Eq. (37) is satisfied, so the criterion Eq. (38) is applicable here and the closed DWs with parameters in this region will finally collapse into PBHs. Region III (beyond line 3) is the parameter space where $r_{1}>T_{1}/T_{c}$ (i.e. $t_{2}>t_{c}$), so we should use the criterion Eq. (39) here. We see that region III is well above the criterion Eq. (39), so the closed DWs with parameters in this region will finally collapse into PBHs. Region II (between line 2 and line 3) where $0.57(T_{1}/T_{c})<r_{1}<T_{1}/T_{c}$ is more subtle. The collapse of closed DWs with parameters in this region will pass through QCD transition, i.e. experience the ‘knee’ of axion mass expression Eq. (2). Since region II satisfies well the criterion of PBH formation from the perspective of both the changing axion mass (Eq. (38)) and the constant axion mass (Eq. (39)), we should expect the closed DWs with parameters in this region will collapse into PBHs 999One may notice that in Fig. 6, the lower three lines (line 1, 2 and the dashed line) intersect with one another at $f_{a}\gtrsim 10^{11}$ GeV and are thus not in good order, which might slightly affect $r_{\rm 1,min}$ in the range $f_{a}\gtrsim 10^{11}$ GeV. However, we may safely ignore the tiny difference since the three lines are very close to each other in this range of $f_{a}$. Also, as we discussed in the main text, the parameter space $f_{a}\gtrsim 10^{11}$ GeV is less interesting since it is almost excluded by observational constraints on $\Omega_{\rm PBH}$. The most interesting part is $f_{a}\sim 10^{9}$ GeV which results in sublunar-mass PBHs, and $r_{\rm 1,min}$ can be well determined for $f_{a}\lesssim 10^{11}$ GeV as we see in Fig. 6..

To conclude, region I, II, and III are all parameter spaces (the shaded region) where closed axion DWs can collapse into PBHs. Thus, the criterion Eq. (38) denoted as line 1 in Fig. 6 is indeed the lowest limit of $r_{1}$ for PBH formation (the tiny difference in the range $f_{a}\gtrsim 10^{11}$ GeV can be ignored as we discussed in footnote 9), which is also plotted in Fig. 1 in the main text. Note that we cannot use Eq. (39) (dashed line) as the final criterion although it is lower than line 1, because the parameter space around the dashed line satisfies the condition Eq. (37) and thus should be checked by the criterion Eq. (38) rather than Eq. (39).