Accelerating Early Massive Galaxy Formation with Primordial Black Holes

Understanding the onset of star and galaxy formation at the end of the cosmic dark ages, a few hundred million years after the Big Bang, is one of the key goals of modern cosmology (e.g., Barkana & Loeb 2001; Bromm & Yoshida 2011). With the successful launch of the James Webb Space Telescope (JWST), this formative period of cosmic history is now becoming accessible to direct observations, thus finally testing the theoretical framework of early structure formation. This framework extends the ΛCDM model, which is highly successful in accounting for galaxy formation and evolution following the epoch of reionization (Springel et al. 2006; Mo et al. 2010) to the first billion years after the Big Bang.

The initial JWST imaging via the Cosmic Evolution Early Release Science (CEERS) survey has revealed a population of surprisingly massive galaxy candidates at z ≳ 10, with inferred stellar masses of ≳10⁹ M_⊙ (Atek et al. 2022; Finkelstein et al. 2022; Harikane et al. 2022; Naidu et al. 2022; Yan et al. 2022). The current record among these photometric detections reaches out to z ≃ 16.7 (Donnan et al. 2022). Such massive sources so early in cosmic history would be difficult to reconcile with the expectation from standard ΛCDM (Boylan-Kolchin 2022; Inayoshi et al. 2022; Lovell et al. 2022), and this includes similarly overmassive (up to ∼10¹¹ M_⊙) galaxy candidates detected at z ≃ 10 (Labbé et al. 2022). An important caveat here is that the early release JWST candidate galaxies are based on photometry only, rendering their redshift and spectral energy distribution (SED) fits uncertain (Steinhardt et al. 2022) until spectroscopic follow-up will become available. Part of the seeming discrepancy may also be alleviated by the absence of dust in sources at the highest redshifts, thus boosting their rest-frame UV luminosities (e.g., Jaacks et al. 2018; Ferrara et al.2022).

It is a long-standing question in cosmology when the first galaxies emerged and how massive they were, going back to the idea that globular clusters formed at the Jeans scale under the conditions immediately following recombination (Peebles & Dicke 1968). The modern view of early structure formation, based on the dominant role of cold dark matter (CDM), posits initial building blocks for galaxy formation that are low mass, of order 10⁶ M_⊙, virializing at z ≃ 20–30 (Couchman & Rees 1986; Haiman et al. 1996).

Numerous models have been proposed to suppress small-scale fluctuations, thus delaying the onset of galaxy formation to later times, in response to empirical hints for the lack or absence of low-mass objects that are otherwise predicted by CDM models (Bullock & Boylan-Kolchin 2017). An extreme scenario here is the fuzzy dark matter (FDM) model, which assumes ultralight axion-like dark matter particles with corresponding de Broglie wavelengths of ∼1 kpc (Hui et al. 2017). Quantum pressure would thus prevent the collapse of structure on the scale of dwarf galaxies (e.g., Sullivan et al. 2018).

The opposite effect, deriving models to accelerate early structure formation beyond the ΛCDM baseline prediction, as may be required to explain the massive JWST galaxy candidates, is much more challenging. One recent study invokes the presence of an early dark energy (EDE) component, resulting in such an accelerated formation of high-redshift structures (Klypin et al. 2021). Here, we explore a different possibility of boosting the emergence of massive galaxies in early cosmic history with primordial black holes (PBHs), considering the isocurvature perturbations from PBHs that increase the power of density fluctuations in addition to the standard ΛCDM adiabatic mode (i.e., the "Poisson" effect; Carr & Silk 2018). In addition, we assess the scenario in which the most massive galaxies reported in Labbé et al. (2022) form in halos seeded by massive PBHs that are very rare and evolve in isolation at high z (i.e., the "seed" effect; Carr & Silk 2018).

For simplicity, we adopt a monochromatic ³ mass function for PBHs, such that a PBH model is specified by the black hole mass m_PBH and fraction f_PBH of dark matter in the form of PBHs. We start with the "Poisson" effect (Carr & Silk 2018) in which PBHs produce isocurvature perturbations in the density field on top of the standard adiabatic mode due to the random distribution of PBHs (at small scales). These isocurvature perturbations only grow in the matter-dominated era. As an extension of the formalism in Afshordi et al. (2003), Kashlinsky (2016), and Cappelluti et al. (2022), where PBHs make up all dark matter (f_PBH ∼ 1), now treating f_PBH as a free parameter, the linear power spectrum (extrapolated to a = 1) of dark matter density fluctuations can be written as (Inman 2019; Liu et al. 2022):

