Primordial Black Holes from the QCD axion

We propose a mechanism to generate Primordial Black Holes (PBHs) which is independent of cosmological inflation and occurs slightly below the QCD phase transition. Our setup relies on the collapse of long-lived string-domain wall networks and is naturally realized in QCD axion models with domain wall number $N_{DW}>1$ and Peccei-Quinn symmetry broken after inflation. In our framework, dark matter is mostly composed of axions in the meV mass range along with a small fraction, $\Omega_{\text{PBH}} \gtrsim 10^{-6} \Omega_{\text{CDM}} $ of heavy $M \sim 10^4-10^7 M_\odot$ PBHs. The latter could play a role in alleviating some of the shortcomings of the $\Lambda$CDM model on sub-galactic scales. The scenario has distinct signatures in ongoing axion searches as well as gravitational wave observatories.


INTRODUCTION
The recent detection of gravitational waves emitted by the merging of relatively heavy black holes (M O(10)M ) [1] has revived interest in the proposal that the DM of the universe comprises Primordial Black Holes (PBHs) [2][3][4][5][6]. Although there are constraints on the abundance of PBHs for almost all viable masses (see e.g. [7]), a small relic abundance of heavy (M 10 5 M ) PBHs may play an important role in the generation of cosmological structures and alleviate shortcomings of the standard CDM scenario on sub-galactic scales [8,9]. Furthermore, such PBHs could shed light on the origin of the super-massive black holes (SMBHs) in the centers of most galaxies, some of which were already in place at very early times [10,11].
Several fundamental physics scenarios may explain the existence of PBHs. Arguably, the most studied proposal relies on the gravitational collapse of density fluctuations generated during inflation (see e.g. [7] and refs. therein). Nevertheless, it is interesting to understand whether PBHs could naturally arise in other contexts.
In this Letter we propose an alternative PBH formation mechanism, independent of inflationary physics, that relies on the collapse of axionic topological defects (see e.g. [12] for an introduction). The generation of PBHs from defects has been investigated in different contexts including PBHs from the collapse of string loops [13][14][15][16], and from domain walls (DWs) during inflation [17]. Here we discuss for the first time the formation of PBHs from long-lived string-DW networks [18,19] (see [20] for a setup closely related to ours) appearing in well-known realizations [21][22][23][24] of the Peccei-Quinn (PQ) solution to the strong CP problem [25][26][27]. These so-called hybrid networks have multiple DWs attached to strings and can arise more generally from sequences of phase transitions in the early Universe.
When the PQ symmetry is broken after inflation, the axion abundance receives comparable contributions from the 1) misalignment mechanism, 2) radiation of quanta from string defects [28][29][30], and 3) annihilation of the string-wall network [31]. We show in this Letter that there can be a fourth small contribution to the relic axion DM abundance in the form of heavy, 10 4−7 M , PBHs. Interestingly, this provides a concrete realization of the proposed role of massive PBHs in the early Universe [9] in the context of QCD axion DM. Moreover, our scenario is not subject to some of the strong constraints arising from µ-distortions in the CMB, which plague PBH formation mechanisms from gaussian inflationary fluctuations (see [9] for a recent discussion).
The hybrid network dynamics is, of course, hard to analyze. However, for our purposes the essential features can be captured by focusing on the closed walls that arise in the network [20].

COLLAPSE OF CLOSED DOMAIN WALLS
Once the Hubble length becomes comparable to the closed wall size R , the DW rapidly shrinks because of its own tension. This occurs at the temperature T defined by R ∼ H −1 g eff (T ) −1/2 M p /T 2 , where M p = (8πG N ) −1/2 and g eff (T ) is the effective number of degrees of freedom at T . The total collapsing mass has two contributions: one induced by the wall tension σ, and another one coming from any possible difference in energy density between the two regions separated by the DW: For closed DWs arising in the network ∆ρ ≥ 0 (see below), thus the wall bounds a region of false vacuum. Another important parameter for the formation of PBHs is the ratio of the Schwarzschild radius R S of the collapsing wall to the initial size R : Intuitively, if p is close to 1 then the DW rapidly enters its Schwarzschild radius and forms a BH. If p 1, however, the wall has to contract significantly before falling inside R S . It is then less likely to form a BH, since asphericities, energy losses and/or angular momentum may severely affect the dynamics of the collapse. We will thus refer to p as the figure of merit for PBH formation from the collapse of DWs.
The temperature behavior of p and M is crucial to our proposal. Whenever the tension terms dominate in (1) and (2), we have M ∼ T −4 and p ∼ T −2 . If, instead, the energy difference terms dominate, we have M ∼ T −6 and p ∼ T −4 .
Therefore, the duration of the hybrid network has a huge impact on the likelihood of forming PBHs, as well as on their masses. The use of long-lived string-wall networks is the essential new idea of our proposal. This requires multiple DWs attached to each string [12]. Interestingly, this can be realized in QCD axion models with domain wall number larger than one. 1

