Novel Conservative Methods for Adaptive Force Softening in Collisionless and Multi-Species N-Body Simulations

TL;DR

The paper develops a fully general, energy- and momentum-conserving framework for adaptive gravitational softening in collisionless N-body simulations, accommodating a broad class of softening rules. It derives generalized equations of motion that include grad-ε corrections and multi-species interaction terms, and introduces several new schemes, notably a tidal-tensor–based softenings, along with consistent timestep criteria and maximum/minimum softening bounds. Through analytical derivations and cosmological tests, the authors demonstrate that tidal softenings preserve substructure better than traditional adaptive schemes while maintaining energy conservation, and they provide practical guidance on symmetry choices, kernel functions, and computational costs. The work culminates with a public implementation in the GIZMO code, offering a flexible toolkit for accurate, conserved adaptive gravity in multi-physics simulations with strong dynamic range.

Abstract

Modeling self-gravity of collisionless fluids (e.g. ensembles of dark matter, stars, black holes, dust, planetary bodies) in simulations is challenging and requires some force softening. It is often desirable to allow softenings to evolve adaptively, in any high-dynamic range simulation, but this poses unique challenges of consistency, conservation, and accuracy, especially in multi-physics simulations where species with different softening laws may interact. We therefore derive a generalized form of the energy-and-momentum conserving gravitational equations of motion, applicable to arbitrary rules used to determine the force softening, together with consistent associated timestep criteria, interaction terms between species with different softening laws, and arbitrary maximum/minimum softenings. We also derive new methods to maintain better accuracy and conservation when symmetrizing forces between particles. We review and extend previously-discussed adaptive softening schemes based on the local neighbor particle density, and present several new schemes for scaling the softening with properties of the gravitational field, i.e. the potential or acceleration or tidal tensor. We show that the tidal softening scheme not only represents a physically-motivated, translation and Galilean invariant and equivalence-principle respecting (and therefore conservative) method, but imposes negligible timestep or other computational penalties, ensures that pairwise two-body scattering is small compared to smooth background forces, and can resolve outstanding challenges in properly capturing tidal disruption of substructures (minimizing artificial destruction) while also avoiding excessive N-body heating. We make all of this public in the GIZMO code.

Figures (13)

