Appendix B: Numerical Stability and Convergence

This appendix documents the numerical stability properties and convergence behavior of the simulations used throughout this work. The goal is not to optimize numerical performance, but to demonstrate that the reported phenomenology is robust under systematic refinement and is not a discretization artifact.

B.1 Discretization Scheme

All simulations are performed on a uniform Cartesian grid with spacing $\Delta x$ and fixed time step $\Delta t$. Spatial derivatives are computed using second-order finite differences. Temporal evolution is carried out with an explicit time-stepping scheme.

No adaptive mesh refinement or implicit solvers are employed. This choice is intentional, ensuring that all observed effects arise from the underlying dynamics rather than numerical sophistication.

B.2 Courant and Stability Constraints

Stability of the explicit scheme requires that $\Delta t$ satisfies a Courant-type condition determined by the characteristic velocity scale:

where $C < 1$ is a safety factor.

All production runs satisfy this constraint with a conservative margin. Reducing $\Delta t$ further does not alter qualitative behavior, only increasing computational cost.

B.3 Damping-Controlled Stability

The damping coefficient $\gamma$ plays a dual role:

• It stabilizes high-frequency velocity fluctuations.
• It controls the inertial lifetime of angular momentum.

Numerical instability is not observed for any $\gamma$ explored in the parameter scans. However, the physical regime transitions (orbital, collapsing, flyby) remain sharply distinct across the entire stable range.

B.4 Grid Convergence

Convergence tests were performed by repeating representative simulations at multiple spatial resolutions:

For fixed physical parameters $(\gamma, \mu, \Delta t)$, the following quantities were monitored:

• Mean absolute angular momentum $\langle |L| \rangle$
• Estimated orbit count $N_{\text{orbit}}$
• Radial oscillation index $\Delta_r$

All three quantities converge monotonically with increasing resolution. No regime transitions are induced by grid refinement.

B.5 Temporal Horizon Scaling

To test long-time stability, simulations were extended to multiple temporal horizons:

with identical initial conditions.

The classification outcome (orbital, collapse, flyby) remains unchanged under horizon extension. While instantaneous diagnostics (e.g. final separation) evolve with time, the existence or non-existence of sustained orbital motion does not.

This demonstrates that the closure condition is not a transient artifact.

B.6 Sensitivity to Initialization

Small perturbations were introduced in the initial density and velocity fields with amplitude:

Across multiple random seeds and spatial shifts, the system consistently returns the same phase classification for a given parameter set. This indicates structural stability rather than fine-tuned behavior.

B.7 Absence of Numerical Forcing

No artificial angular momentum injection is present in the scheme. In particular:

• There is no explicit rotational bias.
• Boundary conditions are symmetric.
• The screened Poisson solver conserves global neutrality.

Observed angular momentum storage therefore arises dynamically from the interaction of inertia, damping, and history-dependent evolution.

B.8 Summary

The numerical evidence supports the following conclusions:

• The simulations are stable across all explored parameters.
• Phase classification is invariant under grid and time refinement.
• Orbital behavior is not a discretization or finite-time artifact.
• The closure-based definition of gravity is numerically robust.

These results justify interpreting the observed phenomena as intrinsic to the dynamical framework rather than numerical pathology.