AXION DARK MATTER FROM STRING-WALL NETWORKS
Let us now embed the basic mechanism illustrated in the previous section in the cosmological history of the QCD axion. Consider a scalar field Φ with a U (1) P Q symmetry broken at some temperature T P Q after inflation. The field acquires a VEV while its phase is identified with the QCD axion, i.e. Φ = ve ia(t,x)/v , and string defects are formed (see e.g. [12]). Below T P Q , the evolution of the axion is: 1 Most of the energy density in the strings dilutes as ρ strings ∼ µ s H 2 , where µ s is the string tension. 2 In addition, the strings radiate axions [28][29][30][31]. Away from the strings, the homogeneous axion field is frozen because of Hubble friction.
2 At T O(GeV) the QCD phase transition occurs. Non-perturbative effects generate a periodic potential for a where N DW is the model dependent color anomaly, also known as DW number. The periodicity of V is given by 2πF ≡ 2πv/N DW . The dependence of the axion mass m(T ) with temperature is usually parametrized as 1 The original DFSZ [23,24] axion has N DW = 6, while the simplest KSVZ realization [21,22] has N DW = 1. However, generalizations of the latter with N DW > 1 can be considered. 2 See however [29,30] for recent claims of small logarithmic deviations from such scaling regime. follows: where n ≈ 7, T 0 100 MeV are numerical parameters which we take from [31] (see [32] for the original computation and the appendix for more details). Here, m 0 0.01 Λ 2 QCD /F is the zero-temperature axion mass, with Λ QCD 400 MeV.
The potential in (3) leads to the existence of DWs, with tension These become relevant once Hubble friction is comparable to the axion mass, i.e. at the temperature For topological reasons, each string gets attached to N DW DWs at T 1 . Thus, a string-wall network is formed, which also contains closed structures. At the same time, the homogeneous component of the axion field starts to oscillate and generates CDM from the misalignment mechanism.
3 Below T 1 , the energy density of the network is quickly dominated by horizon-size DWs. The subsequent evolution crucially depends on the DW number (see e.g. [12]). If N DW = 1 the network is unstable and rapidly decays, generating an additional contribution to the axion DM abundance. If instead N DW > 1, the network is stable because strings are pulled in different directions by the DWs. The network dilutes more slowly than radiation and matter, and one faces the socalled DW problem [33]. To avoid this cosmological catastrophe, a bias term can be added to the axion potential [34], [31,35] of the form: 3 Notice that the periodicity of (6) is different from (3), in such a way that there is only one global minimum per period 2πv. Furthermore, the phase δ represents a generic offset between the bias term and the QCD potential. The addition of (6) to (3) leads to an energy difference between the false and true minima, ∆ρ A 4 B . This generates pressure which competes against the wall tension and renders the network unstable [34]. Balance between the two competing effects is obtained when σ A 4 B H −1 , which is confirmed by detailed numerical analysis of the network evolution in the presence of a bias term [31]. Most of the network disappears at a temperature: where ∼ O(0.1 − 1) is a parameter which increases with N DW and has been numerically determined in [31].
Depending on their initial size, most of the closed DWs in the network will collapse at different temperatures between T 1 and T 2 . In the process, axions are radiated in such a way that the total axion DM abundance today is given by Ω a = Ω mis + Ω strings + Ω nw . Here, we have included contributions from misalignment, strings and the network. The axion DM abundance has been numerically studied in [31] (see however [28][29][30] for more recent estimates of the string contribution) and we review the dependence of Ω a on F and T 2 in the appendix.
There are two crucial points to take from the discussion above: 1) for N DW > 1 there can be a significant separation between T 1 , the temperature of network formation, and T 2 , the temperature at which its annihilation is efficient, and 2) since (6) lifts the degeneracy of the N DW vacua, closed structures surrounding regions with energy A 4 B can exist in the network. Therefore, from now on we assume N DW > 1.
In Fig. 1 we plot the constraints on the F -T 2 plane for the case N DW = 2. The blue-shaded region is excluded because of DM overproduction. The region in parameter space that is in conflict with constraints from supernovae cooling according to the standard analysis in [36] is shown in gray, while the orange-shaded region displays a more conservative recent estimate [37]. The thick black lines signal the largest allowed value of the offset phase δ in (6) that does not spoil the axion solution to the strong CP problem. We thus conclude that a viable region of parameter space exists, around T 2 5 MeV and corresponding to A B 10 −3 Λ QCD , where no tuning of δ is required. This is in contrast with the conclusion reached in [31], where an overconservative bound on θ QCD was assumed. The untuned region of parameter space is slightly reduced as N DW increases. Nonetheless, even for N DW = 6 only a mild tuning δ ∼ 0.1 is required.
In Fig. 1 we also show the relevant would-be BH masses (red lines) and the figure of merit (dashed lines) for closed DWs which collapse at T T 2 . In the most interesting region of parameter space, we find p ∼ 10 −6 , five orders of magnitude larger than for T ∼ T 1 . This shows the advantage of considering N DW > 1. Nevertheless, p remains quite small and at this point it is unclear whether this leads to a significant fraction of PBHs.

