# Constraints on cosmologically coupled black holes from gravitational wave
observations and minimal formation mass

Luca Amendola1, Davi C. Rodrigues1,2,3, Sumit Kumar4, Miguel Quartin3,5,6
1Institut für Theoretische Physik, Universität Heidelberg, Philosophenweg 16,
69120. Heidelberg, Germany  2Departamento de Física & Cosmo-Ufes, Universidade
Federal do Espírito Santo, 29075-910, Vitória - ES, Brazil.  3PPGCosmo,
Universidade Federal do Espírito Santo, 29075-910, Vitória - ES, Brazil  4Max-
Planck-Institut für Gravitationsphysik (Albert-Einstein-Institut), D-30167
Hannover, Germany.
Leibniz Universität Hannover, D-30167 Hannover, Germany  5Instituto de Física,
Universidade Federal do Rio de Janeiro, 21941-972, Rio de Janeiro, RJ, Brazil
6Observatório do Valongo, Universidade Federal do Rio de Janeiro, 20080-090,
Rio de Janeiro, RJ, Brazil

###### Abstract

We test the possibility that the black holes (BHs) detected by LIGO-Virgo-
KAGRA (LVK) may be cosmologically coupled and grow in mass proportionally to
the cosmological scale factor to some power $k$, which may also act as the
dark energy source. This approach was proposed as an extension of Kerr BHs
embedded in cosmological backgrounds and possibly without singularities. In
our analysis, we test these cosmologically coupled BHs either with or without
connection to dark energy. Assuming that the minimum mass of a BH with stellar
progenitor is $2M_{\odot}$, we estimate the probability that at least one BH
among the observed ones had an initial mass below this threshold, thereby
falsifying the hypothesis. We consider either the primary $m_{1}$ or the
secondary $m_{2}$ BHs of 72 confidently detected gravitational wave events and
adopt two different approaches. In the first one, we directly derive the
probability from the observed events, and we obtain a tension with the $k=3$
scenario at the level of 2.6$\sigma$ and 3.05$\sigma$ for the $m_{1}$ or
$m_{2}$ cases respectively. In the second approach, we assume the LVK power-
law-plus-peak (PLPP) mass distribution, which takes into account the
observational bias, and we find tensions at the level of 3.7$\sigma$ and
$4.0\sigma$ for the $m_{1}$ and $m_{2}$ masses respectively. We show that
these bounds can be alleviated by allowing lower $k$ values or faster BHs
mergers (i.e., shorter delay times $t_{\rm d}$). In particular the
$m_{1}$-based results in the following 2$\sigma$ upper bounds: $k\leq 2.4$ for
the direct approach or $k\leq 1.7$ in the PLLP approach. For $k=0.5$, a value
previously studied in the gravitational wave context, we find no relevant
constraints, even increasing the minimum BH mass to $\sim 4M_{\odot}$.
Finally, we show that future observations should quickly strengthen these
bounds.

##  I Introduction

Recently, a new intriguing hypothesis about the origin of the cosmic
acceleration has been put forward in [1, 2], with further developments in [3].
According to this scenario, black holes (BHs) grow in mass due to a form of
cosmological coupling unrelated to local accretion. If this growth is fast
enough, it could compensate the decrease in number density due to the
expansion, and generate a form of effective cosmological constant. These BHs
deviate from the standard Kerr solution since they are expected to model non-
singular solutions that are asymptotically Friedmann-Robertson-Walker, rather
than Minkowski [4, 5, 1, 2]. These non-standard BH solutions could provide an
average pressure that would constitute the entire amount of dark energy needed
to explain the cosmic acceleration if the BHs have the necessary abundance.
Farrah et al. [3] argue that this is the case.

This new BH solution is at the moment a conjecture, but it presents several
theoretical advantages: it does not need a new enigmatic ingredient in the
cosmic recipe as dark energy and it automatically solves the issue of the
“cosmic coincidence”. It would have nevertheless been just an interesting but
very speculative idea were it not for the fact that in a companion paper,
Farrah et al [6] have found strong indication in favor of just such a
cosmological growth of supermassive BHs in elliptical galaxies. This growth
seems very difficult to explain in terms of the standard local growth channels
of accretion.

More recently Lei et al. [7] used JWST data and found a conflict with Farrah
et al. [6] parametrization at redshifts $z\sim 4.5-7$. These high redshift
results are mostly independent of our analysis: the constraints we find come
from lower redshifts, as it will be shown when constraining the maximum
redshift of binary BH formation ($z_{\rm max}$). These results, if confirmed,
have an impact on the dark energy interpretation.

This paper is devoted to testing the cosmologically coupled BHs (CCBH) by
looking at the current and future datasets from LIGO-Virgo-KAGRA (LVK). The
current understanding is that the gravitational waves (GW) detected by LVK
come from the merging of BHs with stellar progenitors. If the CCBH hypothesis
is correct, they must have grown to the observed mass from an initially lower
mass. However, BHs with a stellar progenitor cannot be formed with arbitrarily
low masses. Observationally, there is evidence of a paucity of BH masses
between $2-5M_{\odot}$ [8, 9, 10], while there is no conclusive evidence for a
BH with mass below $2M_{\odot}$ (see also [11]). Neutron star (NS) stability
studies based on the Tolman-Oppenheimer-Volkoff (TOV) equation predict that
nonrotating NS could have masses at least as high as $2.2M_{\odot}$ [12, 13],
which sets a natural lower bound for BHs forming through stellar collapse.
Incidentally, rotating NS have been measured with masses as high as
$2.7M_{\odot}$ [14], and studies of binary systems containing NS find an
empirical upper limit as high as $2.6M_{\odot}$ [15]. Here we adopt a more
conservative lower bound, namely we assume that stellar BHs should not be
formed with a mass lower than $2M_{\odot}$. This is less restrictive than,
e.g., the minimum BH mass of $2.7M_{\odot}$ considered to find the proper
$\Omega_{\Lambda}$ value from CCBHs [3].

In this paper, we quantify the probability of BHs with initial mass below the
$2M_{\odot}$ threshold with two complementary approaches, explore ways to
alleviate the tensions111We use the word “tension” whenever we find a
rejection to a level higher than 2$\sigma$. we find, and briefly discuss
future prospects. We conclude that the CCBH as proposed in [3] is in strong
tension with what we know about stellar progenitor BHs, but there still is an
open parameter space where it can survive the present test. In particular, we
find no relevant tension for the CCBH case studied in [16]. The forthcoming
new GW datasets will soon shed further light on the CCBH conjecture.

The codes we used for this work are available at
https://github.com/itpamendola/CCBH-direct and https://github.com/davi-
rodrigues/CCBH-PLPP.

##  II Cosmologically coupled black holes

Farrah et al. [6] considered three samples of red-sequence elliptical galaxies
at different redshifts and found that the growth of supermassive black holes
(BHs) is significantly larger than the growth of stellar mass, being a factor
of 20 from $z\sim 2$ to $z\sim 0$. This growth is too large to be compatible
with the expected accretion rate [6]. This suggests a different growth
mechanism such that $m_{\rm BH}=m_{*}(1+z)^{-3.5\pm 1.4}$, at $90\%$
confidence level, where $m_{*}$ is the stellar mass of the galaxy and $m_{\rm
BH}$ the supermassive BH mass of the same galaxy.

