Cosmologically Coupled Compact Objects: A Single-parameter Model for LIGO–Virgo Mass and Redshift Distributions

There is overwhelming evidence for the existence of ultrarelativistic compact objects, both as the end stages of stellar collapse and as the supermassive compact objects in the centers of galaxies. These compact objects are interpreted within the framework of General Relativity (GR), which has withstood tests on scales ranging from cosmological to the binary pulsar. A contemporary triumph in this field, following the development of numerical methods for modeling gravitational radiation from compact object mergers (Sperhake 2015), was the detection of such mergers by gravitational-wave observatories (Abbott et al. 2016). The LIGO Scientific Collaboration and the Virgo Collaboration (LVC) have currently published over 50 confident detections of compact binary mergers, the majority being consistent with Kerr black holes (BHs; Abbott et al. 2019, 2021a, 2021c).

The mass spectrum of these objects is less straightforward to interpret. It was expected (e.g., Bethe & Brown 1998) that isolated-binary stellar evolution would be a dominant channel for synthesizing binary black hole mergers visible by LIGO–Virgo. In this scenario, massive progenitor stars in a binary system proceed through one or more phases of stable and/or unstable mass transfer during their evolution that can significantly tighten the orbit of the resultant black hole binary. Simulations investigating the population properties of such systems (e.g., Dominik et al. 2012) predicted a relatively narrow mass spectrum near ∼8 M_⊙ at solar metallicity and a broader spectrum that pushes up to ∼35 M_⊙ for subsolar populations. Though lower metallicities lead to more massive remnants through weakened wind mass loss, the pair-instability process in the cores of massive stars is expected to cause a dearth in the remnant mass spectrum between ∼50 and 120 M_⊙ (e.g., Woosley 2017; Farmer et al. 2019; Marchant et al. 2019). The remnants observed by LIGO–Virgo have a broad distribution that contaminates the pair-instability mass gap, with several remnants consistent with being above 50M_⊙, and a couple with significant support for component masses above 80M_⊙ (LIGO Scientific Collaboration & Virgo Collaboration 2020; Abbott et al. 2021d).

Several explanations have been proposed to maintain consistency between the LIGO–Virgo mass spectrum and Kerr mergers of stellar remnants. Systematic explanations include altered mass priors (Fishbach & Holz 2020; Nitz & Capano 2021). Pre-explosion mechanisms include mass-loss suppression in low-metallicity (Z ≲ 10⁻³) progenitors (e.g., Liu & Bromm 2020; Tanikawa et al. 2021; Belczynski 2020; Spera et al. 2019; Farrell et al. 2021; Kinugawa et al. 2021; Vink et al. 2021), collisions between massive stars in dense clusters (e.g., Di Carlo et al. 2020; Kremer et al. 2020; Renzo et al. 2020), different rate parameterization for ¹²C + α → ¹⁶O + γ (i.e., the initial alpha process; e.g., Belczynski et al. 2010a, 2010b; Farmer et al. 2020; Costa et al. 2021), different overshoot parameter in Population III stars (e.g., Umeda et al. 2020; Tanikawa et al. 2021), or outer layer ejection during the common-envelope phase (e.g., Kruckow et al. 2016; Klencki et al. 2021). Post-explosion mechanisms for populating the high-mass end of the observed black hole mass spectrum include hierarchical merging in regions of high stellar density (e.g., Antonini & Rasio 2016; Fragione et al. 2020; Fragione & Silk 2020; Gayathri et al. 2020; Liu & Lai 2021; Kimball et al. 2021; Anagnostou et al. 2020; Gondán & Kocsis 2021; Tiwari & Fairhurst 2021; Rodriguez et al. 2019), accretion in dense gas clouds (e.g., Natarajan 2021; Safarzadeh & Haiman 2020; Rice & Zhang 2021), accretion and/or hierarchical merging in active galactic nucleus (AGN) disks (e.g., McKernan et al. 2012; Secunda et al. 2019; Yang et al. 2019; Tagawa et al. 2021), or sequential merging in triple systems (Vigna-Gómez et al. 2021). It has been proposed that multiple formation channels drawn from the above may be needed to explain the diversity of systems in the current observational catalog (Zevin et al. 2021), but it is not clear which combinations can reproduce the LIGO–Virgo mass spectrum at both low and high masses. There are also proposed explanations using physics beyond the standard model, including: pre-explosion mass-loss suppression due to core cooling via new light particles (Anastassopoulos et al. 2017; Di Luzio et al. 2020; Croon et al. 2020; Ng et al. 2021), an additional energy source (to fusion) throughout the progenitor (Ziegler & Freese 2021), a varying gravitational constant (Straight et al. 2020), and efficiently accreting primordial BHs (De Luca et al. 2021; Wong et al. 2021; Khalouei et al. 2021).