PBHS FROM LATE COLLAPSES
In the previous section we have seen that most of the axionic string-wall network disappears around T 2 . Crucially, T 2 roughly coincides with the temperature when the vacuum energy contribution to M and p starts dominating over the wall tension, as dictated by (1) and (2). Therefore, for closed DWs collapsing at T < T 2 , p increases more steeply, as T −4 .
Let us focus on the region around F 10 9 GeV and T 2 7 MeV in Fig. 1, which leads to the best case scenario for PBH formation. In Fig. 2 we plot the figure of merit (dashed lines) and PBH masses (red lines) for DWs collapsing at T < T 2 . Fig. 2 shows that DWs that collapse roughly when T ∼ 0.1 T 2 are quite likely to lead to the formation of PBHs. These structures only have to contract by one order of magnitude before entering their Schwarzschild radius. Energy losses via radiation of axions as well as the growth of asphericities can be neglected for such short contractions. Indeed, the radiation of energy from a closed spherical Sine-Gordon DW was studied in [38] and shown to become relevant only once the wall has contracted to a size R ∼ R 2/3 m −1/3 0.1 H −1 . Similarly, in the thin wall approximation asphericities do not spoil the formation of PBHs for large p [39]. Furthermore, we have numerically simulated the collapse of Sine-Gordon non-spherical DWs and checked that they can indeed contract down to r min 0.1 R [40]. 4 The resulting PBHs would Let us now estimate the fraction f ≡ Ω P BH /Ω CDM . After T 2 , the energy density of the network is dominated by the bias contribution. However, at any given T < T 2 only a small fraction P nw of the original network survives. Therefore, Assuming that the formation of PBHs occurs mostly at a single temperature T , the actual fraction f is then given by: where N takes into account the effects of asphericities and angular momentum. One should keep in mind that f might be further suppressed by the probability of finding closed structures in the network. In (9), ρ CDM (T ) ∼ ρ CDM (T 2 )(T 2 /T ) 3 is the energy density of CDM at T . In the most interesting region of parameter space in Fig. 1, ρ CDM (T 2 ) is dominated by the contribution from axion quanta radiated by the network. Hence, ρ CDM (T 2 ) ≈ ρ nw (T 2 ) ∼ ∆ρ. Putting everything together, we find: To estimate the actual value of f requires knowledge of P nw . In this respect, the simulations of [31] show that, for 0.5 expect that for large p the inclusion of the latter should not significantly alter our picture. in (7), only 10% of the original network survives at T 2 , while for 0.3 the percentage is reduced to 1%. We do not know the subsequent evolution of the network. Nevertheless, let us assume for simplicity that the network decay follows a power law beyond T 2 : Fitting (11) to the aforementioned results of [31] gives α ≈ 7.
The final fraction f then depends significantly on N , which we do not expect to be large for p close to one. For example, for N 2 the right hand side of (10) peaks at T ∼ T 2 /30, corresponding to p ∼ 1 and giving f ∼ 10 −6 and M ∼ 10 6 M .
However, this result is sensitive to the precise numerical scaling of P nw after T 2 . In this regard, it is interesting to notice that numerical simulations hint at slight deviations from the scaling regime [31]. The decay of the network can then be slower, resulting in smaller α and larger f .
Observations require that f 10 −5 for PBHs with M ∼ 10 6 M [7] (see also [9] for a more recent discussion). As we have seen, this constraint is easily satisfied in our scenario.
Nevertheless, to confidently estimate the actual fraction requires additional numerical studies, which we leave for future work. Let us remark that if a larger f can be obtained with our mechanism, then (non)observations of PBHs may actually give additional constraints on axion models with N DW > 1.