A possible explanation for the above physics comes from the proposal of
cosmologically coupled BHs (CCBHs) [4, 1, 16]. In this case, BHs would grow
following the parametrization [16]

$m(a)=m(a_{i})\left(\frac{a}{a_{i}}\right)^{k}\,,$ (1)

where $k\geq 0$ is a constant, $a_{i}$ is the cosmological scale factor at the
time of the BH formation and $m(a_{i})$ is its mass at that time.

In [3] the value $k\approx 3$ was advocated, which would both explain the
supermassive BHs growth and provide a source of dark energy capable of
generating the observed $\Omega_{\Lambda}$ value. The latter requires further
assumptions, in particular a proper star formation rate, that all the remnants
with mass $>2.7$$M_{\odot}$ are BHs, and that all the BHs follow the above
mass parametrization. Criticisms of this connection with dark energy have
appeared [17, 18, 19, 20].

If all BHs are cosmologically coupled with $k=3$, Rodriguez [21] pointed out
that this would be in contradiction with globular cluster NGC 3201 data, since
it would imply that one of the BHs would have a mass below $2.2M_{\odot}$ by
the time of its formation, a less conservative minimum BH mass than our
threshold. This mass value is too low for a stellar BH, being instead
compatible with a non-rotating neutron star [13, 15]. A similar test on two
Gaia DR3 stellar-BH systems with reliable age estimation has been carried out
in Ref. [22], resulting in the 2$\sigma$ upper limit $k\leq 3.2$ assuming the
same $2.2M_{\odot}$ limit and fixing the background to $\Lambda$CDM. In Ref.
[23] it was also found that the rate of mergers and their typical masses in a
CCBH scenario would be hardly compatible with LVK observations; they also
point out that CCBHs should exhibit lower spins due to their increase in mass.

Our purpose here is to test if the BHs detected from their coalescence waves
could be cosmologically coupled. Before the results from [3], Croker et al
[16] (see also [5]) developed simulations of merging BHs and considered the
impact of the cosmological coupling on the LVK detected BHs, showing that
$k=0.5$ would be preferred over the standard $k=0$, at least for certain
isolated-binary-evolution model. We use here the most recent data from LVK,
together with more recent delay time expectations. A crucial difference
between this work and the works [16, 23] on CCBH and LVK data is that they
started from a given BBH formation mass, assumed to be realistic, and consider
if they could mimic LVK data from that initial mass distribution. Here we aim
to estimate what is the probability that at least one of the observed BHs via
LVK would be formed with an unrealistically low mass value (thus in part
similar to [21]). A key quantity for modelling the CCBH effects on LVK data is
the delay time $t_{\rm d}$ (i.e. the interval between BBH formation and
merger), which is detailed in Sec. III.

Within the general class of CCBHs, we distinguish two cases. If BHs have dark
energy implications and constitute its only source, as proposed in [2, 3],
then the constant $k$ has a direct connection with the dark energy equation of
state parameter $w$ (assumed constant), with

$k\equiv-3w\,.$ (2)

We call this scenario dark energy BHs, DEBH.

If instead BHs have no impact on dark energy, then the mass increase of BHs
can still be parametrized by a function of $z$, but the microphysics is
supposed to be independent of cosmology. This might occur if the growing BHs
contribution to the total energy density is subdominant, or if there exist
alternative non-cosmological growth channels. In this scenario, we assume no
deviation from $\Lambda$CDM cosmology. This is the case studied in, for
instance, [16]. We call this model growing BHs, GBH.

The two models, DEBH and GBH, have identical background cosmological evolution
for $k=3$ but diverge otherwise. We consider both in the following.

##  III Delay time and the mass correcting factor

The merging of binary black holes (BBH) systems detected by LVK is commonly
considered to be the end of a pair of BHs that orbited together from several
Myrs to several giga-years before the merger [24]. These BHs masses are
consistent with them being remnants of star progenitors, and this constitutes
the standard interpretation (e.g., [25, 26, 27, 28, 9]).

It is expected that most of the BBHs detected by LVK were formed from binary
stellar systems that evolved into a pair of BHs and later merged (e.g., [25,
27, 29] and references therein). The relevant time during which the CCBH
effect (1) is active extends between the BHs formation and their merger. We
call this time the BBH delay time and denote it by $t_{\rm d}$. We note that
another delay-time definition, as the time between the stellar pair formation
and the BBH merger, is also used in the literature. However, they typically
differ by a few Myr [27], hence both definitions are essentially the same.

For the GBH scenario, we consider any value of $k$ in the range $0\leq k\leq
3$, where $k=0$ corresponds to the standard case (uncoupled Kerr BHs). Apart
from the $k=3$ case, a particular value we will consider with some attention
is $k=0.5$ since it was preferred in the analysis of [16].

For the DEBH case, changing the $k$ value changes cosmology. For clarity, this
case will be parameterized with a constant $w$, instead of $k$. We mainly
consider $-0.6\leq w\leq-1$. We do not consider more negative $w$ values since
they can only strengthen our constraints; moreover [16] argues that the
maximum value of $k$ is 3. This case has no limit that leads to standard BHs.

The references [24, 9] suggested that the distribution of delay times can be
approximated by a log-uniform distribution (i.e., $p(t_{d})\propto 1/t_{d}$)
with $0.05<t_{d}(\mbox{Gyrs})<13.5$ for BBH. It is also pointed out that the
formation of the first BBHs is restricted to $z<10$. This picture is in good
agreement with simulations and observational constraints (see e.g., [25, 27,
29]). These simulations are in the context of ordinary (Kerr) BHs. CCBHs, on
the other hand, grow with time and therefore their delay times will also
change: since larger BH masses dissipate energy faster through GW emission,
the delay time of CCBHs should be smaller than for ordinary BHs for the same
initial mass and orbit [2, 3, 23]. Although for given initial conditions the
merger time will be shorter for CCBH, there are other factors that may
increase the delay time of the detected BBH mergers. In particular, as
commented in [3, 23], BBHs that would not merge before $z=0$ in the standard
picture may merge if CCBH is true. Due to such unknowns, we use a log-uniform
distribution for $t_{\rm d}$ varying three parameters that have a direct
impact on the $t_{\rm d}$ values ($t_{\rm min},\,t_{\rm max}$ and $z_{\rm
max}$, as detailed below). Moreover, in Appendix A we explore the possibility
of a steeper PDF for the $t_{\rm d}$ distribution (i.e., smaller delay times
on average), which reduces the tensions we find in the main analysis. Our
results impose constraints on the parameters that describe the delay time
distribution for the CCBH theory. We believe that the parameters we consider
are wide enough to include $t_{d}$ distributions for CCBHs.

We adopt then here the log-uniform $t_{\rm d}$ distribution,

$\log t_{\rm d}\sim U(\log t_{\rm min},\log t_{\rm x}),$ (3)

where $t_{\rm x}$ is the minimum between the maximum delay time $t_{\rm max}$
and the time difference between the merger redshift ($z_{\rm m}$) and the
maximum redshift with BBH formation ($z_{\rm max}$). As reference values, we
consider

$\displaystyle t_{\rm min}=0.05\,{\rm Gyr},$ $\displaystyle t_{\rm
max}=13.0\,{\rm Gyr},$ $\displaystyle z_{\rm max}=10,$ (4) $\displaystyle
w=-1,$ $\displaystyle k=3\,.$

