Cosmological coupling of nonsingular black holes

We show that — in the framework of general relativity (GR) — if black holes (BHs) are singularity-free objects, they couple to the large-scale cosmological dynamics. We find that the leading contribution to the resulting growth of the BH mass (M BH) as a function of the scale factor a stems from the curvature term, yielding M BH ∝ ak , with k = 1. We demonstrate that such a linear scaling is universal for spherically-symmetric objects, and it is the only contribution in the case of regular BHs. For nonsingular horizonless compact objects we instead obtain an additional subleading model-dependent term. We conclude that GR nonsingular BHs/horizonless compact objects, although cosmologically coupled, are unlikely to be the source of dark energy. We test our prediction with astrophysical data by analysing the redshift dependence of the mass growth of supermassive BHs in a sample of elliptical galaxies at redshift z = 0.8–0.9. We also compare our theoretical prediction with higher redshift BH mass measurements obtained with the James Webb Space Telescope (JWST). We find that, while k = 1 is compatible within 1σ with JWST results, the data from elliptical galaxies at z = 0.8–0.9 favour values of k > 1. New samples of BHs covering larger mass and redshift ranges and more precise BH mass measurements are required to settle the issue.


Introduction
Recently, Farrah et al. [1], following previous proposals [2][3][4][5] and building on recent experimental results [6], found observational evidence for cosmologically-coupled mass growth in supermassive black holes (SMBHs) at the center of elliptical galaxies at different redshift.To do so, they employed the mass formula (firstly proposed in Ref. [7]) where a is the scale factor, and a i is the value of the latter evaluated at a particular reference time.Physically, this is identified as the initial time at which the compact object formed and became coupled to the cosmological background.The slope parameter k is a constant, quantifying the strength of the coupling between the small-scale matter inhomogeneities and the large-scale cosmological dynamics.The physically relevant values of k are limited by causality to the interval −3 ≤ k ≤ 3 [5], with k = 0 corresponding to zero coupling.The authors of Ref. [1] found a preference consistent with k ∼ 3 at 90 % confidence level (C.L.), pointing therefore to nonzero coupling.This would also imply that the density of these objects in a cosmological volume a 3 is constant, thus leading to the conclusion of Ref. [1] that BHs could be responsible for the observed accelerated expansion of the universe (see, however, Refs.[8][9][10][11][12][13][14][15][16][17]).
More recently, tests of cosmological BH-mass coupling have been performed using higher redshift Active Galactic Nuclei (AGN) data coming from the JWST survey [14].Strikingly, these tests favour much smaller values of k, with a best-fit value of k ≃ −1.Moreover, a k = 3 cosmological coupling seems to be in tension with the predictions of the LIGO-Virgo-KAGRA network on the rate and typical masses of black hole mergers [18].In view of this intricate situation, an exhaustive theoretical understanding of the cosmological coupling of BHs and other compact astrophysical objects is clearly crucial.
The main purpose of this letter is to investigate the cosmological coupling of nonsingular compact objects, such as regular BHs and horizonless configurations, in an exact general relativity (GR) framework.This is realized by sourcing Einstein's field equations with an anisotropic fluid (see, e.g., Ref. [19], and references therein), which allows to circumvent Penrose's singularity theorem [20].
We start from a parametrization of the spacetime metric that naturally couples the cosmological dynamics with local, spherically-symmetric objects, since BH angular momentum effects are expected to be negligible on cosmological scales.We then show that Einstein's equations allow for solutions describing nonsingular BHs, whose Misner-Sharp (MS) mass is given by the universal form (1.1), but with k = 1.In the case of horizonless objects, instead, we find subleading model-dependent corrections to this formula (not considered in Ref. [1]).We apply our cosmologically-coupled mass expression to some explicit models of regular objects, representing BH mimickers.
In order to test our theoretical prediction with astrophysical data, we perform an analysis of the redshift dependence of the mass growth of SMBHs in elliptical galaxies, along the lines of the study performed in Ref. [6].Our sample consists of 108 SMBHs, at redshift z = 0.8−0.9.We also compare our theoretical result with the outcome of the tests performed using higher redshift AGN data from JWST [14].We find that, while k = 1 is compatible with the latter within a 1σ confidence interval, it is quite in tension with the data from elliptical galaxies, which favour higher values of k instead.