ORIGIN OF THE BIAS TERM
Let us now discuss some ideas related to the bias term. The possibility that (6) arises from an explicit breaking of the PQ symmetry due to Planck-suppressed effective operators has been considered in the literature [41,42]. However, the lore is that gravitational effects would break the PQ symmetry only at the non-perturbative level (i.e. via gravitational instantons, see [43] and [44] for a recent discussion). In this case, the size of the induced potential would be of order A B ∼ M p e −#Mp/F . The corresponding T 2 would then be too small and fall out of the allowed region of parameter space in Fig. 1.
Here, we would like to point out an alternative possibility to generate the bias term. Consider a dark gauge sector, which also breaks the U (1) PQ via anomalies and has DW number N dark DW = 1. The specific matter spectrum and couplings of this hidden sector are not crucial to our discussion, even though cosmological and collider constraints should be checked in any concrete realization. Such a dark sector would then precisely generate a contribution to the axion potential of the form (6), with A B related to the scale of dark gluon condensation, and δ containing the dark sector θ-term.
Interestingly, this naturally allows for a non-trivial temperature dependence of the bias term. Indeed, the scale A B will now have a temperature dependence analogous to the QCD axion potential in (3): with (see also appendix) The natural expectation is that m B will increase as temperature decreases until T 0,B ∼ Λ B , and remain constant afterwards. Here, Λ B is the dark confinement scale and d T , n depend on the dark spectrum. These parameters have an important impact on PBH formation. For instance, for T 2 ∼ 5 MeV and d T , n ∼ 1 the bias term has not yet reached its asymptotic value at T 2 . Therefore, p and M now scale as T −4−n and T −6−n respectively from T 2 to T ∼ T 0,B . A large figure of merit can then be attained in less than one order of magnitude in T , and lighter PBHs may be generated (down to ∼ 10 4 M ). 5 Alternatively, if d T 1 and/or n 6, Λ B roughly coincides with T 2 and we recover the previous case. We leave a more detailed investigation of the dark sector for future work.

CONCLUSIONS
We have discussed a new mechanism to generate PBHs in the context of QCD axion models. It proceeds by the late collapses of closed DWs in a long-lived string-DW network, which arises in QCD axion realizations with N DW > 1 and PQ symmetry broken after inflation.
Lacking accurate knowledge of the network evolution and collapse, we cannot give precise predictions for the fraction and masses of the PBHs. However, under reasonable assumptions, depending on the temperature behavior of the bias term, PBHs with masses in the range M ∼ 10 4 − 10 7 M and representative fraction f 10 −6 can be created. Interestingly, such heavy PBHs can play an important role as seeds for the formation of cosmological structure, alleviating several problems of the CDM scenario on subgalactic scales, and providing an avenue to explain the origin of the SMBHs [8,9].
Our proposal appears to prefer small values of the axion decay constant, F 10 9 GeV, corresponding to axion masses in the meV range. These values are close to the lower bounds from the cooling of supernovae [36,37,45], which are however subject to astrophysical uncertainties and are not universal (see e.g. [46]). On the other hand, small values of F might be observationally interesting. In this respect, it is intriguing that several stellar systems show a mild preference for a non-standard cooling mechanism, which 5 The lines of constant δ in Fig. 1 get modified, but viable regions with δ 0.1 will persist.
can be interpreted in terms of a DFSZ QCD axion [47,48]. In addition, several experiments will be probing this region of QCD axion masses in the near future. In particular: axion helioscopes IAXO [49] and TASTE [50], light-shiningthrough-walls experiments such as ALPS II [51], and the long-range force experiment ARIADNE [52]. Our mechanism might be probed at gravitational wave observatories via the detection of gravitational radiation from: SMBH binaries at LISA [53], the annihilation of the stringwall network [54] at aLIGO (O5) [55], LISA, ET [56] and SKA [57].
Finally, let us mention that similar dynamics to the one described in this Letter can occur for generic Axion-Like-Particles (ALPs), not necessarily coupled to QCD gluons. In this case, it is possible to raise the figure of merit p by considering very light ALPs and thereby delaying the network collapse to temperatures T 100 keV. This scenario may be strongly constrained by the production of extremely heavy PBHs. Office of High Energy Physics, under Award Number No. DE-FG02-91ER40628.

