Clustering of JWST LRDs. Implications for PBHs and other early-forming seeds

abstract · abstract

We use the observed abundance and spatial clustering of JWST little red dots (LRDs) over z = 4--8 to confront two formation scenarios: the standard picture in which LRDs form in virialized dark-matter halos, and one in which they are seeded by a dynamically subdominant population of primordial black holes (PBHs). The halo formation scenario accounts for both observables naturally and without fine-tuning. In the PBH scenario, matching the observed LRD abundance requires a negligibly small PBH mass fraction, and a power-spectrum analysis shows that such a Poisson isocurvature component predicts a clustering signal far below what is observed. The same conclusion holds for any population of dynamically subdominant early-forming seeds, unless their initial clustering is finely tuned or their activity is tied to massive host halos. We conclude that the spatial clustering of JWST LRDs strongly favours the standard halo formation picture and disfavours PBHs, or similar early-forming seeds, as a randomly activated LRD population; PBHs remain viable as seeds only if their activity is governed by massive host halos, in which case clustering alone cannot distinguish them from astrophysical seeds.

S3: III Purely adiabatic fluctuations. Formation at virialized density peaks / block 1 · html

The results of this calculation are shown in Fig. 1 . The two families of curves trace, respectively, the locus of constant lower mass cutoff (varying redshift) and constant redshift (varying M low M_{\rm low} ) in the ( n ¯ , b ¯ ) (\bar{n},\bar{b}) plane. The grey shaded bands mark the observational constraints on number density and bias described in Sect. II .

S3: III Purely adiabatic fluctuations. Formation at virialized density peaks / block 1 · html

Taken together, these results show that the simple halo model naturally accommodates both the observed number density and clustering strength of LRDs over z = 4 z=4 –8, without fine-tuning. A more sophisticated halo occupation distribution (HOD) approach [ 53 ] would allow additional flexibility (scatter in the occupation number, a central-satellite decomposition, etc.), but the current observational precision does not motivate the additional parameters. This analysis concerns the probable formation sites of LRDs; the detailed physical formation mechanisms of the black holes and the conditions giving rise to AGN activity remain active areas of research (see Sect. I and references therein).

S5: V Discussion and Conclusions / block 1 · html

We have used the observed comoving number density and spatial clustering of JWST little red dots over z = 4 z=4 –8 to compare two broad scenarios for their origin. In the standard adiabatic scenario, LRDs form at virialized dark-matter density peaks. In the alternative isocurvature scenario, they are associated with a dynamically subdominant population of primordial black holes. Our main findings are as follows.

S5: V Discussion and Conclusions / block 1 · html

The halo model provides a natural and self-consistent account of both observables. A simple model in which every dark-matter halo above a mass threshold M low M_{\rm low} hosts an LRD with duty cycle f on f_{\rm on} reproduces the observed bias range b ∼ 3 b\sim 3 –8 across z = 4 z=4 –8 for M low ∼ 10 10.5 M_{\rm low}\sim 10^{10.5} – 10 12 ​ h − 1 ​ M ⊙ 10^{12}~h^{-1}{\rm M}_{\odot} , while the observed number density n ¯ LRD ∼ 10 − 5 \bar{n}_{\rm LRD}\sim 10^{-5} – 10 − 4 ​ h 3 ​ Mpc − 3 10^{-4}~{h^{3}~{\rm Mpc}}^{-3} is then matched by an appropriate duty cycle, without fine-tuning. Bias values at the lower end of the observed range can also be accommodated, in the limiting case of passive evolution, by sources formed at earlier epochs. Extending the model to a full HOD framework [ 53 ] is possible but not motivated by current data.

S5: V Discussion and Conclusions / block 1 · html