Besides these reference values, we also explore other combinations. We
anticipate that the CCBH tension that we find here either increases or stays
constant if $t_{\rm max}$, $z_{\rm max}$, or $t_{\rm min}$ are increased. For
the cosmological model we assume $\Omega_{\rm m}=0.32$ and $H_{0}=70$ km s-1
Mpc-1.

Ref. [29] studied a possible correlation between $t_{\rm d}$ and $m_{1}$
considering observational data. It was found a marginal preference for smaller
masses values to have larger $t_{\rm d}$. We do not consider such a
correlation here, but if future analyses confirm this mass delay-time
correlation, it will result in stronger bounds on the CCBH model from GW data.
In any case, the true delay-time distribution and its dependence on the mass
are still uncertain (see for instance [27]).

Let us now consider a BH that is formed at $z_{i}$ and merges after a delay
time $t_{d}$ at $z_{m}$. Then

$t_{d}=\int_{z_{i}}^{z_{m}}\frac{dz}{(1+z)H(z)}\,.$ (5)

This can be inverted for a given $w$ using

$H^{2}=H_{0}^{2}[\Omega_{m}(1+z)^{3}+(1-\Omega_{m})(1+z)^{3+3w}]\,,$ (6)

leading to

$z_{i}=z_{i}(z_{m},t_{d})\,.$ (7)

Then the initial mass $m_{i}$ of BHs will be a function of
$z_{m},\,t_{d},\,k$, and proportional to the mass $m_{m}$ at merging time,

$m_{i}=m_{m}\left(\frac{1+z_{i}(z_{m},t_{d})}{1+z_{m}}\right)^{k}.$ (8)

From a given set of observed BHs masses and a $t_{\rm d}$ distribution (3), we
aim to find the probability that none of the observed BHs was formed with mass
below $2M_{\odot}$ (i.e., $M_{i}<2M_{\odot}$). If this probability is close to
1, then there is no tension between observations and the CCBH model.
Otherwise, there is tension between the observational data, the model and the
given assumptions. We will show that the latter is the case. Alternatively,
one should invoke significant changes in some of the basic assumptions, e.g. a
lower threshold for BH formation, shorter delay times, a late BH formation, or
a different growth scheme.

Table 1:  Selected confident GW events from GWTC-3 catalog [30] that are classified in [9] as BBH or NSBH systems (we require $p_{\rm astro}>0.5$ and FARmin < 1 yr-1), for a total of 72 events. The columns show the event name, the merging redshift and the primary and secondary masses. The full table is provided electronically. Event |  | $z_{m}$ |  |  $m_{1}$ (M⊙) |  |  $m_{2}$ (M⊙)
---|---|---|---|---|---|---
GW150914 |  | $0.090_{-0.030}^{+0.030}$ |  | $35.6_{-3.1}^{+5.}$ |  | $30.6_{-4.}^{+3.0}$
GW151012 |  | $0.21_{-0.09}^{+0.09}$ |  | $23._{-6.}^{+15.}$ |  | $14._{-5.}^{+4.}$
GW151226 |  | $0.09_{-0.04}^{+0.04}$ |  | $13.7_{-3.2}^{+9.}$ |  | $7.7_{-2.5}^{+2.2}$
GW170104 |  | $0.20_{-0.08}^{+0.08}$ |  | $31._{-6.}^{+7.}$ |  | $20._{-5.}^{+5.}$
GW170608 |  | $0.070_{-0.020}^{+0.020}$ |  | $11.0_{-1.7}^{+6.}$ |  | $7.6_{-2.2}^{+1.4}$
… |  | … |  | … |  | …

We estimate this probability in two complementary ways. In the first one, we
use directly the current dataset and find the joint probability that at least
one BH was born with a mass below threshold under the CCBH hypothesis. In this
way, we do not need to assume a mass distribution. Since we do not correct the
observed distribution for the selection effects, we implicitly discard low-
mass BHs from the estimation, and therefore we end up with more conservative,
but possibly more robust, estimates. We call this the direct method. In the
second method, we derive the expected initial mass distribution of BHs taking
into account the selection bias of the detectors. We use the power-law-plus-
peak (PLPP) profile [31, 9] to this end. This method leads to stronger
constraints. We denote this as the PLPP method.

The GWTC-3 data we use are shown in Table 1 and in Fig. 1. These data come
from confident BBH and NSBH events [30, 9] that satisfy $p_{\rm astro}>0.5$
and FARmin < 1 yr-1. In our analysis, we use separately either the primary
$m_{1}$ masses or the secondary $m_{2}$ ones. Therefore, each selected mass
corresponds to an independent history of a compact binary evolution and
merger. Considering all the $m_{1}$ and $m_{2}$ in a single analysis would be
incorrect since binary BHs have the same $t_{\rm d}$ and would therefore not
be independent. Although we consider here results with either primary or
secondary masses, emphasis is given on the results for $m_{1}$ masses, since
these produce more robust and more conservative constraints. For the direct
method, this choice automatically removes BHs that are outliers with
particularly low mass and have a large impact on the statistics used. For the
PLPP method, the $m_{1}$ data is more robust since its distribution depends on
one less parameter, with non-negligible uncertainty, than the $m_{2}$
distribution. In principle, one could consider $m_{2}$ masses for all the BBH
cases, and change to $m_{1}$ masses for the NSBH systems, but the impact on
the results is small since there are only one or two NSBH systems in our
selected sample.

$z_{m}$$m$ ($M_{\odot}$) Figure 1: The distribution of $m_{1}$ (disks) and
$m_{2}$ (squares) as a function of the merging redshift ($z_{m}$). Data from
Table 1.

The selected sample, Table 1, has only two systems classified as NSBH by LVK,
namely: GW200105_162426 and GW190917_114630. This classification depends on
the adopted minimum mass for BHs, and [9] considers $2.5M_{\odot}$. When
considering that the minimum mass is $2M_{\odot}$, there remains a single
NSBH, GW200115_042309. Excluding all events with secondary mass less than
$5M_{\odot}$ as potential NS or outliers, we are left with 69 $m_{2}$ data
points.

##  IV Direct constraints from the observed events

Figure 2: Direct method. Exclusion plots for the GBH and DEBH tensions. The
reference values are always in the upper right corner of each plot. The first
and second column correspond respectively to the GBH and the DEBH cases. The
two last columns show results that are common to both approaches.

Here we discuss the direct method. The formation redshift that a BH of merging
mass $m_{m}$ observed at $z_{m}$ should have to initially form with a given
threshold mass $m_{\rm th}$ is given by Eq. (8) as

$z_{\rm th}=(1+z_{m})\left(\frac{m_{\rm th}}{m_{m}}\right)^{-1/k}-1\,,$ (9)

and the corresponding delay time is

$t_{\rm th}=t_{d}(z_{m},z_{\rm th})=t_{d}(z_{m},m_{m},m_{\rm th})\,.$ (10)

If the delay time is larger than $t_{\rm th}$, the BH would have formed with a
mass below the threshold. If $t_{\rm th}$ plus the merger age $t(z_{m})$ is
larger than the cut-off $t_{\rm max}$, we take $t_{\rm th}=t_{{\rm
max}}-t(z_{m})$. Analogously, if $z_{\rm th}$ is larger than, say, $z_{\rm
max}=10$, we should cut it at $z_{\rm max}$ to prevent formation at an
unrealistically early epoch. Now, given a normalized delay-time distribution
$\Psi(t_{d})$, the probability that a BH has formation mass above $m_{\rm th}$
is