Cosmological coupling of astrophysical compact objects
The anisotropic fluid we employ to source the gravitational field is described by the stressenergy tensor 1  T where u µ and w µ satisfy g µν u µ u ν = −1, g µν w µ w ν = 1, u µ w µ = 0. ρ, p , p ⊥ are, respectively, the fluid density and the pressure components, parallel and perpendicular to the spacelike vector w µ .We consider a general parametrization of the spacetime metric, which is suitable to describe the cosmological coupling of local, spherically-symmetric objects: ds 2 = a 2 (η) −e α(η,r) dt 2 + e β(η,r) dr 2 + r 2 dΩ 2 .
(2.2) α(η, r) and β(η, r) are metric functions depending on both η (the conformal time) and r, while a(η) is the scale factor.It was shown in Ref. [21] that this parametrization of the metric, together with eq.(2.1), can be used to describe different regimes: the dynamics of cosmological inhomogeneities on small scales, the dynamics of large cosmological scales, but also their coupling at intermediate scales.
Appropriate combinations of Einstein's equations give (see Ref. [21] for further details) α β = 0, where the dot means derivation with respect to (w.r.t.) the conformal time.The field equations allow, therefore, for two classes of solutions, characterized either by β = 0 or by α = 0.The first case corresponds to a complete decoupling between the cosmological, timedependent dynamics (that of the scale factor) from the small, r-dependent one, pertaining instead to the local inhomogeneity.This model allows for both standard Friedmann-Lemaître-Robertson-Walker (FLRW) and other solutions, but not for black-hole ones in a cosmological background (for further details, see Ref. [22]).
As we are here interested in a possible cosmological coupling of BHs, we will consider the α = 0 case only.Then, the field equations can be rewritten in the form (see Ref. [21]) where the prime refers to derivation w.r.t.r, while g(r) is an integration function.The remaining equation, stemming from the conservation of the stress-energy tensor, is used to compute p ⊥ .The system (2.3) has to be closed by imposing additional external conditions like, e.g., the equation of state (EoS) of the fluid.The large-scale, cosmological dynamics, determining the scale factor a, can be completely decoupled from the equations.In this limit, our universe becomes homogeneous and isotropic, so that the dependence on the radial coordinate r is negligible.Moreover, we will only consider asymptotically flat sections at constant time, so that, for r → ∞, e α , e −β → 1.In this limit, eq.(2.2) describes a spatially flat, FLRW metric, solution of the system (2.3), which now reduces to the Friedmann equations where ρ 1 , p 1 are the usual energy density and pressure at cosmological scales, respectively.
Let us now consider the local, small-scale limit of eq. ( 2.3), describing compact, sphericallysymmetric objects at a time scale for which both α and β are almost independent of time.The system (2.3) describes now a static, spherically-symmetric spacetime, where eq.(2.3a) reads e −β(r, η 0 ) ≡ e −β 0 (r) = g(r) a rα ′ R , where η 0 represents a generic reference time.Inverting the above yields the function g(r), i.e., ln g = −β 0 − rα ′ ln a R .
The constant a R gives the scale factor at the reference time, which can be fixed as the present time, namely a R = 1.Equation (2.3a), thus, becomes where we defined k(r) ≡ rα ′ (r) . (2.6) At fixed time, the system of equations (2.3) reduces to where a i is the scale factor at a particular, but arbitrary, reference time η i .In order to describe the cosmological coupling, η i is identified as the initial time at which the compact object forms and couples to the cosmological background (see also Croker et al [1,5]).These equations allow for any static GR solutions, like singular or nonsingular BHs and compact objects (see, e.g., Refs.[19,23,24] and references therein).The specific EoS of the fluid sourcing the static solution at fixed time determines the relation between β 0 and α.The explicit form of the latter is, however, not constrained and is determined by the density profile ρ(r).
We now focus on eq.(2.3b), which can be written as (2.8) Using this equation, we can easily compute the comologically-coupled MS mass of the compact object, contained in a cosmological volume a(η) 3 L 3 , where L is the radius of the compact object, measured in the local r coordinate.We have where k L ≡ k(L) comes from eq. (2.6).Here we defined M (a i ) ≡ a i L/2G, which is the mass of the object computed at the coupling epoch.Peculiarly, this expression defines the proper Schwarzschild radius a i L = 2GM (a i ) at the initial time at which the compact object formed and became cosmologically coupled.
Equation (2.9) is our most important result.It shows that the mass of compact astrophysical objects is naturally coupled to the cosmological dynamics.With our formula, we are essentially describing the redshifting of the mass of the object with respect to its formation epoch.
Let us now briefly discuss the different contributions appearing in Equation (2.9).The first term gives a purely cosmological contribution, which depends on the background cosmological energy density ρ 1 and is therefore expected to be relevant only beyond the transition scale to homogeneity and isotropy.It will not play a role at the typical scales L of the compact object and, thus, we will not consider it in the reasoning that follows2 .The second term has the form of a "universal cosmological Schwarzschild mass".Finally, the last term encodes model-dependent corrections to the universal term.At fixed time, asymptotic flatness of the metric functions constrains e −β 0 (L) ≤ 1.This implies that, for k L > 0, the third term in eq.(2.9) gives a subleading correction to M (η) at any redshift z = (1 − a)/a and it decreases as z grows.It is important to stress again that eq.(2.9) is valid for any compact object (either singular or nonsingular) allowed by GR.In the following, we will apply separately our formula to the three possible cases: singular BHs, nonsingular BHs and regular horizonless (star-like) objects.