The PBH scenario faces a fundamental challenge. Matching the observed LRD number density requires f pbh ∼ 10 − 13 f_{\rm pbh}\sim 10^{-13} – 10 − 10 10^{-10} across M pbh = 10 3 M_{\rm pbh}=10^{3} – 10 5 ​ h − 1 ​ M ⊙ 10^{5}~h^{-1}{\rm M}_{\odot} , making PBHs entirely negligible for structure formation. For such a dynamically subdominant Poisson distribution of seeds, the power spectrum analysis of Sect. IV shows directly that b pbh ≃ 1 b_{\rm pbh}\simeq 1 at z ∼ 5 z\sim 5 — far below the observed LRD clustering signal. More broadly, the same conclusion applies to any population of very early-forming seeds with initially unclustered (Poisson) statistics that is dynamically unimportant. If instead such seeds were initially clustered, the Fry equation ( 4 ) governs the subsequent bias evolution: being dynamically subdominant, they follow the adiabatic CDM potential and their bias decays toward unity, so that without fine-tuning of the initial clustering, b ∼ 1 b\sim 1 is again the generic prediction at z ∼ 5 z\sim 5 .

S5: V Discussion and Conclusions / block 1 · html

Non-Gaussian scenarios for PBH formation could in principle provide such initial clustering. Strongly non-Gaussian initial perturbations could both bypass the μ \mu -distortion bounds [ 110 , 111 , 112 ] and produce enhanced initial PBH clustering [ 99 , 100 ] . However, generating b ∼ 5 b\sim 5 –8 at z ∼ 5 z\sim 5 from PBH clustering alone would require a level of primordial non-Gaussianity that is difficult to motivate without fine-tuning.

S5: V Discussion and Conclusions / block 1 · html

Another possible way to reconcile PBH seeds with the observed clustering is an environment-dependent activation probability, f on = f on ​ ( M h ) f_{\rm on}=f_{\rm on}(M_{\rm h}) , such that only PBHs residing in sufficiently massive host halos accrete efficiently enough to appear as LRDs. Such a coupling is physically plausible: the accretion rate is set by the local gas density and temperature, which are strongly enhanced in the centres of massive halos, where dynamical friction also brings the most massive black holes. Conversely, in less massive halos supernova feedback efficiently expels gas from the central regions and stunts black-hole growth, so that efficient accretion sets in only above a critical host-halo mass [ 114 , 115 , 116 , 117 ] . The accretion-driven growth of intermediate-mass PBHs considered in Ref. [ 92 ] relies on an efficient gas supply, which naturally points to such massive host environments. Since PBHs trace the matter distribution, the fraction residing in halos above M low ∼ 10 10.5 M_{\rm low}\sim 10^{10.5} – 10 11.5 ​ h − 1 ​ M ⊙ 10^{11.5}~h^{-1}{\rm M}_{\odot} at z ∼ 5 z\sim 5 , and hence the effective f on f_{\rm on} , is only ∼ 10 − 3 \sim 10^{-3} – 10 − 2 10^{-2} , so that by Eq. ( 5 ) the required PBH abundance increases correspondingly to f pbh ∼ 10 − 11 f_{\rm pbh}\sim 10^{-11} – 10 − 7 10^{-7} , still dynamically negligible. In this case, however, the clustering of LRDs is determined entirely by their host halos, and the scenario becomes observationally indistinguishable from the halo picture of Sect. III : the PBH merely provides the seed, while both the location and the activity of the LRD are governed by the halo. The clustering measurements therefore rule out PBHs (or similar early-forming seeds) as a randomly activated, unbiased LRD population, but cannot by themselves exclude a primordial origin of black holes that subsequently grow in massive halos. Discriminating between primordial and astrophysical seeds in this case requires other observables, such as the seed mass function at z ≳ 10 z\gtrsim 10 , the black-hole occupation fraction of low-mass halos, or gravitational-wave signatures. Gravitational-wave observati

S5: V Discussion and Conclusions / block 2 · html

