Observational Evidence for Cosmological Coupling of Black Holes and its Implications for an Astrophysical Source of Dark Energy

Astrophysical black holes (BHs), with masses spanning a few to several billion solar masses, are found in systems ranging from stellar binaries to supermassive BHs (SMBHs) in the centers of galaxies. Observations of gravitational waves from binary BH mergers (Abbott et al. 2019, 2021) and of supermassive BHs by the Event Horizon Telescope Collaboration et al. (2019) and Akiyama et al. (2022), have shown excellent consistency with the Kerr (1963) solution on timescales from milliseconds to hours, and spatial scales of up to milliparsecs. BHs are thus established as a universal phenomenon, across at least 10 orders of magnitude in mass.

Existing models for astrophysical BHs are necessarily provisional. They feature singularities, horizons, and unrealistic boundary conditions ²⁷ (e.g., Visser 2009). Though singularities and horizons are of theoretical interest (e.g., Harlow 2016), the Kerr solution reduces to flat spacetime at spatial infinity. This is incompatible with our universe, which is in concordance with a perturbed Robertson Walker (RW) cosmology to sub-percent precision (e.g., Aghanim et al. 2020; Dodelson & Schmidt 2020). Thus, regardless of singularities and horizons, Kerr is only appropriate for intervals of time short compared to the reciprocal expansion rate of the universe, and can only be consistently interpreted as an approximation to some more general solution.

Efforts to construct a BH model in general relativity (GR) with realistic RW boundary conditions have been ongoing for nearly a century, but have met with limited success. Early work by McVittie (1933) generalized the Schwarzschild solution to arbitrary RW spacetimes. Nolan (1993) constructed a non-singular interior for this solution, and progress has been made in understanding its horizon/causal structure (e.g., Kaloper et al. 2010; Lake & Abdelqader 2011; Faraoni et al. 2012; da Silva et al. 2013). Faraoni & Jacques (2007) constructed solutions featuring dynamical phenomena such as horizons that comove with the universe's expansion, evolution of interior energy densities and pressures, and time-varying mass. These solutions are significant, because they show how heuristic application of Birkhoff's theorem in cosmological settings can fail in the presence of strong gravity ²⁸ (see Lemaître 1931; Einstein & Straus 1945; Callan et al. 1965; Peebles 1993). Time-varying mass in particular has been studied by Guariento et al. (2012) and Maciel et al. (2015), but its interpretation remains largely unexplored. All of these solutions, however, are incompatible with Kerr on short timescales because they do not spin. A BH solution that satisfies observational constraints at small and large scales simultaneously has yet to be found.

Progress on these problems in GR has become possible with advances that resolve a long-standing ambiguity in Friedmann's equations (e.g., Ellis & Stoeger 1987). In RW cosmology, the metric is position-independent and has no preferred directions in space. In order for Einstein's equations to be consistent, the stress-energy must share these properties. Einstein's equations, however, give no prescription for converting the actual, position-dependent, distribution of stress-energy observed at late times into a position-independent source. Croker & Weiner (2019) resolved this averaging ambiguity, showing how the Einstein–Hilbert action gives the necessary relation between the actual distribution of stress-energy and the source for the RW model.

A consequence of this result is that relativistic material, located anywhere, can become cosmologically coupled to the expansion rate. This has implications for singularity-free BH models, such as those with vacuum energy interiors (e.g., Gliner 1966; Dymnikova 1992; Chapline et al. 2002; Mazur & Mottola 2004; Lobo 2006; Mazur & Mottola 2015; Dymnikova & Galaktionov 2016; Posada 2017; Beltracchi & Gondolo 2019; Posada & Chirenti 2019). The stress-energy within BHs like these, and therefore their gravitating mass, can vary in time with the expansion rate. The effect is analogous to the cosmological photon redshift, but generalized to timelike trajectories.

The presence or absence of cosmologically coupled mass in BHs strongly constrains observationally viable BH solutions. In general, the way in which a BH's mass M changes in time depends on the BH model. Croker et al. (2021) give a parameterization of the effect in terms of the RW scale factor a,