$p_{i}(z_{m},m_{m},m_{\rm th})=\int_{0}^{t_{\rm
th}}\Psi(t_{d})\textrm{d}t_{d}\,.$ (11)

The combined probability of having $N$ BHs within the acceptable formation
mass range $>m_{\rm th}$ is

$P(m_{1}>m_{\rm th})=\prod_{i}^{N}p_{i}$ (12)

and therefore the probability of at least one below-threshold BH is $1-P$. In
order to reject the CCBH hypothesis, we should find a small $P$ for the
currently observed BHs. In other words, the $p$-value for rejecting the CCBH
hypothesis is $p=P(m_{1}>m_{\rm th})$.

Since GW observations pick preferentially high-mass BHs, the constraints we
derive are on the conservative side.

We consider both the DEBH case, in which the BH growth is linked to the dark
energy so that $k=-3w$, and the alternative GBH scenario in which the BH
growth does not influence the cosmological expansion. For simplicity, in this
second case, we fix $w=-1$, i.e. the standard cosmological constant. In each
case, there is then just one BH-cosmological parameter (in addition to the
astrophysical ones, namely $t_{\rm min},\,t_{\rm max},\,z_{\rm max},\,m_{\rm
th}$): either $k$ for the GBH case, or $w$ for the DEBH case.

We illustrate in the corner plot Fig. 2 the exclusion plots for various
combinations of parameters for both scenarios, implicitly fixing all the other
parameters to the reference case. In this and in all subsequent corner plots,
the reference case (for which GHB and DEBH coincide) is always at the upper
right corner; moving beyond this point increases the rejection level of the
CCBH hypothesis.

The main result is that, for DEBH and in the reference case, the probability
of having no BHs below threshold is 0.0083, corresponding to 2.64$\sigma$.
Reducing the threshold mass to $1M_{\odot}$, the rejection level falls to
approximately 2$\sigma$. Using instead the $m_{2}$ masses, and excluding as
potential outliers (perhaps neutron stars) the two compact objects with masses
in the range $2-5M_{\odot}$, we obtain, as expected, a higher rejection level
of 3.05$\sigma$. Decreasing $w$ into the phantom regime $w<-1$ makes the
result stronger. For $w<-1.2$, using $m_{1}$ masses the rejection is at the
3$\sigma$ level (again, fixing all the other parameters to reference).

In Fig. 3 we show the distribution of probabilities for 1000 realizations of
the current data randomly chosen within the $z_{m}$ and $m_{m}$ errorbars
(assuming a bivariate normal uncorrelated distribution for each point222We
converted the 90$\%$ credibility regions of the GWTC-3 catalog into
approximated Gaussian distributions, and imposed a minimum merger mass of
$5M_{\odot}$.) for the reference case. The average and standard deviation is
$(2.655\pm 0.023)\sigma$. This narrow distribution shows that sticking with
the best fit $z_{m},\,m_{m}$ values is an acceptable approximation. We also
notice that the merger redshifts $z_{m}$ are obtained from the luminosity
distance by assuming a $\Lambda$CDM evolution. However, in the DEBH scenario,
the background is $\Lambda$CDM only for $w=-1$ so for any other value of $w$
we should derive a new set of $z_{m}$. This correction is however on average
$\Delta z=0.02$ for $w=-0.6$, and smaller for $w$ closer to $-1$. This is
negligible with respect to the current uncertainty in $z_{m}$, so we neglect
it.

Tension ($\sigma$)Counts Figure 3: Direct method. Distribution of
$P(m_{1}>m_{\rm th})$ for the reference case for 1000 random realizations of
the current GW data. The vertical black dashed line shows the median value of
this tension distribution, 2.655$\sigma$. The red dashed line shows the
tension found for the real data, 2.64$\sigma$.

We can also estimate $P(m_{1}>m_{\rm th})$ for a number of events twice as
large as the current one by simply duplicating each current event. More in
general, one trivially has $P(\alpha N)=P(N)^{\alpha}$, if $\alpha N$ is the
number of events in the forecast. With $\alpha=2$ (a number of events that can
be reached in a few months) the reference scenario would be rejected with
$p\approx 7\cdot 10^{-5}$, i.e. 4$\sigma$.

The dependence on the threshold mass is illustrated in Fig. 4. The level of
rejection of the DEBH increases quite fast with $m_{\rm th}$. For the GBH
scenario, as it can be seen from Fig. 2 and Fig. 5, values of $k<2.4$ are
acceptable at, or better than, the 2$\sigma$ level. For the other parameters,
the range for which the tension is reduced below the 2$\sigma$ level are
$t_{\rm max}<8.7$ Gyr and $z_{\rm max}<4$.

Finally, in Fig. 6 we plot the individual probabilities $p_{i}$ as a function
of the BH mass for the reference case. As expected, small BHs are more likely
to originate from below-threshold masses. However, all values of $p_{i}$ are
relatively close to unity (larger than 0.8), implying that currently observed
BH are more likely to be formed above the threshold than below it. It is the
combined probability, rather that some peculiar outlier, that leads to the
conclusion that at least one BH should have been formed below threshold.
However, one has to be cautious that when many more events will be detected,
it is likely that some outlier might bias the product of probabilities. The
probability $p(m)$ for a BH of observed mass $m$ (in solar mass units) can be
very well approximated in the reference case by the following function

$p(m)=C\log\big{(}A-Bm^{-1/2}\big{)}\,,$ (13)

with $A=26.72,B=30.88,C=0.3096$. The form of this function is suggested by the
analytical integration of Eq. (5) and Eq. (11) for a pure CDM model and a
$1/t_{d}$ distribution; the coefficients are then obtained as a best fit to
the actual $p_{i}$ values.

Minimum BH mass ($M_{\odot}$)Probability

Minimum BH mass ($M_{\odot}$)Probability

Figure 4:  Direct method. Top: plot of $P(m_{1}>m_{\rm th})$ versus minimum BH
mass in the DEBH scenario. The dotted horizontal lines mark the $\sigma$
levels. Bottom: same for the GBH scenario with $k=0.5$. The dashed curve is a
forecast for twice as many data as the current set. $k$Probability  Figure 5:
Direct method. Plot of $P(m_{1}>m_{\rm th})$ versus $k$ in the GBH scenario.
The full line is for current data, the dashed line is a forecast for twice as
many. The horizontal dotted lines mark the 1, 2, 3, 4$\sigma$ levels.

##  V Constraints using the power-law-plus-peak distribution

###  V.1 General procedures

We now move to the PLPP method. The true population of merged BBH is not well
described by the detected BBH since detection bias has an important role. In
particular, it is known that it is easier to detect massive BBH systems than
low-mass systems: many low-mass BBH mergers are expected to happen but are
undetected.

A successful profile for the primary mass ($m_{1}$) of merged BBH is the
power-law-plus-peak (PLPP) one, as proposed in [32] and analysed with current
data in [31, 9]. We consider this profile with their most probable parameter
values [31].