All of the aforementioned studies share a common thread: they have all assumed that the merging objects are Kerr BHs. This is reasonable, as the merging binary waveforms are consistent with numerical relativity simulations of two Kerr objects. Like all BH solutions within GR, however, the Kerr solution is provisional: it features horizons, singularities, and unrealistic boundary conditions (e.g., Wiltshire et al. 2009). The asymptotically flat boundary of Kerr, in particular, is incompatible with established observations. We inhabit a universe that precisely agrees with the predictions of Robertson–Walker (RW) cosmology at large distances (Aghanim et al. 2020). In other words, Kerr can only be consistently interpreted as an approximation to some astrophysically realized solution, appropriate for intervals of time short compared to the Hubble time t_H.

Attempts to construct solutions, which capture such astrophysical aspects, have made some progress. For example, gravitational collapse to the gravastar solution produces a nonsingular, horizon-free, null surface wrapping a stable bubble of dark energy (DE; Mazur & Mottola 2015; Beltracchi & Gondolo 2019). While the gravastar solution does not spin, it has been shown to provide the only known material source to Kerr exterior spacetimes to first (Posada 2017) and second (Beltracchi et al. 2021) order in slow-rotation extensions. With respect to correct boundary conditions, there are many known exact GR solutions describing various objects embedded within cosmologies. For example, the McVittie (1933) analog of a Schwarzschild BH, and the matched Nolan (1998) interior solution, exist within an arbitrary RW cosmology.

The presence of cosmological boundary conditions opens new frontiers for dynamics in time. While the Kerr approximation is excellent for timescales observable to LIGO–Virgo during any individual merger, cosmological evolution may occur over the lifetimes of compact objects. Locally, solutions with dynamical gravitating mass and horizons that comove with the cosmological expansion have been constructed (Faraoni & Jacques 2007). These solutions are significant because they are explicit counterexamples to arguments that local evolution must occur decoupled from cosmological evolution (e.g., Einstein & Straus 1945, 1946; Weinberg 2008; Peebles 2020). Recent global results in GR are consistent with these findings. It has been shown that a population of objects, over which the averaged pressure does not vanish, must couple cosmologically (Croker & Weiner 2019). For example, pure DE objects acquire a dynamical gravitating mass proportional to the RW scale factor a, cubed. Given number densities that diminish ∝ 1/a³, such a population then mimics a cosmological constant (Croker et al. 2020b). The relativistic effect is entirely analogous to the cosmological photon redshift.

In this Letter, we consider the consequences of cosmological coupling within the compact objects observable by LIGO–Virgo. We restrict our attention to objects in excess of the Tolman–Oppenheimer–Volkoff (TOV) limit, as neutron stars (NSs) are governed (in principle) by known physics. We will refer to these objects as cosmologically coupled compact objects, or C3O for short. This work extends preliminary investigations of pure DE objects (Croker et al. 2020a) and considers generic C3Os in the context of current LIGO–Virgo detection sensitivities. We use a contemporary stellar population synthesis code to produce a fiducial compact binary population from the isolated-binary evolution channel. We then track the orbital evolution of each member of the population, determining merger redshift and source-frame masses. Finally, we incorporate semianalytic estimates to the sensitive spacetime volume appropriate for the LIGO–Virgo detector network during its recent observational runs and compare the properties of the C3O population with the current observational sample of compact binary mergers.

