Long-term 3+1 simulations of primordial black hole formation during radiation domination

This paper develops an efficient three-dimensional numerical-relativity framework for simulating primordial black-hole (PBH) formation from superhorizon curvature perturbations in a radiation-dominated Universe. The authors implement flux-conservative relativistic hydrodynamics in the adaptive-mesh-refinement code GRChombo and introduce a cosmologically scaled Gamma-driver that allows the cosmic-time step to grow in proportion to the scale factor. For a representative long-term simulation, the scaled driver preserves the apparent-horizon mass evolution and constraint behavior while reducing the number of coarse-level advances by a factor of approximately 94 relative to the standard driver. They also construct a conformal-time version of the moving-puncture gauge as an independent check. Applying the framework to a spherical Gaussian curvature profile, they find a collapse threshold 0.79578 < µc < 0.79580 and a critical exponent γ ≃ 0.3559, consistent with previous spherically symmetric results. They further fit the late-time PBH mass growth to the Zel’dovich–Novikov accretion law, demonstrating that the code can follow both near-critical collapse and long-term post-formation evolution in three dimensions. The framework provides a foundation for future studies of PBH formation beyond spherical symmetry.

The paper begins by noting that PBHs may form from the gravitational collapse of large primordial curvature perturbations during radiation domination, and that accurate determination of the collapse threshold is important because the predicted PBH abundance is exponentially sensitive to the formation criterion. Most calculations have imposed spherical symmetry, while fully three-dimensional simulations remain rare. PBH formation poses a demanding multiscale problem, requiring a simulation to encompass an initially superhorizon curvature perturbation while resolving the much smaller local length scales associated with collapse after horizon entry. The authors note that adaptive mesh refinement (AMR) addresses the spatial scale separation, but the remaining challenge is temporal, since the coordinate speeds of physical modes decline inversely with the scale factor as a−1, making a fixed time step increasingly conservative relative to the physical Courant–Friedrichs–Lewy (CFL) limit.

The authors implement a flux-conservative relativistic-fluid evolution module in GRChombo and construct a complete framework for simulating PBH formation. They introduce a cosmologically scaled Gamma-driver with b ∝ a−2 and ηGD ∝ a−1, which permits ∆t ∝ a∆x, corresponding to an approximately constant conformal-time step while retaining cosmic time as the evolution coordinate. They also construct a conformal-time version of the moving-puncture gauge as an independent gauge and implementation check. The paper states: "The decreasing coefficients do not imply that the driver becomes ineffective on the expanding background. Since the background proper distance is dl = a dx, the gauge-mode propagation speed measured with respect to this proper distance is avL and remains constant under the scaling in Eq. (36)."

For the numerical setup, the authors use the BSSN formulation for spacetime evolution, with a perfect fluid described by a linear barotropic equation of state p = ωρ with ω = 1/3 for radiation. The initial data is constructed using the long-wavelength, or gradient-expansion, solution, with a spherical Gaussian curvature profile ζ(r) = µ exp(−k2r2)W(r), where µ is the amplitude and W(r) is a window function. The fiducial initial-data parameters are Hi = 50, L = 1, rW = 0.8L, with rm = √6/10, ti = 0.01, and tH = MH = 1.5. The base grid contains 643 cells on the octant, with five additional AMR levels and a refinement ratio of two, yielding the finest grid spacing ∆xmin = L/2048.

The gauge performance tests show that the scaled and standard drivers give nearly identical histories for subcritical and supercritical configurations. The apparent-horizon masses are robust across all gauge choices. The computational advantage is demonstrated by the step-count comparison: reaching t/tH ≃ 266.8 requires 1092 coarse-level advances with the scaled driver and 102436 with the standard driver, a reduction by a factor of approximately 94. The wall-clock time is also reduced by a factor of approximately 92.

The paper also compares the cosmologically scaled cosmic-time gauge with the conformal-time gauge, finding that the apparent-horizon mass histories agree well through the evolution. For the representative µ = 0.825 pair, the two prescriptions have comparable costs, with the cosmic-time and conformal-time runs requiring 1092 and 1090 coarse-level advances, respectively.

Boundary-condition tests show that the reflection conditions used in the fiducial octant simulations reproduce the full periodic solution for the spherical configuration over their common valid interval. The authors introduce an excision-like interior regularization procedure to extend the full periodic-box evolution, which multiplies the right-hand side of every evolved variable by a window function inside the apparent horizon. The regularized full-box runs reach t/tH = 266.8, and their apparent-horizon masses all agree with the octant result.