The PLPP is a combination of a power law, described by $\beta(m_{1})$, a
Gaussian peak given by $G(m_{1})$ and a smoothing function $S(m_{1})$ that
smooths the minimum mass probability transition. The PLPP depends on seven
parameters to describe the $m_{1}$ distribution: the power $\alpha$, the
minimum and maximum masses ($m_{\rm min},m_{\rm max}$), the Gaussian mean and
standard deviation ($\mu,\sigma$), the smoothing parameter $\delta_{m}$ and
the $\lambda$ parameter that adjusts the relative importance of the peak and
the Gaussian. The peak is interpreted as a consequence of pair-instability
supernovae [32]. The smoothing function $S$ is introduced since the most
probable $m_{1}$ values are not expected to be at the minimum $m_{\rm min}$:
expectations from X-ray binaries and simulations [32] suggest a smoother
transition. Explicitly, the PDF reads,

$\displaystyle\pi(m_{1})\propto(1-\lambda)\beta(m_{1})S(m_{1})+\lambda
G(m_{1})S(m_{1})\,,$ $\displaystyle\beta(m_{1})=\frac{\alpha-1}{m_{\rm
min}^{1-\alpha}-m_{\rm max}^{1-\alpha}}m_{1}^{-\alpha},$ $\displaystyle
G(m_{1})=\frac{1}{\sqrt{2\pi}\,\sigma}\exp\left(-\frac{(m_{1}-\mu)^{2}}{2\sigma^{2}}\right)\,,$
(14) $\displaystyle
S(m_{1})=\begin{cases}\left[1+\exp\left(\frac{\delta_{m}}{\delta
m_{1}}-\frac{\delta_{m}}{\delta
m_{1}-\delta_{m}}\right)\right]^{-1}&\\!\\!\\!\\!\\!\\!,\,\delta
m_{1}<\delta_{m}\\\\[14.22636pt] 1&\\!\\!\\!\\!\\!\\!,\,\delta
m_{1}>\delta_{m}\end{cases}\,,$

where $\delta m_{1}\equiv m_{1}-m_{\rm min}$. It is also imposed that
$\pi(m_{1})=0$ for $m_{1}<m_{\rm min}$ or $m_{1}>m_{\rm max}$.

Mass ($M_{\odot}$)Probability Figure 6:  Direct method. Individual
probabilities for each BH to be formed with a mass above threshold, as a
function of its observed mass at merging for the reference case. The cyan
curve is the best fit given by Eq. (13). The small dispersion about the curve
is due to differences in the observed redshift. The orange dots represent the
probabilities for $k=1$; they are of course much closer to unity.
ObservedInfluenced by detector bias$m_{1}$ (M⊙)PDFExpected: mergingStandard
PLPP profile$m_{1}$ (M⊙)Expected: formation ($k=3$)Modified PLPP
profile$m_{1}$ (M⊙)PDFExpected: formation ($k=0$)Standard PLPP profile$m_{1}$
(M⊙)PDF Figure 7:  PLPP method. Illustration of the relation between three
$m_{1}$ distribution contexts: the detected distribution, the expected
distribution of all the $m_{1}$ masses from BBHs that merge (considering
observational bias), and the expected $m_{1}$ distribution when it was formed.
In the standard picture ($k=0,w=-1$), the distributions for formation and
merging are the same. For CCBH, between formation and merger, BHs increase
their mass, hence the formation distribution will favour lower masses than the
standard picture.

The parameters are found from GW observational data and considering the
detector bias, through a hierarchical Bayesian approach [9]. The PLPP model
represents the source frame mass distribution which is corrected for the
selection effects. The minimum and maximum masses of the binary system are
also parameters in this model and their recovery depends on the prior
distribution used in population inference. In GWTC-3 analysis, a uniform prior
on minimum mass was used $m_{1,min}\in[2,10]M_{\odot}$. For the GWTC-3 data,
the eight parameters that describe the $m_{1}$ and $m_{2}$ distributions are
[33] (90% credible intervals):

$\displaystyle\alpha=3.40_{-0.49}^{+0.58},$
$\displaystyle\delta_{m}=4.8_{-3.2}^{+3.3}\,{\rm M}_{\odot},$ $\displaystyle
m_{\rm min}=5.08_{-1.5}^{+0.87}\,{\rm M}_{\odot},\;\;$ $\displaystyle m_{\rm
max}=86.9_{-9.4}^{+11.}\,{\rm M}_{\odot},$ (15)
$\displaystyle\mu=33.7_{-3.8}^{+2.3}\,{\rm M}_{\odot},$
$\displaystyle\sigma=3.6_{-2.1}^{+4.6}\,{\rm M}_{\odot}\,,$
$\displaystyle\lambda=0.039_{-0.026}^{+0.058}\,,$
$\displaystyle\beta_{q}=1.1_{-1.3}^{+1.8}\,.$

The $\beta_{q}$ parameter, only needed to describe the $m_{2}$ distribution,
will be discussed later on. The central values above are the medians of the
posteriors, while the uncertainties represent the $5\%$ and $90\%$ quantiles
of the posteriors. From [33] we also infer the maximum likelihood values,
which read,

$\displaystyle\alpha=3.55,\;m_{\rm min}=4.82{\rm M}_{\odot},\;m_{\rm
max}=83.14{\rm M}_{\odot}\,,$ $\displaystyle\delta_{m}=5.45{\rm
M}_{\odot},\;\mu=34.47{\rm M}_{\odot},\;\sigma=1.87{\rm M}_{\odot}\,,$ (16)
$\displaystyle\lambda=0.019,\;\beta_{q}=0.76\,.$

Since there are correlations among the parameters, the PLPP maximum likelihood
values need not to be close to the central values of each parameter. Our main
resuts use the values in Eq. (16). Variations of these values will also be
considered.

The PLPP profile for the primary mass provides the PDF of the true $m_{1}$
distribution at any merging redshift $z_{m}$ (current data do not favour a
PLPP profile with redshift dependence [9]). The detected $m_{1}$ distribution
should not match the true one since detection bias is not negligible. In
particular, low mass BHs ($\sim 5-10M_{\odot}$) are less likely to be detected
by LVK than higher masses BHs at the same redshift. The PLPP profile depends
on the proper modelling of the detector bias and on the expected BHs degrees
of freedom. Although CCBHs are not identical to standard (Kerr) black holes,
upon merging CCBHs are expected to behave as ordinary BHs [16] and hence the
PLPP should properly take into account detector bias in this scenario.

To test the CCBH hypothesis, we are interested in finding not the merged BBH
mass profile, but its expected version at the time of the BBH formation. By
finding the product distribution between the PLPP distribution with the mass
factor (8) we find the expected $m_{1}$ distribution at the formation time.
These processes are illustrated in Fig. 7.

More precisely, let $M_{1,m}$ be a random realization of the PLPP
distribution, where the index $m$ stands for merger time, and let $F_{z_{m}}$
be a random realization of the mass factor correction of eq. (8) at given
$z_{m}$. Then the $m_{1}$ distribution at BBH formation time is the
distribution of the random variate $M_{1,i}$, with

$M_{1,i}=F_{z_{m}}M_{1,m}\,.$ (17)

Our results are found using at least $10^{5}$ realizations of each random
variable.

A caveat of the above procedure is that the mass factor distribution depends
on $z_{m}$, hence the $m_{1}$ distribution by formation time is $z_{m}$
dependent. This dependence on $z_{m}$ can be either considered by using the
$z_{m}$ values from observations or by using a mean $z_{m}$ value as an
approximation. The dependence on $z_{m}$ is weak (see also Fig. 6), thus
simply using an average $z_{m}$ value is commonly sufficient (the quality of
the approximation can be directly verified by computing any probability at the
extremal $z_{m}$ values). The mass factor ($m_{i}/m_{m}$) distribution and its
dependence on $z_{m}$ are shown in Fig. 8. In particular, for a given merger
mass value $m_{m}$, the CDF coincides with the probability that the formation
mass $m_{i}$ was smaller than a given treshold. For any $z_{m}$, one sees that
the probability of $m_{i}/m_{m}<0.05$ is $\gtrsim\\!4\%$. Hence, for a given
observed BH with a merger mass of $40M_{\odot}$, the probability that the
initial mass was less than $2M_{\odot}$ is about or larger than 4% (using the
reference values).

