Poisson Clustering and Seed Effect of PBHs: Implications for Microlensing, Gravitational Waves, Little Red Dots and Ultra-Faint Dwarf Galaxies

abstract · abstract

The existence of primordial black holes (PBHs) could have significant effects on the formation of cosmic structures through both the seed effect and the Poisson effect. For an extended PBH mass function, this should lead to dark virialised bound clusters of black holes with modest mass, together with central black holes of larger mass. We consider the combination of three effects: the inevitable gravitational binding of small-scale inhomogeneities due to Poisson and seed fluctuations; the dilution through dynamical heating of all clusters smaller than around $10^{6}\,M_{\odot}$; the survival of clusters in the outer Galactic halo (i.e. towards the Magellanic Clouds) but their destruction in the central halo region through tidal disruption and collisions. We examine the implications of this for the constraints on the fraction of dark matter in PBHs from microlensing surveys and gravitational-wave events. Although the clustering of PBHs with the expected mass and radius is unlikely to relax the microlensing limits, the decline of the Galactic rotation curve might do so, allowing much of the dark matter to be in solar-mass PBHs formed at the quantum chromodynamics transition. We also emphasise that binary merger rates must be revised: binaries formed in PBH clusters can dominate the merger rate, which may be compatible with the LIGO-Virgo-KAGRA gravitational-wave limits if PBHs provide the dark matter. More positively, Ultra-Faint Dwarf Galaxies and Little Red Dots might arise through a combination of the seed and Poisson effects. Our conclusions could be modified for PBHs produced from primordial fluctuations with exponential non-Gaussian tails, since the clustering would then exceed the Poissonian prediction, opening up new possibilities to explain dark matter with solar-mass PBHs.

S2: II  Formation of PBH Clusters from Seed and Poisson Effects / block 1 · html

Some of the discussion in this section appeared in our earlier review [ 6 ] but it has been updated to cover an extended PBH mass function and we correct some minor errors in our earlier analysis. Carr and Silk [24] have clarified the relationship between these two effects and we start by presenting their heuristic analysis before discussing the seed and Poisson effects in more detail. If the PBHs have a single mass m m , the initial fluctuation in the matter density on a scale M M is [ 24 ]

S5: V  PBH Clustering from Non-Gaussian Effects / block 1 · html

Our analysis of the fate of PBH clusters and PBH binaries presented here is essentially independent of the cluster origin. One could try to apply our results to determine the size and mass of more compact clusters that would be predicted in scenarios with additional clustering from non-Gaussian fluctuations. However, it is difficult to predict the initial size and mass of clusters in the non-Gaussian case and numerical simulations are typically needed for an accurate treatment. Choi et al. [73] have recently used such simulations to study the influence of spatial correlations on clustering properties of PBHs with a non-Gaussian exponential tail. Their analysis includes the non-sphericities of the emerging clusters, as well as the PBH mass function. These effects arise because any PBH-producing primordial power spectrum implies, for a range of scales, an increased likelihood for the occurrence of peaks near black holes. The associated dynamics take a finite time, possibly several e-folds (cf. Ref. [ 105 ] ). Therefore, a Hubble patch at the formation of an apparent horizon contains numerous initial Hubble patches (connected to the time of horizon entry) and therefore numerous other peaks which will merge and yield increased black hole masses.

S6.SS1: VI.1 Efficiency Analysis / block 1 · html

A recent paper [ 112 ] , reporting the latest results from the Optical Gravitational Lensing Experiment (OGLE) microlensing survey, illustrates the problem with the reliability of constraints on the MACHO content of the halo. The OGLE collaboration started their microlensing survey of the LMC in 1996 and have developed it through successive phases up to the present. The original idea was to replicate or confirm the MACHO collaboration’s microlensing detections, but the results did not turn out as expected. From the outset, the OGLE detection rate of microlensing events was very low, which they claimed to be consistent with known stellar populations in the Milky Way and self-lensing from stars in the LMC. The combined OGLE-III/OGLE-IV analyses covered 300 deg 2 and included some 62 million stars. The final result was the detection of 13 microlensing candidates, which is less than the number expected from known stellar populations. The SMC results have also been reported [ 113 ] but are notoriously difficult to analyse as one is looking at a small tube-like structure end on, where the expected number of Galactic halo events is very small and typically dominated by Poisson statistics. The probability of self-lensing is then very high and hard to quantify due to the orientation of the SMC relative to the Galaxy. For this reason the SMC has generally been excluded from the microlensing analysis.

