Further Mystifying the Hubble Tension: Invalidating the Black Hole Spacetime Sink Model, Painlevé–Gullstrand Congruences as a Coarse-Graining Tool, and N-Cell Buchert Backreaction Insignificance

Abstract

We investigate whether black holes, interpreted through the Painlevé–Gullstrand (PG) “river model” of spacetime, can provide a mechanistic picture of cosmic expansion and generate cosmologically significant backreaction. In Phase 1, we show that the PG infall velocity at the Wigner–Seitz cell boundary of a lattice of identical masses exactly reproduces the Friedmann equation: H² = 8πGρ/3. A critical review confirms this is a mathematical identity via Birkhoff’s theorem—not new physics—and the model is reframed from “BH sinks drive expansion” to “PG congruences as a coarse-graining tool.” In Phase 2, we test whether the observed BH mass function (spanning ~10 to ~10⁹ M_☉) generates non-trivial Buchert backreaction using an N-cell Voronoi–Buchert averaging scheme. The BH mass spectrum contributes only ~10% of the total expansion variance; spatial clustering dominates (~89%). After shear cancellation, the net backreaction is consistent with |β_net| ~ 10⁻⁵, far below the Hubble tension threshold. We conclude that BH demographics cannot explain the Hubble tension through Buchert backreaction, but the PG/Wigner–Seitz framework is a valid coarse-graining tool for cosmology.

Phase 1: The River Model of Gravity

In 2004, Hamilton & Lisle formalized the “river model” interpretation of general relativity: in Painlevé–Gullstrand coordinates, spacetime around a black hole can be described as flat space flowing radially inward toward the singularity. This is not a mere analogy—it is an exact coordinate transformation of the Schwarzschild solution. The velocity of this inflow is:

v(r) = −√(2GM/r)

At r = r_s (the Schwarzschild radius), space flows inward at the speed of light. Inside the horizon, |v| > c, making escape impossible—not because of a “force,” but because the river of space itself exceeds lightspeed.

However, this description is gauge-dependent—it is a property of the PG coordinate choice, not an invariant physical statement. In Schwarzschild coordinates, the exterior is static; there is no “flow.” The “consumption of space” is coordinate ontology, not physical ontology. What is physically invariant is that the event horizon removes causal accessibility: signals and matter that cross it cannot return to external observers.

The Wigner–Seitz Cell Model

Consider N identical Schwarzschild black holes, each of mass M, uniformly distributed throughout a volume V. We approximate each BH’s domain of influence as a spherical Wigner–Seitz cell of radius R = (3V/(4πN))^1/3. For any physically realistic BH population, the cell radius is vastly larger than the Schwarzschild radius—for a 10 M_☉ BH at the critical density, R/r_s ~ 5 × 10¹⁶, making the weak-field approximation spectacularly good.

At the cell boundary r = R, the effective Hubble parameter is defined as the velocity-to-distance ratio:

H_eff = |v(R)|/R = √(2GM/R³)

Squaring and substituting M = ρ · (4/3)πR³:

H_eff² = 8πGρ/3

Velocity Field Superposition

For BHs on a cubic lattice, the total PG velocity field is the vector sum of individual contributions. In the weak-field regime (v/c ~ 10⁻¹⁶ at the cell boundary), linear superposition is an excellent approximation with corrections of order (r_s/r)² ~ 10⁻³³.

At the center of each cell, cubic symmetry forces v = 0 (the “comoving observer” position). At the cell boundary, the nearest BH dominates and the velocity is directed radially inward—from the perspective of the cell center, this appears as outward flow, reproducing the Hubble flow. The superposed field is irrotational (∇ × v = 0), consistent with the FLRW requirement of curl-free expansion.

Observed BH Demographics

How much of the observed expansion can BHs actually account for? We compare the required and observed BH number densities across five mass classes:

BH Class	M (M☉)	nreq (Mpc−3)	nobs (Mpc−3)	Gap Factor
Stellar	10	1.26 × 1010	~5 × 106	2.5 × 103
Heavy stellar	30	4.20 × 109	~106	4.2 × 103
Intermediate	103	1.26 × 108	~0.1	1.3 × 109
SMBH	106	1.26 × 105	~0.02	6.3 × 106
Massive SMBH	109	126	~10−5	1.3 × 107

The total observed BH mass density is ρ_BH,obs ≈ 8 × 10⁷ M_☉ Mpc⁻³, a factor of ~1600 below ρ_crit ≈ 1.26 × 10¹¹ M_☉ Mpc⁻³. Observed BHs contribute H_BH = 1.70 km/s/Mpc—only 2.5% of the observed H₀ = 67.4 km/s/Mpc. Even placing all matter (baryons + dark matter) in BHs yields only H = 37.8 km/s/Mpc (56% of H₀), leaving no room for dark energy.