ons are particularly promising in this respect. The LISA mission will detect mergers of massive black holes with masses ∼ 10 4 \sim 10^{4} – 10 7 ​ M ⊙ 10^{7}\,M_{\odot} out to z ∼ 20 z\sim 20 , and the inferred mass and redshift distributions of these events are sensitive to the seeding channel [ 118 , 119 , 120 , 121 , 122 ] . Proposed mid-band atom interferometers, such as AEDGE and AION, would extend this sensitivity to intermediate-mass black-hole binaries in the ∼ 10 3 \sim 10^{3} – 10 5 ​ M ⊙ 10^{5}\,M_{\odot} range directly relevant for LRD seeds [ 123 , 124 , 125 ] , while the nanohertz gravitational-wave background reported by pulsar timing arrays provides complementary constraints on the growth history of the most massive black holes, including those observed by JWST [ 126 , 127 , 128 ] . Third-generation ground-based detectors, such as the Einstein Telescope, could in addition detect black-hole mergers at z ≳ 30 z\gtrsim 30 , before the onset of significant star formation, which would provide a strong indication of a primordial origin for at least part of the black-hole population [ 129 , 130 , 131 , 132 ] ; more generally, future detectors can probe PBH binaries over a broad mass range, building on detailed modelling of their formation and merger rates [ 133 , 134 , 135 , 136 , 137 ] .

S5: V Discussion and Conclusions / block 1 · html

Observational evidence for a close connection between LRDs and overdense large-scale environments provides additional qualitative support for the halo picture. LRDs appear preferentially associated with bright companion galaxies and enhanced galaxy overdensities [ 45 ] , consistent with the expectation that black-hole growth is facilitated by the dense gas reservoirs and elevated Lyman-Werner radiation backgrounds characterising biased nodes of the cosmic web. Strikingly, the companions of the brightest LRDs are found to be both more numerous and bluer than those of control galaxies at > 3 ​ σ >3\sigma significance [ 138 ] , consistent with UV-bright neighbours supplying the Lyman-Werner flux required for direct collapse. This environmental correlation supports the direct-collapse black-hole scenario, in which nearby UV sources suppress molecular hydrogen cooling in metal-poor atomic-cooling halos, enabling the formation of massive seeds without fragmentation.

S3: 3 Formation of UHZ1-type Systems / block 1 · html

We can now tie this discussion to Fig. 2 which indicates the halo mass range that can collapse at any z z , for various f PBH ​ M PBH f_{\rm PBH}M_{\rm PBH} and ( 1 , 3 ) ​ σ (1,3)\sigma density excursions. The figure shows that at the 3 σ \sigma level, halos of M H ≃ ( 1.5 × 10 5 , 3 × 10 5 , 7 × 10 5 , 2.5 × 10 6 , 10 7 ) ​ M ⊙ M_{\rm H}\simeq(1.5\times 10^{5},3\times 10^{5},7\times 10^{5},2.5\times 10^{6},10^{7})M_{\odot} , can collapse at z = 100 , 70 , 50 , 30 , 20 z=100,70,50,30,20 , if f PBH ​ M PBH = 20 ​ M ⊙ f_{\rm PBH}M_{\rm PBH}=20M_{\odot} , decreasing to M H ≃ ( ≤ 10 5 , 1.7 × 10 5 , 4 × 10 5 , 1.3 × 10 6 , 5 × 10 6 ) M ⊙ M_{\rm H}\simeq(\leq 10^{5},1.7\times 10^{5},4\times 10^{5},1.3\times 10^{6},5\times 10^{6})M_{\odot} at z = 100 , 70 , 50 , 30 , 20 z=100,70,50,30,20 , if f PBH ​ M PBH = 10 ​ M ⊙ f_{\rm PBH}M_{\rm PBH}=10M_{\odot} . The Press-Schechter ( 66 ) formalism suggests that a sufficient fraction of such halos would exist at ( 3 − 5 ) ​ σ (3-5)\sigma levels to produce - for suitable PBH mass functions - at least a few central objects of 10 7 − 10 8 ​ M ⊙ 10^{7}-10^{8}M_{\odot} by z = 10 z=10 , as observed for UHZ1. E.g., at 3 ​ σ 3\sigma we can expect ≃ 0.3 % \simeq 0.3\% of the density field to have the required parameters. Currently, theoretically-motivated PBH mass function proposals contain 10 to 20% by mass in PBHs of 10 ​ M ⊙ 10M_{\odot} or higher ( 12 ) , i.e., those subject to dynamical friction, and 10-20% of them will have the required low angular momentum, λ ≤ 0.1 − 0.2 \lambda\leq 0.1-0.2 . Dynamical friction would therefore drive ∼ 1 − 2 % \sim 1-2\% of the total mass of the PBH component to the central 1 pc in ( 0.5 − 1 ) × t cosm (0.5-1)\times t_{\rm cosm} after halo collapse. These clusters of PBHs at the centers of low mass halos will merge into larger units as hierarchical clustering proceeds. At z ≃ 15 z\simeq 15 , halos of 10 9 ​ M ⊙ 10^{9}M_{\odot} , corresponding to peaks of 3 − 4 ​ σ 3-4\sigma , could have collapsed and given rise to concentrations of up to ∼ 10 7 ​ M ⊙ \sim 10^{7}M_{\odot} in the central parsec by z = 10 z=10 , depending on the fraction, f PBH f_{\rm PBH} , of DM consisting of PBHs. Baryonic accret