S6.SS1: VI.1 Efficiency Analysis / block 1 · html

Perhaps the most convincing reason for distrusting the OGLE results is the status of their recent 13 detections of LMC microlensing events. Quite apart from the question of whether these events can be attributed to PBHs, their number falls far short of the roughly 40 detections expected from known stellar populations. In earlier surveys [ 116 ] , with very small numbers of detections, it could be claimed that they were consistent with known stellar populations. However, in the new much larger sample, the number of OGLE detections [ 112 ] is inconsistent with that expected for known stellar populations at a high significance level. This failure to detect expected microlensing events motivates a sensitivity analysis of the results from the OGLE survey, so any proposed constraints are questionable and should not be considered as firm. It should also be noted that Sugiyama et al. [117] claim to have found evidence for microlensing towards M31 by planetary mass black holes.

S6.SS2: VI.2 Effect of Falling Rotation Curve / block 1 · html

The conclusion that the observed microlensing rate in the Galactic halo is inconsistent with a dark matter halo composed of stellar-mass compact bodies rests fundamentally on the assumption in the analysis of the surveys [ 3 , 118 , 119 ] of a heavy Milky Way halo with a flat rotation curve. The results of the MACHO microlensing survey were published some 20 years ago, when very little was known of the shape of the rotation curve of the Milky Way, as direct measurement was challenging. At that time it was generally believed that the rotation curve was flat, as observed for a number of nearby spiral galaxies. This implied a massive halo, and a large optical depth to microlensing to the Magellanic Clouds, inconsistent with the observations if the halo was made up of compact bodies. When it became possible to make direct measurements of the Galactic rotation curve, it was clear that it declined out to the distance of the Magellanic Clouds [ 120 , 121 , 122 ] . This decline was confirmed by subsequent measurements [ 123 ] and put beyond reasonable doubt by data from Gaia Data Release 3 [ 124 ] . This implies a halo mass which is smaller by an order of magnitude than assumed in the earlier analyses, with a corresponding reduction in the optical depth to microlensing from MACHO dark matter. This would be consistent with the results of Alcock et al. [3] for an all-MACHO halo [ 111 ] .

S6.SS4: VI.4 Effect of Clustering on Microlensing Constraints / block 1 · html

where D L D_{\mathrm{L}} is the distance of the PBH cluster and D S D_{\mathrm{S}} and D LS D_{\rm LS} are the observer-source and lens-source distances, respectively, all of these being comparable. The resulting event may be missed if its timescale lies outside the survey window or if finite-source effects are important. Even if r c r_{\mathrm{c}} exceeds the value given by Eq. ( VI.7 ), one no longer detects standard microlensing by individual PBHs if the Einstein radius r E r_{\mathrm{E}} of the PBH lens is more than the projected separation r sep r_{\text{sep}} of PBHs in the cluster, i.e. if

S2: II  Formation of PBH Clusters from Seed and Poisson Effects / block 1 · html

We mainly emphasise the Poisson effect in this paper but the seed effect may also be important at redshifts above f PBH ​ z eq f_{\text{PBH}}\hskip 1.42262ptz_{\text{eq}} and this may be relevant in the discussion of UFDGs and LRDs. We note that the Poisson effect really leads to clusters with a range of masses at a given redshift. So even if the PBHs have a monochromatic mass function, the clusters would have a range of masses. We derive the associated mass function later.

S7.SS3: VII.3 Two-Body PBH Binary Formation in Clusters / block 1 · html

