4. Numerical Construction and Simulation Protocol

4.1 Discretization and Grid Geometry

All numerical experiments are performed on a two-dimensional uniform Cartesian grid. Space is discretized as

with equal spacing $\Delta x$ in both directions.

The fields $\rho(\mathbf{x},t)$, $\mathbf{v}(\mathbf{x},t)$, and $\Phi(\mathbf{x},t)$ are represented as grid-based arrays evaluated at cell centers. Spatial derivatives are computed using second-order centered finite differences, ensuring consistency between gradient and divergence operators.

Time evolution is performed with an explicit second-order Runge–Kutta (Heun) scheme. The time step $\Delta t$ is held constant throughout each simulation. No adaptive stepping is employed, allowing direct comparison across runs and parameter scans.

4.2 Boundary Conditions and Symmetry Breaking

Unless otherwise stated, simulations employ periodic boundary conditions on all fields:

for each spatial direction $\alpha$.

Periodic boundaries eliminate external forces and artificial confinement, ensuring that all observed structure arises internally from the dynamics. Importantly, no central potential, reflecting walls, or imposed symmetries are introduced.

Symmetry breaking is achieved solely through initial conditions. Small spatial offsets between density concentrations and optional stochastic perturbations are sufficient to prevent trivial head-on collapse and to probe the full dynamical phase space.

4.3 Initialization, Noise, and Reproducibility

Initial conditions consist of two localized density concentrations embedded in an otherwise uniform background. The initial velocity field is set to zero everywhere:

To test robustness and avoid fine-tuned trajectories, two controlled perturbations are optionally applied:

• A discrete spatial shift of one or more grid cells applied asymmetrically to one density component.
• Additive Gaussian noise of small amplitude $\epsilon \ll 1$ applied to the initial density field.

All stochastic elements are governed by an explicit random seed. Each simulation records the seed value, ensuring exact reproducibility of every run. Repeated simulations with identical parameters but different seeds are used to assess phase robustness rather than trajectory-level coincidence.

4.4 Control Parameters ($\gamma$, $\mu$, $\Delta t$)

The model is governed by a small set of dimensionless control parameters:

• $\gamma$ — the damping coefficient, controlling the transition between overdamped, inertial, and weakly damped regimes.
• $\mu$ — the screening parameter in the Poisson equation, setting the effective interaction range.
• $\Delta t$ — the fixed numerical time step.

Parameter scans are performed by varying $\gamma$ while holding all other parameters fixed, isolating its role in enabling or suppressing orbital dynamics. Horizon tests are conducted by extending the total integration time while keeping $\Delta t$ constant, allowing verification of long-term stability versus transient behavior.

No parameters are tuned to enforce orbital motion. The appearance, disappearance, or re-emergence of bound trajectories is treated as an empirical outcome of the numerical evolution.

Together, these numerical choices ensure that observed gravitational behavior is not a numerical artifact, but a genuine dynamical phase supported by the underlying equations.