1 Introduction↩︎

The James Webb Space Telescope (JWST) has opened an unprecedented observational window into the early Universe, enabling detailed investigations of the birth of the first galaxies and stars. JWST have uncovered a surprisingly abundant population of massive galaxies at high redshifts (\(z \gtrsim 10\)), significantly exceeding the predictions of standard galaxy formation model [1]–[5]. In particular, the observed excess at the bright end of the ultraviolet luminosity function (UV LF) presents a notable challenge. [6]–[14] Reconciling these observations with standard models would require an unrealistically high star formation efficiency (SFE) [15], [16], or alternatively, explaining the data within the \(\Lambda\)CDM framework by allowing the SFE to evolve with redshift [17].

It is intriguing to consider the potential role of primordial black holes (PBHs) [18]–[36], which have been extensively studied as viable candidates for dark matter, especially supermassive PBHs (SMPBHs) with \(M_{\rm PBH}\gtrsim 10^6M_\odot\). PBHs, especially SMPBHs, might be explanations for the supermassive galaxies[37]–[39] and "Little Red Dots" (LRDs)[40] recently discovered by the JWST. It has been showed that accreting PBHs can contribute additional ultraviolet luminosity, boosting the number density of bright galaxies [41], [42], while the Poisson fluctuations in spatial distribution of PBHs can modulate the halo mass function (HMF), enhancing the abundance of collapsed structures at high redshift [43], see also the impact of initial clustering of SMPBHs [39].

These SMPBHs could help to bridge the gap between the standard \(\Lambda\)CDM model and JWST observations. The mass range of PBHs can extend from \(10^{-18}\rm~M_\odot\) to \(10^{16}\rm~M_\odot\) ⁴, although current observations indicate that SMBHs in galactic nuclei typically have masses extending up to only nearly \(10^{11}\rm~M_\odot\) [47]. In the standard scenario that PBHs were sourced by the collapse of large density fluctuations [48]–[56], massive PBHs with \(M\gtrsim 10^5M_\odot\) is strictly constrained by CMB \(\mu\)-distortion observations [57]–[60]. However, it can be expected that significant non-Gaussianity of primordial perturbations [61]–[69] might be helpful for overcoming this difficulty, or alternative mechanisms, such as PBHs that the supercritical bubbles nucleated during inflation developed into (which are not constrained by \(\mu\) distortion, e.g.[64]), must be considered. In the latter case, the resulting SMPBHs can have a (multiple) peak-like mass spectrum [70], [71]. Therefor, it is significant to further investigate the role of SMPBHs in shaping the early Universe.

This work focuses on how SMPBHs can influence the UV LF of high-redshift galaxies through both Poisson effects and accretion-induced emission, potentially offering an explanation for the recent JWST observations. The values of cosmological parameters adopted in this paper are the bestfit values of Planck based on the \(\Lambda\)CDM model [72]⁵.

2 Influence of SMPBHs on UV luminosity↩︎

It has been showed in Ref.[70] that in string landscape scenario with multiple metastable vacua (usually modelled as a multi-dimensional potential), when the inflaton slowly passes by a neighboring vacuum, the nucleating rate of supercritical bubbles would inevitably present a peak, so that the resulting PBHs have a non-negligible distinctive peak-like mass distribution [80]: \[\label{eq:psi95m1} \psi_{1}(m) = e^{-\sigma^2/8} \cdot \sqrt{\frac{M_{PBH}}{2\pi \sigma^2 m^3}} \cdot \exp\left[ -\frac{\ln ^2 (m / M_{PBH})}{2\sigma^2} \right],\tag{1}\] where \(M_{PBH}\) is the characteristic mass of PBHs, and \(\sigma\) is the width of mass peak. This power-law term causes the distribution to go down more rapidly at larger masses. The corresponding mean mass of the distribution is: \(\langle m \rangle = M_{PBH} e^{-\sigma^2}\)