Critical Review and Reframing

A rigorous theoretical review identified six fundamental problems with the original “BH sinks drive expansion” claim. The reviewer was correct on all major points:

Birkhoff’s theorem: A black hole and a dust ball of equal mass M produce identical exterior gravitational fields. Our simulation confirms H_eff^BH/H_eff^dust = 1.000000000000000 to machine precision. The Friedmann equivalence depends on ρ (mean density), not on whether that mass is concentrated in event horizons.

Gauge dependence: The PG “river” description is a coordinate artifact. In Schwarzschild coordinates, the exterior is static—there is no flow. The gauge-invariant content is the tidal tensor (Weyl curvature) and the expansion scalar of specific congruences.

No sink term: The Raychaudhuri equation governing the expansion scalar contains no “sink term.” For the PG congruence in Schwarzschild vacuum, R_μν = 0 (Ricci-flat), and expansion dynamics are purely kinematic. The v1 “space flow continuity equation” was retracted—there is no conserved quantity “ρ_space” in general relativity.

Superposition limits: GR is nonlinear. The superposed PG velocity field is a Newtonian-order approximation, not an exact GR solution. The correct exact treatment is the Einstein–Straus Swiss-cheese model, where Schwarzschild vacuoles are embedded in FLRW with junction conditions.

Backreaction magnitude: The literature consensus is |Q_D|/H² ≲ 10⁻⁵ (Green & Wald 2014) to ~10⁻² (optimistic Buchert estimates)—far too small to explain dark energy.

Kerr vorticity: Angular momentum conservation requires vorticity to decay as a⁻² in an expanding universe, eliminating any cosmologically significant accumulation.

Phase 2: Does the BH Mass Spectrum Generate Backreaction?

Phase 1 showed that identical BHs on a lattice exactly reproduce Friedmann with zero backreaction. The critical review recommended testing whether the heterogeneous BH mass function—spanning eight orders of magnitude from stellar remnants to supermassive black holes—produces additional Buchert backreaction beyond what standard density perturbations generate. This is the Phase 2 question.

The Buchert Backreaction

The full Buchert backreaction is the difference between expansion variance and shear:

Q_D = ⅔(⟨θ²⟩_D − ⟨θ⟩_D²) − 2⟨σ²⟩_D

where θ = 3H is the expansion scalar and σ² is the shear scalar. We compute only the first term (expansion variance):

Q_var = 6 ∑_i F_i (H_i − H_D)²

and define β_var = Q_var/H_D².

Algorithm: N-Cell Voronoi–Buchert Averaging

The algorithm proceeds in six steps: (1) Generate N BH positions in a periodic box with masses drawn from observed demographics. (2) Mass normalization: scale total mass to Ω_mρ_critV_box—the standard N-body cosmology approach where each “particle” represents a chunk of cosmological matter, not a literal black hole. Mass ratios (the spectrum shape) are preserved. (3) Compute a Voronoi tessellation with periodic boundary conditions. (4) For each cell i:

H_i = √(2GM_i/R_i³) = √(8πGρ_i/3)

(5) Domain-average with volume weighting (F_i = V_i/V_total). (6) Report β_var = Q_var/H_D².

Mass and Spatial Models

Four mass models and three spatial models are tested to isolate the contribution of each physical effect:

Mass Model	Description
uniform	All masses equal (10 M☉, pre-normalization)
lognormal	Log-normal around 10 M☉, σ = 0.5 dex
bimodal	99.9% stellar (~10 M☉) + 0.1% SMBH (~3 × 106 M☉)
observed	5-class Phase 1 demographics, number-weighted

Spatial Model	Description	Purpose
lattice	Regular cubic grid	Control (QD = 0 for uniform mass)
uniform	Poisson random	Isolate spatial variance
clustered	Nested Poisson, rcl ~ 1 Mpc	Realistic clustering

Population and Voronoi Tessellation

The primary experiment uses N = 1000 BHs with observed demographics and clustered positions in a 50 Mpc box. After mass normalization, the total density matches Ω_mρ_crit = 3.97 × 10¹⁰ M_☉ Mpc⁻³. The Voronoi tessellation produces cells with volumes spanning four orders of magnitude (0.024 to 3910 Mpc³). Dense cluster cores have very small cells (high ρ_i, high H_i), while inter-cluster voids have very large cells (low ρ_i, low H_i). Since volume-weighting favors the large, underdense cells, H_D = 16.4 km/s/Mpc ≪ H_Friedmann = 37.8 km/s/Mpc (a consequence of Jensen’s inequality, not backreaction).

Buchert Averaging Diagnostics

