New mechanism for primordial black hole formation from the QCD axion

We present a new mechanism for the primordial black hole (PBH) production within the QCD axion framework. We take the case where the Peccei-Quinn symmetry breaks during inflation, resulting in a $N_{\rm DW}=1$ string-wall network that re-enters horizon sufficiently late. Therefore, closed axion domain walls naturally arising in the network are sufficiently large to collapse into PBHs. Our numerical simulation shows that $\sim 0.3\%$ of the total wall area is in the form of closed walls. In addition, the relic abundance of dark matter is dominantly accounted for by free axions from the collapse of open walls bounded by strings. In this framework, the abundance of PBH within dark matter is calculated to be $\sim 0.9\%$. This fraction remains unaffected by axion parameters or the re-entering horizon temperature, as it is determined by the fixed proportion of closed walls in the network, governed by the principles of percolation theory. The resultant PBHs uniformly share the same mass, which spans from about $10^{-9}$ to $1$ solar mass, corresponding to the classical QCD axion mass window $10^{-5}-10^{-2}$~eV and the re-entering horizon temperature $300-1$~MeV. Intriguingly, PBHs in this mechanism can naturally account for the ultrashort-timescale gravitational microlensing events observed by the OGLE collaboration.

Introduction-QCD axion is an elegant solution to the strong CP problem via the Peccei-Quinn (PQ) mechanism [1][2][3][4][5][6][7][8], which is one of the most promising directions for physics beyond the Standard Model.This mechanism introduces a complex scalar field Φ with a global U (1) symmetry in a potential The symmetry spontaneously breaks when the cosmological temperature falls below the PQ energy scale v a .Axion cosmology is usually studied in two scenarios, depending on the relation between v a and the inflationary Hubble scale H I ; see e.g., Refs.[9][10][11] for a review.If v a ≲ H I /2π, PQ symmetry breaks after inflation.Axion strings first form as the Goldstone field θ ≡ a/v a of Φ winds around from 0 to 2π at boundary, where the angular variable θ is the rescaled dimensionless version of the axion field a. Later, when temperature drops to T 1 ∼ O(GeV), axion mass becomes comparable to Hubble rate, m a (T 1 ) ∼ H(T 1 ), and the axion effective potential from non-perturbative QCD explicitly breaks PQ symmetry down to a discrete Z(N DW ) symmetry with f a = v a /N DW .Here, N DW is the domain wall number.Consequently, axion domain walls form with the pre-existing strings as boundaries.However, for v a ≳ H I /2π, PQ symmetry breaks before inflation, and the formed strings are blown away out of the current Hubble horizon H −1 0 by inflation.Thus, the axion field gets homogenized on the same scale, and the axion topological defects have no observable effects.
In addition to the above post-and pre-inflationary scenarios, PQ symmetry can also break during inflation.This can happen naturally if the PQ breaking is driven by inflation and can be achieved in various ways (e.g., Refs.[12,13]; see also Ref. [14]).For example, one can couple the PQ field Φ to the inflaton field ϕ via cϕ 2 Φ † Φ where c is the coupling.Then, during inflation, the excursion of ϕ from large field induces PQ breaking when ϕ rolls down to ϕ = λ/2cv a .This example works for v a > H I /2π where the finite temperature correction due to H I will not affect the breaking scale.Also, for v a much larger than H I , the fluctuation δθ ∝ H I /v a during inflation is greatly suppressed [13], which will not disturb the correlation length of topological defects.
In the during-inflationary scenario, the string-wall network re-enters horizon at temperature T en .Supposing that PQ breaks after N PQ e-foldings of visible inflation, we have T en ≃ T 0 e NPQ+4 [13] where T 0 is the current cosmological temperature.T en should be higher than ∼ 1 MeV to avoid possible violations of BBN processes.On the other hand, T en should be smaller than T 1 ∼ 1 GeV, such that walls form when strings are still super-horizon.Otherwise, it makes no difference with the post-inflationary scenario.Since the network is superhorizon at T 1 , it will not collapse until T en .Dominated by wall energy, the network collapses right after T en as strings eat walls rapidly under the wall tension [12].This holds for the case N DW = 1 which naturally avoids the domain wall problem.Since such collapse can be much later than that in the post-inflation scenario, it opens new parameter space for the resultant free axions accounting for dark matter [12,13].
In addition to walls bounded by strings, a small portion of the network is in the form of closed walls [15,16].Closed walls naturally arise due to the same Kibble-Zurek mechanism [17,18] when one vacuum in space is fully surrounded by another topologically different vacuum.Compared with the walls bounded by strings, closed walls received less attention in the literature.However, closed walls can play an important role in the early universe by collapsing into primordial black holes (PBHs) [19,20].
In this study, we will show that in the duringinflationary scenario, the closed walls collapse into PBHs naturally for most of the parameter space after re-entering horizon at T en .Closed walls are big enough because T en is sufficiently late.It overcomes the difficulty in Refs.[19,20] of finding large closed walls for forming PBHs.Moreover, this mechanism differs fundamentally from the following ones, including Ref. [21] which requires the presence of false vacuum induced by a bias potential to form PBHs and requires N DW ≥ 2, Ref. [22] which needs closed walls generated by fluctuations of axion field, and Refs.[23][24][25][26] which need the generation of overdensity fluctuations of particle axions to form PBHs.
In our scenario, PBH abundance in dark matter is 0.92% ± 0.05% and the remaining dark matter is accounted for by free axions from walls bounded by strings.This proportion is fixed, independent of axion parameters or T en , because the proportion of closed walls in the network is fixed, which can be determined by numerical simulations.The resultant PBH mass is monochromatic, which can be 10 −9 − 1 M ⊙ for the classical QCD axion mass window 10 −5 − 10 −2 eV.We want to emphasize that apart from assuming that PQ symmetry breaks during inflation, we do not make any changes or fine-tunings to the standard QCD axion cosmology, which is remarkable that PBHs can form in this simple picture.Furthermore, such PBHs can potentially explain the gravitational microlensing events observed by the OGLE collaboration [27,28] and also serve as a candidate for the intriguing Planet 9 in the solar system if it is interpreted as a PBH [29,30].
Closed axion domain walls-We are going to simulate the string-wall network formation for the case N DW = 1, particularly focusing on the proportion of closed walls.First of all, we stress that a notoriously difficult simulation problem in the post-inflationary scenario does not arise here.The problem is essentially a multi-scale problem that strings and walls are formed at drastically different times, f a ≫ T 1 , making it impossible for computers to fully simulate the long-time string evolution before T 1 unless taking some compromises [31][32][33][34][35][36][37][38][39][40][41][42].In our scenario, strings are frozen in the super-horizon scale without dynamics.Therefore, it is equivalent to treating walls and strings as forming simultaneously in the sense that string evolution is trivial before wall formation.
One easy way to obtain the network structure is as follows [15,16].Basically, the vacuum circle θ ∈ [0, 2π] is divided into three equal-length parts: part I, centering at 0; part II, centering at 2π/3; and part III, centering at 4π/3.After PQ symmetry breaking, each correlated patch randomly chooses one of three parts with equal probability.A string forms at the center if the surrounding correlated patches are assigned part numbers in clockwise (or anti-clockwise) order.A wall forms at the boundary of two correlated patches if they are assigned with II and III respectively.Consequently, the network contains two types of structures: open walls bounded by strings (each string connecting to a single piece of wall, N DW = 1) and closed walls.We apply a large N 3 cubic lattice simulation in the realistic 3D case, and each cube represents a correlated Hubble patch.The simulation is purely a mathematical problem related to percolation theory, independent of specific axion parameters.For illustrative purposes, we show the simulation results (part) in Fig. 1.Although simple, this simulation captures the key features: the network is dominated by open walls bounded by strings and a small portion is closed walls ("bubbles").
Then, we generalize the above simple model to a more realistic case where continuous variables are used in generating initial conditions; see e.g., Refs.[31,40,[43][44][45].The initial values of the PQ field Φ = (ϕ 1 + iϕ 2 )/ √ 2, can be assigned based on the correlation functions in momentum space during inflation.By applying the Fourier transformation to ϕ i in momentum space, one can obtain their values on each lattice site in the position space.The lattice spacing is taken as the correlation length, ξ = 1/H I .θ = arg(Φ) is continuous and domain walls occur at θ = π.We focus on two parameters that can be derived from the statics of simulation: the area parameter A which is the number of walls per lattice cell on average, and γ, the ratio of total closed wall area to the total wall area.As we increase the simulation size N 3 , A and γ a become stable.We adopt N 3 = 512 3 , and we obtain A ≈ 37% ± 0.007%, γ = 0.285% ± 0.002%.(1) with relatively small uncertainties.Also, most of the bubbles have a volume of one cell.Two-cell or larger bubbles are exponentially rare [46] (even rarer with the strings present), which we neglect for now.We refer the readers to the supplemental material (SUPP) for more simulation details.
PBH formation-We are going to demonstrate the criterion for a closed wall to collapse into a PBH.The stringwall network, including walls bounded by strings and bubble walls, re-enters horizon at T en .The diameter of a bubble wall is 2R 0 ≃ H −1 (T en ) identical to the distance between strings.A simple argument for PBH formation is that the bubble radius is smaller than its Schwarzchild radius R S = 2GM [19,20,47].Bubbles lose energy by radiating free axions during collapse, which only becomes significant when their size R ≲ m −1 a [19,20,48].Therefore, the only criterion we need for PBH formation is that, when the bubble shrinks to We have used M = 4πR 2 0 σ w with the wall tension σ w = 9.23f 2 a m a [49,50].Also, we have used the fixed relation between axion mass m a and PQ energy scale , where χ 0 ≈ (75.6 MeV) 4 is the (zerotemperature) QCD topological susceptibility [51].The derived condition ( 2) is quite simple and independent of axion parameters.It only demands the network (therefore, the bubbles) re-enter horizon later than 89. 4 MeV.Note that bubble collapse indeed happens in this picture rather than expands because the gravitational repulsion of walls is negligible compared to the energy of enclosed cosmic plasma [52].
Eq. ( 2) is a rather conservative estimate.A much more accurate approach is to numerically solve the equation of motion of the closed axion wall field a [19,20].We define E(r, t) as the energy enclosed within a radius r at a time t during collapse.Then, a PBH can form if 2GE(r, t) ≳ r is satisfied for some r and t during the collapse.Furthermore, we define S(r, t) ≡ 2E(r, t)/r, and S max ≡ max (r,t) (2E(r, t)/r) is the maximum S achieved during collapse.Therefore, PBH formation criterion can be equivalently expressed as S max ≳ m 2 P where m P is the Planck mass.Our numerical results show that S max is a linear function of the initial bubble radius R 0 in the logarithmic scale, This is in the piecewise form because QCD axion mass is a function of T which can be approximated as [51], One of the most updated calculations [51] indicates that T c ≈ 150MeV near QCD transition and n ≃ 4. Initially,  axion mass increases rapidly as temperature drops.After T c , it settles to constant m a .Our result ( 3) is consistent with that in Refs.[19,20], except a small discrepancy between the prefactor in the first branch and that in Ref. [19] because we have included the Universe's expansion in solving the bubble dynamics.Plugging Eq. ( 3) into the criterion S max ≳ m 2 P , we get the parameter space of (T en , f a ) that allows PBH formation, shown in Fig. 2. Also, the resultant PBH mass can be estimated as M PBH ≃ 4πR 2 0 • σ w (T en ).From Fig. 2, PBH mass range in this scenario is The later the network re-enters horizon, the heavier the resultant PBH will be.The lower bound in Eq. ( 5) corresponds to T en ∼ 300 MeV, below which PBH can form, while the higher bound corresponds to the lowest temperature of re-entering horizon, T en ∼ 1 MeV, which is required by not disturbing BBN observations 1 .This mass range ( 5) is drawn for the case that free axions from the decay of open walls bounded by strings explain the dark matter relic abundance (see Eq. ( 7) below).Note that PBHs obtained here are light, so their mass accretion after formation is negligible (see e.g., Ref. [52]).
If free axions only account for a fraction of dark matter, lighter PBHs can form, and the PBH formation criterion can be extended up to T en ∼ 400 MeV where it reaches the constraint f a ≳ 10 8 GeV from Supernova 1987A [56,57], as can be seen in Fig. 2.
PBH abundance-In addition to closed walls, walls bounded by strings collapse immediately after T en , since strings are pulled toward each other by the wall tension.This is because wall energy dominates the system, σH −1 (T en )/f 2 a ∼ m a H −1 (T en ) ≫ 1.Such collapse of walls bounded by strings releases energy mainly in the form of free axions 2 , with the energy density estimated as the wall energy density, ρ a (T en ) ≃ ρ w (T en ) ≃ Aσ(T en )H(T en ).Furthermore, numerical simulations show that the resultant free axions are mildly relativistic, with a Lorentz factor γ a ≈ 3.23 ± 0.18 [32,33].Thus, the present-day energy density of free axions is a is the scale factor.In deriving the last step, we have used the relations σ w = 9.23f 2 a m a [49,50] and χ 0 = m 2 a f 2 a for QCD axion.Then, the present-day ratio of the axion energy density to the critical ρ cri can be calculated as f a 10 9 GeV 20 MeV T en .
(7) g * (T en ) is the effective numbers of relativistic degrees of freedom at T en .Compared with Ref. [12], we have included γ a and also the area parameter A. Axions can fully account for the dark matter if Ω a h 2 = 0.12, corresponding to the red band in Fig. 2. In comparison, the misalignment mechanism contributes in this scenario with the number density n a,mis (T 1 ) ∼ m a (T 1 )f 2 a θ 2 , which is negligible since m a (T en )n a,mis (T en )/ρ a (T en ) ∼ 0.1T en /T 1 ≪ 1.
We can easily express the PBH abundance as If axions from open walls can fully explain dark matter, Ω a = Ω DM , then PBH abundance is fixed as f PBH = γγ a = 0.92% ± 0.05%.We emphasize that γ is insensitive to axion parameters or T en .It is determined by the ratio of closed walls to total walls, which is a constant predicted by percolation theory.The possible uncertainty from γ a simulations [32,33] has been included.
2 A tiny portion of energy is released as gravitational waves (GWs).This generates GW spectra drastically different from that in the network scaling regime [58].For certain parameter space of axion-like particles, such GW spectra can possibly explain the reported nano-Hertz stochastic GW background and can be tested by various GW interferometry experiments [58].GeV.The shaded elliptic region represents the PBH explanation to the OGLE events [27,28].The existing PBH constraints in the shaded gray region mainly come from the lensing [27,28,59,60] under the null hypothesis that there is no PBH.The projected sensitivities are from the future ET by observing PBH binaries [61].
Our PBH mass spectrum is monochromatic, since closed walls have a uniform volume, enclosing a single cube of correlated patch.We show f PBH in Fig. 3, together with the constraints on f a , T en , and also the PBH formation criterion mapped from Fig. 2.
In comparison, we show the existing constraints on f PBH based on gravitational microlensing observations [27,28,59,60].Besides, GWs generated from PBH binaries can be probed in principle by GW detectors, such as LIGO-Virgo [62], Tianqin [63], Taiji [64] and LISA [64], etc.For our planetary-mass PBHs, they can be potentially probed by the prospective Einstein Telescope (ET) via searching for GWs from binaries in our galaxy and in the solar system vicinity [61].
More interestingly, OGLE collaboration has observed 6 ultrashort-timescale microlensing events [28] which can be well reproduced by PBHs [27].As we can see from Fig. 3, PBHs required to explain such events partly overlap with the PBHs generated in our mechanism.Hence, from a positive perspective, observations by OGLE [27,28] might have potentially revealed supporting evidence for the existence of PBHs proposed in this work.
Another intriguing possibility is that Planet 9 in the solar system can be a PBH with mass ∼ 10 −5 M ⊙ , which explains the growing observational anomalies of trans-Neptunian objects' orbits [29].Such a PBH can be naturally generated in our mechanism.This idea is testable via searching for signals from PBH interactions with the surrounding environment [29,65] or a wilder proposal of sending a small laser-launched spacecraft to approach the PBH [30].
Lastly, it should be noted that the numerical values γ a are different given by different groups [12,16,32,33].The string core's size is f −1 a while the wall thickness is ∼ m −1 a .This multi-scale (m a /f a ∼ 10 −25 ) makes it difficult for lattice simulations of the dynamical collapse of open walls bounded by strings [16].One of the most updated simulations regarding γ a [32,33] gives γ a ≈ 3.23 ± 0.18 by extrapolating the simulation results at the artificial value of Λ QCD /f a = 0.275 ∼ 0.4 to the physical value of 10 −12 (Λ QCD ∼ 400 MeV).We have adopted this γ a in making Figs. 2 and 3, but one can simply shift the results if using other values of γ a .For example, an old simulation [16] obtains γ a ≃ 60 by extrapolating the simulation results at the artificial ln(f a /m a ) ∼ 4.6 to the physical ln(f a /m a ) ∼ 60.This would give f PBH = γγ a ∼ 17% and the PBH mass range as ∼ 4 × 10 −8 − 21M ⊙ .One caveat we would like to point out is that Refs.[16,32,33] carried out simulations for the traditional post-inflationary scenario, while our during-inflation scenario is different where the network re-enters horizon and collapse at a much later time, m a ≫ H(T en ).Therefore, it is possible that the level of uncertainty in calculating the PBH abundance caused by γ a has been underestimated.Besides, Ref. [12] simply assumes the typical momentum of free axions to be the Hubble scale, H(T en ), and therefore γ a ≃ 1.This would give f PBH = γγ a ∼ 0.3% and the PBH mass range as ∼ 2 × 10 −10 − 0.35M ⊙ .Simulation of γ a is well beyond the topic of PBH formation in this work, since it is related to open walls bounded by strings while PBHs are from closed walls.
Conclusions-A new possible mechanism is introduced for producing PBHs within the QCD axion framework.In this mechanism, PQ symmetry breaks during inflation, so the N DW = 1 axion topological defects re-enter horizon and collapse at a temperature sufficiently later than the post-inflationary case.Closed axion domain walls are sufficiently large to collapse into PBHs.Numerical lattice simulation shows that ∼ 0.3% of walls are closed walls.This outcome remains independent of axion parameters and the re-entering horizon temperature T en , with its determination rooted in the mathematical principles of percolation theory.
Dark matter is accounted for by free axions from the collapse of open walls bounded by strings.If we adopt the Lorentz factor of free axions simulated by Refs [32,33], our PBHs constitute ∼ 0.9% of dark matter, independent of T en .The resultant PBH mass is monochromatic, which is between ∼ 10 −9 − 1M ⊙ for the classical QCD axion window m a ∼ 10 −5 −10 −2 eV, depending on T en ∼ 300− 1 MeV.If axions only partially account for dark matter, more parameter space is released, allowing lighter PBHs to form and extending the PBH formation criterion to T en ≲ 400 MeV.
We have also discussed some cosmological and astronomical implications.Without fine-tunings, PBHs produced in this mechanism can naturally explain the microlensing events observed by OGLE [27,28].They could also be a candidate for the PBH-origin Planet 9 in our solar system [29,30].Additionally, binaries of such PBHs can be potentially probed by GW detectors such as the prospective Einstein Telescope [61].
Acknowledgment-We would like to thank Haipeng An, Heling Deng, Bin Guo, Huai-Ke Guo, Yong Tang, and Lian-Tao Wang for useful discussions.

Supplemental Material
This Supplemental Material provides more details of the model and simulations.

A. Equations of Motion
As described in the main text, the Lagrangian of the model can be written as (9) c is the coupling between the PQ field Φ and the inflaton ϕ.The inflation drives the PQ symmetry breaking when ϕ rolls down to the value of λ/2cv a .
After the PQ symmetry breaks during inflation, the initial conditions will be stretched out of the horizon by inflation and become homogenized over enormous distances.The system is kept frozen until the network reenters the horizon.This is valid as long as the spatial derivative is ignorable.This can be seen from the equation of motion of the PQ field Φ = (ϕ a(t) is the scale factor and ∇ is the derivative with respect to the co-moving spatial coordinates.At the PQ symmetry breaking, the correlation length is ∼ H I , so we can take the initial values in a Hubble radius to be uniform, and they will be further stretched out of horizon exponentially due to inflation.Due to this homogeneity, the spatial derivative suppressed by a(t) −2 can be neglected during the evolution out of horizon.Later on, the axion mass m a effectively turns induced by non-perturbative QCD effects which becomes m a ∼ H at ∼ 1 GeV.This overcomes the Hubble friction proportional to 3H, and the axion field in one homogeneous patch (still out of horizon) starts to roll down (misalignment mechanism) either to 0 or 2π determined by the initial values in our simulations.Our final results may be modified by the following effects: the closed walls are not perfectly spherical, the closed walls re-enters the horizon not exactly at the same time, etc.However, these effects are subdominant and will not destroy the robust steps of forming PBHs described in the main text.
Next, we comment on the generality of the model.For the PQ symmetry breaking to be driven by inflation, the assumption here is that Φ couples to the ϕ.One simple realization is via cϕ 2 Φ † Φ as shown in Eq. ( 9).Such coupling has been widely discussed in the literature.It dates back to early literature of axion cosmology; see e.g., Ref. [66] (Section 4.3), which indicated that this coupling is more realistic to occur than the case of no coupling.Ref. [66] also discussed the case that the PQ symmetry breaks during inflation at ϕ = ϕ c = λ/2cv a , which is exactly what we need for our setting 3 .It implies that v a is significantly larger than H I to suppress the thermal correction, which ensures that the symmetry breaking is driven by the inflaton field at ϕ c .In addition, we naturally benefit from the fact that the subsequent axion fluctuation ∝ H/v a is suppressed.Therefore, we make no non-trivial assumptions beyond the classical setting in Ref. [66].In addition to Ref. [66], the couplings between the PQ field and the inflaton field, represented by cϕ 2 Φ † Φ and its variants, have also been commonly explored in many different contexts and for many different purposes; see e.g., Refs.[13,14,[67][68][69][70][71][72][73].Our novelty lies in the PBH formation within this basic framework.
Additionally, our model can be realized not necessarily by the coupling above.For example, the inflationary Hubble scale H I may also drive the PQ symmetry breaking if it varies significantly during inflation (see e.g., Ref. [13]).This may be realized via the coupling between Φ and the Ricci scalar, −ξRΦ † Φ = 2H 2 I Φ † Φ, where the Ricci scalar R = −12H 2 I and a natural choice of ξ is 1/6 with conformal symmetry imposed.
We extend the simulation with three discrete variables [15,16] to a more realistic one with continuous variables, by adopting a well-studied lattice simulation of strings and domain walls based on, for example, Refs.[40,43,44].In our scenario, strings and initial values are stretched and blown out of the horizon during the inflation.As a result, the spatial dynamics become less significant, as shown in the section above.As a result, the general two-stage simulations (PQ era and QCD era) become one-stage, where only the initial values of the field configuration in the PQ era play an important role.
We write the PQ field as Φ = (ϕ 1 + iϕ 2 )/ √ 2. During inflation, the amplitude distribution of the real scalars, ϕ 1 and ϕ 2 , in momentum space at the finite temperature T I = H I 2π can be derived as, where k is the four-momentum and a(t) is the scale factor.m eff = 2cϕ 2 − λv 2 a is the effective mass of the scalars ϕ 1 and ϕ 2 .When the inflaton field ϕ rolls down to the value of λ/2cv a , PQ symmetry breaks.We have omitted the finite-temperature correction of m 2 eff which is λT 2 I /3.This is because T I ≪ v a such that the inflation drives the symmetry breaking while the temperature is less important [12].n k is the Bose-Einstein distribution at T I .The correlation function of the PQ fields ϕ i in the momentum space can be written as In the discrete space, this can be expressed as [43,44] where N 3 is the lattice size and ∆x phy is the physical lattice spacing.We then can generate the Gaussian random distribution of ϕ i (k) in momentum space with the variance as ⟨|ϕ i (k))| 2 ⟩.Physically, the momentum k suffers from the IR and UV truncation based on the lattice, which requires that |k| ∈ [ 2π N •∆x phy , √ 3π ∆x phy ] [44].Besides, one also needs to generate two random phases that follow uniform distribution in the range [0, 2π) (For more details, one can refer to Sec. 7.1 in Ref. [44]).Finally, by applying the discrete Fourier transform to the field in momentum space, one can obtain the initial values of the scalar fields ϕ i in the real space.
With the initial values of ϕ 1 and ϕ 2 at each lattice site known, one can calculate the angular variable θ (axion field rescaled by v a ) at the site as ϕ1 , ϕ 1 > 0 and ϕ 2 > 0; arctan ϕ2 ϕ1 + π, ϕ 1 < 0; arctan ϕ2 ϕ1 + 2π, ϕ 1 > 0 and ϕ 2 < 0. ( 14) Domain walls will form at the positions of θ = π.Equivalently speaking, a domain wall intersects the link between two sites if ϕ 1 < 0 and ϕ 2 has different signs at the two sites (see e.g., Refs.[31,43]).A closed wall can form if a site and all of its neighboring sites satisfy this condition.Following this method, one can count the number of walls and closed walls, from which one can further obtain the area parameter A and the ratio of the total closed wall area to the total wall area, γ.
There are some differences between our simulations and those in the post-inflationary scenario (PQ symmetry breaks after inflation).In the post-inflationary scenario, the simulations usually set the power spectrum of ϕ i as [45]: In comparison, in our during-inflation scenario, the power spectrum Eq. ( 11) has an extra term (H from the superhorizon effect which becomes important during inflation.In fact, this effect is very common in dealing with the initial fluctuations during inflation (see e.g., Chapter 8 in Ref. [74]).For subhorizon, k > H I , Eq. ( 11) becomes which is the same as Eq. ( 15).While for superhorizon, k < H I , Eq. ( 11) becomes The physical meaning is clear, especially in the position space.For the case within the horizon, the correlation function, ⟨Φ(x)Φ(y)⟩, decreases with the increase of their distance, as if in the Minkowski space.When the distance exceeds the horizon length, the correlation function is determined by this Hubble length scale.For the traditional post-inflationary scenario, the typical momentum of the fluctuation is k ∼ T , which is much larger than the Hubble constant H ∝ T 2 /M P .As a result, the superhorizon term can be ignored.While for our case where T I ∼ H I , the distance between two points will be persistently stretched out of horizon much more efficiently by the e-folding expansion of inflation, so the superhorizon term cannot be neglected.As the inflaton field slowly rolls down, PQ symmetry breaks as m eff = 2cϕ 2 − λv 2 a approaches zero.Therefore, at the phase transition, m eff ≲ H I and the correlation is determined by H −1 I rather than m −1 eff , and the lattice spacing can be properly set as H −1 I .We perform the simulations in different lattice sizes starting from N 3 = 100 3 .For each size, we repeat the simulation ten times to get the uncertainties of A and γ.The results are shown in Tab.I and Fig. 4 with error bars.We can see that both the values of A and γ become stable as the lattice size increases, with the uncertainties dropping progressively.Finally, we adopt the values at N 3 = 512 3 which are A ≈ 37% and γ ≈ 0.285% with relatively small uncertainties.
A(%) γ(%) 100 3    Next, we comment on the size distribution of closed walls.The percolation predicts that the number of clusters (closed walls) decreases exponentially with the size increasing [46].Large clusters are even rarer with the strings present.The probability for a lattice cell to be occupied by a closed wall is ∼ γ.Then, the probability for a two-cell closed wall to occur can be simply estimated as ∼ γ 2 which is much smaller.Therefore, the number density of closed walls larger than one cell is tiny, which has been neglected in our further discussion.The PBH mass spectrum is monochromatic because the lattice cell is uniform in volume, which is an artificial effect.Indeed, in the realistic case, the spectrum is not perfectly monochromatic.We expect the volume of closed walls to follow a Gaussian-like distribution centering at one cell.The dispersion (variance) of the distribution should be small due to the following reasons.On the one hand, the initial sizes of closed walls bigger than one cell are rare, and the network dynamics during inflation are trivial which does not affect the initial conditions.On the other hand, after re-entering horizon, the collapse of open walls bounded by strings would not significantly alter the dynamics of closed walls, as these two kinds of objects are separate and there are no strings intersecting with closed walls.Also, the collapse is very quick, which can be completed in about one Hubble time, H(T en ) −1 , which has less effect on the initial bubble size distribution compared to the post-inflationary scenario with a scaling-regime evolution.In this sense, the PBH mass in Fig. 2 and Fig. 3 in the main text should be regarded as the average mass, although the dispersion is narrow which is close to a monochromatic distribution.

Figure 1 .
Figure 1.Simulation result of the string-wall network, based on the simple method in Ref. [15, 16].We show a part where a closed wall emerges.A closed wall is highlighted in red.Gray translucent slices and orange lines represent open domain walls and strings, respectively.

10 -11 10 - 1 10 - 3 10 - 2 10 - 1 MFigure 3 .
Figure 3. Plot of PBH abundance and its monochromatic mass.The solid red band denotes fPBH ≃ 0.92% ± 0.05%, when the dark matter relic abundance is accounted for by free axions from open walls bounded by strings.The corresponding mass range is given in Eq. (5).The three labeled dashed red lines are mapped from the requirements shown in Fig. 2, which are: ① PBH formation criterion; ② Ten ≳ 1 MeV; and ③ fa ≳ 10 8 GeV.The shaded elliptic region represents the

Figure 4 .
Figure 4. Simulation results of A and γ.The lattice size is 512 3 .

1 2
Figure 2. Parameter space for PBH formation.The shaded brown region is the parameter space where PBH cannot form.Additionally, the red band represents that free axions from the decay of walls bounded by strings fully account for dark matter.The shaded gray region is excluded because of dark matter overproduction.The blue lines are the PBH mass contours.
The work of SG is supported by the National Natural Science Foundation of China (NSFC) under Grant No. 12247147, the International Postdoctoral Exchange Fellowship Program, and the Boya Postdoctoral Fellowship of Peking University.The work of JL is supported by NSFC under Grant No. 12075005, 12235001.