Reference [ 137 ] illustrates the dependence of R cl R_{\text{cl}} on f PBH f_{\text{PBH}} for several primordial spectra, using the NFW profile and cluster-evolution prescription adopted there. Its results are sensitive to the lightest surviving haloes. In applying Eq. ( VII.10 ) here, the lower limit M min ​ ( z ) M_{\text{min}}(z) must therefore be calculated consistently with the assumed expansion and disruption history; the NFW results of that reference should not be identified with the uniform-density illustration below.

S8: VIII  Thermal History Model and Multimodal Mass Functions / block 1 · html

The previous discussion motivates scenarios in which a large fraction of the dark matter consists of PBHs with an extended mass distribution. In this context, one must distinguish between a simple power-law function and a multimodal function [ 138 ] . The power-law case was discussed in Ref. [ 24 ] . Of particular interest in the multimodal case until now are distributions with a dominant solar-mass population and a subdominant massive tail. This can yield distributions in which the density peaks at a solar mass, while ⟨ m ​ f ⟩ \langle\mspace{1mu}m\mspace{1mu}f\rangle could lie in a much broader range of masses.

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.

S7: VII Conclusion / block 1 · html

We investigate the origin of the high-redshift SMBHs in LRDs observed by JWST from small-scale clustering of PBHs followed by runaway mergers. The derived mass function of LRDs under the small-scale PBH clustering matches the JWST observations for 1 ​ – ​ 10 % 1\textup{--}10\% accretion duty cycle. Other properties of LRDs—including their overmassive nature and compactness—are also explained within our framework when the seed effect and the halo spin distribution are taken into account.

article / block 1 · html

Supermassive black holes (SMBHs) in “little red dots” (LRDs) discovered the James Webb Space Telescope (JWST) may result from runaway mergers of primordial black holes (PBHs) in clusters—through long–short mode coupling on small scales in the early Universe. In this framework, we derive the SMBH mass function, together with the compactness and overmassive features of LRDs. We also estimate that the dense gas residing in PBH clusters is consistent with LRD observations. In addition, SMBHs formed from PBH clusters can help accelerate galaxy formation at high redshifts, thus alleviating tension with Λ \Lambda CDM cosmology.

S2.SS2: II.2 Star Formation and Stellar Feedback / block 1 · html

Table 2: Summary of supernova (SN) feedback parameters from stellar particles. We show values for a single stellar particle of m ⋆ ≃ 600 ​ M ⊙ m_{\star}\simeq 600~\mathrm{M}_{\odot} , based on the average over the stellar IMF. More specifically, we assume that all Pop II stellar particles are identical, whereas Pop III particles contain individual stars sampled from the input IMF [ 47 , 46 ] . Here, E SN E_{\rm SN} is the total SN energy, M Z M_{Z} the mass of the metals released, t ⋆ t_{\star} the typical life time of massive stars, and r final r_{\rm final} the final radius of the SN bubble. One can refer to table 1 in Liu and Bromm [62] for a more detailed description.

S3: III Results / block 1 · html

With the stellar feedback recipes described in Section II , we now analyze the resulting PBH simulation runs, comparing to previous work 4 4 4 Also see the appendix in [ 73 ] . that did not include the full feedback physics. To study the synergy between stellar and BH feedback, metal enrichment and gas outflow, we first track the star formation history (SFH) and the assembly of galaxies in Section III.1 , followed by a discussion of the metal enrichment history in Section III.2 .

S0.SS0.SSS0.Px1: Introduction and Summary.— / block 1 · html

The James Webb Space Telescope (JWST) has recently discovered a class of compact, red-hued galaxies at high redshifts ( 4 ≲ z ≲ 11 4\lesssim z\lesssim 11 ), dubbed “little red dots” (LRDs) [ 1 , 2 , 3 , 4 , 5 ] , for which the origin remains open. Moreover, around 80% of these galaxies display broad Balmer emission lines, suggesting that they host active galactic nuclei (AGNs) with central supermassive black holes (SMBHs) of 10 5 ​ – ​ 8 ​ M ⊙ 10^{5\textup{--}8}{\rm\,M}_{\odot} [ 6 , 7 , 8 , 9 , 10 , 11 , 12 , 13 , 14 , 15 ] .