$$\begin{eqnarray}\begin{array}{rcl}P(k) & = & {P}_{\mathrm{ad}}(k)+{P}_{\mathrm{iso}}(k),\\ {P}_{\mathrm{iso}}(k) & \simeq & {[{f}_{\mathrm{PBH}}{D}_{0}]}^{2}/{\bar{n}}_{\mathrm{PBH}},\end{array}\end{eqnarray} \tag{ 1 }$$

where P_ad(k) is the standard adiabatic mode in ΛCDM cosmology, ⁴ ${\bar{n}}_{\mathrm{PBH}}={f}_{\mathrm{PBH}}\tfrac{3{H}_{0}^{2}}{8\pi G}({{\rm{\Omega }}}_{{\rm{m}}}-{{\rm{\Omega }}}_{{\rm{b}}})/{m}_{\mathrm{PBH}}$ is the cosmic (comoving) number density of PBHs, and D₀ is the growth factor of isocurvature perturbations evaluated at a = 1, given by Inman (2019):

$$\begin{eqnarray}\begin{array}{rcl}D(a) & \simeq & {\left(1+\displaystyle \frac{3\gamma }{2{a}_{-}}s\right)}^{{a}_{-}}-1,\quad s=\displaystyle \frac{a}{{a}_{\mathrm{eq}}},\\ \gamma & = & \displaystyle \frac{{{\rm{\Omega }}}_{{\rm{m}}}-{{\rm{\Omega }}}_{{\rm{b}}}}{{{\rm{\Omega }}}_{{\rm{m}}}},\quad {a}_{-}=\displaystyle \frac{1}{4}\left(\sqrt{1+24\gamma }-1\right),\end{array}\end{eqnarray} \tag{ 2 }$$

where a_eq = 1/(1 + z_eq) is the scale factor at matter–radiation equality with z_eq ≃ 3400. So far, we have ignored the higher-order (nonlinear) "seed" effect (Carr & Rees 1984; Carr & Silk 2018), as well as mode mixing, ⁵ lacking a self-consistent linear perturbation theory that can take them into account. Such effects tend to enhance the perturbations induced by PBHs at large scales but suppress structure formation at small scales with nonlinear dynamics around PBHs (Liu et al. 2022). Particularly, we expect linear perturbation theory to break down at a certain (mass) scale M_bk, overproducing the abundance of structures below M_bk. We hypothesize M_bk to be between m_PBH and the mass ${M}_{{\rm{B}}}({m}_{\mathrm{PBH}},z)\sim {[(z+1){a}_{\mathrm{eq}}]}^{-1}{m}_{\mathrm{PBH}}$ of overdense (δ ≳ 1) regions bound to individual PBHs, when they evolve in isolation (Carr & Silk 2018). In light of this, as a heuristic approach, we further impose a conservative cutoff to the isocurvature term in Equation (1) to suppress the power at scales smaller than M_bk ∼ m_PBH:

$$\begin{eqnarray}&&{P}_{\mathrm{iso}}(k)=0,\quad k\gt {\left(2{\pi }^{2}{\bar{n}}_{\mathrm{PBH}}/{f}_{\mathrm{PBH}}\right)}^{-1/3}.\end{eqnarray} \tag{ 3 }$$

As examples, Figure 1 shows the power spectra for three models that turn out to be representative in our analysis below (see Section 3 and Table 1).

Zoom In Zoom Out Reset image size

Once P(k) is known, we use the Press–Schechter (PS) formalism (Press & Schechter 1974; Mo et al. 2010) with a Gaussian window function to calculate the halo mass function (HMF), ${dn}(M,z)/{dM}$, including corrections for ellipsoidal dynamics (Sheth & Tormen 1999). Following Boylan-Kolchin (2022), given the star formation efficiency (SFE), ≡ M_⋆/(f_b M_halo), we then derive the (comoving) cumulative stellar mass density contained within galaxies above a certain stellar mass M_⋆ as