Here, a conversion relation is required to connect Eq. (1 ) with the number-density-related mass function of PBHs: \[\label{eq:psi95m3} \psi_{3}(m) = \frac{\langle m\rangle}{m} \psi_{1} = e^{-9\sigma^{2}/8} \sqrt{\frac{M_{PBH}^{3}}{2\pi\sigma^{2}m^{5}}} \exp\left[-\frac{\ln^{2}(m/M_{PBH})}{2\sigma^{2}}\right].\tag{2}\] This distribution differs from the standard log-normal distribution by including an extra \(m^{-3/2}\) factor. The interrelations between various definitions of mass distribution and the corresponding expectation values are summarized in Table 3. Thus we have \[\label{eq:dn95dm95PBH} \frac{dn}{n_0 dm} = \psi_3,\tag{3}\] with the comoving number density of PBHs \(n_{0} = \frac{ \rho_{c}\Omega_{dm}f_{\mathrm{PBH}}}{\langle m\rangle} = \frac{\rho_{c}\Omega_{dm}f_{\mathrm{PBH}}}{M_{PBH}} e^{\sigma^{2}}\), where \(n\) is the number density of PBHs, \(\Omega_{dm}\) is the dark matter density parameter and \(\rho_{c}\) is the critical density of the Universe.

The SMPBHs sourced by such inflationary vacuum bubbles can avoid constraint from spectral distortions in the CMB, and thus can be supermassive, \(M_{PBH}\gtrsim 10^6M_\odot\). This mechanism offers a natural route to SMPBHs coming into being without requiring large curvature perturbations.

It can be expected that SMPBHs would enhance galaxy luminosity through the accretion of matter and subsequent emission of radiation. This process could potentially boost the bright end of the UV LF. Here, we quantify the contribution of SMPBHs to the luminosity using the Eddington ratio, \(\lambda_E\). Specifically, the bolometric luminosity of a PBH is⁶ \(L_{\mathrm{bol}}(m_{\mathrm{PBH}}) = \lambda_E L_E(m_{\mathrm{PBH}})\). Here, \(L_E\) represents the Eddington luminosity, defined as: \[L_E(m_{\mathrm{PBH}}) = \frac{4\pi G m_p c}{\sigma_T} m_{\mathrm{PBH}} = 3.28 \times 10^4 \frac{m_{\mathrm{PBH}}}{M_{\odot}} L_{\odot},\] where \(m_p\) denotes the proton mass, \(c\) is the speed of light, \(\sigma_T\) signifies the Thomson scattering cross-section, \(m_{\mathrm{PBH}}\) is the mass of PBHs, and solar Luminosity is \(L_{\odot} \approx 3.828 \times 10^{33} \, \text{erg s}^{-1}\).

According to the prescription of [81], the bolometric luminosity can be converted into the UV luminosity, \[\label{eq:LUV-PBH} L_{\mathrm{UV, PBH}}(m_{\mathrm{PBH}}) = \frac{L_{\mathrm{bol}}(m_{\mathrm{PBH}})}{c_1 \left( \frac{L_{\mathrm{bol}}(m_{\mathrm{PBH}})}{10^{10} L_{\odot}} \right)^{k_1} + c_2 \left( \frac{L_{\mathrm{bol}}(m_{\mathrm{PBH}})}{10^{10} L_{\odot}} \right)^{k_2}},\tag{4}\] where \[\begin{pmatrix} c_1 \\ k_1 \\ c_2 \\ k_2 \end{pmatrix} = \begin{pmatrix} 1.862 \\ -0.361 \\ 4.870 \\ -0.0063 \end{pmatrix}.\] The UV luminosity \(L_{\mathrm{UV,\,PBH}}\) here refers to the band UV luminosity of AB magnitude system measured in a top-hat filter centered at rest-frame 1450 Åwith a bandwidth of 100 Å, or equivalently, the luminosity \(\nu L_\nu\) at 1450.

3 Confronting SMPBH model with JWST observations↩︎

3.1 Total UV luminosity with stars and PBHs↩︎

The UV LF: \[\label{eq:UV95LF} \Phi_{\rm UV} \equiv \frac{dn}{dM_{\rm UV}}(z,M_{\rm UV})\tag{5}\] quantifies the number density of galaxies per luminosity or absolute UV magnitude interval. It can be related to the dark matter HMF by using \(\Phi_{\rm UV}= {\left|\frac{dM_{\rm UV}}{dM_{\rm halo}}\right|}^{-1}\frac{dn}{dM_{\rm halo}}\), where the HMF \(\frac{dn}{dM_{\rm halo}}\) can be calculated with the Sheth-Tormen algorithm [82].