$$\begin{eqnarray}&&M(a)=M({a}_{i}){\left(\displaystyle \frac{a}{{a}_{i}}\right)}^{k}\qquad a\geqslant {a}_{i},\end{eqnarray} \tag{ 1 }$$

where a_i is the scale factor at which the object becomes cosmologically coupled and k ≥ 0 is the cosmological coupling strength. ²⁹ The Nolan (1993) solution can be regarded as cosmologically coupled with k = 0 because its stress-energy evolves such that the mass remains fixed throughout the cosmological expansion. Vacuum energy interior solutions with cosmological boundaries have been predicted to produce k ∼ 3, which is the maximum value for causal material with positive energy density (Croker & Weiner 2019).

Cosmologically coupled mass change allows for experimental distinction between singular and non-singular BHs, complementing constraints from short-timescale data (e.g., Sakai et al. 2014; Cardoso et al. 2016; Uchikata et al. 2016; Yunes et al. 2016; Cardoso & Pani 2017; Chirenti 2018; Konoplya et al. 2019; Maggio et al. 2020). Observing cosmologically coupled mass change, however, is challenging. Between an initial scale factor a_i and a final one a_f , mass evolution only becomes apparent when ${({a}_{f}/{a}_{i})}^{k}\gg 1$. For example, an observational test of cosmological mass change at z ≲ 3 requires a population of BHs whose masses can be tracked across at least a gigayear, and in which accretion or merging can be independently estimated.

In this paper, we perform a direct test of BH mass growth due to cosmological coupling. A recent study by Farrah et al. (2023) compares the BH masses M_BH and host galaxy stellar masses M_* of "red-sequence" elliptical galaxies over 6–9 Gyr, from the current epoch back to z ∼ 2.7. The study finds that the BHs increase in mass over this time period by a factor of 8–20× relative to the stellar mass. The growth factor depends on redshift, with a higher factor at higher redshifts. Because SMBH growth via accretion is expected to be insignificant in red-sequence ellipticals, and because galaxy–galaxy mergers should not on average increase SMBH mass relative to stellar mass, this preferential increase in SMBH mass is challenging to explain via standard galaxy assembly pathways (Farrah et al. 2023, Section 5). We here determine if this mass increase is consistent with cosmological coupling and, if so, the constraints on the coupling strength k.

We consider five high-redshift samples, and one local sample, of elliptical galaxies given by Farrah et al. (2023). For the high-redshift samples we use: two from the WISE survey (one at $\widetilde{z}=0.75$ measured with the Hβ line, and one at $\widetilde{z}=0.85$ measured with the Mg ii line), two from the SDSS (one at $\widetilde{z}=0.75$ and one at $\widetilde{z}=0.85$, with Hβ and Mg ii, respectively), and one from the COSMOS field (at $\widetilde{z}=1.6$). We then determine the value of k needed to align each high-redshift sample with the local sample in the M_BH–M_* plane. If the growth in BH mass is due to cosmological coupling alone, regardless of sample redshift, the same value of k will be recovered.

To compute the posterior distributions in k for each combination, we apply the pipeline developed by Farrah et al. (2023), which we briefly summarize. Realizations of each galaxy sample are drawn from the sample with its reported uncertainties. The likelihood function applies the expected measurement and selection bias corrections to the realizations, as appropriate for each sample. The de-biased, and so best actual estimate, BH mass of each galaxy is then shifted to its mass at z = 0 according to Equation (1) with some value of k. Using the Epps–Singleton test, an entire high-redshift realization is then compared against a realization of the local ellipticals, where BH masses are shifted to z = 0 in the same way. The result is a probability that can be used to reject the hypothesis that the samples are drawn from the same distribution in the M_BH–M_* plane, i.e., that they are cosmologically coupled at this k.

We present posterior distributions in k, for each high-redshift to local comparison, in the top row of Figure 1. The redshift dependence of mass growth translates to the same value k ∼ 3 across all five comparisons, as predicted by growth due to cosmological coupling alone. As further verification, we compute k from a comparison between high-redshift WISE and COSMOS samples. This comparison requires no BH bias corrections. We find a consistent value of $k={2.96}_{-1.46}^{+1.65}$. Combining the results from each local comparison gives

