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

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 if $k\approx 3$. This approach was proposed as an extension of Kerr BHs embedded in cosmological backgrounds and possibly without singularities or horizons. In our analysis, we develop and apply two methods to test these cosmologically coupled BHs (CCBHs) either with or without connection to dark energy. We consider different scenarios for the time between the binary BH formation and its merger, and we find that the standard log-uniform distribution yields weaker constraints than the CCBH-corrected case. 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. We obtain these probabilities either directly from the observed data or by assuming the LVK power-law-plus-peak mass distribution. In the latter case we find, at $2\sigma$ level, that $k<2.1$ for the standard log-uniform distribution, or $k<1.1$ for the CCBH-corrected distribution. Slightly weaker bounds are obtained in the direct method. Considering the uncertainties on the nature of CCBHs, we also find that the required minimum CCBH mass value to eliminate the tensions for $k=3$ should be lower than 0.5 $M_\odot$ (again at 2$\sigma$). Finally, we show that future observations have the potential to decisively confirm 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 [1,2,4,5].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 ∼ 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 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 [8][9][10], while there is no conclusive evidence for a BH with mass below 2M (see also [11]).Neutron star (NS) stability studies based on the Tolman-Oppenheimer-Volkoff (TOV) equation pre-dict that nonrotating NS could have masses at least as high as 2.2M [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 [14], and studies of binary systems containing NS find an empirical upper limit as high as 2.6M [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 .This is less restrictive than, e.g., the minimum BH mass of 2.7M considered to find the proper Ω Λ value from CCBHs [3].
In this paper, we quantify the probability of BHs with initial mass below the 2M threshold with two complementary approaches, explore ways to alleviate the tensions 1 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.

II. COSMOLOGICALLY COUPLED BLACK HOLES
Farrah et al. [6] considered three samples of redsequence 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 ∼ 2 to z ∼ 0. This growth is too large to be compatible with the expected accretion rate [6].This suggests a different growth mechanism such that m BH = m * (1 + z) −3.5±1. 4, at 90% confidence level, where m * is the stellar mass of the galaxy and m 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) [1,4,16].In this case, BHs would grow following the parametrization [16] where k ≥ 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 ≈ 3 was advocated, which would both explain the supermassive BHs growth and provide a source of dark energy capable of generating the observed Ω Λ value.The latter requires further assumptions, in particular a proper star formation rate, that all the remnants with mass > 2.7M 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 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σ upper limit k ≤ 3.2 assuming the same 2.2M limit and fixing the background to Λ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 isolatedbinary-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 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 ≡ −3w . ( 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 Λ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., [9,[25][26][27][28]).
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 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 ≤ k ≤ 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 ≤ w ≤ −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 [9,24] suggested that the distribution of delay times can be approximated by a log-uniform distribution (i.e., p(t d ) ∝ 1/t d ) with 0.05 < t d (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 d varying three parameters that have a direct impact on the t d values (t min , t max and z max , as detailed below).Moreover, in Appendix A we explore the possibility of a steeper PDF for the t 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 d distribution, where t x is the minimum between the maximum delay time t max and the time difference between the merger redshift (z m ) and the maximum redshift with BBH formation (z max ).As reference values, we consider t min = 0.05 Gyr, t max = 13.0Gyr, 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 max , z max , or t min are increased.For the cosmological model we assume Ω m = 0.32 and H 0 = 70 km s −1 Mpc −1 .
Ref. [29] studied a possible correlation between t d and m 1 considering observational data.It was found a marginal preference for smaller masses values to have larger t 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 This can be inverted for a given w using leading to 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, TABLE I. Selected confident GW events from GWTC-3 catalog [30] that are classified in [9] as BBH or NSBH systems (we require pastro > 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 From a given set of observed BHs masses and a t d distribution (3), we aim to find the probability that none of the observed BHs was formed with mass below 2M (i.e., M i < 2M ).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.
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 [9,31] 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 I and in Fig. 1.These data come from confident BBH and NSBH events [9,30] that satisfy p astro > 0.5 and FAR min < 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 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  I.
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.
The selected sample, Table I, 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 .When considering that the minimum mass is 2M , there remains a single NSBH, GW200115_042309.Excluding all events with secondary mass less than 5M as potential NS or outliers, we are left with 69 m 2 data points.

IV. DIRECT CONSTRAINTS FROM THE OBSERVED EVENTS
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 th is given by Eq. ( 8) as and the corresponding delay time is  BH has formation mass above m th is The combined probability of having N BHs within the acceptable formation mass range > m th is and therefore the probability of at least one belowthreshold 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 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 min , t max , z max , m 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σ.Reducing the threshold mass to 1M , the rejection level falls to approximately 2σ.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 , we obtain, as expected, a higher rejection level of 3.05σ.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σ 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 point2 ) for the reference case.The average and standard deviation is (2.655 ± 0.023)σ.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 ΛCDM evolution.However, in the DEBH scenario, the background is Λ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 ∆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.
We can also estimate P (m 1 > m 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 (αN ) = P (N ) α , if αN is the number of events in the forecast.With α = 2 (a number of events that can be reached in a few months) the reference scenario would be rejected with p ≈ 7 • 10 −5 , i.e. 4σ.
The dependence on the threshold mass is illustrated in Fig. 4. The level of rejection of the DEBH increases quite fast with m 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σ level.For the other parameters, the range for which the tension is reduced below the 2σ level are t max < 8.7 Gyr and z 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 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.

A. 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 [9,31].We consider this profile with their most probable parameter values [31].
The PLPP is a combination of a power law, described by β(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 α, the minimum and maximum masses (m min , m max ), the Gaussian mean and standard deviation (µ, σ), the smoothing parameter δ m and the λ parameter that adjusts the relative importance of the peak and the Gaussian.The peak is interpreted as a consequence of pairinstability supernovae [32].The smoothing function S is introduced since the most probable m 1 values are not expected to be at the minimum m min : expectations from X-ray binaries and simulations [32] suggest a smoother transition.Explicitly, the PDF reads, where δm 1 ≡ m 1 −m min .It is also imposed that π(m 1 ) = 0 for m 1 < m min or m 1 > m max .

Mass (M )
Probability FIG. 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.
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 ∈ [2, 10]M .For the GWTC-3 data, the eight parameters that describe the m 1 and m 2 distributions are [33] (90% credible intervals): The β 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, λ = 0.019, β 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.Illustration of the relation between three m1 distribution contexts: the detected distribution, the expected distribution of all the m1 masses from BBHs that merge (considering observational bias), and the expected m1 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 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 (∼ 5 − 10M ) 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 zm 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 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 4%.Hence, for a given observed BH with a merger mass of 40M , the probability that the initial mass was less than 2M is about or larger than 4% (using the reference values).
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 th , denoted by p(m th , z m ).4).Numerically derived from eq. ( 8) using simulations with 10 7 random realizations.Current confident BBH observations are in the range 0.05 ≤ zm ≤ 0.92.Inset plot: the corresponding CDFs for small mi/mm values.
For the various scenarios studied here, this probability for a single event with a mass larger than 2M 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 where 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 z m .The numerical differences are small, being negligible when stating the tension with σ units.
The above analysis applies for m 1 masses.Following [32], the m 2 distribution is given by the following conditional probability, where β q is the single PLPP parameter that only appears in the m 2 distribution, S is the smoothing function as defined in (14) and Θ is the Heaviside theta function.In order to find the PDF π(m 2 ) we marginalize over m 1 , where π(m 1 ) is given in eq. ( 14) and π(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 th ), eq. ( 18), for m 2 masses.

B. 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 d distribution (t max , t min , z max ) are considered.The latter three parameters can be tuned to reduce the probability of a BH formation with mass below 2M , but, as long as this PDF is not negligible for m < 2M , 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 is P ≈ 2 × 10 −4 , thus implying a minimum tension of 3.7σ 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σ.This number assumes that the PLPP is exactly valid.
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, β q .In our events selection, Table I, there are 69 m 2 masses larger than 5M : these are clearly BHs.Using the m 2 distribution, eq. ( 20), we find that the tension becomes 4.0σ.Considering a forecast with double the data, the tension goes up to 5.8σ.Considering only data whose mass central value is larger than 5M , 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σ.
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 , which was considered in [3], the tension is larger than 4σ (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σ if the minimum mass is about or larger than 4.7M .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  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: tmax = 5 Gyr (dashed), tmin = 5 Myr (dotted) and zmax = 2 (dot-dashed).The plots were generated with 10 7 simulated events.Here we show also the forecast result with twice as many events as the current data.
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 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σ to an acceptable level, below 2.0σ.They are: k ≤ 1.7, t max < 6 Gyr, z max < 2. Reducing t 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 ), as a function of k.This is shown for current data and a forecast.

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  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σ.The red dot-dashed line shows the tension found with the reference parameters and the best-fit PLPP parameters (16), which is 3.68σ at face value.The vertical limits of the light-gray background show the 5% and 95% quantiles.
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 .
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σ or more, using only the primary masses (m 1 ) data.More specifically, we obtain 2.6σ for the direct method of Sec.IV, and 3.7σ 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σ, 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 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.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, β.For β ≥ 1.2, the probability for w = −1 decreases below 2σ, bringing the DEBH model into the non-rejection region.Therefore, as far as current data are concerned, such steeper powerlaws might alleviate or solve the tension.

FIG. 2 .
FIG. 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.

FIG. 3 .
FIG.3.Direct method.Distribution of P (m1 > m 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σ.The red dashed line shows the tension found for the real data, 2.64σ.

FIG. 4 .
FIG.4.Direct method.Top: plot of P (m1 > m th ) versus minimum BH mass in the DEBH scenario.The dotted horizontal lines mark the σ 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.

FIG. 5 .
FIG.5.Direct method.Plot of P (m1 > m 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σ levels.

FIG. 7 .
FIG.7.PLPP method.Illustration of the relation between three m1 distribution contexts: the detected distribution, the expected distribution of all the m1 masses from BBHs that merge (considering observational bias), and the expected m1 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.

FIG. 9 .
FIG.9.PLPP method.The m1 (left plot) and m2 (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: tmax = 5 Gyr (dashed), tmin = 5 Myr (dotted) and zmax = 2 (dot-dashed).The plots were generated with 10 7 simulated events.

FIG. 10 .
FIG.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.

FIG. 11 .
FIG.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σ.The red dot-dashed line shows the tension found with the reference parameters and the best-fit PLPP parameters(16), which is 3.68σ at face value.The vertical limits of the light-gray background show the 5% and 95% quantiles.

FIG. 12 .
FIG.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.
table is provided electronically.