The monochromatic UV luminosity and AB magnitude are related by [83] \[\label{eq:MUV-LUV} M_{\rm UV} = -2.5\log\left(\frac{L_{UV}/ \nu}{\mathrm{erg\,s^{-1}\,Hz^{-1}}}\right) + 51.60.\tag{6}\] In this work, we consider the ultraviolet wavelength band at \(\nu = 1450\,\text{\AA}\). It is important to note whether the luminosity used is the monochromatic UV luminosity \(L_\nu\) or the band UV luminosity \(L_{\mathrm{UV}}\), both are approximately related by \(L_{\mathrm{UV}} \approx \nu L_\nu\).

Here, we adopt a purely empirical fit of the \(M_{\rm halo}\)-\(M_{\rm UV}\) relation [42]: \[\label{eq:mh95muv95fit} M_{\rm UV}= p_1 + \left[ p_2 \left(\frac{M_{\rm halo}}{10^{11} M_\odot}\right)^{p_3} + p_4\right] \log \left(\frac{M_{\rm halo}}{10^{11} M_\odot}\right) + p_5 z^{p_6}\tag{7}\] where \(p_i\) are set with a least-squares fit of Eq. (7 ) to the observed UV LF data [84] (see Figure 1), \[\label{eq:fit-results-muv-vs-mh} \begin{pmatrix} p_1 \\ p_2 \\ p_3 \\ p_4 \\ p_5 \\ p_6 \end{pmatrix}_{\rm fit} = \begin{pmatrix} -13.12 \\ -8.94 \\ -0.0355 \\ 5.78 \\ -3.84 \\ 0.341 \end{pmatrix}\tag{8}\]

The total UV luminosity can be contributed by not only stars but also SMPBHs, when we consider the impact of SMPBHs, i.e., \[\label{eq:LUV-tot} L_{\mathrm{UV,\,tot}}(M_{\rm halo},z,m_{\mathrm{PBH}}) = L_{\mathrm{UV}}(M_{\rm halo}, z) + L_{\mathrm{UV,\,PBH}}(m_{\mathrm{PBH}}),\tag{9}\] where the stellar UV magnitude \(M_{\mathrm{UV}}\) (7 ) can be converted to the UV luminosity \(L_{\mathrm{UV}}(M_{\rm halo}, z)\) using Eq. (6 ).

3.2 Method for performing MCMC analysis↩︎

The number distribution of PBHs within dark matter halo follows a Poisson distribution with mean \(\lambda\), \[\label{Poisson} \lambda = \frac{n_0}{n_{\mathrm{halo}}} = \frac{f_{\mathrm{PBH}} \Omega_{dm} \rho_c}{n_{\mathrm{halo}} \langle m_{\mathrm{PBH}} \rangle} = \frac{f_{\mathrm{PBH}} \Omega_{dm} \rho_c}{n_{\mathrm{halo}} M_{\mathrm{PBH}}} e^{\sigma^2},\tag{10}\] where⁷ \[n_{\mathrm{halo}}=\int_{M_{min}(z)}^{\infty} \frac{dn}{dM_h}(M_h, z) dM_h \,\] is the comoving number density of halos.

In light of the mass function 3 of PBHs and the luminosity model 4 , we can have the luminosity distribution of a single PBH whose number density follows a Poisson distribution, see Eq.(10 ). The sum of luminosities of all PBHs follows a compound Poisson distribution. Therefore, the total luminosity function contributed by stars and PBHs is the convolution of the stellar luminosity distribution and the compound Poisson distribution of PBHs.

Here, we will perform the Markov Chain Monte Carlo (MCMC) analysis with observational data. The emcee package [85] was employed to sample the parameters \(\log M_{\mathrm{PBH}}\), \(\log f_{\mathrm{PBH}}\), and \(\sigma\) while fixing the Eddington ratio \(\lambda_E=0.01\), \(0.1\) and \(1\), respectively.

Accounting for the range of the luminosity function, the log-likelihood function is defined as: \[\ln \mathcal{L} = -\frac{1}{2} \sum_{i} \left( \frac{\log \phi_{\mathrm{model},i} - \log \phi_{\mathrm{obs},i}}{\sigma_{\log,i}} \right)^2\] where \(\phi_{\mathrm{model},i}\) is the model prediction at magnitude \(M_{UV,i}\), \(\phi_{\mathrm{obs},i}\) is the observed value, and the logarithmic error is \(\sigma_{\log,i} = \sigma_i / (\phi_i \ln 10)\).

