From Black Hole to Galaxy: Neural Operator: Framework for Accretion and Feedback Dynamics

TL;DR

Modeling SMBH–galaxy co-evolution requires resolving physics from milliparsecs to megaparsecs, which is computationally infeasible with full physics. The authors introduce a neural-operator-based subgrid black hole that learns the small-scale GRMHD dynamics and embeds it within a two-level multi-level solver to supply dynamic boundary conditions and fluxes for larger scales. The method achieves large speedups and stabilizes long-horizon rollouts, capturing intrinsic variability in accretion-driven feedback that standard closures miss. Training on GRMHD data and the use of radial shell embeddings, radial scaling baselines, and physics-informed losses yield a data-driven, transferable subgrid closure applicable to SMBHs and neutron stars in astrophysical simulations. This approach reframes subgrid modeling in computational astrophysics and promises scalable closures for cosmological simulations like FIRE and IllustrisTNG.

Abstract

Modeling how supermassive black holes co-evolve with their host galaxies is notoriously hard because the relevant physics spans nine orders of magnitude in scale-from milliparsecs to megaparsecs--making end-to-end first-principles simulation infeasible. To characterize the feedback from the small scales, existing methods employ a static subgrid scheme or one based on theoretical guesses, which usually struggle to capture the time variability and derive physically faithful results. Neural operators are a class of machine learning models that achieve significant speed-up in simulating complex dynamics. We introduce a neural-operator-based ''subgrid black hole'' that learns the small-scale local dynamics and embeds it within the direct multi-level simulations. Trained on small-domain (general relativistic) magnetohydrodynamic data, the model predicts the unresolved dynamics needed to supply boundary conditions and fluxes at coarser levels across timesteps, enabling stable long-horizon rollouts without hand-crafted closures. Thanks to the great speedup in fine-scale evolution, our approach for the first time captures intrinsic variability in accretion-driven feedback, allowing dynamic coupling between the central black hole and galaxy-scale gas. This work reframes subgrid modeling in computational astrophysics with scale separation and provides a scalable path toward data-driven closures for a broad class of systems with central accretors.

Figures (8)

• Figure 1: Illustration of our method in a two-level setting. The neural operator efficiently replaces the evolution on small scales (zoom in), enabling proper feedback on large scale (zoom out).
• Figure 2: MHD outputs. Simulated ('DNS'), CNN-predicted, and NO-predicted magnetic field component $B_z$, density $\rho$, internal energy $E_{int}$, and velocity component $v_z$ from the MHD simulation on a horizontal slice through the domain center (after 50 steps). Unphysical artifacts are marked. (2) Comparison of the mass density distributions on a 2D sphere, between analytically prescribed, directed simulated, and NO-rollout (after 10 steps). (3) The radial profiles of predicted quantities: mass density $\rho$, temperature $T$, and total mass flux $\dot{M}$ (Left to Right). Profile from simulation after 50 steps is also shown. The grey-shaded region shows the part that will be coupled to the direct simulation (see Appendix \ref{['apdx: insert back']}).
• Figure 3: GRMHD outputs. Simulated and predicted magnetic field components $(B_x,B_y,B_z)$, density $\rho$, internal energy $E_{int}$, and velocity components $(v_x,v_y,v_z)$ from the GRMHD simulations. Jet structure near the polar region and disk structure near the midplane is preserved.
• Figure 4: Ablation study on a central slice at step $t{=}50$. Columns (left$\to$right): $B_x$, $B_y$, $B_z$, mass density $\rho$, internal energy $e$, and velocity components $(v_x,v_y,v_z)$. Values are in the model’s normalized domain (see Methods). Rows: (a) Ground truth; (b) Ours (LocalNO + shell positional encoding + radial baseline/constraints + ROI and $H^1$ terms); (c) CNN baseline (pure convolutional surrogate); (d) Plain $L^2$ (no component weights/ROI/$H^1$); (e) Positional Embedding only (Fourier features, no shells); (f) No Shell Encoding; (g) No Radial Constraint/Baseline. Colored boxes mark typical failure modes observed in ablations: yellow: over-smoothing/texture blur; orange: spurious mass/energy (unphysical deficits or excess). The full model (b) most closely matches the ground truth, preserving small-scale magnetic structure and inner-torus morphology while avoiding excessive smoothing and mass artifacts.
• Figure 5: Ablation study on a central slice at step $t{=}100$. Columns (left$\to$right): $B_x$, $B_y$, $B_z$, mass density $\rho$, internal energy $e$, and velocity components $(v_x,v_y,v_z)$. Values are in the model’s normalized domain (see Methods). Rows: (a) Ours (LocalNO + shell positional encoding + radial baseline/constraints + ROI and $H^1$ terms); (b) CNN baseline (pure convolutional surrogate); (c) Plain $L^2$ (no component weights/ROI/$H^1$); (d) Positional Embedding only (Fourier features, no shells); (e) No Shell Encoding; (f) No Radial Constraint/Baseline. Colored boxes mark typical failure modes observed in ablations: yellow: over-smoothing/texture blur; green: anisotropic diffusion of the torus; orange: spurious mass/energy (unphysical deficits or excess). The full model (a) most closely maintains proper structure, preserving small-scale magnetic structure and inner-torus morphology while avoiding excessive smoothing and mass artifacts.