APPENDIX: QCD AXION DARK MATTER
The aim of this appendix is to review the relevant formulae for the total axion dark matter abundance. The material presented here can be partially found in [31] (and refs. therein), together with more detailed explanations.
The total axion dark matter abundance is given by where the three terms on the right hand side of (14) represent respectively the contribution from: the misalignment mechanism, the radiation from axionic strings and the radiation from the string-wall network. Let us first provide formulae for the axion mass. Following [31], we have The parameters c 0 , c T and n can be determined semianalytically using the Dilute Instanton Gas Approximation (DIGA) [58], valid at high temperatures. This approach gives c 0 ≈ 10 −3 , c T ≈ 10 −7 and n ≈ 7 [32] (see also lattice QCD results which agree [59] or deviate [60] from these values). From (15) one finds T 0 100 MeV.
Let us now move to the relic abundance. Firstly, let us focus on the contribution from the misalignment mechanism Ω mis . The QCD axion starts to oscillate at the temperature T 1 given by 3H(T 1 ) = m(T 1 ). By means of (15), we find where A n = (7.5 · 10 16 c T ) 1 4+n . For c T ≈ 10 −7 and n ≈ 7, this gives T 1 ≈ 3 GeV for F 10 9 GeV. The relic abundance is given by B n 0.8 · 10 9 × (2.2 · 10 10 ) − 6+n 4+n . Using the DIGA values for these parameters, the misalignment contribution saturates the observed dark matter abundance for F 10 11 GeV.
Let us now move on to Ω strings . At the moment, there is some controversy in the literature regarding the magnitude of this contribution (see [28][29][30] for recent estimates with a different take on previous calculations). Without entering into details, we focus on the parametric dependence of Ω strings on F Ω strings h 2 C n F 10 9 GeV 6+n 4+n where C n is a numerical prefactor. In order to produce Fig. 1, we have used the formulae provided in [31], where C n ∼ 10 −3 . In this case the contribution from strings is generically larger or comparable to the contribution from the misalignment angle; and the sum of the two contributions saturates the observed dark matter abundance for F ∼ 2 × 10 10 GeV. The precise behavior of Ω strings is not crucial to our proposal, since we are especially interested in the region of small F , where Ω mis and Ω strings represent a subdominant contribution to the total axion abundance. Let us finally discuss the contribution from the string-wall network, which is especially important for our proposal. The crucial difference with respect to the abundances from the misalignment mechanism and from strings is that the network radiates axions at T 2 < T 1 . Thus, their abundance is less diluted and for small values of F it dominates over the other contributions. Assuming a so-called exact scaling regime for the evolution of the network, i.e. ρ nw ∼ σH, we have Ω nw h 2 0.14 × Our expression for Ω nw differs from the one presented in [31] in that we keep the dependence on T 2 , rather than trading it for A B according to (7). Furthermore, there are numerical prefactors in (18) which we have fixed according to the results of [31]. The sum of (16), (17) and (18) generates the solid blue curve in Fig. 1.
Let us now discuss the solid lines of constant δ in Fig. 1. The addition of the bias term (6) misaligns the axion from the CP conserving minimum determined by the QCD potential (3). In particular, the QCD angle θ ≡ a/F at the minimum is approximately given by: Inverting (7) and using ∆ρ A 4 B [1 − cos (2π/N DW )], we find: In order to preserve the solution to the strong CP problem, we require θ min 10 −10 [36]. At constant δ, (19) corresponds to a line in the logarithmic F − T 2 plane, as shown in Fig. 1. Notice that the position of these lines depends on N DW : in particular, the region of phenomenologically viable parameter space where δ ∼ 1 shrinks as we increase N DW . Nevertheless, even for N DW = 6 there is an allowed region of parameter space where only δ ∼ 0.1 is required.