3.3 Results↩︎

The relevant model parameters are \(f_{\rm PBH}\), \(M_{\rm PBH}\), \(\lambda_{E}\), and \(\sigma\). In Table 2 and Figure 2, we present the posterior distributions of corresponding parameters.

The mass function 2 of PBHs has an additional mass correction factor of \(m^{-1.5}\) compared to the standard log-normal distribution. This results in a faster decay of the probability distribution at the high-mass end, necessitating relatively larger PBH mass parameters \(M_{\mathrm{PBH}}\) (with \(\lambda_E = 0.1\), \(M_{\mathrm{PBH}} = 10^{7.26} M_\odot\)).

In Figure 2, a positive correlation between peak values of \(f_{\mathrm{PBH}}\) and \(M_{\mathrm{PBH}}\) is revealed, which physically originates from the requirement that the number density of PBHs is approximately constant. This conserved number density ensures statistically significant enhancement of the UV LF. In the vicinity of peak, \(M_{\mathrm{PBH}}\) and \(f_{\mathrm{PBH}}\) exhibit a negative correlation, while a corresponding negative correlation exists between peak \(M_{\mathrm{PBH}}\) and \(\lambda_E\): \[\label{M-lambda} M_{\mathrm{PBH}} \lambda_E \approx 10^{6.26} M_\odot.\tag{11}\] This arises because, under relatively constant number density conditions, \(L_{\mathrm{PBH}}\) maintains a positive correlation with both \(M_{\mathrm{PBH}}\) and \(f_{\mathrm{PBH}}\). It can be also seen that lower PBH masses necessitate broader mass distribution widths to offer sufficient massive PBHs that impact the observational data range, notably the width parameter \(\sigma\) shows nearly no changes across different Eddington ratios.

Therefore our result (11 ) suggests that the PBHs must be supermassive for \(\lambda_E<1\). In Figure 3, we compare the UV LF of the \(\Lambda\)CDM+SMPBH model with the standard \(\Lambda\)CDM scenario at \(z=11\). The results show that the inclusion of SMPBHs, especially with accretion emission, leads to a significantly improved fit to the observed data at the bright end [6], [14]. This demonstrates that our SMPBHs model can account for the excess of luminous galaxies detected by JWST.

4 DISCUSSION↩︎

In this work, we explore how SMPBHs in different mass ranges, sourced by inflationary vacuum bubbles, can influence the UV LF of high-redshift galaxies. It is found that the high-redshift galaxy UV LF can be notably boosted by the intrinsic UV emission of SMPBHs (\(M\gtrsim 10^6M_\odot\)), with the bright-end excess primarily driven by sub-Edington accreting SMPBHs, see (11 ). The bright-end enhancement observed in Figure 3 can originate from the high-mass tail of the contribution (2 ) of SMPBHs, which requires the distribution to have a relatively large width parameter \(\sigma\).

However, our discussion has not yet considered the possible mass growth of PBHs with redshift. It seems that if this effect is taken into account, the corresponding initial mass for PBHs would be relatively smaller, however, for SMPBHs the case might be different, the SMPBHs \((\gtrsim 10^6M_\odot)\) can follow the self-similar accretion, so that their mass growth during pregalactic era might be negligible, e.g.[38]. According to our results, SMPBHs sourced by the inflationary vacuum bubbles are very likely to be regarded as SMBHs when the redshift \(z\lesssim 20\), which can accelerate seeding high-redshift galaxies and their baryonic content. It will be interesting to further explore precise constraints on the mass distribution parameters (such as the characteristic mass \(M_{PBH}\) and width \(\sigma\)) of SMPBHs with more upcoming observations from JWST.

5 Definitions of normalized PBH mass distributions in different references↩︎

6 Influence of the PBH mass distribution on the power spectrum↩︎