$m_{i}/m_{m}$ PDF$m_{i}/m_{m}$ CDF Figure 8:  PDFs of the initial-to-merger
mass ratios ($m_{i}/m_{m}$), considering two different merger redshifts and
the reference values, eq. (4). Numerically derived from eq. (8) using
simulations with $10^{7}$ random realizations. Current confident BBH
observations are in the range $0.05\leq z_{m}\leq 0.92$. Inset plot: the
corresponding CDFs for small $m_{i}/m_{m}$ values.

$m_{1}$ ($M_{\odot}$)PDF

$m_{2}$ ($M_{\odot}$)

Figure 9: PLPP method. The $m_{1}$ (left plot) and $m_{2}$ (right plot)
distributions at formation for different parameter values, assuming the PLPP
distribution at merger (see Fig. 7, using PLPP parameters (16)). In blue we
show the $k=0$ case (no cosmological coupling). The green line shows the
$k=0.5$ case, within the GBH approach. The solid orange curve show the case
$k=3$ with the reference values (4); while the other three orange curves show
variations with respect to the latter values: $t_{\rm max}=5$ Gyr (dashed),
$t_{\rm min}=5$ Myr (dotted) and $z_{\rm max}=2$ (dot-dashed). The plots were
generated with $10^{7}$ simulated events.

From the $m_{1}$ mass distribution at the BBH formation, one can compute what
is the probability that one of the merged BBHs could have been formed with
$m_{1}$ larger than a given mass threshold $m_{\rm th}$, denoted by $p(m_{\rm
th},z_{m})$. For the various scenarios studied here, this probability for a
single event with a mass larger than $2M_{\odot}$ is not far from, but clearly
below, unity. Since the total number of merged BBH is larger than the total
number of confidently detected mergers, a minimum bound can be found by using
the confidently detected BHs from gravitational waves (denoted by $N$), which
we will use here. Hence, the probability of a given CCBH realization being
compatible with existing data, similarly to eq. (12), is

$P(m_{1}>m_{\rm th})=\prod_{j}^{N}p_{j}(m_{\rm th},z_{m,j})\approx
p^{N}(m_{\rm th},\overline{z_{m}})\,,$ (18)

where $\overline{z_{m}}$ is an average over all the $z_{m,j}$ values. The
above equation can be computed either using the redshift values of the
observed merged BBH or by using the last approximation, which only depends on
$\overline{z_{m}}$. The numerical differences are small, being negligible when
stating the tension with $\sigma$ units.

The above analysis applies for $m_{1}$ masses. Following [32], the $m_{2}$
distribution is given by the following conditional probability,

$\pi(m_{2}|m_{1})\propto\left(\frac{m_{2}}{m_{1}}\right)^{\beta_{q}}S(m_{2})\Theta(m1-m_{2})\,,$
(19)

where $\beta_{q}$ is the single PLPP parameter that only appears in the
$m_{2}$ distribution, $S$ is the smoothing function as defined in (14) and
$\Theta$ is the Heaviside theta function. In order to find the PDF
$\pi(m_{2})$ we marginalize over $m_{1}$,

$\pi(m_{2})=\int_{m_{\rm min}}^{m_{\rm
max}}\pi(m_{2}|m_{1})\pi(m_{1})\textrm{d}m_{1}\,,$ (20)

where $\pi(m_{1})$ is given in eq. (14) and $\pi(m_{2}|m_{1})$ needs to be
normalized for each $m_{1}$ value. With this result, one can find the initial
mass distribution and compute the probability $P(m_{1}>m_{\rm th})$, eq. (18),
for $m_{2}$ masses.

###  V.2 Results

Using the approach illustrated in Fig. 7, in Fig. 9 we show $m_{1}$ and
$m_{2}$ distributions at formation time for different parmeter values. Several
values of $k$ and of the three parameters related to the $t_{\rm d}$
distribution ($t_{\rm max}$, $t_{\rm min}$, $z_{\rm max}$) are considered. The
latter three parameters can be tuned to reduce the probability of a BH
formation with mass below $2M_{\odot}$, but, as long as this PDF is not
negligible for $m<2M_{\odot}$, the tension between observational data and the
minimum BH mass will increase with the number of detected BHs. This figure
also shows that the $k=0.5$ case, which was specially studied in [16] as a GBH
model, is the safest among the considered cases.

By applying the mass factor correction on the PLPP distribution (17) and using
eq. (18), with $N=72$ for the $m_{1}$ BHs alone, the probability that no
merged BBH was formed with mass smaller than $2M_{\odot}$ is $P\approx 2\times
10^{-4}$, thus implying a minimum tension of 3.7$\sigma$ for the reference
values. The root of this tension is the mass factor distribution, see Fig. 8.
If the number of detected events is doubled, and if the current mass
distribution is confirmed, the tension will increase to a theoretical
$5\sigma$. This number assumes that the PLPP is exactly valid.

Minimum BH mass ($M_{\odot}$)Probability

Minimum BH mass ($M_{\odot}$)Probability

Figure 10: PLPP method. Left plot: Probability that the DEBH model is in
agreement with the minimum BH mass threshold, as a function of the latter, and
for different $w$ values. This case uses $k=-3w$. Right plot: Similar to the
previous plot, but for the GBH model with $k=0.5$. Here we show also the
forecast result with twice as many events as the current data.

Instead of the $m_{1}$ values one may consider the $m_{2}$ masses, which lead
to stronger constraints but depend on an additional PLPP parameter,
$\beta_{q}$. In our events selection, Table 1, there are 69 $m_{2}$ masses
larger than $5M_{\odot}$: these are clearly BHs. Using the $m_{2}$
distribution, eq. (20), we find that the tension becomes 4.0$\sigma$.
Considering a forecast with double the data, the tension goes up to
5.8$\sigma$. Considering only data whose mass central value is larger than
$5M_{\odot}$, thereby ensuring they are BHs, the strongest constraint from the
72 events would come from considering 72 BH masses composed of 69 $m_{2}$ and
3 $m_{1}$ masses. The resulting tension slightly increases, but remains around
4.0$\sigma$.

In Fig. 10 we show how the tension increases with the increase of the BH
minimum mass. In particular, for the minimal mass of $2.7M_{\odot}$, which was
considered in [3], the tension is larger than 4$\sigma$ (using the reference
values). The left plot of this figure also considers changes in $w$, in the
context of the DEBH model. The tension decreases for larger $w$ values, which
implies lower $k$ values. It is important to point out that the case $k=0.5$
is quite safe [16]. The tension only becomes larger than $1\sigma$ if the
minimum mass is about or larger than $4.7M_{\odot}$. Also, a forecast with
double the observed events and the same PLPP distribution does not change
appreciably this picture.