S3: 3 Formation of UHZ1-type Systems / block 1 · html

Figure 3: Dynamical friction collapse time t df t_{\rm df} from solving numerically eqs. 6 , in units of the crossing/orbital time, R H / V 0 R_{\rm H}/V_{0} as a function of λ \lambda and 𝒞 {\cal C} . These results correspond to m ∙ ​ ln ⁡ Λ / M H = 0.01 m_{\bullet}\ln\Lambda/M_{\rm H}=0.01 and can be rescaled to other values since t df ∝ ( m ∙ ​ ln ⁡ Λ / M H ) − 1 t_{\rm df}\propto(m_{\bullet}\ln\Lambda/M_{\rm H})^{-1} . Dashed horizontal lines mark ( 1 , 2 , 3 ) × t cosm ​ ( z ) (1,2,3)\times t_{\rm cosm}(z) , whereas the dotted line indicates the time for an individual halo to remain intact from z = 50 z=50 to z = 10 z=10 . Left: t df t_{\rm df} vs. λ \lambda ; NFW halo profiles with 𝒞 = 10 , 20 , 30 , 40 , 100 {\cal C}=10,20,30,40,100 , are shown with blue, orange, green, red, and purple solid lines, and isothermal halo profiles with the same concentrations are shown with corresponding dash-dotted lines. The isothermal, non-truncated distribution ℱ ⁡ ( E ) {\cal F}(E) was adopted. Right: Same vertical axis vs. the concentration parameter, 𝒞 {\cal C} . Solid lines represent NFW profiles with λ = 0.05 , 0.1 , .0.15 , 0.2 , 0.3 \lambda=0.05,0.1,.0.15,0.2,0.3 in blue, orange, green, red, and purple, respectively, and dash-dotted lines correspond to isothermal profiles.

article / block 1 · html

in terms of the net metal yield y ∼ 𝒪 ⁡ ( 10 − 2 ) y\sim\mathcal{O}(10^{-2}) [ 130 ] . As expected, a larger fragmentation parameter would result in a larger metal enrichment of the environment. Notice that our modeling neglects the contribution of heavy elements, associated with the usual structure formation history at later redshifts [ 130 ] . In this sense, our results have to be intended as an upper bound on the fragmentation fraction f frag f_{\text{\tiny frag}} of QSO1-like LRDs (see gray region in Fig. 3 ).

article / block 1 · html

Conclusions. The discovery of LRDs poses a challenge to standard models of early galaxy formation, given their combination of compact sizes, low stellar masses, and potentially massive BHs. In this work, we have explored whether these BHs could have a primordial origin.

S6: VI Conclusions and Discussions / block 1 · html

We proposed a new origin for the SMBHs associated with JWST LRDs. We considered a general catalyzed FOPT in a subdominant, decoupled dark sector at sub-MeV scales, with cosmic expansion always governed by SM radiation and the latent heat deposited solely into DR. We showed that, when the dark-sector PT is catalyzed by DWs, the formation of SMPBHs with masses up to M PBH ∼ 10 10 ​ M ⊙ M_{\rm PBH}\sim 10^{10}M_{\odot} and an abundance appropriate for LRDs is feasible. We also confirmed that strong cosmological and phenomenological constraints, including CMB spectral distortions, Δ ​ N eff \Delta N_{\rm eff} , PT completion, and DW overclosure, can be consistently satisfied in this scenario. Such SMPBHs can naturally serve as seeds of compact, low-metallicity, BH-dominated galaxies associated with LRDs.