We adopt the following single-parameter model of cosmological coupling:

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

Here m₀ is the active gravitational mass of the input stellar remnant, a is the scale factor, a_i is the scale factor at which the input stellar remnant was formed, and k is a dimensionless constant. This form is motivated by the known cosmological energy shift in photons (k = −1), and the predicted cosmological energy shift in pure DE objects (k = 3). Note that k → 0 is the decoupled limit, k < 0 gives a cosmological energy loss, and k > 0 gives a cosmological energy gain. From global analysis, physically realistic values of k are constrained to

$$\begin{eqnarray}&&-3\leqslant k\leqslant 3,\end{eqnarray} \tag{ 2 }$$

by the requirement that all causal observers perceive causal flux.

For each input binary system, we numerically integrate the orbit through the linear radiative regime. Let R denote semimajor axis, L angular momentum, e eccentricity, M primary mass, and q mass ratio. The typical linear evolution equations are derived assuming that M does not evolve in time. We extend these equations to C3O systems. Because cosmological coupling introduces no preferred spatial directions, angular momentum and eccentricity of the binary are unaffected. ⁷ This means that dL/da and de/da are unaltered from the standard expressions of Peters (1964), apart from substitution of the time-dependent mass. The remaining evolution equation for semimajor axis then follows from the chain rule,

$$\begin{eqnarray}&&\displaystyle \frac{{dR}}{{da}}=\mathop{\underbrace{{\left.\displaystyle \frac{\partial R}{\partial L}\displaystyle \frac{{dL}}{{da}}\right|}_{e,M}+{\left.\displaystyle \frac{\partial R}{\partial e}\displaystyle \frac{{de}}{{da}}\right|}_{L,M}}}\limits_{\mathrm{Radiative}\ \mathrm{decay}}+\mathop{\underbrace{{\left.\displaystyle \frac{\partial R}{\partial M}\displaystyle \frac{{dM}}{{da}}\right|}_{e,L}}}\limits_{\mathrm{Adiabatic}\ \mathrm{decay}}.\end{eqnarray} \tag{ 3 }$$

The radiative decay portion is the standard expression for dR/da, substituted with time-dependent mass. The adiabatic decay portion is computed from conservation of angular momentum in Newtonian binary evolution, combined with the derivative of Equation (1).

Note that we have omitted explicit reference to q in Equation (3). In any realistic binary system, one object will become a C3O before the other. The timescale of this delay, however, is extremely short compared to the Hubble time t_H for stellar progenitors leading to compact binaries observable by the LIGO–Virgo network. We thus neglect this difference and regard both objects as having converted at the same redshift. This implies that q is nondynamical.

The presented model requires an input stellar remnant population distributed over formation redshift. For simplicity, we assume that the full population of compact remnants originates from the canonical isolated-binary evolution channel, which includes hardening of the compact binary progenitor through either stable mass transfer or common-envelope phase. Though a large number of uncertainties exist in isolated-binary evolution that can affect the population properties and merger rates of compact remnants (e.g., winds, mass transfer stability and efficiency, common-envelope onset and evolution, natal kicks, remnant mass prescriptions, star formation history, and metallicity evolution), we simulate a single model with standard assumptions using the open-source binary population synthesis code COSMIC ⁸ (Breivik et al. 2020) and focus on the impact of cosmological coupling on this fiducial population.

Our isolated-binary population uses default parameterizations of COSMIC for the majority of its settings (see Breivik et al. 2020, and reference therein for details). We simulate 16 populations with metallicities equally log-spaced between Z = 0.0001 and Z = 0.03. As cosmological coupling can cause systems with gravitational-wave inspiral times greater than a Hubble time t_H to merge, we determine whether our population is converged using the properties of binary black hole systems that are both merged systems and those that are unmerged after t_H. Critical mass ratios for the onset of unstable mass transfer are implemented following Neijssel et al. (2019), and we assume a "pessimistic" common-envelope survival scenario (e.g., Dominik et al. 2012). We use a maximum neutron star mass (i.e., TOV limit) of 2.5 M_⊙ with a delayed remnant mass prescription (Fryer et al. 2012) as implemented in Zevin et al. (2020). We create a resampled population of 10⁶ binaries from the 16 metallicity runs and populate these systems in redshift using the star formation rate density and mean metallicity redshift dependencies in Madau & Fragos (2017), where we assume a log-normal distribution about the mean metallicity at a given redshift with a dispersion of 0.5 dex. We numerically integrate the evolution equations in the C3O scenario until merger, or stop following systems that have not merged by z = 0.

