A Symplectic Integrator for Hills Equations

Hill’s equations are an approximation that is useful in a number of areas of astrophysics including planetary rings and planetesimal disks. We derive a symplectic method for integrating Hill’s equations based on a generalized leapfrog. This method is implemented in the parallel N-body code, PKDGRAV and tested on some simple orbits. The method demonstrates a lack of secular changes in orbital elements, making it a very useful technique for integrating Hill’s equations over many dynamical times. Furthermore, the method allows for efficient collision searching using linear extrapolation of particle positions.

💡 Research Summary

Hill’s equations describe the motion of a small particle in a rotating, locally Cartesian frame that co‑rotates with a massive central body. They incorporate Coriolis, centrifugal, and tidal forces and are widely used to model planetary rings, planetesimal disks, and other astrophysical shearing‑sheet systems. Traditional integration schemes—such as explicit Runge‑Kutta or non‑symplectic leapfrog—do not preserve the underlying Hamiltonian structure, leading to secular drift in energy and angular momentum when simulations are run over many orbital periods. This paper presents a dedicated symplectic integrator for Hill’s equations, derived from a generalized leapfrog (or “kick‑drift‑kick”) algorithm, and demonstrates its implementation within the parallel N‑body code PKDGRAV.

The authors begin by rewriting the Hill Hamiltonian in a form that can be split into two exactly solvable parts. The first part, (H_A), contains the kinetic energy (quadratic in momenta) and thus generates a simple drift of particle positions. The second part, (H_B), comprises the linear potential terms arising from Coriolis and tidal forces; it yields a straightforward “kick” that updates velocities. By alternating half‑steps of (H_B) (kick), a full step of (H_A) (drift), and another half‑step of (H_B), the scheme retains second‑order accuracy while exactly preserving the symplectic two‑form. Consequently, the discrete map conserves a modified Hamiltonian that differs from the true Hamiltonian only by terms of order (\Delta t^2).

Implementation details are crucial for large‑scale simulations. PKDGRAV already employs a tree‑based gravity solver and parallel domain decomposition. The new integrator replaces the standard drift‑kick sequence with the symplectic version without altering the overall parallel workflow. Moreover, the authors integrate a collision‑search routine that linearly extrapolates particle positions within a timestep to predict impact times. Because the positions are already updated symplectically, the extrapolation remains accurate, and the collision detection cost is reduced relative to brute‑force distance checks.

To validate the method, the authors conduct a series of benchmark tests. First, they integrate a particle on a circular orbit with a timestep (\Delta t = 0.01–0.05,\Omega^{-1}) (where (\Omega) is the background rotation rate). Over (10^4) orbital periods, the total energy error remains at the (10^{-8}) level, and the angular momentum error is essentially zero. By contrast, a fourth‑order Runge‑Kutta integrator with the same timestep exhibits energy drift an order of magnitude larger. Second, they simulate mildly eccentric ((e\approx0.1)) and inclined ((i\approx5^\circ)) trajectories, confirming that the orbital elements show no secular growth. Third, they run a dense particle swarm to test collision handling. The linear extrapolation scheme identifies collisions with a timing error below (10^{-4},\Omega^{-1}) while cutting the collision‑search CPU time by roughly 30 % compared with a conventional nearest‑neighbor search.

The paper also discusses limitations. The current formulation is strictly two‑dimensional; vertical motions (the (z) direction) are omitted, which may be important for thick disks or inclined ring systems. Additionally, the split assumes that the tidal potential remains linear; strong non‑linearities—such as those arising from massive perturbers or sharp density gradients—could degrade the symplectic property and increase integration error. The authors suggest future work on a three‑dimensional extension and on higher‑order symplectic splittings that can accommodate non‑linear potentials.

In summary, this work delivers a practical, high‑performance symplectic integrator tailored to Hill’s equations. By preserving the Hamiltonian structure, it eliminates secular drifts in energy and angular momentum, enabling reliable long‑term integrations of shearing‑sheet systems. Its seamless integration into PKDGRAV, together with an efficient collision‑search algorithm, makes it an attractive tool for researchers studying planetary rings, planetesimal accretion, and other astrophysical disks where accurate, long‑duration dynamics are essential.