$$\begin{eqnarray}\begin{array}{rcl}{\rho }_{\star }(\gt {M}_{\star },z) & = & \epsilon {f}_{{\rm{b}}}\rho (\gt {M}_{\mathrm{halo}},z)\\ & = & \epsilon {f}_{{\rm{b}}}{\displaystyle \int }_{{M}_{\mathrm{halo}}}^{\infty }M\displaystyle \frac{{dn}(M,z)}{{dM}}{dM},\end{array}\end{eqnarray} \tag{ 4 }$$

where f_b = Ω_b/Ω_m is the cosmic average baryon fraction. Exploring whether PBH scenarios can explain the massive galaxy candidates reported by Labbé et al. (2022), in the following sections we compare the ρ_⋆(>M_⋆, z) and number density of massive galaxies predicted by our PBH models with observations and discuss the general signatures of PBHs observable by JWST.

3.1. PBH Models Required to Explain Current JWST Results

Based on 14 galaxy candidates with masses of ∼10⁹–10¹¹ M_⊙ at 7 < z < 11 identified in the JWST CEERS program, Labbé et al. (2022, see their Figure 4) derive the cumulative stellar mass density at z = 8 and 10 for M_⋆ ≳ 10¹⁰ M_⊙. In particular, they find ${\rho }_{\star }(\gtrsim {10}^{10}\ {M}_{\odot })\simeq {1.3}_{-0.6}^{+1.1}\times {10}^{6}\ {M}_{\odot }\ {\mathrm{Mpc}}^{-3}$ and ${\rho }_{\star }(\gtrsim {10}^{10.5}\ {M}_{\odot })\simeq {9}_{-6}^{+11}\times {10}^{5}\ {M}_{\odot }\ {\mathrm{Mpc}}^{-3}$ at z ∼ 10, higher than the maximum achievable in ΛCDM (with = 1) by up to a factor of ∼50. Using the formalism described in the previous section (Equations (1)–(4)), we find the PBH parameters that can reproduce these results for SFE values = 1 and 0.1, comparing with existing observational constraints in the f_PBH–m_PBH space, as shown in Figure 2.

Zoom In Zoom Out Reset image size

It turns out that, to explain the observational data at M_⋆ ∼ 10^{10.5 (10)} M_⊙, we must have m_PBH f_PBH ≳ 1.8 × 10⁵ (2.4 × 10⁴) M_⊙ for < 1 and m_PBH f_PBH ≳ 6.1 × 10⁶ (1.8 × 10⁶) M_⊙ for < 0.1, indicating that we need relatively massive PBHs with m_PBH ≳ 10⁵ M_⊙ since f_PBH ≤ 1. To form such massive PBHs in the standard spherical collapse scenario (Escrivà 2022), very large perturbations are required out of inflation, e.g., from non-Gaussian tails produced by interacting quantum fields (e.g., Frampton 2016; Atal et al. 2019; Panagopoulos & Silverstein 2019), or oscillatory features in the inflationary power spectrum (Carr & Kühnel 2019). Such PBH models are strongly constrained by observations of μ-distortion in the cosmic microwave background (CMB), if primordial density fluctuations are Gaussian (Carr et al. 2021b, see their Figures 10 and 16).

Nevertheless, the constraints can be weaker when the Gaussian assumption is relaxed (e.g., Nakama et al. 2018). For instance, with the long-dashed curve in Figure 2 we show a particular case of the phenomenological non-Gaussian model in Nakama et al. (2016) with p = 0.5, where p is the non-Gaussianity parameter (p = 2 means Gaussian). In this case, PBHs with m_PBH ≳ 10¹¹ M_⊙ can still have large enough f_PBH (≳10⁻⁵) to explain both the JWST observations and the CMB μ-distortion limit (Nakama et al. 2018), as measured by the COBE Far Infrared Absolute Spectrophotometer (FIRAS). However, such models are disfavored by the nondetection of black holes above 10¹¹ M_⊙ (except for Phoenix A; Brockamp et al. 2016). Therefore, an even higher degree of non-Gaussianity is required to explain the JWST results with PBH models of m_PBH ≲ 10¹¹ M_⊙ and high enough f_PBH. Alternatively, the CMB μ-distortion constraint can be evaded if PBHs grow significantly between the μ-distortion epoch (7 × 10⁶ s < t < 3 × 10⁹ s) and matter–radiation equality (Carr et al. 2021b), or if PBHs form in nonstandard scenarios such as inhomogeneous baryogenesis with the modified Affleck–Dine mechanism (Kawasaki & Murai 2019; Kasai et al. 2022).

