Scrutinizing the Primordial Black Holes Interpretation of PTA Gravitational Waves and JWST Early Galaxies

Yann Gouttenoire [email protected] School of Physics and Astronomy, Tel-Aviv University, Tel-Aviv 69978 Sokratis Trifinopoulos [email protected] Center for Theoretical Physics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA Georgios Valogiannis [email protected] Department of Physics, Harvard University, Cambridge, MA, 02138, USA
Miguel Vanvlasselaer [email protected] Theoretische Natuurkunde and IIHE/ELEM, Vrije Universiteit Brussel, & The International Solvay Institutes, Pleinlaan 2, B-1050 Brussels, Belgium

Abstract

Recent observations have granted to us two unique insights into the early universe: the presence of a low-frequency stochastic gravitational wave background detected by the NANOGrav and Pulsar Timing Array (PTA) experiments and the emergence of unusually massive galaxy candidates at high redshifts reported by the James Webb Space Telescope (JWST). In this letter, we consider the possibility that both observations have a common origin, namely primordial black holes (PBHs) in the mass range between $10^{6}~{}M_{\odot}$ and $10^{13}~{}M_{\odot}$ . While superheavy PBHs act as seeds for accelerated galaxy formation capable of explaining the JWST extreme galaxies, they can also form binary mergers that source gravitational waves which can be potentially identified as the PTA signal. The analysis is performed taking into account the constraints on the relevant region of the PBH parameter space including the novel bound imposed by the Ultraviolet Luminosity Function of galaxies observed by the Hubble Space Telescope. We conclude that PTA’s and JWST’s interpretations in terms of PBH binary mergers and Poissonian gas of PBHs, respectively, are strongly excluded.

I Introduction

The NANOGrav collaboration [1, 2] combined with the other Pulsar Timing Array (PTA) Collaborations , EPTA [3, 4], PPTA [5, 6], CPTA [7] and IPTA [8] have recently released further evidence for the Hellings-Downs angular correlation in the common-spectrum process. This points toward the existence of a Gravitational Wave (GW) background in the nHz range permeating the universe, sparking a lot of new physics speculation [9, 10].

A different lens through which we can gain glimpses into the early cosmic evolution is presented by the James Webb Space Telescope (JWST). Initial observations have reported photometric evidence of massive galaxies at unexpectedly high redshifts $7<z<12$ [11, 12, 13, 14]. A subset of them has been recently spectroscopically confirmed [15, 16] and large cosmological hydrodynamical simulation demonstrated compatibility with existing models of galaxy formation [17, 18]. However, the status of extreme galaxy candidates with stellar mass as high as $10^{11}~{}M_{\odot}$ [14] still remains under investigation and if their distance estimates prove accurate, they would pose a significant challenge to $\Lambda$ CDM itself [19, 20, 21]. Also in this case, explanations beyond the standard cosmological paradigm have been postulated [22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32].

Primordial black holes (PBHs) emerge as one of the most long-studied scenarios [33, 34, 35, 36, 37, 38, 39], capable of leaving distinctive imprints on cosmic history. Depending on the fraction of their abundance with respect to the total dark matter (DM) $f_{\rm PBH}=\Omega_{\rm PBH}/\Omega_{\rm DM}$ , the spectrum of possible PBH masses $M_{\rm PBH}$ spans an enormous range which has been tested by various experiments (for comprehensive reviews see refs. [40, 41]). PBHs are expected to form binaries [42, 43, 44, 45, 46, 47, 48, 49, 40], which once decoupled from the universe evolution, continuously emit GWs until a last spectacular burst when they finally merge. Mergers of stellar mass BHs have been observed by terrestrial interferometers [50, 51], with the potential of a primordial origin not ruled out [45, 52, 53, 54, 55, 56].

PBHs heavier than $100~{}M_{\odot}$ are of particular interest due to their influence on the growth and formation of structures. For example, it is a well-known fact that supermassive black holes (SMBHs) with masses above $10^{5}~{}M_{\odot}$ reside within galactic nuclei [57, 58, 59]. Naturally, it has been proposed that PBHs can be their progenitors reaching these masses either by merging and accretion [60, 61, 47, 62, 63, 64] or direct collapse of primordial fluctuations [65, 66]. In the latter case, supermassive PBHs are constrained to be less than $\mathcal{O}(0.1\%)$ of DM and since they would already be present from the dawn of the matter-domination, they can act as cosmic seeds boosting galaxy formation [67, 68]. Moreover, a subdominant component of DM can consist of stupendously large PBHs, which are heavier than $10^{12}~{}M_{\odot}$ [69] and may traverse the intergalactic medium.

In this letter, we focus on supermassive PBHs in the mass range $10^{6}<M_{\rm PBH}/M_{\odot}<10^{13}$ and assess the viability of the following two scenarios: i) they bind in binary systems which leads to late-time merging and radiation of GWs that are detectable by PTA experiments [70, 71], ii) they are responsible for the generation of the early massive galaxies reported by JWST [22, 24]. For the first case, we calculate the GW energy spectrum induced by binary mergers and then perform a Bayesian analysis to determine the posterior distribution compatible with the NANOGrav 15-year [1] and IPTA DR2 [8] signals. For the second case, the accelerated galaxy formation in the presence of PBHs is investigated using both the Poisson and seed effects [67] and we identify the PBH populations that can sufficiently seed a large number of massive galaxies at $z\sim 10$ .

The results of our analysis are presented combined with all relevant cosmological and astrophysical constraints. Notably, we consider for the first time in the PBH literature a constraint inferred by measurements of the Hubble Space Telescope (HST) and encoded in the so-called UV galaxy luminosity function (UV LF) [72, 73]. Observations of UV-bright galaxies collected by the Hubble Space Telescope (HST) [74, 75, 76, 77, 78, 79, 80, 72, 81, 82, 73] trace the universe at redshifts $z=4-10$ , significantly overlapping with the ones performed by the JWST. In the context of our combined analysis, we find that the HST UV LF together with the large-structure constraints (LSS) rule out the PTA’s interpretation in terms of GW from PBH mergers as well as the JWST’s interpretation in terms of Poissonian density fluctuation sourced by PBHs.

II Fitting the PTA signals with gravitational waves from PBH binary mergers

Supermassive PBHs can form in the early universe at redshift $z_{\rm f}$ with a mass $M_{\rm PBH}\simeq\big{(}5\times 10^{12}/(1+z_{\rm f})\big{)}^{2}M_{\odot}$ [40]. Immediately after their formation, PBHs are sparsely distributed in space, with mean separation much larger than the Hubble scale, $l_{\rm PBH}(z_{\rm f})\gg H^{-1}(z_{\rm f})$ . However the mean separation between PBHs, $l_{\rm PBH}(t)\propto t^{1/2}$ , falls below the Hubble distance $H^{-1}\propto t$ before matter-radiation equality and several PBHs can coexist in the same Hubble patch, leading to the possibility of creating pairs of PBHs. In principle, PBHs binaries can form both before or after matter-radiation equality (early [83] and late-time [42, 43, 44] formation), however the early contribution dominates for black holes of primordial origins, e.g. [49, 40]. We present the calculation of the distribution of binaries and of the associated merging rate $\mathcal{R}(z)$ , in Appendix C.¹¹1For the possibility of initial clustering of PBHs and the implication for the PTA signal see ref. [71]. After the creation of the binary system, the two PBHs will orbit around each other and merge approximately within a time [84]

\displaystyle t_{\rm m}

\displaystyle\simeq\frac{3}{170}\frac{1}{G^{3}M^{3}_{\rm PBH}}a^{4}(1-e^{2})^{7/2}~{},

(1)

which depends on the eccentricity $e$ of the orbit of the two PBHs and its major axis $a$ . $G$ is the gravitation constant. Circular binaries with abundance $f_{\rm PBH}\simeq(10^{13}M_{\odot}/M_{\rm PBH})^{5/16}$ will typically merge today $t_{\rm m}\simeq 13.8$ Gyr.

Gravitational waves from PBH binary mergers.

The energy density of the stochastic gravitational waves from PBH binaries reads [85, 86]

h^{2}\Omega_{\rm GW}=\frac{f}{\rho_{c}/h^{2}}\int_{0}^{\infty}dz\frac{\mathcal{R}(z)}{(1+z)H(z)}\frac{dE_{\rm GW}(f^{\prime})}{df^{\prime}}\bigg{|}_{f^{\prime}=(1+z)f}~{},

(2)

where $\rho_{c}/h^{2}=(3.0~{}\rm meV)^{4}$ is the Hubble-rescaled critical density at present time, $dE_{\rm GW}(f^{\prime})/df^{\prime}$ is the GW power emitted by a circular and GW-driven binary, $f^{\prime}$ is the rest frame frequency and $f=f^{\prime}/(1+z)$ the observed one. We find a fit for the GW signal with the following power broken-law

\displaystyle\Omega_{\rm GW}h^{2}\simeq\,\Omega_{\rm peak}S(f)\Theta(2f_{\rm peak}-f)~{},

(3)

where the spectral function is

\displaystyle S(f)=\,\frac{f_{\rm peak}^{b}f^{a}}{\bigg{(}bf^{\frac{a+b}{c}}+af_{\rm peak}^{\frac{a+b}{c}}\bigg{)}^{c}}~{},

(4)

the peak amplitude is

\displaystyle\Omega_{\rm peak}\simeq\,0.05f_{\rm PBH}^{3}\bigg{(}1.5\frac{M_{\rm PBH}}{10^{12}M_{\odot}}\bigg{)}^{-0.3}~{},

(5)

and the peak frequency is $f_{\rm peak}\simeq\,5000M_{\odot}/M_{\rm PBH}$ . The fitting parameters are $a=0.7,\quad b=1.5,\quad c=0.9$ . We refer the reader to Appendix D for more details. The step function $\Theta$ is introduced to cut the GW spectrum at frequencies above the ring-down frequency $f>2f_{\rm peak}$ . On top of this estimate, we consider the following two effects:

Environmental effects.

At very low frequencies (IR), the assumption of GW-driven energy loss breaks down and interactions with the environment generate additional energy losses which change the spectral shape of the GW background [87, 88]. Following [89], we discard the region of frequencies lower than

\displaystyle f_{\rm min}=\bigg{(}\frac{T_{\rm max}}{\delta_{2}}\bigg{)}^{\!-\tfrac{3}{8}}~{}~{},\qquad\delta_{2}=\frac{5}{256\pi^{\!\frac{8}{3}}(GM_{c})^{\!\frac{5}{3}}}~{},

(6)

where the choice $T_{\rm max}=75$ Myr corresponds to the transition from the stellar-scattering dominated phase to the GW dominated one [89].