Singular BHs
For standard Schwarzschild BHs, the sum of the last two terms in eq.(2.9) is identically zero and we are left with the purely cosmological contribution.This is true at every fixed instant of time, as it can be seen from eq. (2.7a), by substituting e −β(r) = 1 − 2GM/r.Thus, Schwarzschild BHs decouple from the large-scale cosmological dynamics.Recent investigations [27] suggest that this is a quite general result, valid for generic, local and singular compact objects.The coupling is, however, present for nonsingular objects.In the following, we will consider this latter possibility.

Nonsingular BHs
When the compact object is a regular BH, L can be naturally identified with the radius of the apparent horizon (trapped surface).The subleading contribution in eq.(2.9) will, thus, always be zero in these cases.This means that the cosmological mass coupling is universal for BHs and our formula (2.9) becomes equivalent to eq. (1.1), with k = 1.

Horizonless compact objects
Let us now discuss the application of our mass formula to the compact, horizonless solutions of a notable class of regular models with EoS p = −ρ, which has been extensively used to construct regular GR solutions [19,[28][29][30][31][32][33].The latter are also of particular interest because, owing to their EoS, they are the GR solutions most similar to dark energy.
Using this EoS into the system (2.7), we have α = −β 0 .In the following, we will focus on two paradigmatic cases: models endowed with a de Sitter core [19] and those with an asymptotically Minkowski core [33].All the latter are characterized by their ADM mass M and by an external length scale ℓ, responsible for the smearing of the classical singularity.Depending on its value w.r.t. the Schwarzschild radius, these solutions can have either two or no horizons or a single degenerate one, representing the extremal configuration (for more details, see Refs.[19,33]).The metric describing all these models can be written in the general form [19] where the function F encodes the information about the MS mass m(r) at the fixed conformal time η i -the conformal time at which the object forms couples to the cosmological backgroundand, thus, realizes the interpolation between the behavior near the core (regularizing the r ∼ 0 singularity) and the asymptotic Schwarzschild solution at r → ∞.With these few ingredients, we can now consider the horizonless configurations of the models (2.10) and describe the behavior of the subleading, model-dependent term in eq.(2.9), µ(a) ≡ e −β 0 (L) a k L = e α(L) a k L , in all generality.
Since we only consider here the case of the absence of horizons, e α(L) is always positive.Additionally, given the asymptotic behavior at r ∼ 0 and r → ∞, it is always bounded by 1 from above, whatever the value of L. Evaluating the latter is anyway essential to assess the sign of k L , which in turn determines the variation of µ with the redshift.Since the density profile of the models considered here goes to zero only at r → ∞, we cannot define a hard surface for these stellar objects.We can, however, define an effective radius, containing ∼ 99 % of the ADM mass, namely m(L)/M = 0.99.The consequence of this prescription is that L is sufficiently big to ensure that eq.(2.10), evaluated at r = L, is well approximated by the Schwarzschild metric.This implies that F ′ (L/ℓ) < 0. k L is, thus, positive for all models.Therefore, µ(z) starts at a maximum value, less than 1, at z = 0 and then decreases with the redshift.
We confirmed these conclusions considering three specific models: a particular case of the Fan & Wang metric [34] (recently revisited in [35,36]), the Gaussian core regular model [32] and the asymptotically Minkowski core model of Ref. [33].We noted that, for all these, the subleading correction µ(z) has an order of magnitude which is comparable with 1 at z = 0.However, while the magnitude of µ stays close to 1 in a wide interval of redshifts for the first and last models (due to the slow fall of their density profiles at infinity), it decreases much faster for the Gaussian core model.

Can nonsingular compact objects contribute to the dark energy?
If the redshift of BH masses is described by eq.(1.1) with k = 3, as the observational evidence with elliptical galaxies found in Ref. [1] seems to suggest, then BHs should give a contribution to the cosmological constant in the Friedmann equations (see Ref. [1]).On the other hand, if k = 1, as we have shown to be the case adopting a GR framework to describe them, BHs contribute to the cosmological energy density as a curvature term, which redshifts as ρ ∼ a −2 .Therefore, the interpretation of vacuum energy in terms of GR BHs is quite problematic, and would imply a nontrivial averaging of inhomogeneities on large, cosmological, scales (see Ref. [21]).

Constraints from the growth of supermassive BHs in elliptical galaxies
Here, we compare our theoretical prediction with the mass measurements of SMBHs in elliptical galaxies at different redshifts.As done in Ref. [6], we select "red-sequence" elliptical galaxy hosts where the SMBH growth via accretion is expected to be negligible, and thus the measured BH mass growth with z can be considered exclusively due to cosmological coupling.We consider a set of low-(0.0016≤ z ≤ 0.19) and high-redshift (0.8 ≤ z ≤ 0.9) SMBHs.The low-redshift sample is directly taken from Table 2 of Ref. [6].The high-redshift sample is built by cross-matching two AGN catalogues, one from the Wide-field Infrared Survey Explorer (WISE) survey [37], and one from the Sloan Digital Sky Survey (SDSS) survey [38].The relevant properties of the host galaxies, such as, the star-forming rate (SFR), AGN reddening (E b−v ) and stellar mass (M * ) are taken from Ref. [37], while the BH masses from Ref. [38].The final high redshift sample M * − M BH used for our analysis is constructed following the procedure described below: 1.The objects from the above-mentioned catalogues are selected cross-matching their relative position (R.A. and Dec.) within 2 ′′ .The resulting sample consists of 160, 005 objects.
2. Only objects within the redshift range 0.8 ≤ z ≤ 0.9 are selected.This ensures that the measurements of the SMBH masses are estimated using the Mg II line.This cut reduces the sample to 605 objects.Here we show the unshifted final high-redshift sample (i.e., assuming zero cosmological couplings) (red circles), as well as the cases when it is shifted with k = 1 (black crosses), and k = 3 (gray dots).In addition, we show the unshifted low-redshift sample from Table 2 of Ref. [6] (blue squares).The vertical dashed line represents the low stellar mass cut of M * > 4 • 10 10 M ⊙ adopted in our data analysis.

5.
A selection on the SFR below a factor 5 from the best-fit SFR-M * is adopted to ensure the quiescence of the selected galaxies.This is done via the relation obtained in Ref. [39] (Equation ( 28)).This cut further reduces the sample size to 211 objects.
6. Finally, in order to select only elliptical galaxies, we reject objects for which the bolometric luminosity of the spiral component is brighter than 5% of that of the elliptical component 3 .
The final high-z {M * − M BH }-sample used in our analysis consists of 108 objects.It is shifted to z = 0 using Equation (1.1), both with k = 1 and k = 3. Figure 1 shows the unshifted final high-redshift sample, i.e., assuming zero cosmological coupling (red circles), as well as the cases when it is shifted with k = 1 (Sample-A, black crosses), and k = 3 (Sample-B, grey dots).In addition, we show the low-redshift sample from Table 2 of Ref. [6] with blue squares (Sample-C).Note that, in order to compare the high-and low-redshift data, also  Sample-C requires to be shifted to z = 0 via Equation (1.1) with k = 1 and k = 3.However, being all objects included in Sample-C at a redshift close to zero, the resulting shift in the M BH -axis of Figure 1 is not visible by eye.We now compute the probability that the two samples (Sample-A and -B) are drawn from the same sample as Sample-C using a two-dimensional Kolmogorov-Smirnov (KS) test [40,41].In order not to bias the fit, we impose a low stellar mass cut of M * > 4 • 10 10 M ⊙ , due to the fact that in the high-redshift sample there are very few objects below this value as in [6].
For a given value of k, we base our statistical analysis on the estimation of the minimum KS distance between the resulting low-and high-redshift samples, each of them shifted to z = 0 with the given k-value.
Per each k spanning the interval [0 − 8], we simulate 10 5 random realizations assuming a 2D Gaussian distribution for each point defined by [M * ± ∆M * , M BH ± ∆M BH ] in the 2D plane shown in Figure 1 4 .This ensures the inclusion of the uncertainty on the mass measurements into the analysis.Per each k-value that we considered, we are thus left with a distribution of 10 5 KS-distances, as shown in the right panel of Figure 2 for the k = 0, 1, 3 cases.Given a value of k, we derive the mean distance, 2σ (95%), and 3σ (99.7%) confidence limits.By iterating this procedure in the range k = [0, 8] we obtain the allowed 2-(dark red patch) and 3-σ (brown patch) bands as a function of k, as reported in the left panel of Figure 2. Figure 2 shows that the preferred value of k is around 3, since it corresponds to the smallest KS-distance.On the contrary, both the cases k = 0 and k = 1, consistent within 2σ with each other, are in tension with the obtained best-fit value of k.
While the preference for k ≃ 3 is consistent with the results of the analysis in Ref.
[1] with a similar data set from WISE/SDSS, our theoretical prediction is consistent with the BH mass-coupling test performed using higher redshift AGN and quasar data from the JWST [14] survey, the result of which is quite in tension with that of Ref. [1].The posterior probability distribution for k obtained in Ref. [14] is centered around a best-fit value k = −0.94± 1.19 at 68% confidence level, deviating from k = 3 at more than 3σ.Our theoretical prediction, k = 1, is compatible within 1σ with the JWST data.Let us note that the results presented in Ref. [14] are obtained with three AGN in early-type host galaxies, namely, CEERS01244, CEERS00397, and CEERS00717, and the conclusions are mainly driven considering the object CEERS01244 at redshift z = 4.8.However, the JWST is promisingly expected to increase the sample of a such measurements in the upcoming years.
The approach adopted in this work relies on the comparison between data sets of BH/early-type galaxies from different redshift regimes.Systematic uncertainties, possibly arising from the different observations and observables used to estimate the mass of the BH and the stellar mass of the host galaxy, may not be completely under control.For example, the JWST results come from a sample of stellar masses smaller than 10 10 M ⊙ .This range of masses is excluded from the analysis of Ref. [1], making our comparison between the two data sets more complicated.

Final remarks
In this letter, we have shown that GR, in its most general settings allowing nonsingular BH solutions, predicts a universal cosmological redshift of their masses, determined by the local curvature term.This consequently implies a linear scaling of the mass with the scale factor a, i.e., M ∝ a.In the case of horizonless regular compact objects, we have found a subleading, model-dependent, correction to the universal curvature term.
In order to check our theoretical prediction, we have performed a statistical analysis (KS-distance test) of the data of the mass growth of SMBHs in elliptical galaxies.We have compared the observed sample at low redshift with a sample at high-redshift, for several values of k ranging from 0 to 8. We have found the smallest KS-distance for k = 3, in accordance with the results of Ref. [1].On the other hand, our theoretical prediction, k = 1, is compatible within 1σ with the results of the statistical analysis performed in Ref. [14], based instead on the data of high-redshift AGNs detected by the JWST.Small values of k are also confirmed by the recent findings obtained with LIGO/Virgo/KAGRA data of binary BH mergers, which indicate a mild preference towards k ∼ 0.5 [5].
In view of this rather intricate situation, we need to consider other astrophysical observations and data, both of SMBHs and binary mergers, in order to assess which value of k is favoured by these samples.Our work will largely benefit from future observations.As underlined in Ref. [14] the JWST will observe a large number of high redshift galaxies and AGNs identifying a larger number of early-type host galaxies with higher stellar mass in the range of masses of the catalogs used in the present analysis.The improved sensitivity of GW observatories, such as LIGO/Virgo/KAGRA, will enable us to directly measure the masses of an increasingly large sample of stellar-mass BHs.Even more strikingly, the next generation of GW detectors, such as the Einstein Telescope and Cosmic Explorer, will enormously increase the number of detections of stellar mass binary BHs (105 ) per year up to higher and higher redshifts.Among these detections, a large fraction will have accurate mass estimates (see Ref. [42] for details).Moreover, the Laser Interferometer Space Antenna (LISA) will enable the access to massive BH coalescences.More precise mass estimates, broader coverage of mass range, and the use of different messengers will be crucial to constrain the cosmological coupling strength while keeping under control observational biases and possible systematics.
A preference for k = 0 from future observations will imply that nonsingular GR BHs are hardly compatible with the observational evidence.If data will, instead, favour k = 1, i.e., a non-zero cosmological coupling, apart from being important per sè, this result could represent striking evidence that GR BHs are indeed nonsingular objects.If this will be the case, we will still have to discriminate between objects endowed with a horizon from horizonless (stellar-like) ones.To this end, we will need a higher level of precision, sufficient to detect the subleading corrections derived in this work.Conversely, a confirmation from the data of either the k = 3 result of Ref. [1] or the k = −1 result of Ref. [14] will, instead, imply that BHs are nonsingular objects, which however do not admit an effective description in terms of solutions of Einstein's field equations.This will give a strong indication of the existence of new gravitational physics beyond GR at astrophysical scales.

3 .
All objects for which M * and/or M BH are undefined are removed, reducing the sample to 554 objects.4.A cut on the AGN reddening parameter E b−v < 0.2 is applied.This ensures to avoid reddened AGN and select unobscured AGN.This cut reduces the sample to 510 objects.

Figure 1 .
Figure 1.Here we show the unshifted final high-redshift sample (i.e., assuming zero cosmological couplings) (red circles), as well as the cases when it is shifted with k = 1 (black crosses), and k = 3 (gray dots).In addition, we show the unshifted low-redshift sample from Table2of Ref.[6] (blue squares).The vertical dashed line represents the low stellar mass cut of M * > 4 • 10 10 M ⊙ adopted in our data analysis.

Figure 2 .
Figure 2. The 95 % (dark red) and 99.7 % (brown) containment of the KS distance estimates for different k values using 10 5 realizations of low-redshift and high-redshift samples shifted to z = 0.The right vertical panel shows the distribution of KS distances simulated using the uncertainty on the mass of the black hole (∆M BH ) and the uncertainty on the mass of the star (∆M * ).