Beyond the μ-distortion constraint, PBHs with m_PBH ≳ 10⁵ M_⊙ are also constrained by X-ray binaries (XB, Inoue & Kusenko 2017), infall of PBHs into the Galactic center by dynamical friction (DF, Carr & Sakellariadou 1999), and large-scale structure statistics (LSS, Carr & Silk 2018), which together require f_PBH ≲ 10⁻⁴–10⁻³ for m_PBH ∼ 10⁵–10¹¹ M_⊙ (see the dashed–dotted curve in Figure 2). Such constraints are generally weaker, allowing a PBH abundance for m_PBH ≳ 10⁹ M_⊙ sufficient to produce the high stellar mass density in massive galaxies at z ∼ 10 inferred by JWST (Labbé et al. 2022). However, in this regime (f_PBH ≲ 10⁻³), the "seed" effect likely dominates at z ≳ 10 according to the criterion that the fraction of mass bound to PBHs in the universe is ≪1, i.e., f_PBH ≪ (1 + z)a_eq (Carr & Silk 2018). Besides, isocurvature perturbations purely from the "Poisson" effect (without the cutoff in Equation (3)) are strongly constrained by high-z observations (e.g., Sekiguchi et al. 2014; Murgia et al. 2019; Tashiro & Kadota 2021), such that PBH models with m_PBH f_PBH ≳ 170 M_⊙ (and f_PBH > 0.05) are ruled out by high-z Lyα forest data (Murgia et al. 2019). This motivates us to further explore the "seed" effect.

We perform an idealized calculation to explore what PBH models are needed to form the most massive galaxies observed by JWST purely with the "seed" effect. Assuming that the observed massive galaxies form in halos seeded by PBHs that grow in isolation as ${M}_{\mathrm{halo}}\sim {M}_{{\rm{B}}}({m}_{\mathrm{PBH}},z)\sim {[(z+1){a}_{\mathrm{eq}}]}^{-1}{m}_{\mathrm{PBH}}$ (Carr & Silk 2018), we require that (i) the average (comoving) number density ${\bar{n}}_{\mathrm{PBH}}$ of PBHs is larger than that of the observed massive galaxies n_g ∼ 2 × 10⁻⁵ Mpc⁻³ (Boylan-Kolchin 2022; Labbé et al. 2022), and (ii) the PBH-seeded halos have enough gas to form M_⋆ ∼ 10¹¹ M_⊙ of stars, i.e., f_b M_B = M_⋆ can be satisfied given the limit of . As shown in Figure 3, we find that the "seed" effect can explain the JWST results with generally less extreme PBH models of m_PBH f_PBH ≳ 3 × 10^{5 (3)} M_⊙ and m_PBH ≳2 × 10^{10 (9)} M_⊙ for ≲ 0.1 (1). Moreover, considering the nonlinear dynamics around PBH-seeded halos, it is found in cosmological simulations that structures smaller than∼M_B (≳10¹² M_⊙ at z ≲ 10 for m_PBH ≳ 10⁹ M_⊙) are significantly less abundant than predicted by the "Poisson" effect (Liu et al. 2022, see their Figure 14). This can weaken/lift the Lyα forest constraint that is sensitive for k ∼ 0.7–10 hMpc⁻¹ (Murgia et al. 2019), corresponding to M_halo ∼ 10⁹–10¹² M_⊙.

Zoom In Zoom Out Reset image size

3.2. Signatures of Representative PBH Models

To further demonstrate the effects of PBHs in high-z galaxy/structure formation potentially observable by JWST, we focus on three models representative of typical regions in the PBH parameter space (Figure 2), as listed in Table 1. Here, M1 satisfies all observational constraints considered above, while M2 and M3 only satisfy the XB+DF+LSS combined constraint. For the "Poisson" effect, M3 is consistent with the recent JWST results in Labbé et al. (2022) with ∼ 0.1–1, while M1 and M2 cannot reproduce the observations even with = 1. For the "seed" effect, M2 can also marginally explain the JWST observations.