S6: VI Conclusions and Discussions / block 1 · html

At the same time, the dark PT generates SGWBs peaking at f peak ∼ 𝒪 ⁡ ( 10 − 2 ​ – ​ 10 − 1 ) ​ nHz f_{\rm peak}\sim\mathcal{O}(10^{-2}\text{--}10^{-1})\,\mathrm{nHz} , partly accessible to current and future PTAs. A joint analysis of the LRD number density, Δ ​ N eff \Delta N_{\rm eff} constraints, and GW data therefore provides a concrete way to test this sub-MeV catalyzed dark-sector PT origin of SMBH seeds. In this sense, the discovery of LRDs may offer an indirect hint of a sub-MeV-scale FOPT and the presence of DWs.

S2.SS1: II.1 Simulation and Initial Conditions / block 1 · html

At z = 1100 z=1100 , when baryons and photons start to decouple, we include an additional uniform baryon matter field with the same resolution as PDM in the initial condition set up, resulting in a total number of N = 2 × 256 3 N=2\times 256^{3} particles with a mass of ∼ 100 ​ ( 18 ) ​ M ⊙ \sim 100(18)~M_{\odot} for PDM (gas) particles. With this setting, we both approximate the initial structure of PDM around the PBH, and let the gas collapse onto this PBH-seeded halo. In running these simulations, the softening length of PDM and gas is set to ϵ PDM = ϵ gas ∼ 0.01 ​ L / N PDM 1 / 3 ≃ 0.01 ​ h − 1 ​ kpc \epsilon_{\rm PDM}=\epsilon_{\rm gas}\sim 0.01\,L/N_{\rm PDM}^{1/3}\simeq 0.01~h^{-1}\rm kpc , and the initial chemical abundances are assigned to the values predicted for the intergalactic medium (IGM) at z = 1100 z=1100 , as summarized in [ 36 ] .

S2.SS1: II.1 Simulation and Initial Conditions / block 1 · html

In addition, we also consider the relative streaming motion between gas and PDM particles by assigning a universal velocity offset. In our earlier work [ 135 ] , we found that streaming can enhance, rather than delay, galaxy formation—contrary to the conclusion generally reached within the standard Λ \Lambda CDM framework [ 106 , 105 , e.g.,] . This enhancement arises from a displacement in the center of mass (away from the PBH) caused by the relative motion between baryons and dark matter, leading to star formation in the dense gas structures that form in the wake of the PBH. These wake-driven overdensities promote gravitational collapse and thus play a catalytic role in early star formation. To systematically explore the impact of streaming velocity on black hole formation, we implement a velocity offset parameterized as a multiple of the root-mean-square streaming velocity σ b ​ χ = 30 ​ km \sigma_{\rm b\chi}=30~\rm km\, s -1 . Specifically, we consider four cases with velocity offsets of 0.4, 0.8, 1.2, and 1.6 σ b ​ χ \sigma_{\rm b\chi} , corresponding to simulation runs labeled as PBH_LW_wstr_fd005 , PBH_LW_str_fd005 , PBH_LW_mstr_fd005 and PBH_LW_sstr_fd005 , respectively.

S5: 5 Conclusions / block 1 · html

Overall, we conclude that halo assembly bias persists into the early Universe, albeit with a reduced amplitude relative to later times. This underscores the importance of accounting for secondary bias when inferring halo masses of high- z z galaxies from clustering measurements. Neglecting assembly bias may lead to systematic offsets in halo mass estimates of up to 1 dex.

S4.SS3: 4.3 Evaluation and limitation of the results / block 1 · html

Overall, we can safely conclude that the SIDM scenario features the lowest halo-mass distributions among the different scenarios considered here. Interestingly, the secondary bias associated with early formation and high concentration largely cancels out.

S3: III Results and Discussions / block 1 · html

