Cosmological black holes are not described by the Thakurta metric: LIGO-Virgo bounds on PBHs remain unchanged

We show that the physical conditions which induce the Thakurta metric, recently studied by B{\oe}hm {\it et al.} in the context of time-dependent black hole masses, correspond to a single accreting black hole in the entire Universe filled with isotropic non-interacting dust. In such a case, the physics of black hole accretion is not local but tied to the properties of the whole Universe. We show that radiation, primordial black holes or particle dark matter cannot produce the specific energy flux required for supporting the mass growth of Thakurta black holes. In particular, this solution does not apply to black hole binaries. We conclude that cosmological black holes and their mass growth cannot be described by the Thakurta metric, and thus existing constraints on the primordial black hole abundance from the LIGO-Virgo and the CMB measurements remain valid.


I. INTRODUCTION
If primordial black holes (PBHs) exist, their binary formation starts before the matter-radiation equality. The merger rate of these binaries and the resulting gravitational wave (GW) signal strength both depend on the PHB mass distribution and abundance [1][2][3][4][5][6][7]. In the standard picture, where PBHs formed spatially with a Poisson distribution, and their accretion rate is low, the PBH merger rate today would exceed the one indicated by the LIGO-Virgo observations [8,9] if a large fraction of dark matter (DM) was in ∼ 1 − 100M PBHs. The LIGO-Virgo observations, therefore, constrain the PBH abundance within such mass range [3-6, 10, 11]. In addition, the scenario where majority of the observed black hole (BH) mergers originated from PBHs is disfavoured while the presence of a PBH subpopulation within the observed BH mergers remains a viable possibility [11][12][13][14][15].
If the PBH masses change significantly during their evolution, e.g. due to PBHs accretion or mergers, GW signals from PBH binaries are also affected. This phenomenon by itself is simple and robust. However, its realization depends on many nontrivial details such as the nature and properties of DM, initial perturbations and structure formation and the highly non-linear accretion dynamics of BHs. The evolving mass will also affect other early Universe constraints on the PBH abundance, such as all the limits due to the PBH accretion [16][17][18][19][20][21][22].
In this paper, we consider the possibility that cosmological BHs are described by the Thakurta metric [23] where a is the scale factor, f (r) = 1 − 2m/r, H =ȧ/a and m is a constant. The Thakurta ansatz was claimed to imply significant modifications to the constraints on the PBH abundance due to GWs from PBH mergers [24] and due to energy injection into the CMB due to accretion [25]. The Thakurta spacetime is a generalized McVittie solution with accretion [26] 2 . The Thakurta metric approximates the FRW metric at r am. At r H −1 , it resembles a Schwarzschild BH, but with an evolving mass: the Misner-Sharp mass enclosed in a sphere of areal radius R ≡ ar is given by where ρ(R) ≡ρ/f (R) andρ denote the energy density of the ambient cosmic fluid at R and at spatial infinity. The first term gives the BH mass, while the last term can be interpreted as the mass of matter surrounding the BH. In the Newtonian limit, the Misner-Sharp mass coincides with the Newtonian mass at the leading order [27]. Thus, at distances M R 1/H, the Thakurta BH is approximately described by a Newtonian point particle with growing mass, The accretion rate of such BHs is proportional to the Hubble rate, i.e.,Ṁ = HM , thus tying together cosmological expansion and the local accretion physics of individual BHs.
In general relativity, the Thakurta metric (1) implies the spherical accretion that the spherical accretion of matter 3 into the BH follows the stress-energy tensor [26,29] T µν = (ρ + P )u µ u ν + P g µν + u (µ q ν) , where u µ u µ = 1, q µ u µ = 0, q µ describes a radial heat flow with q µ q µ = m 2 H 2 /(4πr 2 af ) 2 and P is the pressure of the cosmic fluid. This idealized radial energy flux is completely isotropic and supports the specific accretion of mass given in Eq. (3). As the Thakurta spacetime is constructed by first postulating the metric, the properties of the surrounding matter are implied by the Einstein equations and not by any microscopic model of matter. Most obviously, the metric (1) does not capture the complexities of accretion physics nor the effects that lead to anisotropies and structure formation. The Thakurta metric includes ambient cosmic fluid that is smooth at the scale of the BH horizon and thus can not correspond to PBHs. Thus, all DM can not be in PBHs, but an additional smooth DM component, e.g. , particle or fuzzy DM, is required. We can naively construct a universe filled with Thakurta PBHs by gluing together several Thakurta spacetimes. 4 The energy within a spherical volume V is given by the Misner-Sharp mass (2). If the separation between PBHs is much larger than the horizon, then f ≈ 1 and thus the ambient energy density is roughly constant, ρ ≈ρ. Within this approximation, we can consider non-spherical patches with the energy within a patch of volume V i given by M i ≈ m i a +ρ V i . The average energy density of a universe containing multiple Thakurta patches is then where n PBH is the PBH number density, or equivalently the number density of Thakurta patches, and m a is the average PBH mass. Conservation of PBH number, i.e., n PBH ∝ a −3 implies that ρ PBH ≈ m a n PBH ∝ a −2 .
As this scaling is due to the energy transfer from a smooth DM component to the PBHs, it is possible that their combined energy density still scales as ρ tot ∝ a −3 in which case the DM interpretation would remain valid. In this case 5 , it is relevant to understand the physics behind this energy flow. According to Eq. (3), the PBH mass growth is quite extreme. For example, between the matter-radiation equality and the redshifts relevant for currently detectable BH mergers, the mass of any PBH would have grown ∼ 3400 times. It was argued in Ref. [24] that the mass differences are milder as mass growth should stop once virialized DM halo absorbs the PBH because the Thakurta metric does not describe BHs in such environments. Thus, the postulated effect can be softened further if the PBHs make up most of DM, since then the first virialized structures can form already around matter-radiation equality [5,[30][31][32]. By a similar argument, Eq. (3) may fail already for as small structures as the PBH binaries or even for BHs that simply move with respect to the cosmic fluid, since such BHs are also not described by the Thakurta metric.
In the following, we argue that the authors of Refs. [24,25] postulated a particular mathematical construction, the metric (1), to describe PBHs without checking whether this choice is consistent with our knowledge of cosmology and accretion physics. This allowed them to claim that the constraints on the PBH abundance derived from the LIGO-Virgo observations of BH coalescence are eliminated, opening the possibility that all DM can consist of 30−100 M mass PBHs. This paper aims to point out some phenomenological consequences of the claims made in Ref. [24,25] and to discuss the claimed accretion of PBHs and its effects on the LIGO-Virgo bounds. By studying the consequences of the accretion rate (3), we demonstrate internal inconsistencies of those claims as well as inconsistencies with cosmology, many-body dynamics of PBHs and accretion physics.

II. PHENOMENOLOGY
It is not mathematically forbidden to have a configuration of matter around the PBHs that supports the metric (1). Whether such a configuration of matter is physically meaningful, what are the needed properties of such a matter, and which are the conditions for accreting it at the rate (3) have not been questioned and answered in Ref. [24]. In the following, we attempt to address some of those questions. We will consider general cosmological scenarios with Thakurta-like accretion. Such scenarios contain the Thakurta PBHs, but are not restricted to it and can correspond to, e.g. , perturbations of the exact solution. We conclude that imposing accretion by PBHs, determined by Eq. (3), is unphysical.
We start by considering the case where all DM consists of PBHs, f PBH = 1, where f PBH is the fraction of DM in PBHs. The central claim and the main result of Ref. [24] is that, due to the accretion rate (3), f PBH = 1 is possible for the LIGO BHs in the mass range 30 − 100 M . As there is no particle DM, the PBHs are initially surrounded by the radiation of the Standard Model particles. Around the time of matter-radiation equality, the radiation redshifts and, after the recombination around z ∼ 1100, the Universe is filled with PBHs and a strong wind of baryons moving with the velocity v ∼ 30 km/h created by the baryon acoustic oscillations [33].
The immediate question arises: in this case, what do the PBHs accrete to grow in mass by several orders of magnitude by the present time? The abundance of baryons is insufficient for that, and, in the early Universe, the baryons move with the velocities of the order v ∼ 30 km/h which makes them uncapturable by the PBHs. The only possible answer is that the PBHs accrete themselves. We note that this possibility must already deviate from the exact Thakurta metric as PBH DM is certainly not smooth at scales comparable to the BH horizon, and its influx is not described as a heat flow. Regardless, a sufficiently large PHB merger rate needed to support the mass growth of Eq. (3) is impossible for initially Poisson distributed PBHs. This process must, therefore, be described by hierarchical PBH binary mergers whose rate depends on the initial properties of the PBH population, e.g. , their mass function or their spatial distribution. Although the PBH merger rate can be significantly enhanced for initially clustered PBHs, [3,32,[34][35][36][37] which, by itself, is constrained by the CMB observations [38], it will ultimately be determined by the local parameters shaping the small scale structure of PBHs instead of the specifics of cosmological expansion, which governs Eq. (3).
In general, PBH mergers cannot grow the average PBH mass by several orders of magnitude. However, even if this was the case, an enormous stochastic GW background would be generated, as each generation of mergers, corresponding to an O(2) increase of the average PBH mass, converts O(5%) of the PBH DM into GWs. In either case, the mass growth asserted by Eq. (3) contradicts studies of PBH structure formation or the GW observations.
Next, consider the case f PBH < 1. The neighbourhood of a BH can be modelled using Newtonian physics at distances M r H −1 . Numerical simulations show that early DM haloes form during the radiation dominated epoch when wider and wider shells of DM surrounding the PBHs decouple from expansion [39]. This results in a density profile that scales as ρ(r) ∝ r −9/4 . Such haloes are stable, but can be disrupted by later close encounters with compact objects, which, when f PBH 1, are mostly stars. Again, the relative stability of such DM haloes is in contradiction with the mass growth ma(t). This stability can be understood intuitively by noting that the DM particles will generally not fall into the BH but form a halo due to tidal torque generated by surrounding inhomogeneities. Such orbiting particles are unlikely to fall into the PBH unless there are mechanisms that carry away their energy and angular momentum. The accretion rate (3) is therefore not realized.
Taking the accretion rate seriously, one finds that the microscopic theory of the cosmological fluid must have unrealistic properties also for other backgrounds. For example, in a vacuum energy dominated universe, the growth of a Thakurta BH must due to a heat flow carried by vacuum energy. The energy density of the source for the Thakurta metric ρ(R) ≡ρ/f (R) grows towards the BH and diverges at m = r. Thus, in a radiation dominated universe in which T ∝ ρ 1/4 , the exact Thakurta solution requires an inflow of heat against a temperature gradient in contradiction with the second law of thermodynamics.
Finally, ignoring the problems with physics behind the mass growth (3), we find (see the Appendix) that the radius of BH binaries with an evolving mass (3) shrinks as r ∝ a −3 . This is a significant modification to binary evolution that can enhance the merger rate. Thus the LIGO-Virgo constraints on PBHs, especially from the stochastic GW background, are not obviously avoided. We must stress again, however, that the Thakurta metric does not describe a binary in the first place and, moreover, accretion into BH binaries can be significantly more involved than accretion into a single BH (see e.g. [40][41][42]).

III. ACCRETION PHYSICS
The radical departure of the BH mass growthṀ = HM from the conventional wisdomṀ ≈ 0 is due to the accretion of the surrounding energy into the BH, as described by the metric (1). Two main approaches have been followed in the studies of the highly nontrivial problem of realistically embedding BHs into the cosmological background, mathematically exact solutionsá la McVittie and physical approximationsá la Zel'dovic.
The McVittie metric [43] is an exact solution reducing to the Schwarzschild BH at small radii and to the FLRW at large. It does not, however, allow for accretion [26], and indeed, as most recently shown in Appendix B of Ref. [32], the use of the exact McVittie Ansatz leads to the conclusion that the evolution of the BH mass M can be neglected for practical purposes. This supports the generally accepted view that BHs as local systems are decoupled from the global cosmological expansion, but it does not address the issue of accretion. The generalised McVittie metric of the form (1) seems to describe an accreting BH, but at the price of introducing an imperfect matter source [44,45]. Namely, a radial energy flow in the surrounding fluid is required to incorporate the effect of accretion. It was initially thought that with this amendment, one could alleviate the unphysical properties of the original McVittie solutions, which exhibit spacelike singularities at the horizon and divergent pressures in the matter sector. However, that has been corrected by later calculations, which clarify that the Thakurta geometry does not avoid the unphysical singularities [28,29]. The Thakurta metric with H = 0 does not describe a BH, but an inhomogeneous expanding universe [28,29]. In this light, it is difficult to justify the proposal of Refs. [24,25] that (1) would more realistically describe cosmological BHs.
When considering perfect fluid sources, PBH accretion can be studied approximately. Bondi had considered accretion in a spherically symmetrical star in a steady-state Universe and found it to be proportional to the square of the mass of the star [46]. Zel'dovic and Novikov adapted the result to an expanding Universe by arguing that the system is decoupled from the cosmological expansion [47]. The BH eats matter from within the accretion radius that, in the case of radiation, is thrice the Schwarzschild radius. The model was refined by Carr and Hawking [48], and further generalised by, e.g. , Babichev et al. [49,50], who arrived at the equation (see Section II of Ref. [51] for more references and discussion) where γ is the barotropic index of the perfect fluid and ρ is its cosmological energy density. For dust γ → 1, the accretion radius would be infinite, and for de Sitter γ → 0, so there is no accretion. Thus the qualitative behaviour of Eq. (7) is reasonable in those limits. The most relevant case is radiation γ = 4/3, for which we can easily integrate the equation to obtain where M 0 is the mass of the BH at its formation time t 0 . If the BH is much smaller than the horizon, M 0 t 0 , this implies negligible accretion. However, the accretion becomes quite significant for PBHs of size comparable to the horizon. In fact, the (special and fine-tuned) selfsimilar solution yields M ∼ a 2 m. These special solutions may not be physical, as was argued by Carr and Hawking by taking into account the cosmic expansion (which seems reasonable once M 0 ∼ t 0 ). Further, hydrodynamical simulations suggest that a horizon-sized PBH in a perturbed positive-pressure fluid shrinks in relation to the horizon. In conclusion, theoretical justifications to ignore the mass growth of PBHs due to accretion appears fairly robust, also in the Bondi-Zel'dovic-Novikov-based approach.
Nevertheless, in Section II and Appendix A we considered the implications of the assumption that M = am, concluding that such an assumption would not be phenomenologically viable.

IV. CONCLUSIONS
We conclude that the Thakurta metric does not describe realistic cosmological black holes. Most impor-tantly, the case f PBH = 1 when all the DM is in the form of PBHs, which was the main aim of Ref. [24] to postulate the Thakurta metric, is intrinsically and phenomenologically inconsistent. Therefore, the existing bounds on the PBH abundance arising from the LIGO-Virgo GW observations and from the BH accretion physics, including the CMB measurements, remain valid.
Note added: After the initial version of this paper, a critical comment appeared in Ref. [52]. We replied to their criticism in Ref. [53]. Around the same time, another paper appeared arguing that the Thakurta metric does not describe a cosmological black hole [54]. binary E = T + V can be split into the kinetic and potential energy, which depend on the scale factor as T ∝ a and V ∝ a 2 . In Lagrangian mechanics, time dependence of the energy is by the explicit time-dependence of the Lagrangian, i.e.,Ė = −∂ t L = H(−T + 2V ). If the orbital period is much shorter than the Hubble time, as expected for a decoupled binary, then the energy is approximately conserved within a single oscillation. Thus, by the virial theorem, T = −V /2. This givesĖ = 5HE. As E = −a 2 m 1 m 2 /r a , with r a the semimajor axis, we obtain that the binary shrinks as r a ∝ a −3 . (A2) This behaviour can be generalized to the evolution of the size of any virialized system comprising PBHs obeying M ∝ a. The argument follows along the same lines as above: the relations T ∝ a and V ∝ a 2 and T = − V /2 implyĖ = 5HE. Since the gravitational potential is inversely proportional to the size of the system, V ∝ R −1 , we obtain R ∝ a −3 . Thus, bound systems of PBHs, which can form already at the onset of matter-radiation equality, would shrink approximately as R ∝ a −3 if the constituent PBHs accrete according to M ∝ a.
The coalescence time of constant mass binaries scales as τ ∝ r 4 a . Therefore, Eq. (A2) translates into a rapid decrease of the coalescence time τ ∝ a −12 , and in the scenario proposed in Ref. [24], binaries will shrink due to mass growth until GW emission takes over, after which the binary merges rapidly. For illustration, Fig. 1 shows the evolution of the semimajor axis of a circular 7 binary with constant mass and with M ∝ a. The latter evolves asṙ a = −3Hr a − (16/5)M (t) 3 r −3 a , where the last term comes from GW emission [55]. Therefore, the mass growth postulated in Eq. (3) can enhance the GW emission from PBH binaries instead of suppressing it. Moreover, this enhancement can be sufficiently strong to force most binaries to merge by the present time, thus reducing the PBH merger rate at low redshifts. In this case, constraints on the PBH abundance arise from the non-observation of a stochastic GW background.