To evaluate the dependence of our results on the maximum likelihood PLPP
parameters (16), in Fig. 11 we consider how the CCBH tension, considering
$m_{1}$ data alone, changes by considering $10^{3}$ samples of the PLPP
parameter distribution [33]. Although the tension does depend on the PLPP
parameter values, one sees that the tension is, as expected, always larger
than the one found with the direct approach.

In Fig. 12 we show exclusion plots for the DEBH and GBH models considering
parameter variations within the PLPP approach and with respect to the
reference values (for the 72 observed $m_{1}$ data). All the reference values
changes here considered are such that the tension decreases. This figure
should be compared with Fig. 2. Both figures show qualitatively the same
behaviour, but quantitatively, as expected, the PLPP approach is stronger. The
two most important parameters that may be changed to alleviate the tension are
$k$ and $t_{\rm max}$. For the DEBH model, the tension is reduced by
increasing $w$, which corresponds to decreasing $k$, but even the case
$w=-0.6$, which is far from the standard cosmological model, is not sufficient
to drop the tension to an acceptable level. In the following, using Fig. 12,
we highlight the necessary individual parameter ranges that reduce the
reference tension from 3.7$\sigma$ to an acceptable level, below 2.0$\sigma$.
They are: $k\leq 1.7$, $t_{\rm max}<6$ Gyr, $z_{\rm max}<2$. Reducing $t_{\rm
min}$ from 50 Myr to 5 Myr only marginally improves the picture. It is good to
recall that we have been conservative in our bounds and that further data can
quickly increase these tensions.

In Fig. 13 we show in detail the probability that there is not a single
primary BH whose initial mass was smaller than the threshold ($2M_{\odot}$),
as a function of $k$. This is shown for current data and a forecast.

TensionCounts Figure 11: PLPP method. The histogram shows the distribution of
$10^{3}$ computations of the CCBH tension, assuming the reference values (4)
and randomly selecting different PLPP parameters in accordance with the
distributions of the parameters (including correlations) (15). The vertical
black dashed line shows the median value of this tension distribution, which
is 3.69$\sigma$. The red dot-dashed line shows the tension found with the
reference parameters and the best-fit PLPP parameters (16), which is
3.68$\sigma$ at face value. The vertical limits of the light-gray background
show the 5% and 95% quantiles.  $t_{\rm max}\,{\rm(Gyr)}$$z_{\rm max}$$t_{\rm
min}\,{\rm(Myr)}$$k$$w$$t_{\rm max}\,{\rm(Gyr)}$$z_{\rm max}$ Figure 12: PLPP
method. The GBH and DEBH tensions with observational data, assuming that the
PLPP distribution models the detection bias. The reference values are always
in the upper right corner of each plot. The first and second columns
correspond respectively to the GBH and the DEBH cases. The two last columns
show results that are common to both approaches. $k$Probability  Figure 13:
PLPP method. Plot of $k$ versus $P(m_{1}>m_{\rm th})$ in the GBH scenario. The
solid line is for current data, while the dashed line is a forecast with
double the data.

##  VI Conclusions

According to a recent proposal [1, 3], BHs grow in mass due to a “cosmological
coupling” and might be responsible for the cosmic acceleration. In such a
scenario, dubbed cosmologically coupled BHs, or CCBH, BHs are not of the Kerr
type, and are supposed to match asymptotically the cosmological background.
This bold idea seems supported by recent analyses of the growth of
supermassive BHs in quiescent elliptical galaxies [6].

In this paper we tested this hypothesis by considering the binary BHs with
stellar progenitors observed with gravitational waves by the LIGO-Virgo-KAGRA
(LVK) detectors. If these BHs are cosmologically coupled, they should have
undergone a fast mass growth correlated with the cosmological scale factor
from the moment of their formation until merging, and should have formed
therefore with a mass smaller than the merger one. Since according to the
current understanding there is a mass limit below which stellar BHs cannot
form, the presently observed masses could be in conflict with this mass
threshold. By assuming a distribution of delay times, we estimated the
probability that at least one BHs observed by the GW detectors was formed
below threshold if the CCBH hypothesis is correct. We considered the low
threshold of $2M_{\odot}$.

The main result is that assuming standard astrophysical parameters for the
delay time distribution and the mass threshold, we find a tension at a
confidence level of 2.6$\sigma$ or more, using only the primary masses
($m_{1}$) data. More specifically, we obtain 2.6$\sigma$ for the direct method
of Sec. IV, and 3.7$\sigma$ for the method based on the power-law-plus-peak
distribution of Sec. V. Stronger bounds are obtained by employing the
secondary masses, namely approximately 3.0 and 4.0$\sigma$, respectively. The
result depends on a number of astrophysical parameters, so we explore various
possibilities to ease the tension. Of particular relevance is to notice that,
if the maximum time between the binary BH formation and their merger ($t_{\rm
d}$) is not larger than 8.7 Gyrs (in the direct method) or 6 Gyrs (in the PLPP
method), the tensions essentially disappear. Finally, for the CCBH variation
studied in [16], whose BH masses increase more slowly ($k=0.5$), we found no
relevant tension with minimum BH masses. This conclusion is also valid for the
case $k=1$, which was recently studied in [34].

It is clear that the physics of CCBHs is still to be understood in any detail,
and therefore the assumptions about minimal mass and evolution might need to
be reconsidered in the future. Moreover, if the LVK BHs are of primordial
origin, then a completely different analysis would be needed (see e.g. [23]).

The current O4 run of the LKV detectors is expected to more than double the
number of BH mergers by the end of 2024. With this extended dataset, it will
be possible to rule out a much larger fraction of parameter space.

###### Acknowledgements.

We thank Kevin Croker and Valerio Faraoni for very useful comments and
feedback. We also thank Riccardo Sturani for several discussions about this
project. LA acknowledges support from DFG project 456622116. DCR acknowledges
partial support from CNPq (Brazil) and FAPES (Brazil) (TO 1020/2022, 976/2022,
1081/2022). MQ is supported by the Brazilian research agencies FAPERJ, CNPq
(Conselho Nacional de Desenvolvimento Científico e Tecnológico) and CAPES.
This study was financed in part by the Coordenação de Aperfeiçoamento de
Pessoal de Nível Superior - Brasil (CAPES) - Finance Code 001. We acknowledge
support from the CAPES-DAAD bilateral project “Data Analysis and Model Testing
in the Era of Precision Cosmology”.

##  Appendix A Changing the delay time distribution

The delay time distribution is clearly a crucial assumption, so we discuss
here briefly the impact of changing the $1/t_{d}$ slope. For this appendix we
adopt the more conservative direct method of Sec. IV. Ref. [27] shows a
tendency for low-mass BHs to be formed via the common envelope channel, and
this channel would have a delay time distribution steeper than $1/t_{d}$
(i.e., shorter delay times on average). The full functional form of the
relation delay times versus mass is unknown and it appears difficult to model
exactly. Here we limit ourselves to a preliminary investigation. We modify the
$t_{d}$ distribution by assuming a $t_{d}^{-\beta}$ with $\beta$ slightly
larger than unity (e.g. $\beta=1.1-1.3$), i.e., shorter delay times on
average. Looking at the results of App. B of [27], we see that such a steeper
power-law distribution approximates the predicted behavior for masses below
$30M_{\odot}$. For simplicity, we assume that this power-law extends to all
masses: this has anyway very little impact since large masses are not the main
drivers of our statistics. In Fig. 14 we show the probability contour plot for
$w,\beta$. For $\beta\geq 1.2$, the probability for $w=-1$ decreases below
2$\sigma$, bringing the DEBH model into the non-rejection region. Therefore,
as far as current data are concerned, such steeper power-laws might alleviate
or solve the tension.

