Constraints on Cosmologically Coupled Black Holes from Planck 2018 and Other Cosmological Probes
July 04, 2025
Abstract
Cosmologically coupled black holes (CCBHs) are alternative black hole models whose masses evolve as \(M \propto a^3\) on cosmological scales. This characteristic suggests that CCBHs could contribute to the accelerated expansion of the Universe. In this paper, we consider a CCBH model in which the cosmological constant is effectively induced, while the baryonic mass is conserved within conventional black holes. This model is motivated by the theoretical framework of Schwarzschild–de Sitter black holes. Assuming that the accelerated cosmic expansion is caused by CCBHs, we perform a cosmological parameter estimation using datasets including Planck 2018 CMB, CMB lensing, BAO, and supernovae. The analysis reveals notable shifts in cosmological parameters, such as \(H_0 = 72.24^{+0.34}_{-0.35} \, \mathrm{km/s/Mpc}\), compared to the standard \(\Lambda \mathrm{CDM}\). My \(H_0\) constraint is consistent with the value \(H_0 = 73.04\pm1.04 \, \mathrm{km/s/Mpc}\) reported by SH0ES within \(1 \sigma\). However, the overall fit to the data worsens, with a total \(\chi^2 = 2884.12\) for the CCBH model, compared to \(\chi^2 = 2836.12\) for the \(\Lambda\)CDM model. We show that the effect of cosmological coupling is suppressed by a factor of \(10^{-16}\) at \(\sim\)pc scales, rendering it negligible compared to the standard black hole mass in local astrophysical phenomena, although the CCBH model can explain the accelerated expansion.
Constraints on Cosmologically Coupled Black Holes from Planck 2018 and Other Cosmological Probes
Shintaro K. Hayashi
Graduate School of Science, Division of Particle and Astrophysical Science, Nagoya University, Furocho, Chikusa-ku, Nagoya, Aichi 464-8602, Japan
1 Introduction↩︎
Recently, several observational studies have suggested that the mass evolution of black holes (BHs) follows \(M \propto a^3\) [1], [2], although this interpretation remains debated [3]–[5]. Black holes exhibiting this behavior are referred to as "Cosmologically Coupled Black Holes" (CCBHs). Since the energy density of non-relativistic matter scales as \(\rho \propto a^{-3}\), CCBHs effectively maintain a constant energy density on cosmological scales. This opens up the intriguing possibility that BHs could contribute to the dark-energy component.
In the standard \(\Lambda\)CDM model, dark energy is assumed to be a constant energy density, typically modeled as a cosmological constant (\(\Lambda\)). If the number of BHs is conserved, then CCBHs would collectively behave as a constant energy density component. However, the total energy density of BHs evolves over cosmic time, due to changes in their number density, particularly a sharp increase at high redshift when the BH formation is most active.
Observations of baryon acoustic oscillations (BAO) by the Dark Energy Spectroscopic Instrument (DESI) suggest that the dark energy density may vary over time, as modeled by the \(w_0w_a\)CDM parameterization [6], [7].
In Ref. [8], the authors demonstrated that the CCBH model can effectively reproduce the \(w_0w_a\)CDM evolution in the DESI BAO data. They further extended their analysis using Planck CMB data in Ref. [9]. In these studies, the CCBH framework is interpreted as a type of baryon–dark energy interaction model, with the interaction rate assumed to be proportional to the star formation rate [10]–[13].
From a theoretical perspective, the CCBH model incorporates black hole geometries into an expanding universe described by the Friedmann–Robertson–Walker (FRW) spacetime [14]. This extension is also relevant to efforts to resolve black hole singularities [15]–[17]. Such configurations are related to the Schwarzschild–de Sitter (SdS) solution, a static black hole embedded in a de Sitter universe. While the SdS solution features a constant cosmological term, CCBH models instead allow for a spatially varying, effectively coupled term.
These theoretical considerations suggest that the effects of cosmological coupling are negligible on local observables, such as the BH–bulge mass relation and gravitational wave signals from binary mergers. Furthermore, unlike the models in [8], [9], the baryon energy density is assumed to be conserved throughout cosmic evolution. In studies of the black hole number density, simulations incorporating binary evolution and merger histories have also been conducted [18], [19]. These results provide a more direct estimation of its redshift evolution than extrapolations from the star formation rate, making them particularly well-suited for use in the CCBH framework.
In this paper, we assume the evolution of the CCBH energy density using the simulated BH mass function in Ref. [18], and we test the CCBH scenario using Planck 2018 CMB data, DESI BAO DR2, and Union3 supernova. We also examine possible astrophysical effects of cosmological coupling under the constraints imposed by the combination of the cosmological dataset.
In Sec. 2.1, we introduce the specific CCBH model adopted in this work. Since its energy density depends on the BH density, we describe the mass function and its evolution in Sec. 2.2, following Ref. [18]. The constraints on cosmological parameters obtained via Markov Chain Monte Carlo (MCMC) analysis are presented in Sec. 3. Finally, Sec. 4 provides discussion and concluding remarks.
2 Methods↩︎
2.1 Cosmologically coupled black holes↩︎
In general relativity, black holes are described as vacuum solutions to the Einstein equation. When a cosmological constant \(\Lambda\) is introduced, the corresponding vacuum solution becomes the Schwarzschild–de Sitter (SdS) metric. This metric describes a static black hole embedded in a de Sitter expanding background, where \(\Lambda\) is treated as a constant energy density that is uniform throughout spacetime.
In this paper, we begin with the SdS black hole, whose line element is given by \[ds^2 = - f(r) \, dt^2 + f(r)^{-1} \, dr^2 + r^2 \left( d \theta^2 + \sin^2 \theta \, d \phi^2 \right),\] where \[f(r) = 1 - \frac{2GM}{r} - \frac{\Lambda}{3} r^2, \\None\] and \(M\) is the conventional mass of BH, originating from the matter content involved in the BH formation.
In the cosmologically coupled black hole (CCBH) model, the cosmological constant is not regarded as a universal constant but is instead sourced by the mass distribution of black holes. In this framework, the effective mass profile of a black hole can be written as \[\label{mass95eff} M_\mathrm{eff}(r) = M + C r^3,\tag{1}\] and \[f(r) = 1 - \frac{2GM_\mathrm{eff}\left(r\right)}{r},\] where \(C\) is a constant with dimensions of energy density. Physically, this implies that the cosmological term contributes a position-dependent mass that becomes dominant at large distances.
In the limit \(r \to \infty\), the additional mass term dominates, leading to a constant mass density associated with the black hole. Since the volume enclosed by a comoving radius evolves as \(a^3\), this suggests that the effective mass of a CCBH evolves as \(M_\mathrm{eff} \propto a^3\) on cosmological scales. Therefore, a collection of such black holes can collectively give rise to an energy density that mimics a dark energy component.
Assuming that the cosmological constant arises from a sum of contributions from individual black holes, we write \[\begin{align} \Lambda \; (\text{or } \rho_\mathrm{CCBH}) &= \sum_i C_i \\ &\simeq C_\mathrm{mean} \times N_\mathrm{BH}, \end{align}\] where \(C_i\) denotes the contribution to the effective cosmological constant of the \(i\)-th black hole and \(N_\mathrm{BH}\) is the number of black holes. Although \(C_i\) may in general depend on the properties of individual black holes (e.g., mass), we approximate all contributions by a mean value \(C_\mathrm{mean}\) corresponding to a representative black hole.
The total energy density of CCBHs can then be approximated as \[\rho_\mathrm{CCBH} \simeq C_\mathrm{mean} \, n_\mathrm{BH} \, r^3,\] where \(n_\mathrm{BH}\) is the number density of black holes, and the usual rest-mass contribution is neglected. This number density can be inferred from the black hole mass function obtained in previous studies [18], [19].
To relate this expression to the comoving cosmological framework, we introduce comoving coordinates \(\bar{r} = r/a\) and define \(C' \equiv C \bar{r}^3 / M\), where \(M\) is a typical black hole mass. Then, the CCBH energy density becomes \[\label{rho95ccbh} \rho_\mathrm{CCBH} \simeq C' \, \rho_{\mathrm{tot,BH}} \, a^3,\tag{2}\] where \(\rho_{\mathrm{tot,BH}}\) is the total BH mass density.
Since the region affected by \(C\) is expected to expand at the speed of light, it may induce additional inhomogeneities or perturbations in the dark energy sector. However, in this paper, we neglect such perturbative effects and focus only on the background evolution of the CCBH energy density.
2.2 Time evolution of BH mass density↩︎
In this analysis, we adopt the black hole (BH) mass density derived from simulations based on a binary evolution code, as reported in Ref. [18]. A convenient analytical fitting function, along with its corresponding best-fit parameters, is used to describe the mass function.
As discussed in Ref. [18], the BH mass function is well approximated by a combination of a Schechter function and a Gaussian component [20], given by: \[\label{schechtergaussian} \frac{dN}{dV d \log m_\bullet} \simeq \mathcal{N} \left( \frac{m_\bullet}{\mathcal{M}_\bullet} \right)^{1 - \alpha} e^{- m_\bullet / \mathcal{M}_\bullet} + \mathcal{N}_\mathrm{G} \frac{1}{\sqrt{2\pi \sigma_\mathrm{G}^2}} \exp\left[ -\frac{ \left( \log m_\bullet - \log \mathcal{M}_{\bullet, \mathrm{G}} \right)^2 }{2\sigma_\mathrm{G}^2} \right],\tag{3}\] where \(N\) is the number of black holes, \(V\) is the comoving volume, and \(m_\bullet\) is the black hole mass. The parameters \(\mathcal{N}\), \(\mathcal{M}_\bullet\), \(\alpha\), \(\mathcal{N}_\mathrm{G}\), \(\mathcal{M}_{\bullet,\mathrm{G}}\), and \(\sigma_\mathrm{G}\) are redshift-dependent fitting coefficients.
The best-fit parameters at various redshifts, obtained in Ref. [18], are listed in Table 1. Using these results, the comoving BH mass function at several redshifts is shown in Fig. 1.
| \(z\) | \(\log \mathcal{N} \, (\mathrm{Mpc}^{-3})\) | \(\log \mathcal{M}_\bullet \, (M_\odot)\) | \(\alpha\) | \(\log \mathcal{N}_\mathrm{G} \, (\mathrm{Mpc}^{-3})\) | \(\log \mathcal{M}_{\bullet,\mathrm{G}} \, (M_\odot)\) | \(\sigma_\mathrm{G}\) |
|---|---|---|---|---|---|---|
| \(0\) | \(5.623\) | \(0.607\) | \(-3.781\) | \(2.413\) | \(2.021\) | \(0.052\) |
| \(1\) | \(5.429\) | \(0.609\) | \(-3.859\) | \(2.309\) | \(2.023\) | \(0.051\) |
| \(2\) | \(5.107\) | \(0.612\) | \(-3.914\) | \(2.064\) | \(2.024\) | \(0.051\) |
| \(4\) | \(4.344\) | \(0.634\) | \(-3.902\) | \(1.419\) | \(2.037\) | \(0.049\) |
| \(6\) | \(3.614\) | \(0.659\) | \(-3.866\) | \(0.806\) | \(2.054\) | \(0.045\) |
| \(8\) | \(2.894\) | \(0.676\) | \(-3.868\) | \(0.197\) | \(2.066\) | \(0.043\) |
| \(10\) | \(2.305\) | \(0.680\) | \(-3.884\) | \(-0.344\) | \(2.072\) | \(0.042\) |
The corresponding BH mass density can be computed from the mass function via: \[\label{rho95from95massfunc} \rho_\mathrm{BH}(z) = \int d\log m_\bullet \, m_\bullet \, \frac{dN}{dV d\log m_\bullet}(m_\bullet | z).\tag{4}\] The resulting redshift evolution of the BH mass density is shown as red points in Fig. 2.
For the estimation of the cosmological parameters, a continuous function of \(\rho_\mathrm{BH}(z)\) is required. We therefore perform a polynomial fit of the form: \[\label{rho95polynomial} \log \rho_\mathrm{BH}(z) = a z^3 + b z^2 + c z + d,\tag{5}\] where the current BH mass density is given by \(\rho_{0\mathrm{BH}} = e^d \, [M_\odot\, \mathrm{Mpc}^{-3}]\).
To select the optimal degree of the polynomial in Eq. 5 , we perform cross-validation. The validation scores are: degree 1: 0.070, degree 2: 0.182, degree 3: 0.065, degree 4: 0.135, and degree 5: 0.325. Based on this evaluation, a degree-3 polynomial provides the best balance between bias and variance. The best-fit parameters are listed in Table 2, and the resulting function is shown as the solid blue line in Fig. 2.
| \(a\) | \(b\) | \(c\) | \(d\) |
|---|---|---|---|
| \(0.00658\) | \(-0.104\) | \(-0.348\) | \(18.8\) |
2.3 Implementation into Cobaya↩︎
We use the cobaya package [21], [22]
to estimate the cosmological parameters and implement the evolution of the CCBH energy density into CAMB by modifying the dark energy component at the background level. In this study, we consider only the background evolution and neglect the
perturbations of the dark energy component.
According to Eqs. 2 and 5 , the energy density of the dark energy component is modified as: \[\begin{align} \rho_\mathrm{DE}(z) &= C' \rho_\mathrm{tot,BH}(z) \, a^3 \\ &= C' \times \frac{1}{a^3} \rho_{0\mathrm{BH}} \exp[az^3 + bz^2 + cz] \times a^3 \\ &= \rho_{0\mathrm{DE}} \, \exp[az^3 + bz^2 + cz], \end{align}\] where \(\rho_{0\mathrm{DE}}\) is the current energy density of the dark energy component induced by CCBHs. It is determined by the other cosmological parameters as \[\rho_{0\mathrm{DE}} = \Omega_\Lambda \frac{3H_0^2}{8\pi G} = \left(1 - \Omega_m\right)\frac{3H_0^2}{8\pi G}.\] From Eq. 1 , the total baryon energy density is conserved in this model, since the baryonic mass incorporated into black hole formation is conserved as the conventional black hole mass. Therefore, the number of free parameters in the CCBH model remains the same as in the standard \(\Lambda\)CDM model.
For the Markov Chain Monte Carlo (MCMC) analysis, we adopt flat priors for the cosmological parameters, as listed in Table 3. We use the Planck 2018 low-\(\ell\) and high-\(\ell\) temperature and polarization data [23], [24], and also include CMB lensing data [25], DESI DR2 BAO data [7], and Union3 supernova data [26].
MCMC sampling is terminated when convergence is achieved according to the Gelman–Rubin criterion [27], with condition \(R-1 < 0.01\) applied to all runs.
| Parameter | Prior Range |
|---|---|
| \(\Omega_{\rm b} h^2\) | [0.05, 0.1] |
| \(\Omega_{\rm c} h^2\) | [0.001, 0.99] |
| \(\tau\) | [0.01, 0.8] |
| \(100\Theta_{\rm s}\) | [0.5, 10] |
| \(n_{\rm s}\) | [0.7, 1.2] |
| \(\log_{10} A_{\rm s}\) | [2, 5] |
3 Results and Discussion↩︎
First, we present the triangle plot of the cosmological parameters obtained from the MCMC analysis in Figure 3. A selection of key parameters, \(H_0\), \(\Omega_\mathrm{m}\), \(\Omega_\Lambda\), and \(\Omega_\mathrm{c} h^2\), are highlighted in Figure 4. The CCBH model yields higher \(H_0\), \(\Omega_\Lambda\), and \(\Omega_\mathrm{c} h^2\) and a lower \(\Omega_\mathrm{m}\) compared to the \(\Lambda\)CDM model. This is due to the rapid growth of the dark energy component in the CCBH model around the scale factor \(a \sim 0.2\)–0.4. To compensate for the reduced total energy density at high redshifts, the present total energy density increases, leading to higher values of \(\Omega_\Lambda\), \(H_0\), and \(\Omega_\mathrm{c} h^2\). In contrast, \(\Omega_\mathrm{m}\) decreases because it is given by \(1 - \Omega_\Lambda\).
| Parameters | CCBH | \(\Lambda \mathrm{CDM}\) |
|---|---|---|
| \(H_0\) | \(72.24 \left( 72.32 \right) ^{+0.34}_{-0.35}\) | \(68.33 \left( 68.21 \right) ^{+0.29}_{-0.29}\) |
| \(\Omega_{\Lambda}\) | \(0.7217 \left( 0.7220 \right) ^{+0.0036}_{-0.0036}\) | \(0.6978 \left( 0.6961 \right) ^{+0.0037}_{-0.0037}\) |
| \(\Omega_{\rm m}\) | \(0.2783 \left( 0.2779 \right) ^{+0.0036}_{-0.0036}\) | \(0.3021 \left( 0.3038 \right) ^{+0.0037}_{-0.0037}\) |
| \(100\Omega_{\rm b} h^2\) | \(2.217 \left( 2.226 \right) ^{+0.012}_{-0.012}\) | \(2.253 \left( 2.253 \right) ^{+0.013}_{-0.013}\) |
| \(\Omega_{\rm c} h^2\) | \(0.1224 \left( 0.1225 \right) ^{+0.0006}_{-0.0006}\) | \(0.1178 \left( 0.1182 \right) ^{+0.0006}_{-0.0006}\) |
| \(\tau\) | \(0.0471 \left( 0.0470 \right) ^{+0.0070}_{-0.0062}\) | \(0.0604 \left( 0.0533 \right) ^{+0.0069}_{-0.0076}\) |
| \(100\Theta_{\rm{s}}\) | \(1.0407 \left( 1.0408 \right) ^{+0.0003}_{-0.0003}\) | \(1.0412 \left( 1.0412 \right) ^{+0.0003}_{-0.0003}\) |
| \(n_{\rm{s}}\) | \(0.9598 \left( 0.9591 \right) ^{+0.0032}_{-0.0032}\) | \(0.9703 \left( 0.9701 \right) ^{+0.0033}_{-0.0033}\) |
| \(\mathrm{log}_{10}A_{\rm{s}}\) | \(3.034 \left( 3.034 \right) ^{+0.014}_{-0.013}\) | \(3.053 \left( 3.040 \right) ^{+0.014}_{-0.015}\) |
Next, we examine the impact of the cosmological coupling parameter \(C\). The critical density is calculated as \(\rho_\mathrm{crit} = 3 H_0^2 / 8 \pi G\), which yields \(\rho_\mathrm{crit} \simeq 1.452 \times 10^{11} \, \mathrm{M}_\odot / \mathrm{Mpc}^3\) using the best-fit parameters. The contribution from all CCBHs is \(\Omega_\Lambda \rho_\mathrm{crit} \simeq 1.047 \times 10^{11} \, \mathrm{M}_\odot / \mathrm{Mpc}^3\).
In contrast, the current black hole mass density is given by \(\rho_{0\mathrm{BH}} = e^d\), and with \(d = 18.8\), we obtain \(\rho_{0\mathrm{BH}} \simeq 1.461 \times 10^8 \, \mathrm{M}_\odot / \mathrm{Mpc}^3\). Assuming each black hole has mass \(M = 10 \, \mathrm{M}\odot\), the implied effects of cosmological coupling per black hole is \(10^3 \, \mathrm{M}_\odot\) per \(\mathrm{Mpc}^3\), suggesting that each CCBH contributes far more than its rest mass.
However, since the CCBH contribution arises from a constant density term, its effect becomes negligible at small scales. For example, on \(\sim\)pc scales typical of binary black hole separations, the contribution of cosmological coupling is only \(10^{-15} \, \mathrm{M}_\odot\), effectively undetectable. Thus, gravitational wave observations of binary systems are unlikely to probe the CCBH effect.
Now we turn to the comparison between the CCBH and \(\Lambda\)CDM models. Although the contours in Figure 3 differ, this is expected because the CCBH parameterization is not a smooth extension of \(\Lambda\)CDM (such as \(w_0 w_a\)CDM). Since both models have the same number of free parameters, a direct comparison using total \(\chi^2\) is valid.
Figure 5 and Table 5 present the total and individual \(\chi^2\). The total \(\chi^2\) is decomposed into \(\chi^2_\mathrm{CMB}\), \(\chi^2_\mathrm{BAO}\), and \(\chi^2_\mathrm{SN}\), corresponding to Planck+lensing, DESI DR2 BAO, and supernova data. Although the fits to the CMB are comparable between the two models, the CCBH model yields significantly worse fits for BAO and supernovae data.
| \(\mathrm{CCBH}\) | \(\Lambda \mathrm{CDM}\) | |||
|---|---|---|---|---|
| \(\chi^2\) | \(2884.12 \left( 2871.13 \right) ^{+3.65}_{-7.60}\) | \(2836.12 \left( 2819.00 \right) ^{+1.63}_{-9.53}\) | ||
| \(\chi^2_\mathrm{CMB}\) | \(2796.49 \left( 2783.92 \right) ^{+7.78}_{-3.89}\) | \(2794.59 \left( 2777.11 \right) ^{+1.90}_{-9.69}\) | ||
| \(\chi^2_\mathrm{BAO}\) | \(33.55 \left( 33.01 \right) ^{+0.50}_{-1.22}\) | \(13.23 \left( 13.90 \right) ^{+0.95}_{-3.02}\) | ||
| \(\chi^2_\mathrm{SN}\) | \(54.07 \left( 54.21 \right) ^{+1.90}_{-1.91}\) | \(28.30 \left( 27.99 \right) ^{+0.58}_{-0.66}\) |
Finally, we discuss the implications for the Hubble constant. In the CCBH model, the inferred \(H_0\) is significantly higher than in the \(\Lambda\)CDM model. This aligns with earlier discussions in Ref. [8]. Local measurements of \(H_0\), such as from the SH0ES Collaboration, which reports \(H_0 = 73.04 \pm 1.04 \, \mathrm{km/s/Mpc}\) [28], have persistently exceeded the value inferred from Planck CMB observations under \(\Lambda\)CDM, a discrepancy known as the “Hubble tension” [29], [30]. If local \(H_0\) data are included, the CCBH model may provide a better overall fit than the standard \(\Lambda\)CDM.
4 Summary↩︎
In this paper, we investigated whether cosmologically coupled black holes (CCBHs) can serve as the origin of dark energy by examining their impact on cosmological observations. The cosmological data analysis for the CCBH model was performed using Planck 2018 CMB data, as well as CMB lensing, DESI DR2 BAO, and Union3 supernovae datasets. In this study, the modification of black hole geometry was translated into an effective mass, and the energy density of CCBHs was parameterized in terms of the conventional black hole mass density, as described in Eq.@eq:rho95ccbh . Additionally, the relic black hole mass density estimated from the simulated mass function in Refs.[18], [19] can be used in this context.
According to the MCMC analysis using the aforementioned datasets, several cosmological parameters showed significant deviations from their values in the \(\Lambda\)CDM model. In particular, we found a larger value of the Hubble constant \(H_0\), which may offer a possible hint toward resolving the so-called "Hubble tension." However, the overall \(\chi^2\) values obtained for the CCBH model were worse than those for \(\Lambda\)CDM, indicating a poorer fit to the combined data.
It should be noted that the current results are not yet sufficient for a definitive conclusion, as several approximations were employed throughout the analysis. These include the use of an approximate form for the mass function in Eq.@eq:schechtergaussian and the polynomial fitting of the black hole mass density in Eq.@eq:rho95polynomial .
Finally, while the contribution of black holes as CCBHs is significant on the \(\mathrm{Mpc}^3\) scale, the effect of cosmological coupling is suppressed by a factor of approximately \(10^{-16}\) compared to the conventional mass contribution on the \(\mathrm{pc}^3\) scale. Therefore, it is difficult to probe the CCBH model through astrophysical phenomena at small (galactic or binary) scales.
However, there is another possibility to constrain this cosmological coupling model. This model can be generalized as one in which a singularity generates a cosmological constant, implying that the constant emerges at the initial time of the universe. Hence, methodologies developed to discriminate among inflation models may be applicable to this cosmological coupling, but a detailed investigation is left for future work.
5 Acknowledgement↩︎
Thank you for meaningful discussions, Koichiro Nakashima and Genki Naruse.