To compare the output population with the compact binaries observed by LIGO–Virgo, we estimate the relative sensitive spacetime volume of each simulated merger given its masses and merger redshift. We pre-compute detection probabilities for a grid of systems in primary mass, mass ratio, and redshift space assuming a three-detector network consisting of LIGO-Hanford, LIGO-Livingston, and Virgo with midhighlatelow power spectral densities (Abbott et al. 2020). Using this grid, we utilize a nearest neighbors algorithm to determine the detection probability of each system in our population based on its masses and redshift. The relative detection weighting of each system also accounts for the surveyed spacetime volume, such that the relative weight of each system i is

$$\begin{eqnarray}&&{w}_{i}={p}_{\det }\left({M}_{i},{q}_{i},{z}_{i}\right)\displaystyle \frac{{{dV}}_{c}}{{dz}}\displaystyle \frac{{{dt}}_{{\rm{m}}}}{{{dt}}_{0}},\end{eqnarray} \tag{ 4 }$$

where M_i is the primary mass of system i, q_i ≤ 1 is the mass ratio, dV_c /dz is the differential comoving volume at redshift z_i , and dt_m/dt₀ = (1 + z)⁻¹ is the time dilation between clocks at the merger and clocks on Earth.

To showcase the utility of the model, we present a C3O population (k = 0.5) and a decoupled (k = 0) population in Figure 1. These observable populations begin with identical input stellar remnant populations. The known tensions between the LVC population and BHs, as anticipated from the isolated-binary evolution channel, are apparent. In particular, the lack of systems with a source-frame total mass ≳90 M_⊙ is driven by the pair-instability process limiting the mass of component BHs to ≲45 M_⊙. If evolved as C3O, remnants formed through this same channel show good qualitative agreement. The broad distribution in source-frame total mass of C3O is much closer to the LVC population in slope and peak location. The C3O population also reproduces the LVC population preference for higher-redshift mergers. This is an expected correlation in mass and redshift from observational bias; massive systems are easier to detect, while more systems will be found at the detection horizon where sensitive volume is largest.

Zoom In Zoom Out Reset image size

Dependence of the observable population on the coupling strength k is shown in Figure 2. The trend is toward a smoother distribution with increased merger redshift and progressively larger masses as the coupling strength increases. The horizon redshift levels off since massive mergers that would be visible at this horizon are redshifted out of band.

Zoom In Zoom Out Reset image size

Cosmological coupling in binary systems directly alters their rate of orbital period decay. Electromagnetic observation of a putative BH–pulsar binary over many periods could allow constraint at levels comparable to the Hulse–Taylor system, PSR B1913 + 16. This system has been constrained to ±0.2% of the GR prediction after 35 yr of observation (Weisberg et al. 2010). With the caveat that the analysis of Damour & Taylor (1991, Section IV) assumes flat-space boundary conditions, one can use this constraint to estimate an upper bound on the magnitude of k for NSs:

$$\begin{eqnarray}&&\left|k\right|\lesssim 0.04\qquad (\mathrm{NS}).\end{eqnarray} \tag{ 5 }$$

The upcoming space-based gravitational-wave observatory, LISA, will be able to directly measure the rate of orbital period decay for putative BH–BH binaries in the intermediate-mass regime (Amaro-Seone et al. 2017). The sensitivity of LISA to cosmological coupling is the topic of future work.