Discrete versus continuous signal.

Eq. (2) assumes a smooth distribution of sources. However, above a frequency $f_{\rm max}$ , the number of sources per frequency bins $N(f,\Delta f)$ (we explain in Appendix D how to compute the function $N$ ) can become smaller than one $N(f,\Delta f)<1$ and the assumption of continuity in Eq. (2) breaks down [89]. The stochastic signal is replaced by a pop corn noise which could in principle be subtracted [90]. We follow the approach of ref. [91] to determine when $N(f,\Delta f)<1$ .

III PBH-enhanced Structure formation and the current JWST observations

Assuming a monochromatic mass function, the PBH population can be parameterized by the common mass, $M_{\rm PBH}$ , and the DM fraction $f_{\rm PBH}$ . At scales where the discrete nature of PBHs becomes relevant, there are two different effects that can influence structure formation.

Poisson effect.

In a region of mass $\tilde{M}$ the condition $f_{\rm PBH}>M_{\rm PBH}/\tilde{M}$ implies that one expects to find more than one PBH $\tilde{N}_{\rm PBH}>1$ . For initially Poisson-distributed PBHs the fluctuations in their number $\sqrt{\tilde{N}_{\rm PBH}}$ generate isocurvature density perturbations of magnitude $\delta_{\rm{PBH},i}\approx(f_{\rm PBH}M_{\rm PBH}/\tilde{M})^{1/2}$ (see appendix A in ref. [92] for a detailed derivation). The perturbations evolve linearly right after matter-radiation equality [34] (i.e. $z<z_{\rm eq}\approx 3400$ ). As a result, the matter power spectrum is modified at smaller scales by the addition of an isocurvature component [68]

		$\displaystyle P(k)=P_{\rm ad}(k)+P_{\rm iso}(k)~{},$
		$\displaystyle P_{\rm iso}(k)\simeq\begin{cases}\frac{(f_{\rm PBH}D(0))^{2}}{\bar{n}_{\rm PBH}}~{},&\text{if }k\leq k_{\rm cut}\\ 0~{},&\text{otherwise}\end{cases}~{},$		(7)

where $P_{\rm ad}(k)$ is the adiabatic mode in $\Lambda$ CDM²²2We generate a prediction for the linear matter power spectrum using the Boltzmann solver CAMB [93], for a fiducial $\Lambda$ CDM cosmology corresponding to the following Planck 2018 [94] values : $n_{s}=0.9649$ , $\sigma_{8}=0.8111$ , $\Omega_{b}=0.0493$ , $\Omega_{\rm DM}=0.266$ and $h=0.6736$ ., $\bar{n}_{\rm PBH}=f_{\rm PBH}\rho_{\rm crit,0}\Omega_{\rm DM}/M_{\rm PBH}$ is the co-moving average number density of PBHs, and $D(z)\simeq\left(1+\frac{3\gamma}{2\alpha_{-}}\frac{1+z_{\rm eq}}{1+z}\right)^{a_{-}}$ is the growth factor with $a_{-}=\left(\sqrt{1+24\gamma}-1\right)/4$ and $\gamma=(\Omega_{\rm DM}-\Omega_{\rm b})/(\Omega_{\rm DM}+\Omega_{\rm b})$ . The isocurvature term is truncated at the scale $k_{\rm cut}$ which is the scale where we expect the linear Press-Schechter theory (PS) [95] to break down. Two different cut-off scales $k_{\rm cut}$ have been discussed, i.e. the inverse mean separation between PBHs, $\bar{k}_{\rm PBH}=(2\pi^{2}\bar{n}_{\rm PBH})^{1/3}$ [68, 96, 24], and the approximate scale $k_{\rm NL}\approx(\bar{n}_{\rm PBH}/f_{\rm PBH})^{1/3}$ [22], where non-linear dynamics (see the seed effect below and the case of mode mixing in ref. [97]) starts to dominate. In the absence of a reliable description of the transition between the linear and the non-linear regimes, we perform our analysis reporting both benchmarks, the conservative scenario $k_{\rm cut}=\bar{k}_{\rm PBH}$ being shown in fig. 1 and the more aggressive one being reported to fig. 2 in the Appendix E.

Equipped with the enhanced power spectrum, we utilize the Sheth-Tormen (ST) [98, 99] modification to the PS formalism with a top-hat window function (see Appendix A) and calculate the halo mass function $dn(M_{h},z)/dM_{h}$ . The expected number density of galaxies with stellar mass above the observational threshold, $M_{\star}^{\rm obs}$ , is given by [19]

n_{\rm gal}(M_{\star}\geq M_{\star}^{\rm obs})=\int_{M^{\rm cut}_{h}}^{\infty}\frac{dn(z_{\rm obs},M_{h})}{dM_{h}}dM_{h}~{}.

(8)

Under the assumption that each dark-matter halo contains a single central galaxy, the relation between the halo and total stellar mass is $M_{h}(M_{\ast})=M_{\ast}/(f_{\rm b}\epsilon_{\ast})$ , where $f_{b}=\Omega_{\rm b}/(\Omega_{\rm DM}+\Omega_{\rm b})=0.157$ is the baryon fraction and $\epsilon_{\ast}$ the star formation efficiency. The lower limit of the integral in eq. (8) is defined then as ${M^{\rm cut}_{h}}=M_{h}(M_{\star}^{\rm obs})$ . The JWST signature can be expressed as $n_{\rm gal}(M_{\star}\geq 10^{10.8}~{}M_{\odot})\simeq 10^{-5}\rm Mpc^{-3}$ at $z_{\rm obs}\sim 8$ [14].

Seed effect.

In the opposite limit, $f_{\rm PBH}<M_{\rm PBH}/\tilde{M}$ , the PBHs make up only a small fraction of dark matter and form isolated halos. The density perturbations in this case are $\delta_{\rm{PBH},i}\approx M_{\rm PBH}/\tilde{M}$ , binding gravitationally regions of mass $\tilde{M}(M_{\rm PBH},z)\simeq z_{\rm eq}M_{\rm PBH}/(z+1)$ . Due to its highly non-linear nature, this effect can be examined properly only using simulations [68]. However, following the calculation of refs. [22, 32] we can determine the part of the parameter space that can potentially accommodate accelerated early galaxy formation compatible with the JWST observations solely with the seed effect. The three basic requirements are: (i) $f_{\rm PBH}<M_{\rm PBH}/\tilde{M}=(z+1)/z_{\rm eq}$ , (ii) $\bar{n}_{\rm PBH}>10^{-5}\rm Mpc^{-3}$ , and (iii) $\tilde{M}(M_{\rm PBH},z_{\rm obs}\sim 8)\simeq M_{h}(M_{\ast}^{\rm obs}\sim 10^{11}M_{\odot})$ , assuming also in this case that the whole bound region hosts a single galaxy.

IV Observational constraints

In this section we review the constraints on the PBH populations with masses above $10^{6}~{}M_{\odot}$ and introduce the UV LF bound on the PBH parameter space.

1.

Large-scale structure and dynamical friction. As discussed in Section III, PBHs can significantly contribute to the formation of cosmic structures. In fact, the non-observation of different types of structures can be used to constrain population of PBHs [67]. We comment however that those constraints are approximative as they primarily apply to the bulk of various types of structures. For specific cases (e.g. the JWST early galaxies), the limits may vary within an order of magnitude.

The dynamical friction (DF) limit concerns the accumulation of halo BHs into the galactic nucleus. The effect is induced due to dynamical friction by stellar populations and the subsequent merging of those BHs in the nucleus would result into SMBHs heavier than the currently observed ones unless $f_{\rm PBH}$ is rather small. The relevant derivation can be found in ref. [100].
2.

CMB $\mu$ distortion.

The PBH formation from large-amplitude primordial fluctuations leaves imprints in the CMB spectrum, called $\mu$ distortion [101]. The defining observable $\mu$ is strictly constrained by the COBE/FIRAS measurements to be smaller than $4.7\times 10^{-5}$ [102, 103], which translates into a strong bound on the PBH abundance over the mass range $10^{5}~{}M_{\odot}<M_{\rm PBH}<10^{12}$ . If the primordial fluctuations are Gaussian then practically all PBH models of interest are excluded [66]. ³³3See however the models of refs. [104, 105, 106, 107, 108, 109, 65, 110, 111, 112, 113, 114, 115, 116, 117], which postulate different PBH formation mechanisms and may avoid bounds from CMB spectral distortion altogether. On the other hand, if one allows for significant non-Gaussianities (NG) in the curvature power spectrum the bound can be avoided. In this work, we adopt the phenomenological description of ref. [65] and characterize NGs by the parameter $p\leq 2$ (where $p=2$ corresponds to the Gaussian case). All the relevant expressions can be found in Appendix B.
3.

HST UV luminosity function.

The UV LF of galaxies observed by the Hubble Space Telescope (HST) [74, 75, 76, 77, 78, 79, 80, 72, 81, 82, 73] has emerged as a sensitive early universe probe to constrain not only $\Lambda$ CDM [74, 118, 119, 120, 121, 122, 123], but also a broad portfolio of extensions about it [124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134]. As it was recently pointed out [135], more importantly, any models that attempt to modify the high-redshift halo mass function in relation to the JWST excess of massive galaxies [14] are already strongly limited by the HST UV LF, since the latter traces galaxies of the same redshifts and distances as the ones in the JWST sample. Such is the case of the massive PBHs that we consider in this work, which predict an enhancement of dark matter clustering as discussed in section III, and which have indeed been found theoretically capable of producing a desired excess of massive galaxies in the range $z_{\rm obs}=7-10$ [22, 24].

In this letter, we take advantage of the above connection and use HST UV LF measurements [72, 73] to constrain the properties of PBH populations using this probe. In particular, we proceed as follows: combining eqs. (III), (8) and the ST formalism further explained in Appendix A, we produce theoretical predictions for the galaxy number density above a stellar mass cut-off $M_{\ast}$ in the presence of a PBH population with mass $M_{\rm PBH}$ and abundance $f_{\rm PBH}$ , for $\epsilon_{\star}=0.3$ at $z_{\rm obs}\sim 8$ . These predictions are then contrasted against the maximum deviation from the corresponding $\Lambda$ CDM result that is allowed by the HST UV LF measurements, as it was reported by the likelihood analysis of [135] at the $95\%$ C.L. (see fig. 3 of that work).⁴⁴4We note that the likelihood analysis of ref. [135] considers model-agnostic power spectrum enhancements that were general enough to encompass modifications of the specific form in eq. (III), that we are using, marginalized over variations in the power spectrum amplitude and various astrophysical parameters characterizing the galaxy-halo connection.

