Cosmologically coupled compact objects: a single parameter model for LIGO--Virgo mass and redshift distributions

We demonstrate a single-parameter route for reproducing higher mass objects as observed in the LIGO--Virgo mass distribution, using only the isolated binary stellar evolution channel. This single parameter encodes the cosmological mass growth of compact stellar remnants that exceed the Tolman-Oppenheimer-Volkoff limit. Cosmological mass growth appears in known solutions to General Relativity with cosmological boundary conditions. We consider the possibility of solutions with cosmological boundary conditions, which reduce to Kerr on timescales short compared to the Hubble time. We discuss complementary observational signatures of these solutions that can confirm or invalidate their astrophysical relevance.


INTRODUCTION
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) has currently published over 50 confident detections of compact binary mergers, the majority being consistent with Kerr black holes (BHs) (Abbott et al. 2019(Abbott et al. , 2021a. 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 * NASA Hubble Fellow mergers visible by LIGO-Virgo. In this scenario, massive progenitor stars in a binary system proceed through one of 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 ∼ 8M at solar metallicity and a broader spectrum that pushes up to ∼ 35M 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−120M (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 Observable merger redshift and source frame total mass joint distribution of a decoupled stellar remnant population (k = 0, orange) and the same stellar remnant population evolved as cosmologically coupled objects (k = 0.5, blue). Contours display 1σ (dark) and 2σ (light) confidence. Predicted detectable populations are computed from the same 10 6 input stellar remnant binaries. The decoupled population produced 8417 observable mergers, while the k = 0.5 population produced 167867, a 20× increase in merger rate. The most recent LVC sample of binary mergers (black) is overplotted, with 90% credible intervals in total mass and redshift shown (Abbott et al. 2019(Abbott et al. , 2021a. Marginalized distributions in merger redshift and source frame total mass are displayed on the right and top axes, respectively, with the LVC population given for reference (dashed, black). LVC binary systems with primary or secondary masses below 2.5M have been excluded. Credible intervals at 90% have been truncated for the two extreme LVC systems for clarity of visualization. explanations include altered mass priors (Fishbach & Holz 2020;Nitz & Capano 2021). Pre-explosion mechanisms include mass loss suppression in low-metallicity (Z 10 −3 ) progenitors (e.g. Liu & Bromm 2020;Tanikawa et al. 2020a;Belczynski 2020;Spera et al. 2019;Farrell et al. 2020;Kinugawa et al. 2020;Vink et al. 2020); 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 12 C + α → 16 O + γ (i.e. the initial alpha process) (e.g. Belczynski et al. 2010a,b;Farmer et al. 2020;Costa et al. 2020); different overshoot parameter in Population III stars (e.g. Umeda et al. 2020;Tanikawa et al. 2020b); or outer layer ejection during the common-envelope phase (e.g. Kruckow et al. 2016;Klencki et al. 2020). Postexplosion 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 & Silk 2020;Gayathri et al. 2020;Liu & Lai 2020;Kimball et al. 2020;Anagnostou et al. 2020;Gondán & Kocsis 2020;Tiwari & Fairhurst 2020;Rodriguez et al. 2019); accretion in dense gas clouds (e.g. Natarajan 2020; Safarzadeh & Haiman 2020; Rice & Zhang 2020); accretion and/or hierarchical merging in AGN disks (e.g. McKernan et al. 2012;Secunda et al. 2019;Yang et al. 2019;Tagawa et al. 2020); or sequential merging in triple systems (Vigna-Gómez et al. 2020). 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. 2020); an additional energy source (to fusion) throughout the progenitor (Ziegler & Freese 2020); a varying gravitational constant (Straight et al. 2020); and efficiently accreting primordial BHs Wong et al. 2020;Khalouei et al. 2020).
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 non-singular, 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) analogue 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. 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 counter-examples to arguments that local evolution must occur decoupled from cosmological evolution (e.g. Einstein & Straus 1945, 1946Weinberg 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 3 , 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 paper, 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 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 semi-analytic 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.
Here m 0 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 Dark Energy objects (k = 3). Note that k → 0 is the decoupled limit, k < 0 gives a cosomological energy loss, and k > 0 gives a cosmological energy gain. From global analysis, physically realistic values of k are constrained to 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 semi-major 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. 1 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 semi-major axis then follows from the chain rule, . (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 non-dynamical.

METHODS
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 2 (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 set-1 This can also be shown directly with the theory of adiabatic invariants, see Croker et al. (2020a). tings (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.5M with a delayed remnant mass prescription (Fryer et al. 2012) as implemented in Zevin et al. (2020). We create a resampled population of 10 6 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 rela-tive 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 3detector 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 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 0 = (1 + z) −1 is the time dilation between clocks at the merger and clocks on Earth.

RESULTS
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 90M is driven by the pair instability process limiting the mass of component BHs to 45M . 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.
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.

DISCUSSION
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 years of observation (Weisberg et al. 2010). For comparison, a binary C3O with −0.3 k −0.15 can lead to ∼ 1% alterations in orbital period decay (Croker & Weiner 2019). 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 (Danzmann et al. 2017). The sensitivity of LISA to cosmological coupling is the topic of future work.
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.