The influence of the Poisson effect [95] induced by PBHs with a monochromatic mass distribution on the matter power spectrum is [96]: \[\label{eq:P95tot} P_{\mathrm{tot}} = P_{\Lambda \mathrm{CDM}} + P_{\mathrm{PBH}},\tag{12}\] with \[\label{eq:P95PBH} P_{\mathrm{PBH}} = \frac{f_{\mathrm{PBH}}^2}{n_{\mathrm{PBH}}} D_0^2 \, \Theta(k_{\mathrm{cut}} - k) \qquad for \quad k>k_{eq}\tag{13}\] The PBH contribution depends on the PBH dark matter fraction \(f_{\mathrm{PBH}}\) and the comoving number density \[n_{\mathrm{PBH}} = \frac{f_{\mathrm{PBH}}}{M_{\mathrm{PBH}}} \left( \Omega_{\mathrm{m}} - \Omega_{\mathrm{b}} \right) \frac{3 H_0^2}{8 \pi G}\] where \(M_{\mathrm{PBH}}\) is the mass of PBHs, \(\Omega_{\mathrm{m}}\) represents the total matter density parameter (including both dark and baryonic matter), \(\Omega_{\mathrm{b}}\) denotes the baryonic matter density parameter, \(H_0\) is the Hubble constant.

The linear growth function \(D_0=D(a=1)\) is approximated using the analytic form provided by [96]. \[D_{\mathrm{PBH}}(a) \approx \left(1 + \frac{3 \gamma a}{2 a_{-} a_{\mathrm{eq}}} \right)^{a_{-}}, \quad \gamma = \frac{\Omega_{\mathrm{m}} - \Omega_{\mathrm{b}}}{\Omega_{\mathrm{m}}}\]

\[a_{-} = \frac{1}{4} \left( \sqrt{1 + 24 \gamma} - 1 \right),\] where \(a_{\mathrm{eq}} \approx \frac{1}{3400}\) denotes the scale factor at the epoch of matter-radiation equality. At the epoch of matter-radiation equality, Poisson fluctuations become significant on scales \(k > k_{\rm eq} \approx 0.01\;\mathrm{Mpc}^{-1}\,.\)[72], And a small-scale cutoff is applied with a step function \(\Theta(k_{\mathrm{cut}} - k)\) for \(k > k_{cut} = \left( 2 \pi^2 n_{\text{PBH}}/{f_{\text{PBH}}} \right)^{1/3}\), following [43]. The power spectrum in Eq. (12 ) is normalized by matching to \(\sigma(R = 8\, h^{-1}\mathrm{Mpc}) = \sigma_8\)[72].

In the case with extended mass spectrum, the mass fraction function is correspondingly expressed as: \[f(m) = \frac{m}{ \Omega_{dm} \rho_{c}} \frac{dn}{dm}\] The relation between the mass of PBHs and the wavenumber of the power spectrum can be naturally written as: \[m(k) \equiv 2\pi^{2} \Omega_{dm} \frac{\rho_{c}}{k^{3}}\] The influence of the mass function (3 ) of SMPBHs originating from inflationary vacuum bubbles on the matter power spectrum can be derived from Eq.@eq:eq:P95PBH as follows: \[\label{eq:P95PBH95mf} \begin{align} P_{\mathrm{PBH}}(k) &= \int_{0}^{m(k)} \frac{f(m) m}{\Omega_{dm} \rho_{c}} D^{2} dm \\ &= \frac{D^{2}}{\Omega_{dm}^{2} \rho_{c}^{2}} n_{0} \int_{0}^{m(k)} m^{2} \psi_{3}(m) dm \\ &= \frac{D^{2}}{\Omega_{dm} \rho_{c}} \frac{f_{\mathrm{PBH}}}{M_{c}} e^{\sigma^{2}} \cdot \frac{M_{c}^{2} e^{-\sigma^{2}}}{2}\left[1+\operatorname{erf}\left(\frac{\ln \left(m / M_{c}\right)-\frac{\sigma^{2}}{2}}{\sqrt{2} \sigma}\right)\right] \\ &= \frac{D^{2} f_{\mathrm{PBH}}}{2 \Omega_{dm} \rho_{c}} M_{c}\left[1+\operatorname{erf}\left(\frac{\ln \left(m / M_{c}\right)-\frac{\sigma^{2}}{2}}{\sqrt{2} \sigma}\right)\right]. \end{align}\tag{14}\] Mathematically, it can be rigorously verified that (14 ) reduces to the monochromatic mass spectrum in the limit \(\sigma \to 0\).

References↩︎