V Discussion and conclusions

Refer to caption — Figure 1: Illustration of the PBH parameter space in the mass region $10^{6}~{}M_{\odot}<M_{\rm PBH}/M_{\odot}<10^{13}~{}M_{\odot}$ . The green thick line and trapeze indicate PBHs which can explain JWST’s observation of early massive galaxies assuming a redshift of observation of $z\simeq 8$ and a star formation rate of $\epsilon_{\ast}\simeq=0.3$ . The blue and orange ovals represent the $68\%$ and $90\%$ confidence levels for PTA’s interpretation in terms of GW from PBH binaries. The two dashed brown lines indicate where the GW-driven and smooth assumption used to derive the GW spectrum in eq. (2) break down due to environmental effects, below $f<f_{\rm min}$ defined in eq. (6), and due to the number of emitting sources per frequency bin in eq. (39) is $N(\Delta f)<1$ . The red region shows the constraints on PBHs due to measurements of the galaxy UV luminosity function by the Hubble Space telescope, while the gray regions show the constraints from large-scale structure and dynamical friction. The $\mu$ -distortion bounds are depicted with the dotted lines for two different degrees of primordial non-Gaussianities, $p=0.3,0.5$ , and the Gaussian scenario $p=2$ . The PBH isocurvature power spectrum is suppressed above $k_{\rm cut}=(2\pi^{2}\bar{n}_{\rm PBH})^{1/3}$ . A different choice for $k_{\rm cut}$ is shown in fig. 5.

The collective results of our analysis can be found in fig. 1. This includes the regions on the $M_{\rm PBH}$ - $f_{\rm PBH}$ plane that can account for both JWST and PTA signatures as discussed in sections III and II, in conjunction with the constraints discussed in section IV.

We perform a Bayesian analysis using the software tool ${\tt PTArcade}$ [136] of NANOGrav 15-year [2] and IPTA DR2 signals [8] in terms of GW from PBH mergers. The constraints $f_{\rm min}<f$ and $N(\Delta f)>1$ (see Appendix D for more details), signal to be stochastic, are imposed as priors and displayed as red dashed lines on Fig. 1. We find that PBH binaries with mass $10^{8}\lesssim M_{\rm PBH}/M_{\odot}\lesssim 10^{11}$ and abundance $f_{\rm PBH}\sim 10^{-2}$ can accommodate both PTA signals (68 $\%$ and 95 $\%$ C.L. are blue and orange ovals).

Assuming a Poissonian gas of PBHs and using linear cosmological perturbation theory and the ST approach, we model the enhancement to the galaxy abundance as a function of the grid of PBH model parameters. Two different wavelength cut-offs $k_{\rm cut}$ are used, the conservative one presented in fig. 1 of the main text, and the more aggressive one being reported in SuM V. The observation of massive galaxies at redshift $z\sim 8$ by the JWST is explained along the green thick line. The maximum value of $M_{\rm PBH}$ corresponds to a bound region mass of $\mathcal{O}(10^{12}M_{\odot})$ , which is the typical upper theoretical limit for halo masses [137, 67]. Additionally, using the same assumption for the linear halo evolution, we project the $95\%$ exclusion bound imposed by the UV LF on the parameter space (red shaded region). The constraints from large-scale structure and dynamical friction are also shown (gray regions).

Regarding PBHs evolving in isolation, we display within the green trapeze region of fig. 1, the region that could explain the JWST data with the seed effect. The upper and lower lines correspond to conditions (i) and (ii) discussed in section III (see paragraph on seed effect), respectively. The two vertical horizontal lines constrain the resulting halo masses to be $\mathcal{O}(10^{10}M_{\odot}-10^{13}M_{\odot})$ , allowing for sufficient variation around $M_{h}(M_{\ast}^{\rm obs}\sim 10^{11}M_{\odot})$ in order to compensate for unknown uncertainties due to the simplicity of the calculation.

Finally, it is worth noticing that a common requirement for all scenarios of interest involving PBH populations in the superheavy mass region, is the formation of supermassive PBHs from NG primordial fluctuations in order to circumvent the strict bound from CMB $\mu$ distortion. In particular, in fig. 1 we depict the $\mu$ bound given $p=0.3$ (gray dotted lines), which is marginally compatible with the depicted JWST solution based on the seed effect. This is, in fact, an extremely large degree of NGs. For instance, the benchmark $p=0.5$ , which already corresponds to the infinite limit of more traditional NG measures [66], excludes all of the interesting PBH populations.

Conclusions.

In fig. 1, we demonstrate that the PBH populations needed to source the stochastic GW background observed in PTAs are partly excluded by LSS and decisively excluded by the UV LF constraint derived in this work. Similarly, the PBH population required to explain the JWST anomalous observation with the Poisson effect is excluded due to the same constraint. The novel UV LF bound is more rigorous than the former LSS bound of PBHs, because it has been inferred from an extended likelihood analysis of the HST data at the redshift range contemporaneous with the JWST extreme galaxy sample. The tight constraints we obtain are in agreement with the findings of ref. [135].

On the other hand, the PBH solution of JWST observation based on the seed effect is in principle still viable for $f_{\rm PBH}\lesssim 10^{-3}$ up to the caveat of extreme primordial NGs. A clearer picture of the reach of the seed effect necessitates the usage of dedicated simulations in the $M_{\rm PBH}>10^{7}M_{\odot}$ range, as well as an incorporating of non-linear effects in the derivation of the UV LF. Looking forward, a spectroscopic analysis will provide the final verdict on whether the JWST observations constitute a $\Lambda$ CDM anomaly. Finally, a future increase in observation time of PTAs and in number of detected pulsars might facilitate the resolution of individual sources at larger frequencies and thus enable the more careful examination of new physics interpretations.

Acknowledgements.

We would like to thank Sebastien Clesse, Edoardo Vitagliano, Daniel Eisenstein and Julian Munoz for useful discussions throughout the completion of this work, as well as Omer Katz for sparing computational resources on the TAU cluster during the preparation of this work. YG is grateful to the Azrieli Foundation for the award of an Azrieli Fellowship. ST is supported by the Swiss National Science Foundation - project n. P500PT_203156, and by the Center of Theoretical Physics at MIT (MIT-CTP/5538). GV recognizes partial support by the Department of Energy (DOE) Grant No. DE-SC0020223. MV is supported by the “Excellence of Science - EOS” - be.h project n.30820817, and by the the Strategic Research Program High-Energy Physics of the Vrije Universiteit Brussel.

Appendix A Halo mass function

We consider a spherical dark matter halo at redshift $z$ , containing an average mass $\tilde{M}$ within a region of effective radius $R$ . If, additionally, $\rho_{m}$ is the average matter density of the universe at $z$ , it follows that

R=\left(\frac{3\tilde{M}}{4\pi\rho_{m}}\right)^{\frac{1}{3}}~{}.

(9)

The root-mean square of matter density fluctuations in this region can be further calculated through a convolution of the the linear matter power spectrum, $P(k,z=0)$ , with a spherical top-hat smoothing Kernel of radius $R$ :

W\left(kR\right)=\frac{3\left[\sin(kR)-kR\cos(kR)\right]}{\left(kR\right)^{3}}~{},

(10)

which results in the (Fourier space) integral

\sigma^{2}(\tilde{M})=\int\frac{dkk^{2}}{2\pi^{2}}W^{2}\left(kR\right)P(k,z=0)~{}.

(11)

The PS theory and its generalizations [95, 138] postulate that the abundance of halos in the universe can be evaluated by counting all the regions of radius $R$ (and, equivalently, mass $\tilde{M}$ ) that have gravitationally collapsed by the time of interest $z$ . According to the excursion set approach, these will be the regions with a smoothed density exceeding the critical overdensity for spherical collapse at redshift $z$ , $\delta_{cr}$ . For an Einstein De-Sitter (EDS) cosmology, it is always $\delta_{cr}=1.686$ [139, 138], which turns out to be a very good approximation for $\Lambda$ CDM cosmologies and will be adopted in this work as well. Assuming, finally, density fluctuations that are accurately described by linear perturbation theory and a Gaussian probability distribution, we can derive an analytical prediction for the comoving mean number density of halos $n_{h}(M)$ per logarithmic mass interval $d\ln M$ , given by:

\frac{dn_{h}}{d\ln M}=\frac{\rho_{m}}{M}\nu_{c}(M)f\left(\nu_{c}(M)\right)\frac{d\ln\nu_{c}(M)}{dM}~{},

(12)

where the quantity $\nu_{c}(M)$ , called the peak significance, is defined as

\nu_{c}(M)=\frac{\delta_{cr}}{\sigma(M,z)}=\frac{\delta_{cr}}{D(z)\sigma(M)}~{}.

(13)

