PHYSICAL REVIEW LETTERS 122, 101301 (2019) Primordial Black Holes from the QCD Axion Francesc Ferrer,1,2 Eduard Masso,2,3 Giuliano Panico,2,4 Oriol Pujolas,2 and Fabrizio Rompineve2 1 Department of Physics, McDonnell Center for the Space Sciences, Washington University, St. Louis, Missouri 63130, USA 2 IFAE and BIST, Universitat Autònoma de Barcelona, E-08193 Bellaterra, Barcelona, Spain 3 Departament de Física, Universitat Autònoma de Barcelona, E-08193 Bellaterra, Barcelona, Spain 4 Laboratory of High Energy and Computational Physics, National Institute of Chemical Physics and Biophysics, Rävala pst. 10, 10143 Tallinn, Estonia (Received 24 July 2018; revised manuscript received 31 October 2018; published 13 March 2019) We propose a mechanism to generate primordial black holes (PBHs) that 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, ΩPBH ≳ 10−6 ΩCDM of heavy M ∼ 104 − 107 M ⊙ PBHs. The latter could play a role in alleviating some of the shortcomings of the standard cosmological model on subgalactic scales. The scenario might have distinct signatures in ongoing axion searches as well as gravitational wave observatories. DOI: 10.1103/PhysRevLett.122.101301 Introduction.—The recent detection of gravitational including PBHs from the collapse of string loops waves emitted by the merging of relatively heavy black [14–17], and from domain walls (DWs) during inflation holes [M ≳ Oð10ÞM⊙ ] [1] has revived interest in the [18,19]. Here we discuss for the first time the formation proposal that the dark matter (DM) of the Universe com- of PBHs from long-lived string-DW networks [20,21] prises primordial black holes (PBHs) [2–7]. Although there (see Ref. [22] for a setup closely related to ours and are constraints on the abundance of PBHs for almost all Refs. [23–25] for similar previous work) appearing in viable masses (see, e.g., Ref. [8]), a small relic abundance of well-known realizations [26–29] of the Peccei-Quinn heavy (M ≳ 105 M⊙ ) PBHs may play an important role in (PQ) solution to the strong CP problem [30–32]. These the generation of cosmological structures and alleviate so-called hybrid networks have multiple DWs attached shortcomings of the standard cold dark matter (CDM) to strings [33] and suffer from a DW problem [34] unless scenario on subgalactic scales [9,10]. Furthermore, such the vacua separating different walls are split [35]. The PBHs could shed light on the origin of the supermassive splitting likely requires extra new physics beyond the QCD black holes (BHs) in the centers of most galaxies, some of axion (see below); however this need not interfere with which were already in place at very early times [11,12]. the present mechanism of PBH formation. Several fundamental physics scenarios may explain When the PQ symmetry is broken after inflation, the the existence of PBHs. Arguably, the most studied axion abundance receives comparable contributions from proposal relies on the gravitational collapse of density the (i) misalignment mechanism, (ii) radiation from string fluctuations generated during inflation (see, e.g., Ref. [8]). defects [36–38], and (iii) annihilation of the string-wall Nevertheless, it is interesting to understand whether PBHs network [39] (see also Ref. [40]). We show in this Letter could naturally arise in other contexts. that there can be a fourth small contribution to the axion In this Letter we propose an alternative PBH formation DM abundance in the form of heavy, 104−7 M ⊙ , PBHs. mechanism, independent of inflationary physics, that relies Interestingly, this provides a concrete realization of the on the collapse of axionic topological defects (see, e.g., proposed role of massive PBHs in the early Universe [10] Ref. [13] for an introduction). The generation of PBHs in the context of QCD axion DM. Moreover, our scenario is from defects has been investigated in different contexts not subject to some of the strong constraints arising from CMB μ distortions, which plague PBH formation mechanisms from Gaussian inflationary fluctuations (see Ref. [10]). Published by the American Physical Society under the terms of The hybrid network dynamics is hard to analyze. the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to However, for our purposes the essential features can be the author(s) and the published article’s title, journal citation, captured by focusing on the closed walls that arise in the and DOI. Funded by SCOAP3. network [22]. 0031-9007=19=122(10)=101301(6) 101301-1 Published by the American Physical Society PHYSICAL REVIEW LETTERS 122, 101301 (2019)    Collapse of closed domain walls.—Once the Hubble m2 ðTÞv2 a Vða; TÞ ¼ 1 − cos N DW ; ð3Þ length becomes comparable to the closed wall size R⋆ , N 2DW v the DW rapidly shrinks because of its own tension. This occurs at the temperature T ⋆ defined by R⋆ ∼ H−1 ⋆ ≃ where N DW is the model dependent color anomaly, also −1=2 2 −1=2 geff ðT ⋆ Þ Mp =T ⋆ , where M p ¼ð8πGN Þ and geff ðT ⋆ Þ known as the DW number. The periodicity of V is is the effective number of degrees of freedom at T ⋆ . The 2πF ≡ 2πv=N DW . The dependence of the axion mass total collapsing mass has two contributions: one induced mðTÞ with temperature can be parametrized as by the wall tension σ, and another one coming from any  possible difference in energy density between the two m20 ; if T ≲ T 0 ; 2 m ðTÞ ¼ ð4Þ regions separated by the DW: m20 ðTT0 Þ−n ; if T ≳ T 0 ; 4 4 where n ≈ 7, T 0 ≃ 100 MeV are numerical parameters M⋆ ¼ 4πσR2⋆ þ πΔρR3⋆ ∼ 4πσH−2 −3 ⋆ þ πΔρH ⋆ : ð1Þ 3 3 which we take from Ref. [39] (see Ref. [43] for the original computation [44] for lattice results and the For closed DWs arising in the network Δρ ≥ 0 (see below); thus the wall bounds the region of false vacuum. Supplemental Material [45]). Here, m0 ≃ 0.01Λ2QCD =F is Another important parameter for the formation of PBHs the zero-temperature axion mass, with ΛQCD ≃ 400 MeV. is the ratio of the Schwarzschild radius RS of the collapsing The potential in Eq. (3) leads to the existence of DWs, with wall to the initial size R⋆ : tension σðTÞ ≃ 8mðTÞF2 . These become relevant at the temperature T 1 ∼ GeV defined by 3HðT 1 Þ ¼ mðT 1 Þ (see RS 2GN M⋆ σH−1 ⋆ ΔρH−2 ⋆ also the Supplemental Material [45]). For topological p≡ ∼ ∼ þ : ð2Þ R⋆ H−1 ⋆ M2p 3M2p reasons, each string gets attached to N DW DWs at T 1 . Thus, a string-wall network is formed, which also contains If p is close to 1 then the DW rapidly enters its closed structures. At the same time, the homogeneous Schwarzschild radius and forms a BH. If p ≪ 1, however, component of the axion field starts to oscillate and the wall has to contract significantly before falling inside RS . generates CDM (misalignment mechanism). (iii) Below It is then less likely to form a BH, since asphericities, energy T 1 , the energy density of the network is quickly dominated losses, and/or angular momentum may severely affect the by horizon-size DWs. The subsequent evolution crucially dynamics of the collapse. We will thus refer to p as the figure depends on the DW number (see, e.g., Ref. [13]). If of merit for PBH formation from the collapse of DWs. N DW ¼ 1 the network is unstable and rapidly decays. If The temperature behavior of p and M⋆ is crucial to our instead N DW > 1, the network is stable because strings are proposal. Whenever the tension terms dominate in Eqs. (1) pulled in different directions by the DWs. One thus faces the and (2), we have M⋆ ∼ T −4 −2 ⋆ , p ∼ T ⋆ . If, instead, the energy DW problem [34]. To avoid this catastrophe, a bias term can difference terms dominate, we have M⋆ ∼ T −6 −4 ⋆ , p ∼ T⋆ . be added to the axion potential [35], [39,49] of the form [50] Therefore, the duration of the hybrid network has a huge    a impact on the likelihood of forming PBHs, as well as on V B ðaÞ ¼ A4B 1 − cos þδ : ð5Þ their masses. The use of long-lived string-wall networks v is the essential new idea of our proposal. This requires Notice that the periodicity of Eq. (5) is different from Eq. (3); multiple DWs attached to each string [13]. Interestingly, there is only one global minimum per period 2πv. this can be realized in QCD axion models with domain wall Furthermore, the phase δ represents a generic offset between number larger than one [41]. the bias term and the QCD potential. The addition of Eq. (5) Axion dark matter from string-wall networks.—Let us to Eq. (3) leads to an energy difference between the false embed the basic mechanism illustrated in the previous section in the QCD axion cosmology. Consider a scalar and true minima, Δρ ≃ A4B . This generates pressure, which field Φ with a Uð1ÞPQ symmetry broken at some temper- competes against the wall tension and renders the network unstable [35]. Balance between the two competing effects is ature T PQ after inflation. The field acquires a VEV while its obtained when σ ≃ A4B H−1 , which is confirmed by detailed phase is identified with the QCD axion, i.e., Φ ¼ veiaðt;xÞ=v , numerical analysis of the network evolution in the presence and string defects are formed (see, e.g., Ref. [13]). Below of a bias term [39]. Most of the network disappears at a T PQ , the axion evolution is (i) Most of the energy density in temperature the strings dilutes as ρstrings ∼ μs H 2 , where μs is the string    1=4 tension [42]. In addition, the strings radiate axions [36–39]. Mp Δρ 1=2 90 T2 ¼ ϵ ; ð6Þ Away from the strings, the homogeneous axion field is σ π 2 geff ðT 2 Þ frozen because of Hubble friction. (ii) At T ≲ OðGeVÞ the QCD phase transition occurs. Nonperturbative effects where ϵ ∼ Oð0.1 − 1Þ is a parameter which increases with generate a periodic potential for a N DW and has been numerically determined in Ref. [39]. 101301-2 PHYSICAL REVIEW LETTERS 122, 101301 (2019) Depending on their initial size, most of the closed DWs in the T ⋆ ≃ T 2 . In the most interesting region of parameter space, network will collapse at different temperatures between T 1 we find p ∼ 10−6 , five orders of magnitude larger than and T 2 . In the process, axions are radiated in such a way for T ⋆ ∼ T 1. This shows the advantage of considering that the total axion DM abundance today is given by N DW > 1. Nevertheless, p remains quite small and at this Ωa ¼ Ωmis þ Ωstrings þ Ωnetwork . The axion DM abundance point it is unclear whether this leads to a significant fraction has been numerically studied in Ref. [39] and we review the of PBHs. dependence of Ωa on F and T 2 in the Supplemental Material PBHs from late collapses.—Crucially, for closed DWs [45]. There are two crucial points to take from the discussion collapsing at T ⋆ < T 2 , p increases as T −4 ⋆ , because the above: (i) for N DW > 1 there can be a significant separation vacuum energy contribution dominates over the wall tension, between T 1 , the temperature of network formation, and T 2 , as dictated by Eq. (2). the temperature at which its annihilation is efficient, and The region around F ≲ 109 GeV and T 2 ≃ 7 MeV in (ii) since Eq. (5) lifts the degeneracy of the N DW vacua, closed Fig. 1 leads to the best case scenario for PBH formation. structures surrounding regions with energy A4B can exist in In Fig. 2 we plot the figure of merit and PBH masses for the network. Therefore, from now on we assume N DW > 1. DWs collapsing at T ⋆ < T 2 . In Fig. 1 we plot the constraints on the F − T 2 plane for Figure 2 shows that DWs collapsing roughly when N DW ¼ 2. The blue-shaded region is excluded because of T ⋆ ∼ 0.1T 2 quite likely form PBHs. These structures only DM overproduction. In gray and orange we show the region have to contract by 1 order of magnitude before entering excluded using SN-cooling arguments according to the their Schwarzschild radius. Energy losses via radiation standard analysis [51] and to a more conservative recent of axions as well as the growth of asphericities can be estimate [52], respectively. The thick black lines signal the neglected for such short contractions. Indeed, the radiation largest allowed value of the offset phase δ in Eq. (5) that of energy from a closed spherical sine-Gordon DW was does not spoil the axion solution to the strong CP problem. studied in Ref. [53] and shown to become relevant only We thus conclude that a viable region of parameter once the wall has contracted to a size R ∼ R2=3 ⋆ m −1=3 ≪ space exists, around T 2 ≃ 5 MeV and corresponding to −1 0.1H⋆ . Similarly, in the thin wall approximation aspher- AB ≃ 10−3 ΛQCD , where no tuning of δ is required. This is in icities do not spoil the PBH formation for large p [54]. contrast with the conclusion reached in Ref. [39], where Furthermore, we have numerically simulated the collapse an overconservative bound on θQCD was assumed. The of sine-Gordon nonspherical DWs and checked that they untuned region of parameter space is slightly reduced as can indeed contract down to rmin ≲ 0.1R⋆ [55,56]. N DW increases. Nonetheless, even for N DW ¼ 6 only a mild Let us estimate the fraction f ≡ ΩPBH =ΩCDM . After T 2, tuning δ ∼ 0.1 is required. the energy density of the network is dominated by the In Fig. 1 we also show the relevant would-be BH masses bias contribution. However, at any given T ⋆ < T 2 only a and the figure of merit for closed DWs which collapse at small fraction Pnw of the original network survives. Therefore, FIG. 1. Constraints on F and T 2 from DM overproduction (blue shaded region), and from supernovae cooling (orange shaded). The figure of merit (dashed lines) and the DW masses (red lines) FIG. 2. Constraints on F for T 2 ≃ 7 MeV. The figure of merit are also shown. No tuning of the offset phase is required below (dashed lines) and the masses of the DWs collapsing at T ⋆ < T 2 the line δ ¼ 0.1. (red lines) are also shown. 101301-3 PHYSICAL REVIEW LETTERS 122, 101301 (2019) ρnw ðT ⋆ Þ ∼ Pnw ðT ⋆ ÞΔρ: ð7Þ operators could lead to Eq. (5). These operators were originally investigated in Refs. [61–63], which showed that Assuming that the PBH formation occurs mostly at a single they must come with very small coefficients in order not temperature T ⋆ , the actual fraction f is then given by to spoil the axion solution to the strong CP problem. ρnw ðT ⋆ Þ However, the lore is that gravity affects the PQ symmetry f ∼ pN × ; ð8Þ only at the nonperturbative level [64] (see also ρCDM ðT ⋆ Þ Refs. [65,66]). As a result, the size of the induced potential where N ≳ 1 takes into account the effects of asphericities from gravity can be estimated ∼AB ∼ M p e−#Mp =F , which is and angular momentum. One should keep in mind that f certainly small albeit too small to give a viable cosmology. might be further suppressed by the probability of finding Here, we propose an alternative possibility to generate the closed structures in the network. In Eq. (8), ρCDM ðT ⋆ Þ ∼ bias term. Consider a dark gauge sector, which also breaks ρCDM ðT 2 ÞðT ⋆ =T 2 Þ3 is the energy density of CDM at T ⋆ . In the Uð1ÞPQ via anomalies and has DW number N dark DW ¼ 1. the most interesting region of parameter space in Fig. 1, The specific matter spectrum and couplings of this hidden ρCDM ðT 2 Þ is dominated by the contribution from axions sector are not crucial to our discussion, even though radiated by the network. Hence, ρCDM ðT 2 Þ ≈ ρnw ðT 2 Þ ∼ cosmological and collider constraints should be checked Δρ. Putting everything together, we find in any concrete realization. Such a dark sector would then  3 precisely generate a contribution to the axion potential of T f ∼ pN Pnw ðT ⋆ Þ 2 : ð9Þ the form (5), with AB related to the scale of dark gluon T⋆ condensation, and δ containing the dark sector θ term. To estimate the actual value of f requires knowledge of Interestingly, this naturally allows for the scale AB Pnw . In this respect, the simulations of Ref. [39] show that, to have a T dependence analogous to the QCD axion for N DW ¼ 2, at T 2 defined by Eq. (6) with ϵ ≃ 0.5 only potential 10% of the original network survives. pffiffiffi The percentage is AB ðTÞ4 ¼ m2B ðTÞv2 ¼ m2B ðTÞN 2DW F2 ; ð11Þ further reduced to 1% at T 2 = 2 approximately. We do not know the subsequent evolution of the network. with (see also the Supplemental Material [45]) Nevertheless, let us assume for simplicity that the network   0 2 Λ4B T −n decay follows a power law beyond T 2 [57]: mB ðTÞ ¼ dT 2 ; if T ≳ T 0;B : ð12Þ  −α F ΛB T2 Pnw ðT ⋆ Þ ∼ : ð10Þ The natural expectation is that mB will increase as temper- T⋆ ature decreases until T 0;B ∼ ΛB , and remain constant after- Fitting Eq. (10) to the aforementioned results of Ref. [39] wards. Here, ΛB is the dark confinement scale and dT , n0 gives α ≈ 7. The final fraction f then does not depend depend on the dark spectrum. significantly on N. As long as N > 1, the right-hand side These parameters have an important impact on PBH of Eq. (9) increases as T ⋆ decreases, until p is saturated to formation. For instance, for T 2 ∼ 5 MeV and dT , n0 ∼ 1 one. This occurs at T ⋆ ∼ T 2 =30 (see Fig. 2), which gives the bias term has not yet reached its asymptotic value at T 2 . 0 0 f ∼ 10−6 and M ⋆ ∼ 106 M⊙ . Below this temperature, the Therefore, p and M ⋆ now scale as T −4−n ⋆ and T −6−n ⋆ , fraction falls rapidly. respectively, from T 2 to T ⋆ ∼ T 0;B . A large figure of merit However, this result is sensitive to the precise numerical can then be attained in less than 1 order of magnitude in T, and scaling of Pnw after T 2. In this regard, it is interesting to lighter PBHs may be generated (down to ∼104 M⊙ ) [67]. notice that numerical simulations hint at slight deviations Alternatively, if dT ≪ 1 and/or n0 ≳ 6, ΛB roughly coincides from the scaling regime [39]. The decay of the network can with T 2 and we recover the previous case. In general, the dark then be slower, resulting in smaller α and larger f. sector confinement scale should be 100 keV ≲ ΛB ≲ T 2 in Observations require f ≲ 10−5 for PBHs with M ∼ order for the mechanism presented here to generate an 106 M⊙ [8] (see also Refs. [10,58]). This constraint is interesting and viable fraction of PBHs, as discussed in easily satisfied in our scenario. the previous section. We leave a more detailed investigation Nevertheless, to confidently estimate the actual fraction of the dark sector for future work. requires additional numerical studies, which we leave for Conclusions.—We have discussed a new mechanism future work. Let us remark that if a larger f can be obtained to generate PBHs in the context of QCD axion models. with our mechanism, then (non)observations of PBHs may It proceeds by the late collapses of closed DWs in a long- actually give additional constraints on axion models with lived string-DW network, which arises in QCD axion N DW > 1. realizations with N DW > 1 and PQ symmetry broken after Origin of the bias term.—A minimalistic option to inflation. generate Eq. (5) is from gravity. As pointed out in Lacking accurate knowledge of the network evolution Ref. [59] (see also Ref. [60]), Planck-suppressed effective and collapse, we cannot give precise predictions for the 101301-4 PHYSICAL REVIEW LETTERS 122, 101301 (2019) fraction and masses of the PBHs. However, under [6] M. Sasaki, T. Suyama, T. Tanaka, and S. Yokoyama, Phys. reasonable assumptions, depending on the temperature Rev. Lett. 117, 061101 (2016). behavior of the bias term, PBHs with masses in the range [7] A. Kashlinsky, Astrophys. J. 823, L25 (2016). [8] B. Carr, F. Kuhnel, and M. Sandstad, Phys. Rev. D 94, M ∼ 104 − 107 M ⊙ and representative fraction f ≳ 10−6 083504 (2016). can be created. Interestingly, such heavy PBHs can play an [9] S. Clesse and J. García-Bellido, Phys. Rev. D 92, 023524 important role as seeds for the formation of cosmological (2015). structure, alleviating several problems of the CDM scenario [10] B. Carr and J. Silk, Mon. Not. R. Astron. Soc. 478, 3756 on subgalactic scales, and providing an avenue to explain (2018). the origin of the supermassive BHs [9,10]. [11] X.-B. Wu et al., Nature (London) 518, 512 (2015). Our proposal appears to prefer small values of the axion [12] E. Banados et al., Nature (London) 553, 473 (2018). decay constant, F ≲ 109 GeV, corresponding to axion [13] A. Vilenkin and E. P. S. Shellard, Cosmic Strings and masses in the meV range. These values are close to the Other Topological Defects (Cambridge University Press, lower bounds from the cooling of supernovae [51,52,68], Cambridge, England, 2000). which are however subject to astrophysical uncertainties and [14] A. Vilenkin, Phys. Rev. Lett. 46, 1169 (1981); 46, 1496(E) (1981). are not universal (see, e.g., Ref. [69]). On the other hand, [15] S. W. Hawking, Phys. Lett. B 246, 36 (1990). small values of F might be observationally interesting. In this [16] J. Fort and T. Vachaspati, Phys. Lett. B 311, 41 (1993). respect, it is intriguing that several stellar systems show a [17] J. Garriga and M. Sakellariadou, Phys. Rev. D 48, 2502 mild preference for a nonstandard cooling mechanism, which (1993). can be interpreted in terms of a DFSZ QCD axion [70,71]. [18] M. Yu. Khlopov, S. G. Rubin, and A. S. Sakharov, Astro- In addition, several experiments will be probing this region part. Phys. 23, 265 (2005). of QCD axion masses in the near future. In particular: IAXO [19] H. Deng, J. Garriga, and A. Vilenkin, J. Cosmol. Astropart. [72], TASTE [73], ALPS II [74], and ARIADNE [75]. Phys. 04 (2017) 050. Our mechanism might be probed at gravitational wave [20] T. W. B. Kibble, G. Lazarides, and Q. Shafi, Phys. Rev. D observatories via the detection of gravitational radiation 26, 435 (1982). from: SMBH binaries at LISA [76], the annihilation of [21] A. Vilenkin and A. E. Everett, Phys. Rev. Lett. 48, 1867 the string-wall network [77] at aLIGO (O5) [78], LISA, ET (1982). [22] T. Vachaspati, arXiv:1706.03868. [79], and SKA [80]. [23] A. S. Sakharov and M. Yu. Khlopov, Yad. Fiz. 57, 514 Finally, let us mention that considering very light generic (1994) [Phys. At. Nucl. 57, 485 (1994)]. axionlike particles, the network collapse could be delayed [24] A. S. Sakharov, D. D. Sokoloff, and M. Yu. Khlopov, to T ⋆ ≲ 100 keV raising the figure of merit. The resulting Yad. Fiz. 59, 1050 (1996) [Phys. At. Nucl. 59, 1005 (1996)]. extremely heavy PBHs, however, are strongly constrained. [25] M. Yu. Khlopov, A. S. Sakharov, and D. D. Sokoloff, Nucl. Phys. B, Proc. Suppl. 72, 105 (1999). We thank J. J. Blanco-Pillado, J. Garriga, J. Redondo, [26] J. E. Kim, Phys. Rev. Lett. 43, 103 (1979). K. Saikawa, G. Servant, and T. Vachaspati for useful [27] M. A. Shifman, A. I. Vainshtein, and V. I. Zakharov, Nucl. discussions. We acknowledge support by the Spanish Phys. B166, 493 (1980). Ministry MEC under Grant No. FPA2014-55613-P and [28] M. Dine, W. Fischler, and M. Srednicki, Phys. Lett. 104B, the Severo Ochoa excellence program of MINECO (Grants 199 (1981). No. SO-2012- 0234, No. SEV-2016- 0588), as well as [29] A. R. Zhitnitsky, Yad. Fiz. 31, 497 (1980) [Sov. J. Nucl. by the Generalitat de Catalunya under Grant No. 2014- Phys. 31, 260 (1980)]. SGR-1450. F. F. was also supported in part by the U.S. [30] R. D. Peccei and H. R. Quinn, Phys. Rev. Lett. 38, 1440 Department of Energy, Office of High Energy Physics, (1977). under Awards No. DE-FG02-91ER40628 and No. DE- [31] F. Wilczek, Phys. Rev. Lett. 40, 279 (1978). [32] S. Weinberg, Phys. Rev. Lett. 40, 223 (1978). SC0017987. [33] This situation can arise more generally from sequences of phase transitions in the early Universe. [34] Ya. B. Zeldovich, I. Yu. Kobzarev, and L. B. Okun, Zh. Eksp. Teor. Fiz. 67, 3 (1974) [Sov. Phys. JETP 40, 1 (1974)]. [1] B. P. Abbott et al. (Virgo and LIGO Scientific Collabora- [35] P. Sikivie, Phys. Rev. Lett. 48, 1156 (1982). tions), Phys. Rev. Lett. 116, 061102 (2016). [36] V. B. Klaer and G. D. Moore, J. Cosmol. Astropart. Phys. 11 [2] S. Hawking, Mon. Not. R. Astron. Soc. 152, 75 (1971). (2017) 049. [3] B. J. Carr and S. W. Hawking, Mon. Not. R. Astron. Soc. [37] M. Gorghetto, E. Hardy, and G. Villadoro, J. High Energy 168, 399 (1974). Phys. 07 (2018) 151. [4] S. Bird, I. Cholis, J. B. Muñoz, Y. Ali-Haïmoud, M. [38] M. Kawasaki, T. Sekiguchi, M. Yamaguchi, and J. Yokoyama, Kamionkowski, E. D. Kovetz, A. Raccanelli, and A. G. Prog. Theor. Exp. Phys. (2018), 091E01. Riess, Phys. Rev. Lett. 116, 201301 (2016). [39] M. Kawasaki, K. Saikawa, and T. Sekiguchi, Phys. Rev. D [5] S. Clesse and J. García-Bellido, Phys. Dark Universe 15, 91, 065014 (2015). 142 (2017). [40] A. Ringwald, arXiv:1805.09618. 101301-5 PHYSICAL REVIEW LETTERS 122, 101301 (2019) [41] The original DFSZ [28,29] axion has N DW ¼ 6, while the horizon-size loops scales as H−3 ∼ T −6 , i.e., as a power law simplest KSVZ realization [26,27] has N DW ¼ 1. However, in T. generalizations of the latter with N DW > 1 can be considered. [58] V. Poulin, P. D. Serpico, F. Calore, S. Clesse, and K. Kohri, [42] See, however, Refs. [37,38] for recent claims of small Phys. Rev. D 96, 083524 (2017). logarithmic deviations from such a scaling regime. [59] B. Rai and G. Senjanovic, Phys. Rev. D 49, 2729 (1994). [43] O. Wantz and E. P. S. Shellard, Phys. Rev. D 82, 123508 [60] A. Ringwald and K. Saikawa, Phys. Rev. D 93, 085031 (2010). (2016); 94, 049908(A) (2016). [44] S. Borsanyi et al., Nature (London) 539, 69 (2016). [61] M. Kamionkowski and J. March-Russell, Phys. Lett. B 282, [45] See Supplemental Material at http://link.aps.org/ 137 (1992). supplemental/10.1103/PhysRevLett.122.101301 for a re- [62] R. Holman, S. D. H. Hsu, T. W. Kephart, E. W. Kolb, view on the determination of the axion dark matter mass R. Watkins, and L. M. Widrow, Phys. Lett. B 282, 132 (1992). and relic abundance, which includes Refs. [46–48]. [63] S. M. Barr and D. Seckel, Phys. Rev. D 46, 539 (1992). [46] S. R. Coleman, Subnuclear series 15, 805 (1979). [64] R. Kallosh, A. D. Linde, D. A. Linde, and L. Susskind, [47] S. Borsanyi, M. Dierigl, Z. Fodor, S. D. Katz, S. W. Mages, Phys. Rev. D 52, 912 (1995). D. Nogradi, J. Redondo, A. Ringwald, and K. K. Szabo, [65] G. Dvali, arXiv:hep-th/0507215. Phys. Lett. B 752, 175 (2016). [66] R. Alonso and A. Urbano, arXiv:1706.07415. [48] C. Bonati, M. D’Elia, M. Mariti, G. Martinelli, M. Mesiti, F. [67] The lines of constant δ in Fig. 1 get modified, but viable Negro, F. Sanfilippo, and G. Villadoro, EPJ Web Conf. 137, regions with δ ≳ 0.1 persist. 08004 (2017). [68] K. Hamaguchi, N. Nagata, K. Yanagi, and J. Zheng, Phys. [49] T. Hiramatsu, M. Kawasaki, and K. Saikawa, J. Cosmol. Rev. D 98, 103015 (2018). Astropart. Phys. 08 (2011) 030. [69] L. Di Luzio, F. Mescia, E. Nardi, P. Panci, and R. Ziegler, [50] For the time being, we consider a bias term which switches Phys. Rev. Lett. 120, 261803 (2018). on at T 0 and remains constant thereafter. We discuss more [70] M. Giannotti, I. Irastorza, J. Redondo, and A. Ringwald, about this point later on. J. Cosmol. Astropart. Phys. 05 (2016) 057. [51] M. Tanabashi et al. (Particle Data Group), Phys. Rev. D 98, [71] M. Giannotti, I. G. Irastorza, J. Redondo, A. Ringwald, and 030001 (2018). K. Saikawa, J. Cosmol. Astropart. Phys. 10 (2017) 010. [52] J. H. Chang, R. Essig, and S. D. McDermott, arXiv:1803 [72] E. Armengaud et al., J. Instrum. 9, T05002 (2014). .00993. [73] V. Anastassopoulos et al. (TASTE Collaboration), [53] L. M. Widrow, Phys. Rev. D 40, 1002 (1989). J. Instrum. 12, P11019 (2017). [54] L. M. Widrow, Phys. Rev. D 39, 3576 (1989). [74] R. Bähre et al., J. Instrum. 8, T09001 (2013). [55] F. Ferrer, E. Massó, G. Panico, O. Pujolàs, and F. Rompineve [75] A. Arvanitaki and A. A. Geraci, Phys. Rev. Lett. 113, (to be published). 161801 (2014). [56] We have neglected angular momentum in the numerical [76] A. Klein et al., Phys. Rev. D 93, 024003 (2016). simulation. We expect that for large p the inclusion of the [77] K. Saikawa, Universe 3, 40 (2017). latter should not significantly alter our picture. The resulting [78] J. Aasi et al. (LIGO Scientific Collaboration), Classical PBHs would have masses M ⋆ ∼ 104 − 107 M ⊙ . Quantum Gravity 32, 115012 (2015). [57] In support of this choice, numerical simulations exist [13], [79] M. Punturo et al., Classical Quantum Gravity 27, 194002 which show that at formation the number density of closed (2010). string loops of radius R scales as R−3 . Thus the density of [80] G. Janssen et al., Proc. Sci., AASKA14 (2015) 037. 101301-6