Figure 14: Direct method. Contour plot of $P(m_{1}>m_{\rm th})$ as a function
of $w$ and the slope $\beta$ in the DEBH model.

## References

  * Croker and Weiner [2019] K. S. Croker and J. L. Weiner, Astrophys. J. 882, 19 (2019), arXiv:2107.06643 [gr-qc] .
  * Croker  _et al._ [2020a] K. S. Croker, J. Runburg, and D. Farrah, Astrophys. J. 900, 57 (2020a).
  * Farrah  _et al._ [2023] D. Farrah _et al._ , Astrophys. J. Lett. 944, L31 (2023), arXiv:2302.07878 [astro-ph.CO] .
  * Faraoni and Jacques [2007] V. Faraoni and A. Jacques, Phys. Rev. D 76, 063510 (2007), arXiv:0707.1350 [gr-qc] .
  * Croker  _et al._ [2020b] K. Croker, K. Nishimura, and D. Farrah, Astrophys. J. 889, 115 (2020b), arXiv:1904.03781 [astro-ph.CO] .
  * Farrah  _et al._ [2023] D. Farrah, S. Petty, K. S. Croker, G. Tarlé, M. Zevin, E. Hatziminaoglou, F. Shankar, L. Wang, D. L. Clements, A. Efstathiou, M. Lacy, K. A. Nishimura, J. Afonso, C. Pearson, and L. K. Pitchford, Astrophys. J. 943, 133 (2023), arXiv:2212.06854 [astro-ph.GA] .
  * Lei  _et al._ [2023] L. Lei _et al._ , (2023), arXiv:2305.03408 [astro-ph.CO] .
  * Özel  _et al._ [2010] F. Özel, D. Psaltis, R. Narayan, and J. E. McClintock, Astrophys. J. 725, 1918 (2010), arXiv:1006.2834 [astro-ph.GA] .
  * Abbott  _et al._ [2023] R. Abbott _et al._ (KAGRA, VIRGO, LIGO Scientific), Phys. Rev. X 13, 011048 (2023), arXiv:2111.03634 [astro-ph.HE] .
  * de Sá  _et al._ [2022] L. M. de Sá, A. Bernardo, R. R. A. Bachega, J. E. Horvath, L. S. Rocha, and P. H. R. S. Moraes, Astrophys. J. 941, 130 (2022), arXiv:2211.01447 [astro-ph.HE] .
  * Abbott  _et al._ [2022] R. Abbott _et al._ (LIGO Scientific, VIRGO, KAGRA), (2022), arXiv:2212.01477 [astro-ph.HE] .
  * Ye and Fishbach [2022] C. Ye and M. Fishbach, Astrophys. J. 937, 73 (2022), arXiv:2202.05164 [astro-ph.HE] .
  * Legred  _et al._ [2021] I. Legred, K. Chatziioannou, R. Essick, S. Han, and P. Landry, Phys. Rev. D 104, 063003 (2021), arXiv:2106.05313 [astro-ph.HE] .
  * Romani  _et al._ [2012] R. W. Romani, A. V. Filippenko, J. M. Silverman, S. B. Cenko, J. Greiner, A. Rau, J. Elliott, and H. J. Pletsch, Astrophys. J. Lett. 760, L36 (2012), arXiv:1210.6884 [astro-ph.HE] .
  * Rocha  _et al._ [2021] L. S. Rocha, R. R. A. Bachega, J. E. Horvath, and P. H. R. S. Moraes, (2021), arXiv:2107.08822 [astro-ph.HE] .
  * Croker  _et al._ [2021] K. S. Croker, M. J. Zevin, D. Farrah, K. A. Nishimura, and G. Tarle, Astrophys. J. Lett. 921, L22 (2021), arXiv:2109.08146 [gr-qc] .
  * Sadeghi  _et al._ [2023] J. Sadeghi, S. Noori Gashti, M. R. Alipour, and M. A. S. Afshar, (2023), arXiv:2305.12545 [gr-qc] .
  * Avelino [2023] P. P. Avelino, (2023), arXiv:2303.06630 [gr-qc] .
  * Maes [2023] S. H. Maes,  (2023), 10.31219/osf.io/369pd.
  * Mistele [2023] T. Mistele, Res. Notes AAS 7, 101 (2023), arXiv:2304.09817 [gr-qc] .
  * Rodriguez [2023] C. L. Rodriguez, Astrophys. J. Lett. 947, L12 (2023), arXiv:2302.12386 [astro-ph.CO] .
  * Andrae and El-Badry [2023] R. Andrae and K. El-Badry, Astron. Astrophys. 673, L10 (2023), arXiv:2305.01307 [astro-ph.CO] .
  * Ghodla  _et al._ [2023] S. Ghodla, R. Easther, M. M. Briel, and J. J. Eldridge, (2023), arXiv:2306.08199 [astro-ph.CO] .
  * Abbott  _et al._ [2021a] R. Abbott _et al._ (KAGRA, Virgo, LIGO Scientific), Phys. Rev. D 104, 022004 (2021a), arXiv:2101.12130 [gr-qc] .
  * Belczynski  _et al._ [2016] K. Belczynski, D. E. Holz, T. Bulik, and R. O’Shaughnessy, Nature 534, 512 (2016), arXiv:1602.04531 [astro-ph.HE] .
  * Mapelli [2020] M. Mapelli, Front. Astron. Space Sci. 7, 38 (2020), arXiv:2105.12455 [astro-ph.HE] .
  * van Son  _et al._ [2022] L. A. C. van Son, S. E. de Mink, T. Callister, S. Justham, M. Renzo, T. Wagg, F. S. Broekgaarden, F. Kummer, R. Pakmor, and I. Mandel, Astrophys. J. 931, 17 (2022), arXiv:2110.01634 [astro-ph.HE] .
  * Chen  _et al._ [2022] Z. Chen, Y. Lu, and Y. Zhao, Astrophys. J. 940, 17 (2022), arXiv:2210.09892 [astro-ph.HE] .
  * Fishbach and Kalogera [2021] M. Fishbach and V. Kalogera, Astrophys. J. Lett. 914, L30 (2021), arXiv:2105.06491 [astro-ph.HE] .
  * Abbott  _et al._ [2021b] R. Abbott _et al._ (LIGO Scientific, VIRGO, KAGRA), (2021b), arXiv:2111.03606 [gr-qc] .
  * Abbott  _et al._ [2021c] R. Abbott _et al._ (LIGO Scientific, Virgo), Astrophys. J. Lett. 913, L7 (2021c), arXiv:2010.14533 [astro-ph.HE] .
  * Talbot and Thrane [2018] C. Talbot and E. Thrane, Astrophys. J. 856, 173 (2018), arXiv:1801.02699 [astro-ph.HE] .
  * LIGO Scientific Collaboration  _et al._ [2023] LIGO Scientific Collaboration, Virgo Collaboration, and KAGRA Collaboration, Zenodo v2 (2023), 10.5281/zenodo.7843926.
  * Cadoni  _et al._ [2023] M. Cadoni, A. P. Sanna, M. Pitzalis, B. Banerjee, R. Murgia, N. Hazra, and M. Branchesi, (2023), arXiv:2306.11588 [gr-qc] .