In eq. (13), $\sigma(M,z)$ denotes the density variance at redshift $z$ , evaluated through the linear evolution of $\sigma(M)$ (obtained from eq. (11) to the time of collapse $z$ , using the linear growth factor $D(z)$ . For the latter we adopt the standard normalization $D(z=0)=D_{0}=1$ . We note again that for the critical overdensity to collapse we use the EDS solution $\delta_{cr}=1.686$ .

In the original PS theory, the multiplicity function, $f\left(\nu_{c}(M)\right)$ , had the simple form:

\nu_{c}f\left(\nu_{c}\right)=\sqrt{\frac{2}{\pi}}\nu_{c}e^{-\frac{\nu_{c}^{2}}{2}}~{},

(14)

which is exact in an EDS universe described by a power-law power spectrum. While this prescription, often referred to as the universal mass function, has been used to describe the halo mass function for a broad range of cosmologies, it lacks the necessary accuracy for predictions in the era of precision cosmology. As a result, ST [98] later introduced an alternative function:

\nu_{c}f\left(\nu_{c}\right)=\sqrt{\frac{2}{\pi}}A(p)\left[1+\frac{1}{(q\nu_{c}^{2})^{p}}\right]\sqrt{q}\nu_{c}e^{-\frac{q\nu_{c}^{2}}{2}}~{},

(15)

with $A(p)=\left[1+\pi^{-\frac{1}{2}}2^{-p}\Gamma(0.5-p)\right]^{-1}$ and where $q,p$ are free parameters to be fitted over N-body simulations, reducing to the vanilla PS function for $q=1,p=0$ . The best fit pair was initially proposed to be $(q,p)=(0.707,0.3)$ which was later updated to $(q,p)=(0.75,0.3)$ , commonly considered to be the “standard” ST parameters [98, 99]. When focusing on higher redshift probes like in our work, however, a slightly different choice of values, $(q,p)=(0.85,0.3)$ , has been found to provide a more accurate fit [140], which is the one we adopt in order to match the analysis of [135] that we compare against. We also comment on the fact that we experimented with both choices and found our results to change minimally across the tests, demonstrating the robustness of our analysis against the modeling details of the halo mass function, as was also reported by [135].

To summarize, given a specific choice of the multiplicity function (14), eqs. (9)-(13) provide a prediction for the halo mass function for a given cosmology, that we subsequently use to obtain our galaxy observables of interest from eq. (8) of the main text (MT). Its sensitivity to the underlying cosmological model manifests through the linear matter power spectrum entering eq. (11), which we obtain from the Boltzmann solver CAMB [93] in the $\Lambda$ CDM case and additionally using MT-eq. (7) to account for the PBH-induced effects.

Appendix B PBHs from non-Gaussian perturbations and $\mu$ distortion

According to the standard hypothesis of PBH formation large-amplitude primordial fluctuations undergo spherical gravitational collapse upon horizon reentry [34]. Fluctuations that dissipate via Silk damping [141] during the photon diffusion scale, i.e. $5\times 10^{4}<z<2\times 10^{6}$ , inject energy in the photon bath and modify the number of photons at different frequencies w.r.t the black-body equilibrium.

Let us start with the probability density function for the primordial curvature perturbations, which can be parameterized as [65]

P(\zeta)=\frac{1}{2\sqrt{2}\tilde{\sigma}\Gamma\left(1+1/p\right)}\exp\left[-\left(\frac{|\zeta|}{\sqrt{2}\tilde{\sigma}}\right)^{p}\right]~{},

(16)

where $p=2$ corresponds to the the Gaussian limit of ref. [36]. The variance is

\sigma^{2}=\int_{-\infty}^{\infty}\zeta^{2}P(\zeta)d\zeta=\frac{2\Gamma(1+3/p)}{3\Gamma(1+1/p)}\tilde{\sigma}^{2}~{},

(17)

and the abundance of PBH is

\beta=\int_{\zeta_{c}}^{\infty}P(\zeta)d\zeta=\frac{\Gamma(1/p,2^{-p/2}(\zeta_{c}/\tilde{\sigma})^{p})}{2p\Gamma(1+1/p)}~{},

(18)

where $\zeta_{c}\simeq 0.67$ is the threshold for PBH formation [142], $\Gamma(a)$ is the gamma function and $\Gamma(a,z)$ the incomplete gamma function.

It is shown in ref. [101] that a modification to the curvature power spectrum of the form

{\Delta\cal P}_{\zeta}=2\pi^{2}\sigma^{2}k^{-2}\delta(k-k_{\delta})~{},

(19)

which exhibits an extremely sharp feature at some scale $k_{\delta}$ , generates the following $\mu$ distortion

\mu\simeq 2.2\sigma^{2}\left[\exp\left(-\frac{\hat{k}_{\delta}}{5400}\right)-\exp\left(-\left[\frac{\hat{k}_{\delta}}{31.6}\right]^{2}\right)\right]~{}.

(20)

where $\hat{k}_{\delta}=k_{\delta}\rm Mpc$ .

For a PBH population of mass $M_{\rm PBH}$ ⁵⁵5Eq. (19) does not necessarily correspond to a monochromatic PBH mass function since critical collapse may broaden the mass spectrum [143]. the scale $k_{\delta}$ is given by [144]

k_{\delta}\simeq 13\mathrm{Mpc}^{-1}\gamma^{1/2}\left(\frac{g}{10.75}\right)^{-1/12}\left(\frac{M_{\rm PBH}}{10^{11}M_{\odot}}\right)^{-1/2}~{},

(21)

where $\gamma$ gives the size of the PBH in units of the horizon mass at formation (we take it to be $\gamma=1$ ) and $g$ is the number of degrees of freedom of relativistic particles. The initial abundance $\beta$ is related to $f_{\rm PBH}$ [145, 144], via

\beta\simeq 6\times 10^{-4}\left(\frac{g}{10.75}\right)^{1/4}\left(\frac{\Omega_{\mathrm{DM}}}{0.27}\right)^{-1}\left(\frac{M_{\rm PBH}}{10^{11}M_{\odot}}\right)^{1/2}f_{\rm PBH}~{}.

(22)

Replacing this expression in eq. (18), we can solve for $\tilde{\sigma}$ and substitute in eq. (17) to get $\sigma$ and finally calculate $\mu$ using eq. (20).

We notice also that it is not straightforward to connect the expression of eq. (16) to known models of inflation [146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156]. In more realistic scenarios that exhibit local NGs, one can examine also modifications to the shape of the overdensity and thus the critical threshold for collapse [149, 157, 158, 159]. However, NGs corresponding to values $p<0.5$ cannot be achieved by traditional quadratic $f_{\rm NL}$ and cubic $g_{\rm NL}$ parameters [66].

Appendix C PBH binary distribution during radiation domination

In this appendix, we present the computation of the distribution of binaries as a function of eccentricity and major axis. We follow closely [40] (see also refs. [45, 160] and ref. [83] for the earlier proposals). We need first to compute the moment during radiation domination when the PBH pair decouples from the expansion of the universe. This occurs when the energy density enclosed by the two PBHs becomes larger than the radiation energy density in the same sphere

\displaystyle 1+z_{\rm dec}=(1+z_{\rm eq})\bigg{(}\frac{x_{\rm max}}{x}\bigg{)}^{3}~{},

(23)

where $x$ is the comoving separation between the two PBHs and $x_{\rm max}=f_{\rm PBH}^{1/3}l_{\rm PBH}(z=0)$ . The closest third PBH, at comoving distance $y$ , exerts a torque on the binary and prevents head-on collision and transmits angular momentum to the now rotating binary. The distribution of binaries, parameterized by the semi-major axis $a$ and the eccentricity $e$ , can be described by

	$\displaystyle dP\simeq$	$\displaystyle\,\frac{4\pi x^{2}dx}{n_{\rm PBH}^{-1}}\frac{4\pi y^{2}dy}{n_{\rm PBH}^{-1}}\Theta(y-x)\Theta(y_{\rm max}-y)~{},$
	$\displaystyle=$	$\displaystyle\,\frac{4\pi^{2}}{3}n_{\rm PBH}^{1/2}f_{\rm PBH}^{3/2}(1+z_{\rm eq})^{3/2}a^{1/2}e(1-e^{2})^{-3/2}de\,da~{},$		(24)

with $a=\frac{\rho_{c,0}\Omega_{\rm DM}x^{4}}{(1+z_{\rm eq})M_{\rm PBH}}$ , $e=\sqrt{1-(x/y)^{6}}$ , and $y_{\rm max}=\big{(}\frac{4\pi}{3}n_{\rm PBH}\big{)}^{-1/3}$ . This leads to a maximum on $a$ and $e$ ,

	$\displaystyle a_{\rm max}=$	$\displaystyle\,\frac{x_{\rm max}}{1+z_{\rm eq}}~{},$		(25)
	$\displaystyle e_{\rm max}=$	$\displaystyle\,1-\bigg{(}\frac{4\pi n_{\rm PBH}}{3}\bigg{)}^{2}\bigg{(}\frac{(1+z_{\rm eq})M_{\rm PBH}}{\rho_{c,0}\Omega_{\rm DM}}a\bigg{)}^{3/2}~{}.$		(26)

Levels of approximation.

The analysis of the binary formation and the derivation of the GW background we presented and used in this letter follows the prescription of ref. [45]. Firstly, it relies on the three body approximation, and several effects might bring corrections to the rate we estimated. However, it was shown in ref. [161] that those effects are $\mathcal{O}(1)$ corrections and would not substantially alter the results of our analysis. Moreover, for low enough $f_{\rm PBH}\ll 1$ , the disruption of binaries formed by close encounter with a third PBH is expected to be small [49, 96].

The analyse derived in ref. [45] which we follow analyzes mergers from PBH binary with masses $M_{\rm PBH}\sim 10^{2}M_{\odot}$ . Since the relevant physics in the radiation-dominated universe is mass independent, we don’t expect any special difference between the properties of PBH binaries with $M_{\rm PBH}\sim 10^{2}M_{\odot}$ and the ones in the superheavy regime $M_{\rm PBH}>10^{6}M_{\odot}$ , such as for instance the eccentricities and the major axis (see also ref. [162, 70]).

We find that the late-time contribution to binary formation, using the formula in ref. [45], is subdominant to the early-time contribution for $M_{\rm PBH}\lesssim 10^{9}~{}M_{\odot}$ , assuming $f_{\rm PBH}\sim 10^{-2}$ [45]. The latter accounts for the binary formation due to encounter within a dark matter halo, and energy loss through GW emission. At larger masses, $M_{\rm PBH}\gtrsim 10^{9}~{}M_{\odot}$ , the formula breaks down since there are less than 2 PBHs per dark matter halo, assuming a typical halo mass $M_{\rm halo}\sim 10^{12}~{}M_{\odot}$ . Dedicated studies supported by $N$ -body simulations are required for a more precise determination of binary formation at late time in this mass range, considering contributions induced by halo mergers and energy loss through dynamical friction [163, 164, 165]. We notice that the latter effect leads to constraints derived in ref. [100], excluding PBHs of mass $M_{\rm PBH}\sim 10^{9}~{}M_{\odot}$ unless $f_{\rm PBH}\lesssim 10^{-3}$ .

Finally, PBHs can grow in mass due to accretion [166]. The effect depends on parameters whose values are generically uncertain for all SMBHs, see e.g. [167]. We notice that data seems to suggest that accretion becomes inefficient at large masses $M_{\rm PBH}\gtrsim 10^{9}~{}\rm M_{\odot}$ , e.g. see Fig. 1 in [156]. For those two reasons, we neglect accretion in this analysis.

Appendix D Derivation of gravitational wave signal and PTA observations

After their formation during radiation era, the two PBHs pertaining to the binary will orbit around each other and lose energy mostly via gravitational radiation. We can estimate their time of merging, assuming $e\to 1$ , by [168]

t_{\rm m}(a,e)\simeq\frac{3}{170}\frac{1}{G^{3}M^{3}_{\rm PBH}}a^{4}(1-e^{2})^{7/2}~{}.

(27)

This expression of the time of merging $t_{\rm m}(a,e)$ can be inverted to eliminate $a(t_{\rm m},e)$ in eq. (24) to express the probability distribution as a function of $t_{\rm m}$ and $e$ . This leads to the expression [40]

\frac{dP}{dt_{\rm m}}\simeq\frac{3}{58}\bigg{(}\frac{t}{T}\bigg{)}^{\frac{3}{8}}\frac{1}{t_{\rm m}}\bigg{[}\frac{1}{(1-e_{\rm up})^{\frac{29}{16}}}-1\bigg{]}~{},

(28)

where the typical time of merging is

T=\frac{3}{170}\frac{1}{G^{3}M^{3}_{\rm PBH}}\bigg{(}\frac{3y_{\rm max}}{4\pi f_{\rm PBH}(1+z_{\rm eq})}\bigg{)}^{4}~{},

(29)

and the maximal eccentricity reads

e_{\rm up}=\,\begin{cases}\sqrt{1-(t_{\rm m}/T)^{6/37}}\qquad\text{for }t_{\rm m}<t_{c}~{},\\ \sqrt{1-(4\pi f_{\rm PBH}/3)^{2}(t_{\rm m}/t_{c})^{2/7}}\qquad\text{for }t_{\rm m}>t_{c}\,.\end{cases}

(30)

The time of the transition between the two regimes is of the form $t_{c}=T(4\pi f_{\rm PBH}/3)^{37/3}$ . The final expression for the rate of merging of PBHs can be obtained by multiplying this probability with the PBHs density, $n_{\rm PBH}$ , and we obtain

\displaystyle\mathcal{R}(t_{\rm m})\simeq

\displaystyle\,\frac{3n_{\rm PBH}}{58}\bigg{(}\frac{t}{T}\bigg{)}^{\frac{3}{8}}\frac{1}{t}\bigg{[}\frac{1}{(1-e_{\rm up})^{\frac{29}{16}}}-1\bigg{]}~{}.

(31)

The gravitational wave amplitude emitted from the merger of two PBHs is then [86]

h^{2}\Omega_{\rm GW}=\frac{f}{\rho_{c}/h^{2}}\int_{0}^{\infty}dz\frac{\mathcal{R}(z)}{(1+z)H(z)}\frac{dE_{\rm GW}(f^{\prime})}{df^{\prime}}\bigg{|}_{f^{\prime}=(1+z)f}\,,

(32)

where $dE_{\rm GW}(f^{\prime})/df^{\prime}$ is the GW power emitted by an individual binary [169, 170]

	$\displaystyle\frac{dE_{\rm GW}(f)}{df}=$	$\displaystyle\,\frac{(G\pi)^{2/3}M_{c}^{5/3}}{3}$
	$\displaystyle\times$	$\displaystyle\,\begin{cases}f^{-1/3}~{},\quad f<f_{1}\\ f^{2/3}f_{1}^{-1}~{},\quad f_{1}<f<f_{2}\\ f_{2}^{-4/3}f_{1}^{-1}\bigg{(}\frac{f}{1+4(\frac{(f-f_{2})}{\sigma})^{2}}\bigg{)}^{2}~{},\quad f_{1}<f<f_{2}\end{cases}$		(33)

where $M_{c}$ is the chirp mass, $M_{c}=M_{\rm PBH}/2^{1/5}$ . We approximate the GW power spectrum in eq. (32) by the fitting function power broken-law

\displaystyle\Omega_{\rm GW}h^{2}\simeq\,\Omega_{\rm peak}S(f)\Theta(2f_{\rm peak}-f)~{},

(34)

where the spectral function is

	$\displaystyle S(f)=\,\frac{f_{\rm peak}^{b}f^{a}}{\bigg{(}bf^{\frac{a+b}{c}}+af_{\rm peak}^{\frac{a+b}{c}}\bigg{)}^{c}}~{},$
	$\displaystyle a=0.7,\quad b=1.5,\quad c=0.9$		(35)

the peak amplitude is

\displaystyle\Omega_{\rm peak}\simeq\,0.05f_{\rm PBH}^{3}\bigg{(}1.5\frac{M_{\rm PBH}}{10^{12}M_{\odot}}\bigg{)}^{-0.3}~{},

(36)

and the peak frequency is $f_{\rm peak}\simeq\,5000M_{\odot}/M_{\rm PBH}$ . The goodness of the approximation is shown in fig. 2.

Continuous signal.

The signal observed at PTAs is a continuous signal. In this paragraph, we explain how to determine if the merger of PBH would procude a continuous signal or a pop-corn like signal. The merging rate $\mathcal{R}$ can be related to the comoving number of binaries emitting in the logarithmic interval, in the interval of redshift $[z,z+dz]$ , via

\displaystyle\frac{dN_{\rm mergers}}{dzd\log f}=-\mathcal{R}(t)\frac{dV_{c}}{dz}\frac{dt}{d\log f}~{}.

(37)

We can further define $\tau_{r}=dt/d\log{f}$ , known as the residence time [171], as the amount of time that each binary spends emitting in a logarithmic frequency interval, which takes the following form [172]

\tau_{r}=\frac{5}{96\pi^{8/3}}\left(\frac{1}{GM_{c}}\right)^{5/3}f^{-8/3}~{},

(38)

where $M_{c}=M_{\rm PBH}/2^{1/5}$ in the chirp mass. We can see that $dV_{c}/dz=4\pi d^{2}_{c}(z)/H(z)$ , where $d_{c}$ is the comoving distance, and thus the number of mergers can be obtained by integrating over $f$ and $z$ ,

\displaystyle N(f,\Delta f)=\ \int^{f_{0}+\Delta f/2}_{f_{0}-\Delta f/2}\frac{df}{f}\ \int^{\infty}_{0}\ dz\frac{4\pi d^{2}_{c}(z)}{H(z)}R(z)\tau_{r},~{}

(39)

where $f_{0}\simeq 2~{}\textrm{nHz}$ is the frequency of the first bin and $\Delta f=14f_{0}$ the number of bins considered. We notice that the computation depends only very weakly on $\Delta f$ .

Details of the Bayesian analysis.

In this letter we performed a Bayesian analysis of the PBH mergers interpretation of the PTA data GW signal using eq. (34), which we truncate below the IR cut-off frequency $f_{\rm min}$ when orbital energy is controlled by star slingshot. We impose the condition to have more than one emitting source per frequency bins, $N(\Delta f)>1$ (see eq. (39)) as a prior in our analysis. To search for a GW power spectrum from supermassive PBH binaries in the timing residual data, we used the wrapper ${\tt PTArcade}$ [136] which relies on the software tools ${\tt enterprise}$ [173] and ${\tt enterprise\_extensions}$ [174] to calculate the likelihood accounting for all the different sources of noise. We included the first 14 frequency bins for NANOGrav 15-year dataset [2] and the first 13 frequency bins for IPTA DR 2 dataset [8]. The chains are generated with the parallel-tempering Markov Chain Monte-Carlo sampler ${\tt PTMCMC}$ [175]. Finally, we use the software ${\tt GetDist}$ [176] to visualize the posterior distribution, shown with blue and orange ellipses in MT-fig. 1 and fig. 5.

Case of SGWB from combined PBHs and ABHs.

The favored interpretation of the PTA signal is attributed to SMBH binary of astrophysical origins [177, 178]. Those are supposed to form at the center of halos in the late universe with a mass correlated to the halo mass [166]. For completeness, we perform a Bayesian analysis of the combined GW spectrum from primordial and astrophysical black holes (PBHs and ABHs). The typical strain spectrum for a set of coalescing, GW-driven SMBH binaries follows a red-tilted power-law [179]

h_{c}(f)=A_{\text{ABH}}\left(\frac{f}{1~{}\text{yr}^{-1}}\right)^{-2/3},

(40)

with $A_{\text{ABH}}$ denoting the strain amplitude at a frequency of $1~{}\text{yr}^{-1}\simeq 3.2\times 10^{-8}~{}\rm s^{-1}$ . When expressed as a fraction of energy density, this relates to a blue-tilted power-law:

\Omega_{\text{ABH}}(f)=\frac{2\pi^{2}}{3H_{0}^{2}}f^{2}h_{c}^{2}(f)\sim f^{2/3}.

(41)

We perform a joint analysis for gravitational waves emanating from both PBH and ABH binary systems. We show the credible intervals in fig. 3.

Appendix E PBH parameter space with a higher cut-off

In order to highlight the complementarity between the UV LF and other standard cosmological probes, we plot, in fig. 4, the modification to the $\Lambda$ CDM power spectrum that is predicted by a supermassive PBH population, together with galaxy observations by the SDSS collaboration [180] and the Planck 2018 [94] CMB measurements⁶⁶6In principle, Lyman- $\alpha$ measurements do access wavemodes in the range $0.1\lesssim k/\rm{Mpc}^{-1}\lesssim 10$ [181, 182, 183]. However, it is unclear how to extrapolate the Lyman- $\alpha$ bound into the region $M_{\rm PBH}>10^{5}M_{\odot}$ [69, 41]. As a result, we do not include them in Fig. 3.. Even though the observational landscape is soon expected to significantly change with the advent of Stage-IV surveys such as DESI [184], Euclid [185] and the Rubin Observatory LSST [186], we note that such deviations are currently not yet ruled out by standard constraints, highlighting the importance of the UV LF for this study.

We then proceed to discuss the results of the main analysis using the same cut-off $k_{\rm cut}\approx(\bar{n}_{\rm PBH}/f_{\rm PBH})^{1/3}$ on the isocurvature component in MT-eq. (7) as it is assumed in the heuristic approach of ref. [22]. To illustrate the difference from the cut-off $k_{\rm cut}=\bar{k}_{\rm PBH}$ used in the MT, we plot in fig. 4 the modification to the power spectrum predicted by both cut-offs. As it is evident, this choice allows for higher wavenumbers to contribute and thus the boost to the halo mass function is significantly increased. As a result, the curve in MT-fig. 1 where the JWST observations are satisfied as well as the bound derived from the UV LF are both displaced towards lower $f_{\rm PBH}$ values, as clearly seen in fig. 5.

It is noteworthy that in this case the solution based on the seed effect is partially subjected to the UV LF bound. The underlying assumption is that at intermediate scales between the linear- and non-linear regimes the two effects co-exist and even though the seed effect dominates, the Poisson may still contribute at a smaller percentage (see also fig. 6 in ref. [68], but notice that the authors employ linear perturbation theory only above $k_{\rm cut}=\bar{k}_{\rm PBH}$ ).

References

[1] NANOGrav Collaboration, Z. Arzoumanian et al. Astrophys. J. Lett. 905 (2020), no. 2 L34, [arXiv:2009.04496].
[2] NANOGrav Collaboration, G. Agazie et al. Astrophys. J. Lett. 951 (2023), no. 1 L8, [arXiv:2306.16213].
[3] S. Chen et al. Mon. Not. Roy. Astron. Soc. 508 (2021), no. 4 4970–4993, [arXiv:2110.13184].
[4] J. Antoniadis et al. arXiv:2306.16214.
[5] B. Goncharov et al. Astrophys. J. Lett. 917 (2021), no. 2 L19, [arXiv:2107.12112].
[6] D. J. Reardon et al. Astrophys. J. Lett. 951 (2023), no. 1 L6, [arXiv:2306.16215].
[7] H. Xu et al. Res. Astron. Astrophys. 23 (2023), no. 7 075024, [arXiv:2306.16216].
[8] J. Antoniadis et al. Mon. Not. Roy. Astron. Soc. 510 (2022), no. 4 4873–4887, [arXiv:2201.03980].
[9] NANOGrav Collaboration, A. Afzal et al. Astrophys. J. Lett. 951 (2023), no. 1 L11, [arXiv:2306.16219].
[10] J. Antoniadis et al. arXiv:2306.16227.
[11] N. J. Adams, C. J. Conselice, L. Ferreira, D. Austin, J. A. A. Trussler, I. Juodž balis, S. M. Wilkins, J. Caruana, P. Dayal, A. Verma, and A. P. Vijayan Monthly Notices of the Royal Astronomical Society 518 (nov, 2022) 4755–4766.
[12] R. P. Naidu, P. A. Oesch, P. van Dokkum, E. J. Nelson, K. A. Suess, G. Brammer, K. E. Whitaker, G. Illingworth, R. Bouwens, S. Tacchella, J. Matthee, N. Allen, R. Bezanson, C. Conroy, I. Labbe, J. Leja, E. Leonova, D. Magee, S. H. Price, D. J. Setton, V. Strait, M. Stefanon, S. Toft, J. R. Weaver, and A. Weibel The Astrophysical Journal Letters 940 (nov, 2022) L14.
[13] S. L. Finkelstein, M. B. Bagley, H. C. Ferguson, S. M. Wilkins, J. S. Kartaltepe, C. Papovich, L. Y. A. Yung, P. A. Haro, P. Behroozi, M. Dickinson, D. D. Kocevski, A. M. Koekemoer, R. L. Larson, A. L. Bail, A. M. Morales, P. G. Pérez-González, D. Burgarella, R. Davé, M. Hirschmann, R. S. Somerville, S. Wuyts, V. Bromm, C. M. Casey, A. Fontana, S. Fujimoto, J. P. Gardner, M. Giavalisco, A. Grazian, N. A. Grogin, N. P. Hathi, T. A. Hutchison, S. W. Jha, S. Jogee, L. J. Kewley, A. Kirkpatrick, A. S. Long, J. M. Lotz, L. Pentericci, J. D. R. Pierel, N. Pirzkal, S. Ravindranath, R. E. Ryan, J. R. Trump, G. Yang, R. Bhatawdekar, L. Bisigello, V. Buat, A. Calabrò, M. Castellano, N. J. Cleri, M. C. Cooper, D. Croton, E. Daddi, A. Dekel, D. Elbaz, M. Franco, E. Gawiser, B. W. Holwerda, M. Huertas-Company, A. E. Jaskot, G. C. K. Leung, R. A. Lucas, B. Mobasher, V. Pandya, S. Tacchella, B. J. Weiner, and J. A. Zavala The Astrophysical Journal Letters 946 (mar, 2023) L13.
[14] I. Labbé, P. van Dokkum, E. Nelson, R. Bezanson, K. A. Suess, J. Leja, G. Brammer, K. Whitaker, E. Mathews, M. Stefanon, and B. Wang Nature (London) 616 (Apr., 2023) 266–269, [arXiv:2207.12446].
[15] E. Curtis-Lake, S. Carniani, A. Cameron, S. Charlot, P. Jakobsen, R. Maiolino, A. Bunker, J. Witstok, R. Smit, J. Chevallard, C. Willott, P. Ferruit, S. Arribas, N. Bonaventura, M. Curti, F. D’Eugenio, M. Franx, G. Giardino, T. J. Looser, N. Lützgendorf, M. V. Maseda, T. Rawle, H.-W. Rix, B. R. del Pino, H. Übler, M. Sirianni, A. Dressler, E. Egami, D. J. Eisenstein, R. Endsley, K. Hainline, R. Hausen, B. D. Johnson, M. Rieke, B. Robertson, I. Shivaei, D. P. Stark, S. Tacchella, C. C. Williams, C. N. A. Willmer, R. Bhatawdekar, R. Bowler, K. Boyett, Z. Chen, A. de Graaff, J. M. Helton, R. E. Hviding, G. C. Jones, N. Kumari, J. Lyu, E. Nelson, M. Perna, L. Sandles, A. Saxena, K. A. Suess, F. Sun, M. W. Topping, I. E. B. Wallace, and L. Whitler, Spectroscopic confirmation of four metal-poor galaxies at z=10.3-13.2, 2023.
[16] B. E. Robertson et al. Nature Astron. 7 (2023), no. 5 611–621, [arXiv:2212.04480].
[17] B. W. Keller, F. Munshi, M. Trebitsch, and M. Tremmel Astrophys. J. Lett. 943 (2023), no. 2 L28, [arXiv:2212.12804].
[18] J. McCaffrey, S. Hardin, J. Wise, and J. Regan arXiv:2304.13755.
[19] M. Boylan-Kolchin Nature Astron. 7 (2023), no. 6 731–735, [arXiv:2208.01611].
[20] C. C. Lovell, I. Harrison, Y. Harikane, S. Tacchella, and S. M. Wilkins Mon. Not. Roy. Astron. Soc. 518 (2022), no. 2 2511–2520, [arXiv:2208.10479].
[21] M. Haslbauer, P. Kroupa, A. H. Zonoozi, and H. Haghi Astrophys. J. Lett. 939 (2022), no. 2 L31, [arXiv:2210.14915].
[22] B. Liu and V. Bromm Astrophys. J. Lett. 937 (2022), no. 2 L30, [arXiv:2208.13178].
[23] Y. Gong, B. Yue, Y. Cao, and X. Chen Astrophys. J. 947 (2023), no. 1 28, [arXiv:2209.13757].
[24] G. Hütsi, M. Raidal, J. Urrutia, V. Vaskonen, and H. Veermäe Phys. Rev. D 107 (2023), no. 4 043502, [arXiv:2211.02651].
[25] N. Menci, M. Castellano, P. Santini, E. Merlin, A. Fontana, and F. Shankar Astrophys. J. Lett. 938 (2022), no. 1 L5, [arXiv:2208.11471].
[26] M. Biagetti, G. Franciolini, and A. Riotto Astrophys. J. 944 (2023), no. 2 113, [arXiv:2210.04812].
[27] G.-W. Yuan, L. Lei, Y.-Z. Wang, B. Wang, Y.-Y. Wang, C. Chen, Z.-Q. Shen, Y.-F. Cai, and Y.-Z. Fan arXiv:2303.09391.
[28] C. Ilie, J. Paulin, and K. Freese arXiv:2304.01173.
[29] P. Parashari and R. Laha arXiv:2305.00999.
[30] H. Jiao, R. Brandenberger, and A. Refregier arXiv:2304.06429.
[31] S.-Y. Guo, M. Khlopov, X. Liu, L. Wu, Y. Wu, and B. Zhu arXiv:2306.17022.
[32] B.-Y. Su, N. Li, and L. Feng arXiv:2306.05364.
[33] S. Hawking Mon. Not. Roy. Astron. Soc. 152 (1971) 75.
[34] B. J. Carr and S. W. Hawking Mon. Not. Roy. Astron. Soc. 168 (1974) 399–415.
[35] P. Meszaros Astron. Astrophys. 37 (1974) 225–228.
[36] B. J. Carr Astrophys. J. 201 (1975) 1–19.
[37] G. F. Chapline Nature 253 (1975), no. 5489 251–252.
[38] B. Carr, F. Kuhnel, and M. Sandstad Phys. Rev. D 94 (2016), no. 8 083504, [arXiv:1607.06077].
[39] A. M. Green and B. J. Kavanagh J. Phys. G 48 (2021), no. 4 043001, [arXiv:2007.10722].
[40] M. Sasaki, T. Suyama, T. Tanaka, and S. Yokoyama Class. Quant. Grav. 35 (2018), no. 6 063001, [arXiv:1801.05235].
[41] B. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama Rept. Prog. Phys. 84 (2021), no. 11 116902, [arXiv:2002.12778].
[42] G. D. Quinlan and S. L. Shapiro, STAR CLUSTER COLLAPSE TO A SUPERMASSIVE BLACK HOLE: BINARIES AND GRAVITATIONAL RADIATION, 1987.
[43] H. Mouri and Y. Taniguchi Astrophys. J. Lett. 566 (2002) L17–L20, [astro-ph/0201102].
[44] S. Bird, I. Cholis, J. B. Muñoz, Y. Ali-Haïmoud, M. Kamionkowski, E. D. Kovetz, A. Raccanelli, and A. G. Riess Phys. Rev. Lett. 116 (2016), no. 20 201301, [arXiv:1603.00464].
[45] M. Sasaki, T. Suyama, T. Tanaka, and S. Yokoyama Phys. Rev. Lett. 117 (2016), no. 6 061101, [arXiv:1603.08338]. [Erratum: Phys.Rev.Lett. 121, 059901 (2018)].
[46] A. Kashlinsky Astrophys. J. Lett. 823 (2016), no. 2 L25, [arXiv:1605.04023].
[47] Y. Ali-Haïmoud and M. Kamionkowski Phys. Rev. D 95 (2017), no. 4 043534, [arXiv:1612.05644].
[48] M. Raidal, V. Vaskonen, and H. Veermäe JCAP 09 (2017) 037, [arXiv:1707.01480].
[49] M. Raidal, C. Spethmann, V. Vaskonen, and H. Veermäe JCAP 02 (2019) 018, [arXiv:1812.01930].
[50] LIGO Scientific, Virgo Collaboration, B. P. Abbott et al. Phys. Rev. X 6 (2016), no. 4 041015, [arXiv:1606.04856]. [Erratum: Phys.Rev.X 8, 039903 (2018)].
[51] LIGO Scientific, Virgo Collaboration, R. Abbott et al. Phys. Rev. X 11 (2021) 021053, [arXiv:2010.14527].
[52] A. Hall, A. D. Gow, and C. T. Byrnes Phys. Rev. D 102 (2020) 123524, [arXiv:2008.13704].
[53] K. Jedamzik JCAP 09 (2020) 022, [arXiv:2006.11172].
[54] V. De Luca, G. Franciolini, P. Pani, and A. Riotto JCAP 06 (2020) 044, [arXiv:2005.05641].
[55] G. Franciolini, V. Baibhav, V. De Luca, K. K. Y. Ng, K. W. K. Wong, E. Berti, P. Pani, A. Riotto, and S. Vitale Phys. Rev. D 105 (2022), no. 8 083526, [arXiv:2105.03349].
[56] A. Romero-Rodriguez, M. Martinez, O. Pujolàs, M. Sakellariadou, and V. Vaskonen Phys. Rev. Lett. 128 (2022), no. 5 051301, [arXiv:2107.11660].
[57] L. Ferrarese and H. Ford Space Sci. Rev. 116 (2005) 523–624, [astro-ph/0411247].
[58] K. Gultekin et al. Astrophys. J. 698 (2009) 198–221, [arXiv:0903.4897].
[59] J. Kormendy and L. C. Ho Ann. Rev. Astron. Astrophys. 51 (2013) 511–653, [arXiv:1304.7762].
[60] R. Bean and J. Magueijo Phys. Rev. D 66 (2002) 063505, [astro-ph/0204486].
[61] M. Kawasaki, A. Kusenko, and T. T. Yanagida Phys. Lett. B 711 (2012) 1–5, [arXiv:1202.3848].
[62] S. Clesse and J. García-Bellido Phys. Dark Univ. 15 (2017) 142–147, [arXiv:1603.05234].
[63] S. Clesse and J. García-Bellido Phys. Dark Univ. 22 (2018) 137–146, [arXiv:1711.10458].
[64] P. D. Serpico, V. Poulin, D. Inman, and K. Kohri Phys. Rev. Res. 2 (2020), no. 2 023204, [arXiv:2002.10771].
[65] T. Nakama, T. Suyama, and J. Yokoyama Phys. Rev. D 94 (2016), no. 10 103522, [arXiv:1609.02245].
[66] T. Nakama, B. Carr, and J. Silk Phys. Rev. D 97 (2018), no. 4 043525, [arXiv:1710.06945].
[67] B. Carr and J. Silk Mon. Not. Roy. Astron. Soc. 478 (2018), no. 3 3756–3775, [arXiv:1801.00672].
[68] D. Inman and Y. Ali-Haïmoud Phys. Rev. D 100 (2019), no. 8 083528, [arXiv:1907.08129].
[69] B. Carr, F. Kuhnel, and L. Visinelli Mon. Not. Roy. Astron. Soc. 501 (2021), no. 2 2029–2043, [arXiv:2008.08077].
[70] V. Atal, A. Sanglas, and N. Triantafyllou JCAP 06 (2021) 022, [arXiv:2012.14721].
[71] P. F. Depta, K. Schmidt-Hoberg, and C. Tasillo arXiv:2306.17836.
[72] P. A. Oesch, R. J. Bouwens, G. D. Illingworth, I. Labbé, and M. Stefanon The Astrophysical Journal 855 (Mar., 2018) 105, [arXiv:1710.11131].
[73] R. J. Bouwens, P. A. Oesch, M. Stefanon, G. Illingworth, I. Labbé, N. Reddy, H. Atek, M. Montes, R. Naidu, T. Nanayakkara, E. Nelson, and S. Wilkins The Astrophysical Journal 162 (Aug., 2021) 47, [arXiv:2102.07775].
[74] R. J. Bouwens et al. Astrophys. J. 803 (2015), no. 1 34, [arXiv:1403.4295].
[75] S. L. Finkelstein, J. Ryan, Russell E., C. Papovich, M. Dickinson, M. Song, R. S. Somerville, H. C. Ferguson, B. Salmon, M. Giavalisco, A. M. Koekemoer, M. L. N. Ashby, P. Behroozi, M. Castellano, J. S. Dunlop, S. M. Faber, G. G. Fazio, A. Fontana, N. A. Grogin, N. Hathi, J. Jaacks, D. D. Kocevski, R. Livermore, R. J. McLure, E. Merlin, B. Mobasher, J. A. Newman, M. Rafelski, V. Tilvi, and S. P. Willner The Astrophysical Journal 810 (Sept., 2015) 71, [arXiv:1410.5439].
[76] H. Atek et al. Astrophys. J. 814 (2015), no. 1 69, [arXiv:1509.06764].
[77] R. C. Livermore, S. L. Finkelstein, and J. M. Lotz Astrophys. J. 835 (2017), no. 2 113, [arXiv:1604.06799].
[78] R. J. Bouwens, P. A. Oesch, G. D. Illingworth, R. S. Ellis, and M. Stefanon The Astrophysical Journal 843 (July, 2017) 129, [arXiv:1610.00283].
[79] V. Mehta, C. Scarlata, M. Rafelski, T. Gburek, H. I. Teplitz, A. Alavi, M. Boylan-Kolchin, S. Finkelstein, J. P. Gardner, N. Grogin, A. Koekemoer, P. Kurczynski, B. Siana, A. Codoreanu, D. F. de Mello, K.-S. Lee, and E. Soto The Astrophysical Journal 838 (Mar., 2017) 29, [arXiv:1702.06953].
[80] M. Ishigaki, R. Kawamata, M. Ouchi, M. Oguri, K. Shimasaku, and Y. Ono The Astrophysical Journal 854 (Feb., 2018) 73, [arXiv:1702.04867].
[81] H. Atek, J. Richard, J.-P. Kneib, and D. Schaerer Mon. Not. Roy. Astron. Soc. 479 (2018), no. 4 5184–5195, [arXiv:1803.09747].
[82] S. Rojas-Ruiz, S. L. Finkelstein, M. B. Bagley, M. Stevans, K. D. Finkelstein, R. Larson, M. Mechtley, and J. Diekmann The Astrophysical Journal 891 (Mar., 2020) 146, [arXiv:2002.06209].
[83] T. Nakamura, M. Sasaki, T. Tanaka, and K. S. Thorne Astrophys. J. Lett. 487 (1997) L139–L142, [astro-ph/9708060].
[84] P. C. Peters and J. Mathews Phys. Rev. 131 (1963) 435–439.
[85] LIGO Scientific, Virgo Collaboration, B. P. Abbott et al. Phys. Rev. Lett. 116 (2016), no. 13 131102, [arXiv:1602.03847].
[86] LIGO Scientific Collaboration and Virgo Collaboration Collaboration Phys. Rev. Lett. 116 (Mar, 2016) 131102.
[87] L. Z. Kelley, L. Blecha, L. Hernquist, A. Sesana, and S. R. Taylor Mon. Not. Roy. Astron. Soc. 471 (2017), no. 4 4508–4526, [arXiv:1702.02180].
[88] S. Burke-Spolaor et al. Astron. Astrophys. Rev. 27 (2019), no. 1 5, [arXiv:1811.08826].
[89] P. A. Rosado, A. Sesana, and J. Gair Mon. Not. Roy. Astron. Soc. 451 (2015), no. 3 2417–2433, [arXiv:1503.04803].
[90] S. R. Taylor, R. van Haasteren, and A. Sesana Phys. Rev. D 102 (2020), no. 8 084039, [arXiv:2006.04810].
[91] A. Sesana, A. Vecchio, and C. N. Colacino Mon. Not. Roy. Astron. Soc. 390 (2008) 192, [arXiv:0804.4476].
[92] T. Papanikolaou, V. Vennin, and D. Langlois JCAP 03 (2021) 053, [arXiv:2010.11573].
[93] A. Lewis, A. Challinor, and A. Lasenby Astrophys. J. 538 (2000) 473–476, [astro-ph/9911177].
[94] Planck Collaboration, N. Aghanim et al. Astron. Astrophys. 641 (2020) A6, [arXiv:1807.06209]. [Erratum: Astron.Astrophys. 652, C4 (2021)].
[95] W. H. Press and P. Schechter Astrophys. J. 187 (1974) 425–438.
[96] V. De Luca, V. Desjacques, G. Franciolini, and A. Riotto JCAP 11 (2020) 028, [arXiv:2009.04731].
[97] B. Liu, S. Zhang, and V. Bromm Mon. Not. Roy. Astron. Soc. 514 (2022), no. 2 2376–2396, [arXiv:2204.06330].
[98] R. K. Sheth and G. Tormen Mon. Not. Roy. Astron. Soc. 308 (1999) 119, [astro-ph/9901122].
[99] R. K. Sheth, H. J. Mo, and G. Tormen Mon. Not. Roy. Astron. Soc. 323 (2001) 1, [astro-ph/9907024].
[100] B. J. Carr and M. Sakellariadou Astrophys. J. 516 (1999) 195–220.
[101] J. Chluba, A. L. Erickcek, and I. Ben-Dayan Astrophys. J. 758 (2012) 76, [arXiv:1203.2681].
[102] D. J. Fixsen, E. S. Cheng, J. M. Gales, J. C. Mather, R. A. Shafer, and E. L. Wright Astrophys. J. 473 (1996) 576, [astro-ph/9605054].
[103] F. Bianchini and G. Fabbian Phys. Rev. D 106 (2022), no. 6 063527, [arXiv:2206.02762].
[104] S. W. Hawking Phys. Lett. B 231 (1989) 237–239.
[105] J. Garriga and A. Vilenkin Phys. Rev. D 47 (1993) 3265–3274, [hep-ph/9208212].
[106] J. Garriga, A. Vilenkin, and J. Zhang JCAP 02 (2016) 064, [arXiv:1512.01819].
[107] H. Deng, J. Garriga, and A. Vilenkin JCAP 04 (2017) 050, [arXiv:1612.03753].
[108] H. Deng and A. Vilenkin JCAP 12 (2017) 044, [arXiv:1710.02865].
[109] H. Deng, A. Vilenkin, and M. Yamada JCAP 07 (2018) 059, [arXiv:1804.10059].
[110] M. Kawasaki and K. Murai Phys. Rev. D 100 (2019), no. 10 103521, [arXiv:1907.02273].
[111] E. Cotner, A. Kusenko, M. Sasaki, and V. Takhistov JCAP 10 (2019) 077, [arXiv:1907.10613].
[112] N. Kitajima and F. Takahashi JCAP 11 (2020) 060, [arXiv:2006.13137].
[113] K. Kawana and K.-P. Xie Phys. Lett. B 824 (2022) 136791, [arXiv:2106.00111].
[114] K. Kasai, M. Kawasaki, and K. Murai JCAP 10 (2022) 048, [arXiv:2205.10148].
[115] J. H. Chang, D. Egana-Ugrinovic, R. Essig, and C. Kouvaris JCAP 03 (2019) 036, [arXiv:1812.07000].
[116] P. Lu, K. Kawana, and A. Kusenko Phys. Rev. D 107 (2023), no. 10 103037, [arXiv:2210.16462].
[117] G. Domènech, D. Inman, A. Kusenko, and M. Sasaki arXiv:2304.13053.
[118] G. Sun and S. R. Furlanetto Monthly Notices of the Royal Astronomical Society 460 (July, 2016) 417–433, [arXiv:1512.06219].
[119] R. P. Naidu, S. Tacchella, C. A. Mason, S. Bose, P. A. Oesch, and C. Conroy arXiv:1907.13130.
[120] C. Mason, M. Trenti, and T. Treu Astrophys. J. 813 (2015), no. 1 21, [arXiv:1508.01204].
[121] M. Sahlén and E. Zackrisson arXiv:2105.05098.
[122] N. Sabti, J. B. Muñoz, and D. Blas Astrophys. J. Lett. 928 (2022), no. 2 L20, [arXiv:2110.13161].
[123] N. Sabti, J. B. Muñoz, and D. Blas Phys. Rev. D 105 (2022), no. 4 043518, [arXiv:2110.13168].
[124] N. Menci et al. Astrophys. J. 900 (2020), no. 2 108, [arXiv:2007.12453].
[125] B. Bozek, D. J. E. Marsh, J. Silk, and R. F. G. Wyse Mon. Not. Roy. Astron. Soc. 450 (2015), no. 1 209–222, [arXiv:1409.3544].
[126] C. Schultz, J. Oñorbe, K. N. Abazajian, and J. S. Bullock Mon. Not. Roy. Astron. Soc. 442 (2014), no. 2 1597–1609, [arXiv:1401.3769].
[127] P. Dayal, A. Mesinger, and F. Pacucci Astrophys. J. 806 (2015), no. 1 67, [arXiv:1408.1102].
[128] P. S. Corasaniti, S. Agarwal, D. J. E. Marsh, and S. Das Phys. Rev. D 95 (2017), no. 8 083512, [arXiv:1611.05892].
[129] N. Menci, A. Merle, M. Totzauer, A. Schneider, A. Grazian, M. Castellano, and N. G. Sanchez Astrophys. J. 836 (2017), no. 1 61, [arXiv:1701.01339].
[130] N. Menci, A. Grazian, A. Lamastra, F. Calura, M. Castellano, and P. Santini Astrophys. J. 854 (2018), no. 1 1, [arXiv:1801.03697].
[131] A. Rudakovskyi, A. Mesinger, D. Savchenko, and N. Gillet Mon. Not. Roy. Astron. Soc. 507 (2021), no. 2 3046–3056, [arXiv:2104.04481].
[132] S. Yoshiura, M. Oguri, K. Takahashi, and T. Takahashi Phys. Rev. D 102 (2020), no. 8 083515, [arXiv:2007.14695].
[133] J. Chevallard, J. Silk, T. Nishimichi, M. Habouzit, G. A. Mamon, and S. Peirani Mon. Not. Roy. Astron. Soc. 446 (2015) 3235–3252, [arXiv:1410.7768].
[134] N. Sabti, J. B. Muñoz, and D. Blas JCAP 01 (2021) 010, [arXiv:2009.01245].
[135] N. Sabti, J. B. Muñoz, and M. Kamionkowski arXiv:2305.07049.
[136] A. Mitridate, D. Wright, R. von Eckardstein, T. Schröder, J. Nay, K. Olum, K. Schmitz, and T. Trickle arXiv:2306.16377.
[137] J. Silk Astrophys. J. 211 (1977) 638–648.
[138] J. R. Bond, S. Cole, G. Efstathiou, and N. Kaiser Astrophys. J. 379 (1991) 440.
[139] J. M. Bardeen, J. R. Bond, N. Kaiser, and A. S. Szalay Astrophys. J. 304 (1986) 15–61.
[140] A. Schneider, S. K. Giri, and J. Mirocha Phys. Rev. D 103 (2021), no. 8 083025, [arXiv:2011.12308].
[141] J. Silk Astrophys. J. 151 (1968) 459–471.
[142] T. Harada, C.-M. Yoo, K. Kohri, and K.-I. Nakao Phys. Rev. D 96 (2017), no. 8 083517, [arXiv:1707.03595]. [Erratum: Phys.Rev.D 99, 069904 (2019)].
[143] B. J. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama Phys. Rev. D 94 (2016), no. 4 044029, [arXiv:1604.05349].
[144] T. Nakama, J. Silk, and M. Kamionkowski Phys. Rev. D 95 (2017), no. 4 043511, [arXiv:1612.06264].
[145] B. J. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama Phys. Rev. D 81 (2010) 104019, [arXiv:0912.5297].
[146] C. T. Byrnes, E. J. Copeland, and A. M. Green Phys. Rev. D 86 (2012) 043512, [arXiv:1206.4188].
[147] S. Young and C. T. Byrnes JCAP 08 (2013) 052, [arXiv:1307.4995].
[148] G. Franciolini, A. Kehagias, S. Matarrese, and A. Riotto JCAP 03 (2018) 016, [arXiv:1801.09415].
[149] V. Atal, J. Garriga, and A. Marcos-Caballero JCAP 09 (2019) 073, [arXiv:1905.13202].
[150] M. W. Davies, P. Carrilho, and D. J. Mulryne JCAP 06 (2022), no. 06 019, [arXiv:2110.08189].
[151] Y.-F. Cai, X.-H. Ma, M. Sasaki, D.-G. Wang, and Z. Zhou JCAP 12 (2022) 034, [arXiv:2207.11910].
[152] G. Ferrante, G. Franciolini, A. Iovino, Junior., and A. Urbano Phys. Rev. D 107 (2023), no. 4 043520, [arXiv:2211.01728].
[153] S. Pi and M. Sasaki Phys. Rev. Lett. 131 (2023), no. 1 011002, [arXiv:2211.13932].
[154] L. Iacconi and D. J. Mulryne arXiv:2304.14260.
[155] H.-L. Huang, Y. Cai, J.-Q. Jiang, J. Zhang, and Y.-S. Piao arXiv:2306.17577.
[156] D. Hooper, A. Ireland, G. Krnjaic, and A. Stebbins arXiv:2308.00756.
[157] C.-M. Yoo, J.-O. Gong, and S. Yokoyama JCAP 09 (2019) 033, [arXiv:1906.06790].
[158] A. Kehagias, I. Musco, and A. Riotto JCAP 12 (2019) 029, [arXiv:1906.07135].
[159] V. Atal, J. Cid, A. Escrivà, and J. Garriga JCAP 05 (2020) 022, [arXiv:1908.11357].
[160] A. Escrivà, F. Kuhnel, and Y. Tada arXiv:2211.05767.
[161] K. Ioka, T. Chiba, T. Tanaka, and T. Nakamura Phys. Rev. D 58 (1998) 063003, [astro-ph/9807018].
[162] H. Deng JCAP 03 (2022), no. 03 037, [arXiv:2110.02460].
[163] M. C. Begelman, R. D. Blandford, and M. J. Rees Nature 287 (1980) 307–309.
[164] F. Dosopoulou and F. Antonini Astrophys. J. 840 (2017), no. 1 31, [arXiv:1611.06573].
[165] J. Ellis, M. Fairbairn, G. Hütsi, M. Raidal, J. Urrutia, V. Vaskonen, and H. Veermäe Astron. Astrophys. 676 (2023) A38, [arXiv:2301.13854].
[166] M. Volonteri Astron. Astrophys. Rev. 18 (2010) 279–315, [arXiv:1003.4404].
[167] N. J. Outmezguine, O. Slone, W. Tangarife, L. Ubaldi, and T. Volansky JHEP 11 (2018) 005, [arXiv:1807.04750].
[168] P. C. Peters Phys. Rev. 136 (Nov, 1964) B1224–B1232.
[169] P. Ajith et al. Phys. Rev. D 77 (2008) 104017, [arXiv:0710.2335]. [Erratum: Phys.Rev.D 79, 129901 (2009)].
[170] P. Ajith et al. Phys. Rev. Lett. 106 (2011) 241101, [arXiv:0909.2867].
[171] S. R. Taylor arXiv:2105.13270.
[172] M. Maggiore, Gravitational Waves. Vol. 2: Astrophysics and Cosmology. Oxford University Press, 3, 2018.
[173] J. A. Ellis, M. Vallisneri, S. R. Taylor, and P. T. Baker, “Enterprise: Enhanced numerical toolbox enabling a robust pulsar inference suite.” Zenodo, Sept., 2020.
[174] S. R. Taylor, P. T. Baker, J. S. Hazboun, J. Simon, and S. J. Vigeland, enterprise_extensions, 2021. v2.3.3.
[175] J. Ellis and R. van Haasteren, jellis18/ptmcmcsampler: Official release, Oct., 2017.
[176] A. Lewis arXiv:1910.13970.
[177] NANOGrav Collaboration, G. Agazie et al. Astrophys. J. Lett. 952 (2023), no. 2 L37, [arXiv:2306.16220].
[178] EPTA Collaboration, J. Antoniadis et al. arXiv:2306.16227.
[179] E. S. Phinney astro-ph/0108028.
[180] B. A. Reid, W. J. Percival, D. J. Eisenstein, L. Verde, D. N. Spergel, R. A. Skibba, N. A. Bahcall, T. Budavari, J. A. Frieman, M. Fukugita, J. R. Gott, J. E. Gunn, Ž. Ivezić, G. R. Knapp, R. G. Kron, R. H. Lupton, T. A. McKay, A. Meiksin, R. C. Nichol, A. C. Pope, D. J. Schlegel, D. P. Schneider, C. Stoughton, M. A. Strauss, A. S. Szalay, M. Tegmark, M. S. Vogeley, D. H. Weinberg, D. G. York, and I. Zehavi MNRAS 404 (May, 2010) 60–85, [arXiv:0907.1659].
[181] S. Chabanier, M. Millea, and N. Palanque-Delabrouille Mon. Not. Roy. Astron. Soc. 489 (2019), no. 2 2247–2253, [arXiv:1905.08103].
[182] R. Murgia, G. Scelfo, M. Viel, and A. Raccanelli Phys. Rev. Lett. 123 (2019), no. 7 071102, [arXiv:1903.10509].
[183] D. Blas and A. C. Jenkins Phys. Rev. Lett. 128 (2022), no. 10 101103, [arXiv:2107.04601].
[184] M. Levi, C. Bebek, T. Beers, R. Blum, R. Cahn, D. Eisenstein, B. Flaugher, K. Honscheid, R. Kron, O. Lahav, P. McDonald, N. Roe, D. Schlegel, and representing the DESI collaboration arXiv e-prints (Aug, 2013) arXiv:1308.0847, [arXiv:1308.0847].
[185] EUCLID Collaboration Collaboration, R. Laureijs et al. arXiv:1110.3193.
[186] LSST Dark Energy Science Collaboration, A. Abate et al. arXiv:1211.0310.