$$\begin{eqnarray}&&k={3.11}_{-1.33}^{+1.19}\qquad (90 \% \ \mathrm{confidence}),\end{eqnarray} \tag{ 2 }$$

which excludes k = 0 at 99.98% confidence, equivalent to >3.9σ observational exclusion of the singular Kerr interior solution.

Zoom In Zoom Out Reset image size

3.1. Implications

3.2. Ω_Λ from the Cosmic Star Formation History

If k ∼ 3 BHs contribute as a cosmological dark energy species, a natural question is whether they can contribute all of the observed Ω_Λ. We test this by assuming that: (1) BHs couple with k = 3, consistent with our measured value; (2) BHs are the only source for Ω_Λ, and (3) BHs are made solely from the deaths of massive stars. Under these assumptions, the total BH mass from the cosmic history of star formation (and subsequent cosmological mass growth) should be consistent with Ω_Λ = 0.68.

In Appendix A we construct a simple model of the cosmic star formation rate density (SFRD) that allows exploring combinations of stellar production rate, stellar initial mass function (IMF), and accretion history. Figure 2 displays models that produce the Planck measured value of Ω_Λ = 0.68 (Aghanim et al. 2020) with the indicated IMF and Eddington ratio λ_e . Any monotonically decreasing path inside the filled region produces Ω_Λ = 0.68 for some observationally viable IMF and accretion history, consuming at most 3% of baryons. This baryon consumption is compatible with the results of Macquart et al. (2020), who find that Ω_b at low redshift agrees with Ω_b inferred from the big bang to within 50%.

Zoom In Zoom Out Reset image size

It follows from Equation (1) that cosmological coupling in BHs with k = 3 will produce a BH population with masses >10² M_⊙. If these BHs are distributed in galactic halos, they will form a population of MAssive Compact Halo Objects (MACHOs). In Appendix B, we consider the consistency of SFRDs in Figure 2 with MACHO constraints from wide halo binaries, microlensing of objects in the Large Magellanic Cloud, and the existence of ultra-faint dwarfs (UFDs). We conclude that non-singular k = 3 BHs are in harmony with MACHO constraints while producing Ω_Λ = 0.68, driving late-time accelerating expansion.

Further tests of cosmological coupling in BHs are essential to confirm or refute our proposal. We list some examples below, highlighting possible tensions.

4.1. Signatures in the Cosmic Microwave Background

4.2. Strong Lensing of γ-Ray Bursts

4.3. Stellar Mass BH–BH Merger Rates

4.4. Direct Measurement from Orbital Period Decay

4.5. Stellar Mass BHs and Stellar Evolution

4.6. Theoretical Considerations

4.7. Validation of Preferential Growth of SMBHs

4.8. Quasars at z > 6 and the SMBH Mass Function

Realistic astrophysical BH models must become cosmological at large distance from the BH. Non-singular cosmological BH models can couple to the expansion of the universe, gaining mass proportional to the scale factor raised to some power k. A recent study of SMBHs within elliptical galaxies across ∼7 Gyr finds redshift-dependent 8–20× preferential BH growth, relative to galaxy stellar mass. We show that this growth excludes decoupled (k = 0) BH models at 99.98% confidence. Our measured value of $k={3.11}_{-1.33}^{+1.19}$ at 90% confidence is consistent with vacuum energy interior BH models that have been studied for over half a century. Cosmological conservation of stress-energy implies that k = 3 BHs contribute as a dark energy species. We show that k = 3 stellar remnant BHs produce the measured value of Ω_Λ within a wide range of observationally viable cosmic star formation histories, stellar IMFs, and remnant accretion. They remain consistent with constraints on halo compact objects and they naturally explain the "coincidence problem," because dark energy domination can only occur after cosmic dawn. Taken together, we propose that stellar remnant k = 3 BHs are the astrophysical origin for the late-time accelerating expansion of the universe.

We thank the referees for very helpful reports. We thank M. Valluri (U. Michigan) for suggesting the bias-free high-redshift crosscheck. We thank the David C. and Marzia C. Schainker Family for their financial support of required computations. M.Z. is supported by NASA through the NASA Hubble Fellowship grant HST-HF2-51474.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. G.T. acknowledges support through DoE Award DE-SC009193. V.F. is supported by the Natural Sciences & Engineering Research Council of Canada (grant 2016-03803). J.A. acknowledges financial support from the Science and Technology Foundation (FCT, Portugal) through research grants PTDC/FIS-AST/29245/2017, UIDB/04434/2020 and UIDP/04434/2020.