In Figure 4, we present the results for ρ_⋆(>M_⋆) in standard ΛCDM (solid), M1 (dashed), M2 (dashed–dotted), and M3 (dotted) for = 1 (thin) and 0.1 (thick) at z = 10, 15, and 20. We also plot the lower mass limits of galaxies detectable by JWST CEERS NIRCam LW imaging (for an exposure time of ∼2000 s) in Figure 4, considering a magnitude limit of m_F444W ∼ 28 for the filter F444W with a signal-to-noise ratio (SNR) ≳ 5, derived from the JWST Exposure Time Calculator. ⁶ We calculate the stellar mass limit as the total mass of stars formed in a star formation event whose peak luminosity matches the magnitude limit according to the stellar population synthesis code yggdrasil (Zackrisson et al. 2011). For Population III (Pop III) stars we adopt their (instantaneous-burst) Pop III.1 model with an extremely top-heavy Salpeter initial mass function (IMF) in the range of 50–500 M_⊙, based on Schaerer (2002), while for Population II (Pop II) stars we use their Z = 0.0004 model from Starburst99 (Leitherer et al. 1999; Vázquez & Leitherer 2005) with a universal Kroupa (2001) IMF in the interval 0.1–100 M_⊙. For both Pop III and II, we consider a nebula covering fraction of f_cov = 0.5 and no Lyα transmission. Note that these limits will be reduced by about one order of magnitude for future deeper surveys with longer (∼100 h) exposures (Zackrisson et al. 2012).

Zoom In Zoom Out Reset image size

As shown in Figure 4, ρ_⋆(>M_⋆) can be significantly increased by PBHs, and the effects are the strongest for halos of M_halo ∼ m_PBH−10m_PBH, and stronger at higher z. In the extreme model M3, the stellar mass budget at z ≲ 20 will be dominated by galaxies more massive than 10⁸ M_⊙ given ≳ 0.1. Fully taking into account the "seed" effect may further increase the mass budget of massive halos (M_halo ∼ m_PBH−M_B) with respect to that of the smaller ones (Liu et al. 2022). Even in M1, ρ_⋆(>M_⋆) is higher by a factor of ≳10 than in ΛCDM for M_⋆ ∼ 10⁵ M_⊙ at z ≳ 15. Considering the area ($\sim 40\ {\mathrm{arcmin}}^{2}$) of CEERS, corresponding to the number count limit in Figure 4 (i.e., the dashed–dotted–dotted line), we conclude that the fact that JWST has detected galaxies above 10⁸ M_⊙ at z ≳ 15 (e.g., Atek et al. 2022; Donnan et al. 2022; Harikane et al. 2022; Naidu et al. 2022; Yan et al. 2022) already requires a very high SFE, ≳ 0.1, for Pop II star formation in ΛCDM, consistent with the UV luminosity function analysis by Inayoshi et al. (2022). That said, for structure formation accelerated by PBHs, such galaxies can also form in more massive halos with a lower SFE (down to 0.01 for M3). Finally, our calculation indicates that at z ≳ 20 detection of any Pop II galaxy by a survey like CEERS is impossible in ΛCDM, and even Pop III galaxies require ≳ 0.1 to be detected with stellar masses above 10⁷ M_⊙, which is inconsistent with the theoretical and observational upper limits, M_{Pop III} ≲ 10⁶ M_⊙, on the total mass of (active) Pop III stars a halo can host (Yajima & Khochfar 2017; Bhatawdekar & Conselice 2021). However, with massive PBHs like those in M3, it is possible to detect Pop II galaxies as massive as ≳4 × 10⁸ M_⊙ with ≳ 0.01.

The recent detection of surprisingly massive (∼10¹⁰–10¹¹ M_⊙) galaxy candidates at z ∼ 10 by JWST (Labbé et al. 2022), if confirmed, brings new challenges to ΛCDM (and a broad range of dynamical dark energy models) that cannot provide enough baryonic matter in collapsed structures for such early massive galaxy formation (Boylan-Kolchin 2022; Lovell et al. 2022; Menci et al. 2022). Note that the source properties derived from photometry are sensitive to the underlying SED fitting templates, such that the stellar masses can be lower by up to 1.6 dex than reported in Labbé et al. (2022), if other templates with improved physical justifications are adopted, removing the tension with ΛCDM (Steinhardt et al. 2022; see their Figure 3). Currently no template can well match the observed photometry self-consistently, and follow-up spectroscopic observations are needed to robustly pin down the nature and properties of these galaxy candidates (Furlanetto & Mirocha 2022).