• Figure 1: Top: Examples of different softening "rules" from § \ref{['sec:ex']}, in a collisionless Plummer sphere ($\rho=3\,M/(4\pi\,a_{\rm Plummer}^{3}\,[1+r^{2}/a_{\rm Plummer}^{2}]^{5/2})$) realized with $\sim 10^{5}$ equal-mass ($m$) $N$-body particles and evolved to $10$ dynamical times ($t_{\rm dyn} = (G\,M/a_{\rm Plummer}^{3})^{-1/2})$. We plot the median (thick line) and $90\%$ inclusion interval of the calculated force softening $\varepsilon$ (defined for the cubic spline softening) at each radius $r$. The normalization of each rule is somewhat arbitrary, but all of the adaptive rules produce a broadly similar behavior of larger $\varepsilon$ at $r\gg a_{\rm Plummer}$, with $\varepsilon$ (by construction) similar to the constant-softening (§ \ref{['sec:ex.fixed']}) around the effective radius. For $\tilde{\rho}$ (§ \ref{['sec:ex.rho']}) and $n_{\rm NGB}$ (§ \ref{['sec:ex.n']}), identical here because particle masses are equal, we adopt $N_{\rm ngb}=32$; for the ${\bf T}$-based model (§ \ref{['sec:ex.T']}), $\xi=3.5$. But for the $\hat{\Phi}$ (§ \ref{['sec:ex.phi']}) and $\hat{\bf a}$ (§ \ref{['sec:ex.a']}) models we must adopt somewhat more "ad hoc" normalizations to obtain reasonable $\varepsilon$, with $\xi \sim 3000$ and $\xi \sim 0.01$ respectively. Thin dotted lines show the analytically expected $\varepsilon$ for a perfect realization of the system. Bottom: Same, for a hernquist:profile profile ($\rho = M\,a_{\rm Hernquist}/(2\pi\,r\,[r + a_{\rm Hernquist}]^{3})$). The qualitative differences between models are similar to the Plummer case. In all cases the system is stable in time and conserves energy and momentum to integration accuracy.
• Figure 2: Errors in the reconstruction of the potential ( top) and acceleration ( bottom) of the Plummer sphere test from Fig. \ref{['fig:plummer.h']}, where $\Phi$ and ${\bf a}$ are the potential and acceleration of each $N$-body particle and $\Phi_{\rm model}$, ${\bf a}_{\rm model}$ are the analytic Plummer values. We confirm the well-known result that force softening actually has a quite weak effect on these errors. For ${\bf a}$, the error is dominated by discreteness noise determined by the number of particles: we can accurately approximate the mean $|{\bf a}-{\bf a}_{\rm model}| / |{\bf a}_{\rm model}|$ as $\sim 2/\sqrt{N(<r)}$ in terms of the number of enclosed particles $N(<r)$. The result is very similar if we compare errors at $t=0$ (when the particle configuration is identical).
• Figure 3: Numerical timesteps in the Plummer test from Fig. \ref{['fig:plummer.h']}. Note the code uses a discretized power-of-two timestep hierarchy (the reason for discrete "levels") and adopts a standard universal timestep criterion for integrating gravity alongside the adaptive-softening-specific additional criterion in Eq. \ref{['eqn:courant']} (with $C=0.25$). The fixed-$\varepsilon$ and potential/acceleration/tidal criteria have timesteps primarily determined by the standard gravity criterion (Eq. \ref{['eqn:courant']} is already met and imposes no significant additional burden). The neighbor-based ($\tilde{\rho}$ and $n_{\rm NGB}$) models require much shorter timesteps to accurately integrate through collisionless particle-particle crossings owing to the sensitivity of $\varepsilon$ to the strictly-local particle distribution (see § \ref{['sec:timesteps']}).
• Figure 4: Toy illustration of the reason for the timestep difference in Fig. \ref{['fig:plummer.dt']}, following discussion in § \ref{['sec:timesteps']}. Consider two streams or sheets ($1$ & $2$) of collisionless particles, intersecting at an arbitrary angle in an external potential (negligible self-gravity), with local densities $\tilde{\rho}_{1}$, $\tilde{\rho}_{2}$ (estimated as § \ref{['sec:ex.rho']}) and velocities ${\bf v}_{1}$, ${\bf v}_{2}$, respectively ( top). There are no local interactions, so this crossing has no effect on the dynamics, and in the fixed-$\varepsilon$ or gravity ($\hat{\Phi}$/$\hat{\bf a}$/${\bf T}$)-based schemes, negligible effect on the timesteps $\Delta t \approx \Delta t_{\rm global}$, which remain constant at their "global" value ($\sim C\,L_{\rm global}/|{\bf v}|$, where $L_{\rm global}$ is some global potential scale-length) set by the gravity integrator. But in the neighbor ($\tilde{\rho}/n_{\rm NGB}$)-based schemes ( bottom), the change in $\tilde{\rho}$/$n_{\rm NGB}$ (equivalently, the local particle/mass arrangement) as a test particle moving with stream $1$ intersects (or moves out of) stream $2$ leads by definition to a change in $\varepsilon$ over a scale length $\Delta x \sim \varepsilon$, or timescale $\delta t \sim \varepsilon/|{\bf v}| \sim 1/| \tilde{\nabla} \cdot {\bf v}|$ (where $\tilde{\nabla}$ is an "SPH-like" divergence estimator based on the in-kernel neighbors). So maintaining energy conservation in the $\tilde{\rho}/n_{\rm NGB}$ schemes necessarily requires a timestep $\Delta t \lesssim C\,\delta t$, a factor of $\sim (\epsilon/L_{\rm global})$ smaller than the orbital/gravitational dynamics timestep $\Delta t_{\rm global}$.
• Figure 5: Behavior of the force in a close particle encounter following different symmetrization rules for $\tilde{\phi}_{bc}=\tilde{\phi}_{cb}$ discussed in § \ref{['sec:symm']}, as a function of the separation between the particles in units of the compact support radius $\tilde{\varepsilon}$ of the (non-Keplerian) softening kernel. We consider an encounter of two particles $a$, $b$ of masses $m_{a}$, $m_{b}$ with unequal softenings $\varepsilon_{a}\ne \varepsilon_{b}$ (ratio of $\varepsilon_{a}$ and $\varepsilon_{b}$ taken to be $=5$ or $=10$ at top and bottom), and plot the magnitude of the total force along the line connecting their centers-of-mass. We compare the un-softened limit (treating both as point masses) and the softened result (all adopt a cubic spline $K$, see § \ref{['sec:kernels']}) for the "average the force" rule ($\tilde{\phi} = (1/2)\,(\phi^{0}(\varepsilon_{a}) + \phi^{0}(\varepsilon_{b}))$, Eq. \ref{['eqn:phisymm.avg']}), "average the softening" ($\tilde{\phi}=\phi^{0}(\tilde{\varepsilon})$ with $\tilde{\varepsilon}=(\varepsilon_{a}+\varepsilon_{b})/2$, Eq. \ref{['eqn:hsymm.avg']}) and "maximum softening" ($\tilde{\varepsilon}={\rm MAX}(\varepsilon_{a},\,\varepsilon_{b})$) rules. We also compare the exact result if we treat each particle as an extended mass distribution with size parameter $\varepsilon_{a}$ obeying $\nabla^{2} \phi_{a}(\varepsilon_{a}) = 4\pi\,G\,\rho$ (§ \ref{['sec:kernels']}). By definition all models become Keplerian at large distance. In close encounters, the "${\rm MAX}$" rule gives nearly identical results to the exact extended-mass calculation (and becomes exact for highly-unequal softenings), while the "average $\tilde{\varepsilon}$" rule deviates systematically and the "average the force" rule remains nearly-Keplerian (un-softened) until both particles are well inside the smaller of the two force softenings.