The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.

Here we construct a simple model to establish whether plausible SFRDs lead to Ω_Λ ≃ 0.7, if massive stellar collapse produces k = 3 BHs. For the overall form of the SFRD, we adopt the model of Madau & Fragos (2017) at z ≤ 4. Between z = 4 and cosmic dawn at z = 25 we adopt a power law for the SFRD with exponent α < 0, matched continuously to the z ≤ 4 SFRD. In units of M_⊙ Mpc⁻³ yr⁻¹,

$$\begin{eqnarray}{\dot{\rho }}_{* }(z,\alpha ):= \left\{\begin{array}{ll}0.01\displaystyle \frac{{\left(1+z\right)}^{2.6}}{1+{\left[(1+z)/3.2\right]}^{6.2}} & z\leqslant 4\\ {\left(z/4\right)}^{\alpha }{\dot{\rho }}_{* }(\alpha ) & 4\lt z\leqslant 25\\ 0 & z\gt 25\end{array}\right..\end{eqnarray} \tag{ A1 }$$

where ${\dot{\rho }}_{* }$ is the time rate of production of stellar mass. The choice of varying the power-law slope at z > 4 is motivated by the disparities in observed ${\dot{\rho }}_{* }$ over 4 < z < 10. Estimates from IR (Rowan-Robinson et al. 2018) and GRBs (Kistler et al. 2009) can be 1–2 dex higher than estimates from UV-dropouts (Harikane et al. 2022a; Bouwens et al. 2022; Donnan et al. 2022). Varying α at z > 4 allows us to encompass SFRDs consistent with all extant data.

To account for variation in the stellar IMF with redshift, we divide the SFRD into a Population III epoch at z > 7 (e.g., Inoue et al. 2014, but see Liu & Bromm 2020) and a standard epoch, characterized by a Kroupa IMF at z ≤ 7. We consider scenarios where the Population III epoch is characterized by either a Kroupa IMF or the top-heavy Harikane et al. (2022a) IMF ∝M^−2.35 for 50 ≤ M ≤ 500, and zero elsewhere. Finally, we use the redshift-dependent, mean-metallicity model of Madau & Fragos (2017), truncated at Z_⊙, to enable use of standard metallicity-dependent zero-age main sequence (ZAMS) mass to remnant mass, delayed, core-collapse supernovae prescriptions by Fryer et al. (2012).

For simplicity, we apply post-remnant formation accretion, with a fixed Eddington ratio λ_e , to all Population III BHs over a duration t_e ,

$$\begin{eqnarray}&&M({t}_{i}+{t}_{e})\propto M({t}_{i})\exp \left({\lambda }_{e}{t}_{e}/\epsilon \right),\end{eqnarray} \tag{ A2 }$$

brief enough that the impact on instantaneous rate from cosmological coupling of accreted mass can be neglected. Because t_e , λ_e , and are degenerate in Equation (A2), there is effectively only one parameter. Farrah et al. (2023) use simulation studies of galaxy mergers to estimate t_e = 44 ± 22 Myr at λ_e > 1 from observations of 21 IR-luminous galaxy mergers at z < 0.3. This value is consistent with ∼10 Myr estimates by Madau et al. (2014) for high-z BH seeds, so we adopt t_e = 20 Myr. The uncertainty in t_e is large, e.g., Safarzadeh & Haiman (2020) find Population III stellar BH accretion across timescales ∼0.5 Myr. Tortosa et al. (2022) measure λ_e ≃ 472 for a Seyfert 1 galaxy, while simulations by Inayoshi et al. (2016) find stable accretion onto isolated Population III BHs at λ_e ≃ 5000. As we apply this accretion to all Population III BHs in our model, we consider models within a more conservative λ_e ≤ 30. Given the large range of plausible t_e and λ_e , we fix the efficiency = 0.16, roughly 30% of the Kerr BH limit (e.g., Bardeen 1970).

To compute viable SFRDs, we use a Monte Carlo approach. Because high-mass remnants formed at high redshift can easily dominate Ω_Λ, care must be taken to sufficiently sample these tails of the distributions. We divide the domains for IMF and SFRD distributions into distinct windows such that, in these windows, the distributions are guaranteed to give Poisson errors ≤1% in the lowest probability bin. The total collection of draws is eventually reweighted by relative areas under the windowed IMFs and Equation (A1). In each window, we draw 10⁶ ZAMS stellar masses from the appropriate IMFs and 10⁶ redshift-dependent metallicities. Given the birth masses and redshifts, we approximate a redshift of stellar death using ${t}_{\mathrm{life}}\sim {10}^{10}{(M/{M}_{\odot })}^{-\tau }\,\mathrm{yr}$, where τ = 2.5 (e.g., Harwit 2006). We discard draws that live beyond z = 0, and then determine collapsed remnant masses. The ZAMS mass to remnant mass model we adopt assumes that all stars drawn are single stars. Consistent with Kalogera & Baym (1996), we regard any remnants with mass >2.7M_⊙ as BHs. To convert the resulting population of BHs into a cosmological density, we regard this drawn population as residing within one co-moving Mpc³. We first reweight all draws by the probability of having come from their respective windows. We then convert from this reweighted mass to a predicted mass by rescaling ${\dot{\rho }}_{* }(z=0)$ by the total mass in long-lived stars, neutron stars, and initial BH remnant mass, as described by Madau & Fragos (2017). We compute ${\dot{\rho }}_{* }(z=0)$ by integrating Equation (A1) in time, and scaling by an IMF-appropriate gas return fraction 1 − R (e.g., R = 0.39 given by Madau & Fragos 2017 for Kroupa, and 1 for Harikane). We then divide the summed predicted mass in BHs at z = 0 by the critical density today $3{H}_{0}^{2}/8\pi G$ to get Ω_Λ.

Here we establish that a k ∼ 3 BH population sufficient to produce Ω_Λ is consistent with constraints on MACHOs. In Figure 3, we display the contribution to halo density from k = 3 BH masses at z = 0 computed in the explicit SFRD models shown in Figure 2. Consistent with the discussion in Section 4.1, we have assumed a uniformly dispersed population in the computation of Figure 3. We adopt a uniform Milky Way (MW) halo density equal to the median dark matter density at the solar system 8.8 ± 0.5 × 10⁻³ M_⊙/pc³ as measured by Cautun et al. (2020). This is conservative relative to a mean halo density 1.7 × 10⁻² M_⊙/pc³, inferred from 1.37 × 10¹¹ M_⊙ within the MW at <20 kpc as measured by Posti & Helmi (2019).

Zoom In Zoom Out Reset image size

In Figure 3 we display limits on MACHOs from microlensing (e.g., Blaineau et al. 2022), wide halo binary disruption (e.g., Monroy-Rodríguez & Allen 2014; Tyler et al. 2022), and UFD disruption (e.g., Brandt 2016). We approximate the UFD constraint from Brandt (2016) as ∝1/m following Binney & Tremaine (2011, Equation (7.104)) and convert from the UFD halo density assumed by Brandt (2016, Figure 4) to our MW halo density scale by a multiplicative factor 0.3/0.0088.

All BH density contributions lie below the microlensing and wide binary limits by several orders of magnitude. High accretion models with a Kroupa IMF show some tension with UFD constraints, which begin to have discriminatory power in our scenario. Comparing the two top-heavy IMF models, the effect of accretion is to redistribute mass density from the IMBH range into the SMBH range. The top-heavy IMF with no accretion model shows how a top-heavy IMF acts to populate the IMBH region.

The aforementioned constraints do not account for comoving BH mass density decrease with increasing redshift (Equation (1)) for k > 0. The constraints as presented are thus overly stringent because the probability of compact object interactions is reduced at earlier times. The UFD galaxy and wide halo binary constraints are most impacted, as they consider the stability of astrophysical phenomena over gigayear timescales. A thorough analysis that incorporates decreasing comoving density of compact objects further back in time is the subject of future work.