Assuming that the results in Labbé et al. (2022) are true and considering only the "Poisson" effect of PBHs, we use a simple analytical model based on linear perturbation theory and the PS formalism to show that such early formation of massive galaxies is possible if PBHs make up part of dark matter with m_PBH f_PBH ≳ 6 × 10⁶ (2 × 10⁵) M_⊙ for SFE < 0.1 (1). The "seed" effect, however, requires m_PBH f_PBH ≳ 3 × 10^{5 (3)} M_⊙ and m_PBH ≳ 2 × 10^{10 (9)} M_⊙ for ≲ 0.1 (1). Such massive PBHs are mostly ruled out by observations of the CMB μ-distortion, although this strong constraint relies on the assumptions that primordial density fluctuations are Gaussian and that PBHs (formed in the standard scenario of primordial density fluctuations) hardly grow during the radiation-dominated era. Besides, strong isocurvature perturbations purely from the "Poisson" effect of PBHs with m_PBH f_PBH ≳ 170 M_⊙ at scales of k ∼ 0.7–10 hMpc⁻¹ are ruled out by high-z Lyα forest data (Murgia et al. 2019). Nevertheless, considering the nonlinear dynamics around massive (m_PBH ≳ 10⁹ M_⊙) PBHs, the abundance of halos at such scales can be much lower than predicted by the "Poisson" effect (Liu et al. 2022), lifting the Lyα forest constraint. The μ-distortion constraint can also be evaded by relaxing the underlying assumptions or considering nonstandard PBH formation mechanisms (e.g., Kawasaki & Murai 2019; Kasai et al. 2022). In this way, the other constraints from X-ray binaries (Inayoshi et al. 2022), dynamical friction (Carr & Sakellariadou 1999), and large-scale structures (Carr & Silk 2018) allow for a substantial region in the PBH parameter space with m_PBH ∼ 10⁹–10¹¹ M_⊙ and f_PBH ∼ 10⁻⁶–10⁻³ that is consistent with the recent JWST observations (and no detection of black holes above 10¹¹ M_⊙).

However, considering that PBHs can grow by up to 2 orders of magnitude through the acquisition of dark matter halos by z ∼ 10 with optimistic spherical ⁷ accretion (e.g., Mack et al. 2007; De Luca et al. 2020), and the nondetection of black holes above ∼10¹¹ M_⊙, the allowed region will further shrink. That said, if the stellar masses measured by Labbé et al. (2022) are overestimated (or some of the galaxy candidates are actually at lower z) such that the stellar mass density at z ∼ 10 is lower in reality, less extreme PBH models would be able to explain them for the same value of , and the same PBH model could allow lower values of . For instance, if the inferred stellar masses were indeed reduced by 1.6 dex (Steinhardt et al. 2022) in follow-up observations, we would only need m_PBH f_PBH ≳ 2 × 10⁵ M_⊙ to form the observed galaxies with < 0.025 from the "Poisson" effect, while the "seed" effect requires m_PBH f_PBH ≳ 200 (2) M_⊙ and m_PBH ≳ 5 × 10^{8 (7)} M_⊙ for ≲ 0.1 (1).