Cosmological coupling may be visible within astrophysical systems besides isolated binaries. For example, Equation (1) predicts that a C3O born at z_i = 25 with a coupling of k = 0.5 would double in mass by the time it is observed as a quasar at z = 6. A maximal coupling of k = 3 would increase the mass by 51×. These objects presumably become supermassive BHs (SMBHs) within the cores of present-day galaxies. Because cosmological mass growth acts perpetually, if we neglect subsequent accretion, then the nonobservation of very massive objects within the local universe can be used to bound k. Following Jiang et al. (2016, Figure 10), assume that quasars with mass ∼10⁹ M_⊙ have an abundance of ∼2 Gpc⁻³ comoving at z = 6 and that uncertainties in SMBH masses are ∼0.4 dex (e.g., Dalla Bontà et al. 2020, Section 5). Then nonobservation of an SMBH with mass >10¹¹ M_⊙ within z ≲ 0.2, combined with Equation (1), conservatively gives

$$\begin{eqnarray}&&k\lesssim 3\qquad (\mathrm{BH}),\end{eqnarray} \tag{ 6 }$$

which is consistent with causal flux constraints. Comparison of the high-z SMBH mass function with the present-day SMBH mass function can be expected to refine such bounds. We also note that cosmological coupling may impact constraint of the primordial BH mass spectrum. Less massive C3O, formed deep within the radiation-dominated epoch, could grow into present-day SMBHs or their seeds. We leave further investigations in these directions to future work.

The model given in Equation (1) is appropriate for a population of identical objects, where the object-averaged energy density and principal pressures do not (independently) change much over the objects' lifetimes. Under these assumptions, k is the negative of the flat-space object-averaged principal pressures summed, divided by the flat-space object-averaged energy density. ⁹ Relativistic material within objects below the TOV limit, such as NSs, will behave as k < 0 because the dominant principal pressures are positive. Although it is suggested by known remnant models above the TOV limit, our decision to consider C3O with k ≥ 0 is pragmatic: the source-frame masses observed by LIGO–Virgo are unexpectedly large. Our exclusion of k ≥ 1 in this study follows from Figure 2. Given the input population considered, such values for k are not supported by current data.

In this work, we only consider a single astrophysical model of the remnant population in order to succinctly demonstrate the impact of cosmological coupling. Our input model may, of course, not be an accurate representation of the birth parameters and redshift evolution of binary BH systems. Contributions from other formation scenarios (e.g., Zevin et al. 2021) and uncertainties in the isolated evolution channel itself (e.g., Belczynski et al. 2021) can strongly affect this input population. One of the most qualitatively interesting features of the C3O scenario is its ability to populate the high end of the mass spectrum, where a dearth of BHs is expected in the isolated evolution channel due to the pair-instability process. Both individual events and analysis of the full binary BH population hint at structure beyond a sharp high-mass cutoff in the BH mass spectrum (Abbott et al. 2021d). Though the presence of a mass gap in the isolated evolution paradigm is a robust prediction, the exact location and width of the gap is uncertain (e.g., Farmer et al. 2019). Other channels are also predicted to be more efficient at generating merging BHs that occupy this gap, such as stellar mergers or hierarchical BH mergers in dense stellar environments. Thus, values of the coupling parameter k that match the observational population of merging BHs will be sensitive to the input birth population of binary BHs. What we have shown is that C3O solutions may act as a facade for the high-mass end of the BH birth mass spectrum. Occurrence of such solutions in nature can ease the tension between the observed BH population and astrophysical models of remnant formation.

The LIGO–Virgo population of merging compact remnants provides an unprecedented window into compact object physics. Compact objects with realistic, cosmological, boundary conditions can experience new dynamics in time. We have shown that a single-parameter model of cosmological mass growth, subsequent to compact remnant formation through the isolated-binary evolution channel, provides good qualitative agreement with the source-frame total masses and merger redshifts of the observed population. This signature of cosmological coupling can be directly measured from the rate of orbital period decay, and may be accessible to next-generation gravitational-wave observatories.

K. Croker thanks N. Warrington and the University of Washington Institute for Nuclear Theory (INT) for hospitality during preparation of this work. M. Zevin 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. Tarlé acknowledges support through DoE Award DE-SC009193.