S0.SS0.SSS0.Px1: Introduction and Summary.— / block 1 · html

Fig. 1 summarizes the parameter space in terms of the modulation strength | η | |\eta| and the comoving wavenumber ratio k l / k s k_{l}/k_{s} , where PBH clusters can form and potentially explain the origin of LRDs (light-gray region) via runaway mergers. The remaining regions represent parameter ranges where SMBHs at high redshift cannot form through mergers, as will be discussed in the following sections.

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: 5 Discussion / block 1 · html

Recently, JWST has discovered a new type of compact galaxies, LRDs, hosting SMBHs of ∼ 10 5 ​ – ​ 8 ​ M ⊙ \sim 10^{5\textup{--}8}{\rm\,M}_{\odot} at redshifts z ∼ 4 ​ – ​ 8 z\sim 4\textup{--}8 Harikane et al. (2023) ; Maiolino and others (2024) ; Kocevski et al. (2023) ; Kokorev et al. (2023) ; Killi et al. (2024) ; Kokorev et al. (2024) ; Wang et al. (2025a) ; Durodola et al. (2025) ; Ananna et al. (2024) ; Kocevski et al. (2025) ; Furtak et al. (2024) ; Maiolino and others (2025) ; Rusakov et al. (2026) . On the other hand, a few unexpectedly large-mass galaxies at z ∼ 10 z\sim 10 are also reported Bogdan and others (2024) ; Kovacs and others (2024) . All of these call for an explanation beyond conventional astrophysics. In particular, the “overmassive” feature of these SMBHs Jones et al. (2025) suggests a heavy-seed origin Zhang et al. (2026b) . Although heavy-seed PBH of ≳ 10 4 ​ M ⊙ \gtrsim 10^{4}{\rm\,M}_{\odot} can just provide a resolution, either some non-Gaussianity in the primordial power spectrum is required or massive black holes (not quite primordial) form only after redshift z ≃ 1100 z\simeq 1100 Qin et al. (2025) ; Cyr (2025) in order to circumvent the CMB μ \mu -distortion constraint De Luca et al. (2025) . The non-Gaussianity can also produce the clustering of PBHs while satisfying all the observational constraints Zhang et al. (2025) ; Belotsky et al. (2019) . If these early galaxies are seeded by PBH clusters, then our runaway-merger model could also account for the rapid emergence of SMBHs at high redshifts.

S5: 5 Discussion / block 1 · html

The high spin of emergent SMBHs is another key feature of the runaway merger of PBH clustering. In Ref. Zhang et al. (2025) , we demonstrate that a PBH cluster can acquire spin via tidal interactions among clusters. The resulting spin can be much larger than the limit of Kerr black holes if all PBHs in the cluster were merged into an SMBH. We expect that the emergent SMBH will inherit most of the spin of the cluster 2 2 2 Even if the original cluster has zero net spin, it has been shown in full general relativity that an initial black-hole cluster of 25 equal-mass, nonspinning black holes produces a most-massive remnant with spin ∼ 0.3 \sim 0.3 times the Kerr limit Bamber et al. (2025) . , and the high-angular-momentum small-mass PBHs remain unmerged, surrounding the central SMBH. The EMRI system resulting from PBH cluster could potentially explain some distinct features of LRDs. For example, the V-shaped spectral energy density can be explained by a central SMBH surrounded by stellar-mass black holes Kritos and Silk (2025) ; Wang et al. (2025b) ; Wang et al. (2025c) ; Chen et al. (2026) . In addition, X-ray quasi-periodic eruptions at the center of galactic nuclei can be interpreted as arising from the interaction between an EMRI and the accretion disk of central SMBHs Kara and García (2025) .

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.

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.