The diagnostics reveal two important features of the N-cell calculation: the backreaction variance is diffuse—764 cells (76%) contribute 90% of Q_var, so no single cell dominates—and the estimated shear 2⟨σ²⟩ exceeds the expansion variance, implying Q_net < 0. This near-cancellation is the standard result in the backreaction literature.

Control Case Validation

A critical test is recovery of known limits. Einstein & Straus (1945) proved that spherical vacuoles embedded in FLRW produce zero backreaction. Our uniform-mass lattice control (Case A) gives β = 5.92 × 10⁻³¹—effectively machine zero, confirming the Swiss-cheese limit.

Case	Mass Model	Spatial Model	βvar	Status
A	uniform	lattice	5.92 × 10−31	PASS (= 0 to machine precision)
B	lognormal	lattice	2.35 ± 0.23	PASS (mass variance only)
C	uniform	uniform	0.27 ± 0.03	PASS (position variance only)
D	observed	clustered	24.9 ± 1.5	Primary result

Primary Result and Decomposition

The primary result for observed demographics with clustered positions is:

This decomposes into three sources:

Source	βvar Contribution	Fraction
Poisson spatial variance	0.27	1%
BH mass spectrum	2.35	10%
Spatial clustering	~22.3	89%
Total	24.9	100%

The BH mass spectrum contributes only ~10% of the total expansion variance. Spatial clustering dominates. This immediately rules out the BH mass function as a special source of backreaction.

Parameter Sweep

A systematic sweep across mass models, spatial models, particle number, and box size reveals several patterns: the bimodal model gives β ~ 800 due to the extreme mass ratio (10 M_☉ vs. 3 × 10⁶ M_☉); clustering amplifies β by ~3–8× over Poisson positions; and β stabilizes for N ≥ 1000.

Mass Model	Spatial Model	NBH	βvar
uniform	lattice	1000	5.9 × 10−31 ± 0
uniform	uniform	1000	0.27 ± 0.03
lognormal	lattice	1000	2.35 ± 0.23
lognormal	uniform	1000	2.72 ± 0.26
bimodal	lattice	1000	802 ± 371
bimodal	uniform	1000	767 ± 383
bimodal	clustered	1000	6161 ± 3648
observed	uniform	1000	1.43 ± 0.08
observed	clustered	1000	24.9 ± 1.5

Shear Estimate and Net Backreaction

The estimated inter-cell tidal shear from nearest-neighbor differential tidal fields gives ⟨σ²⟩ ≈ 1.43 × 10⁻³⁴ s⁻². For comparison, H_D² = 2.84 × 10⁻³⁷ s⁻², so ⟨σ²⟩/H_D² ~ 500. The shear contribution 2⟨σ²⟩ greatly exceeds the expansion variance Q_var, implying:

|β_net| = |Q_var − 2⟨σ²⟩| / H_D² ≪ β_var

This is consistent with the standard result in the backreaction literature: the expansion variance and shear are individually large but nearly cancel, leaving |β_net| ~ 10⁻⁵.

Time Evolution

Time evolution from z = 100 to z = 0 uses linear perturbation theory with 45 active mass bins. Volume fractions use the exponential form F_b ∝ exp(−δ_bD) to ensure positivity even in the nonlinear regime. At z = 0, 16 of the 45 bins have |δ_b| > 1 (deeply nonlinear), so the linear evolution is approximate—it captures the qualitative growth of Q_D with structure formation but underestimates the z = 0 value relative to the full static Voronoi calculation.

Conclusion

Phase 1 established that the PG infall velocity at the Wigner–Seitz cell boundary exactly reproduces the Friedmann equation—a mathematical identity via Birkhoff’s theorem, not new physics. Review determined that this cannot replace dark energy, that the “space flowing into BHs” picture is gauge-dependent, and that no “sink term” exists in the Raychaudhuri equation. The model was reframed as a coarse-graining tool.

Phase 2 asked whether the heterogeneous BH mass function generates non-trivial Buchert backreaction. The answer is no: the BH mass spectrum contributes only ~10% of the total expansion variance; ~89% comes from spatial clustering. After shear cancellation, the net backreaction is consistent with |β_net| ~ 10⁻⁵, far below the Hubble tension threshold.

Limitations: (1) Expansion variance only—not the net Q_D. (2) Newtonian order (v/c ~ 10⁻¹⁶ for cosmological cells). (3) Gauge-dependent velocity field (PG coordinates, not invariant). (4) Birkhoff’s theorem: BHs are gravitationally identical to dust of equal mass. (5) No lightcone (static snapshot at z = 0). (6) Approximate shear (order-of-magnitude estimate). (7) BH mass growth is too slow to mimic dark energy (w > −1/3 for any realistic growth model).