NANOGrav signal as mergers of Stupendously Large Primordial Black Holes
NNANOGrav signal as mergers of Stupendously Large Primordial Black Holes
Vicente Atal, ∗ Albert Sanglas, † and Nikolaos Triantafyllou ‡ Departament de F´ısica Qu`antica i Astrof´ısica, i Institut de Ci`encies del Cosmos,Universitat de Barcelona, Mart´ı i Franqu`es 1, 08028 Barcelona, Spain.
We give an explanation for the signal detected by NANOGrav as the stochastic gravitationalwave background from binary mergers of primordial “Stupendously Large Black Holes”(SLABs) ofmass M ∼ (10 − ) M (cid:12) , and corresponding to roughly 0 .
1% of the dark matter. We show thatthe stringent bounds coming from µ distortions of the CMB can be surpassed if the perturbationsresulting in these BHs arise from the non-Gaussian distribution of fluctuations expected in singlefield models of inflation generating a spike in the power spectrum. While the tail of the stochasticbackground coming from binaries with M (cid:46) M (cid:12) could both fit NANOGrav and respect µ distortions limits, they become excluded from large scale structure constraints. I. INTRODUCTION
The NANOGrav collaboration has recently reportedevidence for a signal consistent with a gravitational wavebackground of frequency ν ∈ (2 . × − , 1 . × − ) Hzand amplitude at 1- σ confidence level Ω GW ∈ (3 × − ,2 × − ). If the signal is modelled as Ω GW ∝ ( ν/ν (cid:63) ) ξ ,the tilt ξ ∈ ( − .
5, 0 .
5) at ν (cid:63) = 5 . σ confidence level [1].The nature of the signal has still to be confirmed withfurther observations and analysis (e.g. whether it is reallya gravitational wave (GW) and whether it is of stochasticorigin or not), but since a potential GW detection mighthave tremendous consequences for our understanding ofthe Universe, it is important to investigate the potentialimplications of such discovery. A suggestion of particularinterest entails the signal being related to a population ofprimordial black holes (PBHs) [2], which might be an im-portant constituent of our Universe (for recent reviews,see e.g [3, 4]). In this line, several possible interpretationsof the NANOGrav signal being the gravitational back-ground induced by large scalar perturbations responsiblefor PBH formation has been proposed [5–12]. In theseworks it has been shown that the signal could be accom-modated by a population of sublunar, solar, or slightlysupersolar PBHs.A quite different possibility, although historically oneof the first to be conceived for explaining a signal at suchfrequencies [13], is that it corresponds to the stochasticbackground resulting from the past mergers of large blackholes, with masses M > M (cid:12) . This possibility wasstudied in [1, 14] where the distributions of the binarieswas inspired from astrophysical models.Here we show that BH-inspirals accounting forNANOGrav could be of primordial origin. For this tohappen, their mass should be ∼ − M (cid:12) , and sothey enter in the class of so-called “Stupendously LargeBlack Holes” (SLABs) [15]. Their abundance with re- ∗ [email protected] † [email protected] ‡ [email protected] spect to dark matter, f ≡ Ω PBH / Ω DM , should be at the O (0 . M > M (cid:12) , is obtained from dynamicalconstraints of large scale structure (LSS) (for a recentreview on this topic, see [15]). II. MERGERS OF SUPERMASSIVE BLACKHOLES AND THEIR STOCHASTICGRAVITATIONAL WAVE SIGNAL
The energy released by BH binaries contribute to astochastic background of gravitational waves, whose en-ergy density Ω GW is given by (we use c = G = 1) [16]Ω GW ≡ ρ c dρ GW d log ν = νρ c H (cid:90) z ∗ R BH ( z )(1 + z ) E ( z ) dE GW dν s ( ν s ) dz (1)where ρ GW is the energy density at a given frequency ν , ρ c is the critical density, d E GW / d ν s is the GW energyspectrum of the merger, R BH is their merger rate at red-shift z and ν s is the frequency in the source frame, relatedto the observed frequency as ν s = (1 + z ) ν . The function E ( z ) ≡ H ( z ) /H = [Ω r (1 + z ) + Ω m (1 + z ) + Ω Λ ] / ,and the GW energy spectrum of the merger is [17] dE GW dν s ( ν s ) = π M c ν − / s ν s < ν ω ν / s ν ≤ ν s < ν ω σ ν s ( σ +4( ν s − ν ) ) ν ≤ ν s < ν ν ≤ ν s (2)where ν i ≡ ( ν , ν , σ, ν ) = ( a i η + b i η + c i ) / ( πM t ), M t = m + m is the total mass of binary system, M c is thechirp mass ( M / c = m m M − / t ), and η = m m M − t is the symmetric mass ratio. The parameters a i , b i and a r X i v : . [ a s t r o - ph . C O ] D ec c i can be found in [17], and ( ω , ω ) are chosen such thatthe spectrum is continuous.The merger rate R BH ( z ) depends on the formationchannel of the binary system. If these are formed fromprimordial fluctuations, then it will further depend onthe their initial distribution, which might be Poissonian(if perturbations are Gaussian [18–20]), or clustered (ifdepartures from Gaussianities are large [21–26]).While we will deal with non-Gaussian primordial fluc-tuations, the effect on their merger rate is going to benegligible since in the minimal model of inflation thatwe will consider here there is no large variance at largescales that can source a spatial modulation in the num-ber density of PBHs. We consider thus a merger rateas coming from a Poissonian distribution of PBHs and amonochromatic mass function ( M ≡ m = m ), that is[27] R BH ( t ) = 8 π ¯ n (cid:90) x max x min x y (cid:12)(cid:12)(cid:12) d y d t ( t, x ) (cid:12)(cid:12)(cid:12) dx , (3)where ¯ n is the comoving number density of BHs, y = x ( x/x min ) , where x min and x max are the minimum andmaximum distances of the binaries such that they mergeat time t [28], x max (cid:39) (cid:18) M t ρ eq (cid:19) and x min = (cid:18) η M t t x (cid:19) x max , (4)with ρ eq the background energy density at matter-radiation equality . In Fig. 1 we show how theNANOGrav signal can be fitted with the stochastic GWsignal given by Eq. (1). It can be attributed to the peakof the signal if PBHs are of masses M ∼ (5 × − ) M (cid:12) and make a fraction f ∼ (2 − × − of thetotal DM. The infrared tail of the stochastic backgroundhas a spectral index ξ = 2 /
3, and is thus within the2- σ interval determined by observations. Thus BHs of M < × M (cid:12) and f > × − could also provide agood fit. However, as we will see in the next section, thesebecome in conflict with LSS constraints. In Fig. 1 weshow the smallest mass that provides a fit to NANOGravthat is not in tension with LSS constraints. The combina-tion of both constraints imply M ∼ (2 × − ) M (cid:12) and f ∼ (2 − × − .While such large BHs have yet to be observed in nature(the largest reported BH has a mass 7 × M (cid:12) [31]), In Eq. (3) we have neglected an exponential factor that is irrel-evant in this case but can be very important when large scalecorrelations are present [26]. Recently BHs described by the so-called Thakurta metric hasbeen studied [29]. In that description, and due to a constantenergy flow towards the BHs, their masses are time dependent,which alters their merger rate. Since in the cosmological setupthat we are interested there is no fluid sustaining this mass grow,these BHs are not the ones that we are interested in (see also[30]). − − − − ν [Hz]10 − − − Ω G W f = × − f = × − M = 1012 M (cid:12) M = 5 × M (cid:12) M = 2 × M (cid:12) NanoGrav
FIG. 1. The stochastic background of binary mergers. In greywe show the 1- σ region consistent with NANOGrav signal. the presence of more massive BHs is not ruled out . Inthe following we discuss limits on the fraction of PBHs ofthis range of masses, most notably the spectral distortionand large scale structure bounds. III. SPECTRAL DISTORTIONS,NON-GAUSSIANITIES AND LSS
Spectral distortions provide the most stringent boundfor the presence of PBHs in the mass range M =(10 − ) M (cid:12) . The large fluctuations necessary forthe production of PBHs dissipate through Silk dampingand might leave a large imprint in the CMB as depar-tures from a black-body spectrum . In particular, theseprimordial inhomogeneities are constrained by a non-detection of µ distortions, which from COBE/FIRAS isbounded to be smaller than 9 × − [35].The amplitude of the µ distortions depends on boththe scale and the variance of the power spectrum of thescalar perturbations [36]. Both of these are related tothe abundance of PBH, and so constraints on µ distor-tions can be used to constraint the presence of PBHs.From here it has been deduced that PBHs in the massrange (6 × − × ) M (cid:12) are excluded [37] (for earlierapplication of this idea, see [38, 39]).These constraints can be relaxed if NG are considered[40, 41]. The abundance of PBH is most sensitive to theratio ν ≡ ζ c /σ , and NG might change both σ , the vari-ance of the perturbations, and ζ c , the critical thresholdfor collapse. In general the role that NG plays in de-termining the threshold for collapse has been neglected, For for astrophysical considerations on their presence, see [32]. A natural exception are models in which the PBHs are not asso-ciated to an enhancement of the power spectrum (like [33, 34]).These models are thus essentially unconstrained by spectral dis-tortions measurements. and only the effect coming from changes in the PDF hasbeen considered [12, 40, 41]. However, the threshold forcollapse depends on the shape of the overdensity [42–48],and non-Gaussianities modify its shape [49–52].For example, if we consider local models of NG, forwhich the curvature perturbations ζ are related to aGaussian variable ζ G with a local function ζ = F ( ζ G ),the profile of an overdensity is simply F ( ζ G ( r )), where ζ G ( r ) is its shape as given by a Gaussian random field[53]. Then the treshold for collapse can be simply deter-mined numerically (with public codes [47]) and/or ana-lytically [48].Note that since local models of NG can be written interms of an underlying Gaussian field, it is not necessaryto find how the PDF changes by the local transformationfor computing the abundances. It is actually sufficient tocount the regions for which the image of the underlyingGaussian field is above the threshold of collapse. That is,if ζ c N G is the critical value for collapse of the NG profile,we can define its counterpart in the underlying Gaussianfield, µ c , given by µ c ≡ F − ( ζ c N G ). Then, if the PBHabundance for a Gaussian field is β G ∼ e − ( ζ c G ) /σ G with ζ cG the treshold for the Gaussian field, then the abun-dance of PBHs in the local NG theory is simply given by β NG ∼ e − µ c /σ G .It is thus fortunate that in the case of single-field infla-tionary models producing a spike in the power spectrum(necessary for a controled production of PBHs), the NGsare of the local type. Its precise shape and amplitudewere established in [49], and are given by ζ = − µ ∗ log (cid:18) − ζ g µ ∗ (cid:19) (5)where µ ∗ is related to the curvature of the local maximathat the inflaton field traverses, as [49, 54] µ ∗ = − (cid:112) − η V | max (6)where η V ≡ V (cid:48)(cid:48) /V , with V ( φ ) the inflaton potential, andwhere (cid:48) denotes derivatives with respect to inflaton field φ . The parameter η V in Eq. (6) is to be evaluated at thelocal maxima of the potential. As the field overcomesthe local maxima it enters a stage of constant-roll, per-turbations are largely amplified and the abundance ofPBHs increases. The local transformation Eq. (5) is anon-perturative completion of the well known perturba-tive version F ( ζ G ) = ζ G + (3 / f NL ζ [52]. Thus, we canrelate the parameter µ ∗ to f NL , as µ ∗ = 56 f NL . (7)The local transformation Eq. (5) induces a change in theshape of the overdensity field with respect to a Gaussianfield. In [52] it was determined how the threshold for col-lapse changes as the parameter µ ∗ (or equivalently f NL )varies. For a peaked power spectrum growing as k , as expected in single-field inflationary models [55], we find µ c = (cid:40) . − . f NL + 0 . f f NL < / (6 f NL ) f NL > F ( ζ G ) = ζ G + (3 / f NL ζ [52]. Also, since thethreshold is found here in terms of ζ , we avoid com-plications related to additional NGs from the non-linearrelation between curvature and density perturbation.For determining the spectral distortions it is necessaryto calculate the variance of the non-Gaussian field . Thisis given by σ = 1 (cid:112) πσ (cid:90) µ ∗ −∞ log (cid:18) − ζ g µ ∗ (cid:19) e − ζ g σ dζ g − (cid:104) ζ (cid:105) . (9)The amplitude of the µ distortion is related to the am-plitude and scale of a perturbation as [36] µ (cid:39) . σ (cid:34) exp (cid:18) − k (cid:63) (cid:19) − exp (cid:32) − (cid:18) k (cid:63) . (cid:19) (cid:33)(cid:35) . (10)Here the wavenumber k (cid:63) is in M − . Assuming amonochromatic mass function -a good approximation forpeaked power spectra- k (cid:63) is related to the mass of thePBH as [56] k (cid:63) (cid:39) . × γ (cid:16) g . (cid:17) − (cid:18) M M (cid:12) (cid:19) − . (11)Here γ is the ratio between the mass of the BH and themass enclosed within the horizon H = ak − , and g isthe number of relativistic species (we fix g = 10 . γ = 1, knowing that we expect somedepartures from unity from the fact that the BH collapseis a critical phenomena [57], and because the size of theoverdensity undergoing the collapse is typically some fac-tors larger than ak − [44, 45].The fraction of PBH to DM, f , is related to their pri-mordial abundance as [56] β (cid:39) − γ − (cid:16) g . (cid:17) (cid:18) Ω DM . (cid:19) − (cid:18) M M (cid:12) (cid:19) f . (12) Actually there are no large differences at the level of the thresholdif we had considered a power spectrum given by a δ function. Fordetails, see [52]. This turns out to be very close to the Gaussian variance, sincethe logarithm in Eq. (5) only affects profiles close to µ ∗ , whichare very rare with respect to the typical profiles of ζ . M PBH /M (cid:12) − − − − − f f NL = 0 f NL = 3 f NL = 10LSSNanoGrav FIG. 2. Limits on the abundance of PBHs coming fromthe non observations of µ distortions in the CMB and fromlarge scale structure considerations. We can see that for f NL (cid:38)
3, PBHs with the mass and abundance appropiateto fit NANOGrav are not constrained by µ distortions. How-ever, for these abundances, the formation of clusters imposes M > × M (cid:12) . The region in black corresponds to the fitsshown in Fig. 1 . The abundance can be calculated using the standardpeak theory [53], since we mention previously that wecan always refer to the underlying Gaussian field . Forthe case of a monochromatic power spectrum, it is simplygiven by [44, 45] β (cid:39) µ c σ exp (cid:18) − µ c σ (cid:19) . (13)Here µ c is the threshold of the underlying Gaussian field,and σ its variance. In fig. 2 we show the constrains fromthe FIRAS limit on µ , and how they vary with increas-ing non-Gaussianity parameter f NL . For f NL > f NL > k , we should also worry about theconstraints at scales larger than the peak of the powerspectrum, since those are also constrained by the boundson µ distortion. In this respect, we should note that byvirtue of a duality between the background before andafter the peak in the power spectrum, the amplitude ofthe non-Gaussianity remains constant during the transi-tion [54, 59]. Thus, the most critical test concerning µ distortions is at the scale of the peak, k (cid:63) , irrespective ofthe steepness of the power spectrum. Note also that insimilar inflationary scenarios than the one we describe It has been shown that computing the abundances using the com-paction function gives a similar result for peaked power spectra[58]. here there might be an additional source of µ distor-tions, coming from shock waves of an expanding bubbleof false vacuum in the thermalized fluid [60]. In our case,when the BH corresponds to regions in the false vacuum( f NL > µ distortions are not the only ones con-straining PBHs of these masses. The effects of gas ac-cretion in the thermal history of the Universe might alsobe important. For this range of masses, their effect wasstudied in [61]. In the simplified model considered there,for an efficiency (cid:15) < . f ∼ − is allowed. Wenote that these constraints should be updated consider-ing more recent studies, that however mostly focused in M < M (cid:12) (see e.g. [62, 63]).We can also consider dynamical constraints. MassivePBHs could destroy galaxies in clusters [64] and/or trig-ger an undesirably early formation of cosmic structure[65]. In Fig. 2 we show the limits coming from the for-mation of dwarf galaxies up to clusters of galaxies. Forexample, in the mass range M ∼ (3 × − ) M (cid:12) ,the strongest limit comes from the formation of clustersof galaxies, and reads [65] f < M M (cid:12) for 3 × M (cid:12) < M < M (cid:12) . (14)These constraints are fully consistent with fittingNANOGrav with the peak of the stochastic background,but put limits on the fits coming from its tail, fix-ing M > × M (cid:12) . All in all, these constraintsimply that the NANOGrav signal can be explained if M ∼ (2 × − ) M (cid:12) , f ∼ (2 − × − and f NL >
3. Note that for these parameters the signal isno longer well described by a power law in the relevantrange of frequencies, and thus motivates an extension ofthe templates used in [1, 14] to asses the significance ofthe fits.
IV. CONCLUSIONS
We have provided an explanation of the NANOGravsignal as coming from the binary mergers of Stupen-dously Large Black Holes, of masses in the range M ∼ (2 × − ) M (cid:12) . For a mild non-Gaussianity re-sulting from the dynamics of single field inflation, thespectral distortions constraints can be avoided. For theabundance needed for explaining NANOGrav, f ∼ − ,the dynamical constraints fix their mass to be M > × M (cid:12) .These results provide yet another motivation for study-ing the presence of SLABs, and illustrate the importanceof non-Gaussianities and their connection to the forma-tion mechanism of PBHs for determining bounds on theirpresence.Finally, we should indicate that in this work we haveneglected the mass evolution of the PBHs. Since mostof the signal from the stochastic background comes fromlate mergers, the mass evolution could be important (seee.g. [66]). Then it would be necessary to invoke PBHsof smaller masses for having a late time signal at the fre-quencies of NANOGrav. While large values of f NL caneasily be obtained to surpass the µ distortion limits [49],bounds on LSS are more difficult to overcome. Consid-ering that at the moment bounds from LSS and gas ac-cretion are only order of magnitude constraints for suchmasses, a deeper examination of them should be pursuedfor a more robust model selection. ACKNOWLEDGMENTS
We thank Jaume Garriga for comments on themanuscript. This work has been partially sup- ported by FPA2016-76005-C2-2-P, MDM-2014-0369 ofICCUB (Unidad de Excelencia Maria de Maeztu), andAGAUR2017SGR-754 (V.A), by an APIF grant fromUniversitat de Barcelona (A.S), and by an INPhINITgrant from laCaixa Foundation (ID 100010434) codeLCF/BQ/IN17/11620034 (N.T). This project has also re-ceived funding from the European Unions Horizon 2020research and innovation programme under the MarieSklodowska-Curie grant agreement No. 713673. [1] Z. Arzoumanian et al. [NANOGrav], [arXiv:2009.04496[astro-ph.HE]].[2] S. Hawking, Mon. Not. Roy. Astron. Soc. (1971), 75[3] M. Sasaki, T. Suyama, T. Tanaka and S. Yokoyama,Class. Quant. Grav. (2018) no.6, 063001doi:10.1088/1361-6382/aaa7b4 [arXiv:1801.05235 [astro-ph.CO]].[4] A. M. Green and B. J. Kavanagh, [arXiv:2007.10722[astro-ph.CO]].[5] V. Vaskonen and H. Veerm¨ae, [arXiv:2009.07832 [astro-ph.CO]].[6] V. De Luca, G. Franciolini and A. Riotto,[arXiv:2009.08268 [astro-ph.CO]].[7] K. Kohri and T. Terada, [arXiv:2009.11853 [astro-ph.CO]].[8] L. Bian, J. Liu and R. Zhou, [arXiv:2009.13893 [astro-ph.CO]].[9] S. Sugiyama, V. Takhistov, E. Vitagliano, A. Kusenko,M. Sasaki and M. Takada, [arXiv:2010.02189 [astro-ph.CO]].[10] G. Dom`enech and S. Pi, [arXiv:2010.03976 [astro-ph.CO]].[11] S. Bhattacharya, S. Mohanty and P. Parashari,[arXiv:2010.05071 [astro-ph.CO]].[12] K. Inomata, M. Kawasaki, K. Mukaida andT. T. Yanagida, [arXiv:2011.01270 [astro-ph.CO]].[13] M. Rajagopal and R. W. Romani, Astrophys. J. (1995), 543-549 doi:10.1086/175813 [arXiv:astro-ph/9412038 [astro-ph]].[14] H. Middleton, A. Sesana, S. Chen, A. Vecchio, W. DelPozzo and P. A. Rosado, [arXiv:2011.01246 [astro-ph.HE]].[15] B. Carr, F. Kuhnel and L. Visinelli,doi:10.1093/mnras/staa3651 [arXiv:2008.08077 [astro-ph.CO]].[16] T. Regimbau, Res. Astron. Astrophys. (2011), 369-390 doi:10.1088/1674-4527/11/4/001 [arXiv:1101.2762[astro-ph.CO]].[17] P. Ajith et al. , Phys. Rev. D (2008) 104017 Erratum: [Phys. Rev. D (2009) 129901]doi:10.1103/PhysRevD.79.129901, 10.1103/Phys-RevD.77.104017 [arXiv:0710.2335 [gr-qc]].[18] Y. Ali-Ha¨ımoud, Phys. Rev. Lett. (2018)no.8, 081304 doi:10.1103/PhysRevLett.121.081304[arXiv:1805.05912 [astro-ph.CO]].[19] V. Desjacques and A. Riotto, Phys. Rev. D (2018) no.12, 123533 doi:10.1103/PhysRevD.98.123533[arXiv:1806.10414 [astro-ph.CO]].[20] G. Ballesteros, P. D. Serpico and M. Taoso, JCAP (2018) 043 doi:10.1088/1475-7516/2018/10/043[arXiv:1807.02084 [astro-ph.CO]].[21] Y. Tada and S. Yokoyama, Phys. Rev. D (2015) no.12, 123534 doi:10.1103/PhysRevD.91.123534[arXiv:1502.01124 [astro-ph.CO]].[22] S. Young and C. T. Byrnes, JCAP (2015) 034doi:10.1088/1475-7516/2015/04/034 [arXiv:1503.01505[astro-ph.CO]].[23] K. M. Belotsky et al. , Eur. Phys. J. C (2019) no.3, 246 doi:10.1140/epjc/s10052-019-6741-4[arXiv:1807.06590 [astro-ph.CO]].[24] T. Suyama and S. Yokoyama, PTEP (2019) no.10,103E02 doi:10.1093/ptep/ptz105 [arXiv:1906.04958[astro-ph.CO]].[25] S. Young and C. T. Byrnes, JCAP (2020)no.03, 004 doi:10.1088/1475-7516/2020/03/004[arXiv:1910.06077 [astro-ph.CO]].[26] V. Atal, A. Sanglas and N. Triantafyllou, JCAP (2020), 036 doi:10.1088/1475-7516/2020/11/036[arXiv:2007.07212 [astro-ph.CO]].[27] T. Nakamura, M. Sasaki, T. Tanaka and K. S. Thorne,Astrophys. J. Lett. (1997), L139-L142doi:10.1086/310886 [arXiv:astro-ph/9708060 [astro-ph]].[28] P. C. Peters, Phys. Rev. (1964), B1224-B1232doi:10.1103/PhysRev.136.B1224[29] C. Boehm, A. Kobakhidze, C. A. J. O’Hare,Z. S. C. Picker and M. Sakellariadou, [arXiv:2008.10743[astro-ph.CO]]. [30] V. De Luca, V. Desjacques, G. Franciolini andA. Riotto, JCAP (2020), 028 doi:10.1088/1475-7516/2020/11/028 [arXiv:2009.04731 [astro-ph.CO]].[31] O. Shemmer, H. Netzer, R. Maiolino, E. Oliva,S. Croom, E. Corbett and L. di Fabrizio, Astrophys.J. (2004), 547-557 doi:10.1086/423607 [arXiv:astro-ph/0406559 [astro-ph]].[32] P. Natarajan and E. Treister, Mon. Not. Roy.Astron. Soc. (2009), 838 doi:10.1111/j.1365-2966.2008.13864.x [arXiv:0808.2813 [astro-ph]].[33] J. Garriga, A. Vilenkin and J. Zhang, JCAP (2016), 064 doi:10.1088/1475-7516/2016/02/064[arXiv:1512.01819 [hep-th]].[34] H. Deng, J. Garriga and A. Vilenkin, JCAP (2017), 050 doi:10.1088/1475-7516/2017/04/050[arXiv:1612.03753 [gr-qc]].[35] D. J. Fixsen, E. S. Cheng, J. M. Gales, J. C. Mather,R. A. Shafer and E. L. Wright, Astrophys. J. (1996),576 doi:10.1086/178173 [arXiv:astro-ph/9605054 [astro-ph]].[36] J. Chluba, A. L. Erickcek and I. Ben-Dayan, Astro-phys. J. (2012), 76 doi:10.1088/0004-637X/758/2/76[arXiv:1203.2681 [astro-ph.CO]].[37] K. Kohri, T. Nakama and T. Suyama, Phys. Rev. D (2014) no.8, 083514 doi:10.1103/PhysRevD.90.083514[arXiv:1405.5999 [astro-ph.CO]].[38] B. J. Carr and J. E. Lidsey, Phys. Rev. D (1993),543-553 doi:10.1103/PhysRevD.48.543[39] B. J. Carr, J. H. Gilbert and J. E. Lidsey, Phys. Rev.D (1994), 4853-4867 doi:10.1103/PhysRevD.50.4853[arXiv:astro-ph/9405027 [astro-ph]].[40] T. Nakama, B. Carr and J. Silk, Phys. Rev. D (2018) no.4, 043525 doi:10.1103/PhysRevD.97.043525[arXiv:1710.06945 [astro-ph.CO]].[41] C. Unal, E. D. Kovetz and S. P. Patil, [arXiv:2008.11184[astro-ph.CO]].[42] M. Shibata and M. Sasaki, Phys. Rev. D (1999) 084002 doi:10.1103/PhysRevD.60.084002 [gr-qc/9905064].[43] T. Harada, C. M. Yoo, T. Nakama andY. Koga, Phys. Rev. D (2015) no.8, 084057doi:10.1103/PhysRevD.91.084057 [arXiv:1503.03934[gr-qc]].[44] C. M. Yoo, T. Harada, J. Garriga and K. Kohri, PTEP (2018) no.12, 123E01 doi:10.1093/ptep/pty120[arXiv:1805.03946 [astro-ph.CO]].[45] C. Germani and I. Musco, Phys. Rev. Lett. (2019)no.14, 141302 doi:10.1103/PhysRevLett.122.141302[arXiv:1805.04087 [astro-ph.CO]].[46] I. Musco, Phys. Rev. D (2019) no.12, 123524doi:10.1103/PhysRevD.100.123524 [arXiv:1809.02127[gr-qc]].[47] A. Escriv`a, Phys. Dark Univ. (2020), 100466doi:10.1016/j.dark.2020.100466 [arXiv:1907.13065 [gr-qc]].[48] A. Escriv`a, C. Germani and R. K. Sheth, Phys. Rev. D (2020) no.4, 044022doi:10.1103/PhysRevD.101.044022 [arXiv:1907.13311[gr-qc]].[49] V. Atal, J. Garriga and A. Marcos-Caballero, JCAP (2019), 073 doi:10.1088/1475-7516/2019/09/073[arXiv:1905.13202 [astro-ph.CO]].[50] C. M. Yoo, J. O. Gong and S. Yokoyama, JCAP (2019), 033 doi:10.1088/1475-7516/2019/09/033[arXiv:1906.06790 [astro-ph.CO]].[51] A. Kehagias, I. Musco and A. Riotto, JCAP (2019), 029 doi:10.1088/1475-7516/2019/12/029[arXiv:1906.07135 [astro-ph.CO]].[52] V. Atal, J. Cid, A. Escriv`a and J. Garriga, JCAP (2020), 022 doi:10.1088/1475-7516/2020/05/022[arXiv:1908.11357 [astro-ph.CO]].[53] J. M. Bardeen, J. R. Bond, N. Kaiser and A. S. Szalay,Astrophys. J. (1986) 15. doi:10.1086/164143[54] V. Atal and C. Germani, Phys. Dark Univ. (2019), 100275 doi:10.1016/j.dark.2019.100275[arXiv:1811.07857 [astro-ph.CO]].[55] C. T. Byrnes, P. S. Cole and S. P. Patil, JCAP (2019), 028 doi:10.1088/1475-7516/2019/06/028[arXiv:1811.11158 [astro-ph.CO]].[56] T. Nakama, J. Silk and M. Kamionkowski, Phys. Rev. D (2017) no.4, 043511 doi:10.1103/PhysRevD.95.043511[arXiv:1612.06264 [astro-ph.CO]].[57] J. C. Niemeyer and K. Jedamzik, Phys. Rev. D (1999),124013 doi:10.1103/PhysRevD.59.124013 [arXiv:astro-ph/9901292 [astro-ph]].[58] C. Germani and R. K. Sheth, Phys. Rev. D (2020) no.6, 063520 doi:10.1103/PhysRevD.101.063520[arXiv:1912.07072 [astro-ph.CO]].[59] W. H. Kinney, Phys. Rev. D (2005), 023515doi:10.1103/PhysRevD.72.023515 [arXiv:gr-qc/0503017[gr-qc]].[60] H. Deng, A. Vilenkin and M. Yamada, JCAP (2018), 059 doi:10.1088/1475-7516/2018/07/059[arXiv:1804.10059 [gr-qc]].[61] B. Carr, Mon. Not. Roy. Astron. Soc. (1981), 639doi:10.1093/mnras/194.3.639[62] M. Ricotti, J. P. Ostriker and K. J. Mack, Astrophys.J. (2008), 829 doi:10.1086/587831 [arXiv:0709.0524[astro-ph]].[63] Y. Ali-Ha¨ımoud and M. Kamionkowski, Phys. Rev. D (2017) no.4, 043534 doi:10.1103/PhysRevD.95.043534[arXiv:1612.05644 [astro-ph.CO]].[64] B. J. Carr and M. Sakellariadou, Astrophys. J. (1999), 195-220 doi:10.1086/307071[65] B. Carr and J. Silk, Mon. Not. Roy. Astron. Soc. (2018) no.3, 3756-3775 doi:10.1093/mnras/sty1204[arXiv:1801.00672 [astro-ph.CO]].[66] J. Garcia-Bellido, M. Peloso and C. Unal, JCAP09