We also find that the effects of PBHs are stronger at higher z, implying that stronger signatures of PBHs may be found in future wider and deeper surveys by JWST. Actually, if the object CEERS-1749 reported in Naidu et al. (2022) is a galaxy at z ∼ 17, its large mass (∼5 × 10⁹ M_⊙) is also in >3σ tension with ΛCDM (Lovell et al. 2022; see their Figure 6). Even if follow-up studies do not find such an excess of very massive galaxies at higher z, this may not necessarily rule out our PBH models, since the SFE can evolve rapidly with redshift at Cosmic Dawn (z ∼ 6–30) due to metal enrichment, the build-up of radiation backgrounds, and the transition in dominant stellar population (see, e.g., Fialkov & Barkana 2019; Mirocha & Furlanetto 2019; Schauer et al. 2019; Liu et al. 2020). Besides, if the average SFE is not significantly lower with PBHs, the accelerated structure formation leads to accelerated star formation that can also facilitate cosmic reionization. Interestingly, this may explain the recent observations of Lyα emitting galaxies that support a double-reionization scenario with the first full ionization event happening at z ∼ 10 (Salvador-Solé et al. 2017, 2022). Very massive (∼100–10³ M_⊙) Pop III stars are required to produce this double-reionization feature in ΛCDM (Salvador-Solé et al. 2017), while less extreme stellar populations may be sufficient with accelerated star formation by PBHs given proper values of SFE.

Note that the PBH masses and the enhancement of density fluctuations by PBHs required to explain the JWST observations in our case (m_PBH ≳ 10⁹ M_⊙ and m_PBH f_PBH ≳ 3 × 10³ M_⊙) are significantly higher (and working at larger scales) than those considered by previous studies for stellar-mass PBHs (m_PBH ≲ 100 M_⊙ and m_PBH f_PBH ≲ 30 M_⊙; see, e.g., Kashlinsky 2016; Gong & Kitajima 2017; Inman 2019; Cappelluti et al. 2022; Liu et al. 2022). In these studies, the standard picture of first star formation in molecular-cooling minihalos remains unchanged, but with earlier onset and increased abundances of star-forming halos. However, in our extreme PBH models, star formation may first occur in atomic-cooling halos (M_halo ≳ 10⁸ M_⊙), seeded by PBHs or formed in a top-down fashion by fragmentation of massive (≳10¹¹ M_⊙) structures around PBHs. The accretion feedback from black holes can also significantly affect nearby formation of stars and direct collapse black holes (e.g., Pandey & Mangalam 2018; Aykutalp et al. 2020; Liu et al. 2022). Considering the gravitational, hydrodynamic, and radiative effects, early star formation in the presence of massive PBHs can leave strong imprints in the CMB, 21 cm signal, reionization, cosmic infrared, and X-ray backgrounds (see, e.g., Sekiguchi et al. 2014; Kashlinsky 2016; Gong & Kitajima 2017; Murgia et al. 2019; Hasinger 2020; Tashiro & Kadota 2021; Cappelluti et al. 2022; Minoda et al. 2022). These important aspects involving nonlinear dynamics and baryonic physics are beyond the scope of our exploratory work. Follow-up studies with more detailed modeling of the interactions between PBHs, baryons, and (non-PBH) dark matter (in cosmological simulations) are required to fully understand the roles played by massive PBHs in structure/galaxy/star formation and evaluate the viability of such PBH models.

In general, the recent discovery of potentially very early massive galaxy formation by JWST (Labbé et al. 2022) hints at faster structure formation in the high-z universe, above the ΛCDM baseline (Boylan-Kolchin 2022; Lovell et al. 2022), which can be achieved if a small fraction (∼10⁻⁶–10⁻³) of dark matter is composed of massive (≳10⁹ M_⊙) PBHs, although more work needs to be done to check whether such fast structure formation is consistent with other observations of the high-z Universe. A similar trend is also seen in observations of (proto-) galaxy clusters that show an excess of strong-lensing sources (Meneghetti et al. 2020, 2022) and of star formation (Remus et al. 2022) that are difficult to explain in ΛCDM. Strikingly, this trend for accelerated structure formation at high z goes in the opposite direction to that of invoking the suppression of fluctuations to account for the well-known small-scale problems of ΛCDM. Together with other hints for PBHs (see, e.g., Clesse & García-Bellido 2018), these findings imply that the nature of dark matter may be more complex than our standard expectation, e.g., involving important subcomponents such as PBHs. This challenge calls for the thorough theoretical exploration of the alternatives to ΛCDM, in conjunction with advanced observational campaigns to probe the high-redshift universe. Here, we are entering an exciting period of discovery with frontier facilities such as JWST, Euclid, and the Square Kilometre Array, as well as the Einstein Telescope and Laser Interferometer Space Antenna gravitational wave observatories.

We would like to thank Antonio Riotto and the anonymous referee for their helpful comments.