For the collapse threshold, the authors use an additional refinement level and halve the AMR refinement thresholds, giving ∆xmin = L/4096. An amplitude scan finds that the threshold lies between µ = 0.79578 and 0.79580. This value refines the earlier three-dimensional result 0.795 < µc < 0.805 obtained with the COSMOS code and is very close to the spherically symmetric Misner–Sharp result µc ≃ 0.79579 ± 1 × 10−5 for the same curvature profile. Near the collapse threshold, the PBH mass obeys the scaling relation MBH = MH K(µ − µc)γ. A joint fit of (K, µc, γ) gives K ≃ 3.624, µc ≃ 0.7957813, and γ ≃ 0.3559, consistent with the universal radiation-fluid value γ ≃ 0.3558.

For the post-formation mass growth, the authors follow the apparent horizon long after formation and measure the declining accretion rate. They compare the numerical mass evolution with the late-time Zel’dovich–Novikov prescription, fitting the effective accretion rate ṀBH = 4πF RBH2 ρb(t) = (3F MBH2)/(2t2), where F is a dimensionless effective efficiency. The fitted efficiencies span 2.958 ≲ F ≲ 3.613, with the three larger amplitudes clustering in the narrower range 3.532 ≲ F ≲ 3.613, consistent with the range 3.5 ≲ F ≲ 3.75 found by Escrivà in the Misner–Sharp simulations. The near-threshold µ = 0.805 run gives the smaller value F = 2.958. The authors also observe a common crossing of the Ψ curves at t/tH = 50.6, where Ψ ≃ 0.62.

The paper concludes that this work presents the first determination of both the critical mass scaling relation and the post-formation accretion law from three-dimensional PBH formation simulations. The framework can be extended to other initial-data families and barotropic perfect fluids with other constant values of ω, including nonspherical profiles, multi-peak curvature perturbations, and other cosmological scenarios such as massless scalar-field collapse, kination-dominated cosmologies, and quartic scalar-field collapse.

Improvement 1: Adaptive Temporal Gauge Control for Multiscale Physical Simulations

The AI system can implement a self-adjusting time-stepping mechanism that scales the gauge-driver coefficients (e.g., b ∝ a-2, η GD ∝ a-1) based on the expanding background metric, allowing the cosmic time step to grow proportionally with the scale factor. This enables the system to maintain accuracy (e.g., apparent-horizon mass evolution, constraint behavior) while reducing computational cost by 94× in long-term cosmological or astrophysical simulations with large dynamic ranges.

Improvement 2: Autonomous Gauge Selection and Validation via Cross-Checks

The AI can automatically construct and compare multiple gauge choices (e.g., cosmologically scaled cosmic-time vs. conformal-time moving-puncture) for the same physical setup, validating results (e.g., mass histories, collapse thresholds) and selecting the most computationally efficient yet accurate gauge for a given problem. This reduces user intervention and improves robustness in general-relativistic simulations.

Improvement 3: High-Precision Critical Collapse Threshold Detection

The AI can perform automated amplitude scans with adaptive mesh refinement (e.g., increasing resolution and lowering refinement thresholds) to bracket critical collapse thresholds (e.g., 0.79578 < μ c < 0.79580) and extract critical exponents (e.g., γ ≃ 0.3559) via joint fitting of scaling relations. This allows the system to predict PBH abundance with exponential sensitivity accurately, even in 3D without spherical symmetry.

Improvement 4: Long-Term Post-Formation Accretion Modeling

The AI can track apparent horizons over extended cosmic times, measure declining accretion rates, and fit them to analytic laws (e.g., Zel’dovich–Novikov Ṁ BH = 3F M BH2 / (2t2)), automatically extracting effective efficiencies (e.g., F in range 2.958–3.613). This enables the system to predict mass growth of compact objects (PBHs, black holes) after formation in expanding universes, including near-critical cases.

Improvement 5: Excision-Like Regularization for Boundary Robustness

The AI can implement interior regularization (e.g., multiplying evolution equations by a window function inside the apparent horizon) to extend simulations in full periodic boxes without spurious boundary effects, while preserving agreement with symmetry-reduced runs. This allows the system to handle non-symmetric or multi-peak initial data reliably.

Improved AI System Capabilities

• Simulates 3D gravitational collapse and black-hole/PBH formation from arbitrary curvature perturbations in radiation-dominated or other barotropic-fluid cosmologies (e.g., kination, scalar-field collapse) with high efficiency and accuracy.

• Automatically balances temporal and spatial resolution using adaptive gauge scaling, reducing wall-clock time by up to two orders of magnitude for long-term runs.

• Provides validated, high-precision predictions of collapse thresholds, critical exponents, and post-formation mass accretion laws, directly usable for cosmological abundance calculations.

• Extends beyond spherical symmetry to study realistic, non-spherical, multi-peak initial conditions, enabling new insights into PBH formation from inflationary models.

• Offers a self-checking framework (via multiple gauges and boundary treatments) that ensures reliability in extreme multiscale regimes, applicable to other fields like neutron-star mergers or core-collapse supernovae.