Based on the initial conditions and numerical setup described above, we now present the simulation results and discuss their physical implications. To highlight the environmental dependence of PBH-seeded evolution, we organize this section around four connected aspects. We first examine how the large-scale environment drives distinct BH–host evolutionary pathways in Section III.1 . Building on these trends, we identify the main evolutionary phases and possible transitions toward normal galaxy–AGN configurations in Section III.2 . We then connect these regimes to stellar growth and metal enrichment in Section III.3 before translating the resulting evolutionary stages into observable signatures for current and future surveys in Section III.4 .

S3.SS2: III.2 Road to Normalcy / block 1 · html

The results of Section III.1 show that the cosmological environment regulates the gas supply and co-evolution of PBH-hosting systems. We now examine how the initial PBH mass and its placement within the assembling overdense region modify this evolution. Here, “normalcy” is defined in an evolutionary, trend-driven sense. It does not require exact agreement with the local BH–galaxy scaling relations. Instead, it refers to a progressive reduction in the relative dominance of the initial seed, such that stellar and halo growth outpace BH growth and drive both M ∙ / M ⋆ M_{\bullet}/M_{\star} and M ∙ / M h , PBH M_{\bullet}/M_{\rm h,PBH} toward the ranges occupied by more conventional high-redshift galaxy–AGN systems.

S2.SS4: II.4 Black hole feedback / block 1 · html

We implement BH feedback through thermal energy injection, resulting from BH accretion. The energy injected to the nearby cold gas at each timestep δ ​ t \delta t is determined as

S3: III Results / block 1 · html

We compare our model outputs with various LRD observations, specifically considering two classes of BH populations: heavy DCBH seeds, and light stellar-remnant seeds that are accreting at Eddington/super-Eddington rates. For the heavy seed models, we consider 4 cases - the Eddington-limited case, super-Eddington limited growth (up to 1.5 times Eddington), forced super-Eddington growth (where all DCBHs accrete at super-Eddington rates when gas is available with f edd = 1.5 f_{\rm edd}=1.5 ), and the case imposing a higher LW flux threshold for DCBH formation (up to J crit = 3000 − 10000 J_{\rm crit}=3000-10000 ). The default heavy DCBH case used for our analysis is the Eddington-limited model, unless stated otherwise. For the light seeds, we consider 2 cases - super-Eddington limited (up to 1.5 times) and forced super-Eddington growth ( f edd = 1.5 f_{\rm edd}=1.5 ). The default model for light seeds is the super-Eddington limited case unless stated otherwise.

S4.SSx4: Parameter space scanning / block 1 · html

We systematically scan 6.5 < log ⁡ ( M vir / M ⊙ ) < 9.5 6.5<\log(M_{\rm vir}/\mathrm{M}_{\odot})<9.5 and 6 < c < 15 6<c<15 to determine how the BH-seed formation and DM-accretion times depend on halo properties. The results are summarized in Fig. 5 .

S4.SSx4: Parameter space scanning / block 1 · html

We also test the sensitivity of these results to the cross section. For σ 0 = 5 ​ cm 2 ​ g − 1 \sigma_{0}=5\,\mathrm{cm^{2}\,g^{-1}} and w = 100 ​ km / s w=100\,\,{\rm km/s} , the gravothermal phases (Stages 1 and 2) last longer (Fig. 6 ). This longer evolution, however, allows the LMFP envelope to reach a higher density 4 4 4 Before the LMFP-to-SMFP transition, the LMFP density continues to increase. A smaller cross section delays this transition, after which the SMFP region continues to contract while the LMFP envelope remains nearly unchanged. . Consequently, after the BH seed forms, its Stage 5 accretion rate is higher. The total evolution time (Stages 1–5) is close to that in the fiducial model ( σ 0 = 20 ​ cm 2 ​ g − 1 \sigma_{0}=20\,\mathrm{cm^{2}\,g^{-1}} and w = 100 ​ km / s w=100\,\,{\rm km/s} ). These results indicate that, over σ 0 = 5 \sigma_{0}=5 – 20 ​ cm 2 ​ g − 1 20\,\mathrm{cm^{2}\,g^{-1}} , the predicted high-redshift SMBH number density is insensitive to the cross section.