Preserving Angular Momentum in Symplectic Numerical Integration

The discrete Noether theorem provides the foundational answer: when a discrete Lagrangian inherits the rotational symmetry of its continuous counterpart, the associated momentum map is exactly conserved at each timestep. This elegant result from variational integrators theory, formalized by Marsden and West (2001), shows that angular momentum preservation is not an add-on but emerges naturally from proper discretization of the variational principle. However, this theoretical guarantee comes with practical constraints—the Ge-Marsden impossibility theorem (1988) proves that no integrator can simultaneously preserve symplecticity, energy, and all momentum maps for general non-integrable systems, forcing practitioners to make fundamental design choices.

The practical landscape divides into three approaches: variational integrators that preserve momentum through symmetry, constrained algorithms like RATTLE that enforce conservation via Lagrange multipliers, and energy-momentum methods that sacrifice symplecticity for exact physical invariants. Each serves distinct application domains, from celestial mechanics requiring billion-year orbital stability to computational solid mechanics demanding exact energy balance.

Mathematical foundations through discrete variational principles

The continuous Noether theorem establishes that rotational invariance of a Lagrangian implies angular momentum conservation. For numerical methods, the critical insight is to discretize Hamilton's principle directly rather than the resulting equations of motion. Given a continuous action S[q] = ∫L(q,q̇)dt, one constructs a discrete action Sₐ = ΣₖLₐ(qₖ,qₖ₊₁) where Lₐ approximates the action integral over one timestep.

The discrete Euler-Lagrange equations emerge from stationarity: D₂Lₐ(qₖ₋₁,qₖ) + D₁Lₐ(qₖ,qₖ₊₁) = 0. When the discrete Lagrangian satisfies G-invariance—meaning Lₐ(g·q₀, g·q₁) = Lₐ(q₀,q₁) for all group elements g in the symmetry group G—the discrete momentum map is exactly preserved between steps. For angular momentum specifically, if SO(3) acts on configuration space and Lₐ respects this action, then J(qₖ,qₖ₊₁) = qₖ × pₖ remains constant.

The momentum map J: M → g* (where g* is the dual of the Lie algebra) captures the conserved quantity. For the rotation group acting on ℝ³, the three components are precisely the angular momenta. The discrete momentum map definitions J⁺ and J⁻, evaluated at consecutive configurations, satisfy J⁺(qₖ₋₁,qₖ) = J⁻(qₖ,qₖ₊₁)—exact conservation without numerical drift.

This framework requires regularity: the discrete Legendre transforms must be local diffeomorphisms. The most common construction uses the midpoint discrete Lagrangian Lₐ(q₀,q₁) = h·L((q₀+q₁)/2, (q₁-q₀)/h), which produces the implicit midpoint method and preserves both linear and angular momentum when the continuous Lagrangian has corresponding symmetries.

Projection methods trade symplecticity for exact constraints

Projection methods offer the most direct approach to enforcing angular momentum conservation: take any integration step, then project the result onto the constraint manifold J⁻¹(μ₀) where μ₀ is the initial angular momentum. The mathematical formulation solves a constrained optimization problem—minimize ‖y - ỹₙ₊₁‖² subject to J(y) = μ₀, where ỹₙ₊₁ comes from the base integrator.

Orthogonal projection uses Newton iteration: y^(k+1) = y^(k) - ∇J·λ^(k) where λ solves J(y^(k) - ∇J·λ) = μ₀. This typically converges in 1-3 iterations. Oblique projection instead moves along a chosen direction v, solving for the scalar λ in yₙ₊₁ = ỹₙ₊₁ + λ·v.

The critical limitation is that generic projection destroys symplecticity. The projected method exactly satisfies the constraint but loses the near-energy conservation property that makes symplectic methods valuable for long-time integration. The symplectic 2-form ω is typically not tangent to the constraint manifold J⁻¹(μ), so projecting off the manifold moves the solution to a region where the original symplectic structure no longer applies.

Symmetric projection methods preserve time-reversibility by projecting both the initial and final points, then combining them symmetrically. This maintains the geometric properties of symmetric integrators while enforcing constraints, though symplecticity remains sacrificed.

When projection is applied along the Hamiltonian vector field of the constraint function—rather than orthogonally—symplecticity can sometimes be preserved, forming the basis of the Simo-Tarnow-Wong energy-momentum methods discussed below.

RATTLE maintains symplecticity on constraint manifolds

The RATTLE algorithm (Andersen, 1983), extending SHAKE (Ryckaert et al., 1977), enforces constraints using Lagrange multipliers while maintaining symplecticity on the constraint manifold. Originally developed for molecular bond constraints, the framework extends naturally to angular momentum.

The algorithm operates in three stages. First, a half-step velocity update: v_{n+1/2} = vₙ + (Δt/2)M⁻¹(F(qₙ) - ∇g(qₙ)ᵀλₙ). Second, position update: qₙ₊₁ = qₙ + Δt·v_{n+1/2}, solving g(qₙ₊₁) = 0 for the Lagrange multiplier λₙ. Third, full velocity update: vₙ₊₁ = v_{n+1/2} + (Δt/2)M⁻¹(F(qₙ₊₁) - ∇g(qₙ₊₁)ᵀμₙ₊₁), solving ∇g·vₙ₊₁ = 0 for μₙ₊₁.

For angular momentum constraints specifically, g(q,p) = q × p - L₀ = 0, with the hidden constraint ġ = q̇ × p + q × ṗ = 0 enforced at the velocity level. The Lagrange multipliers λ and μ represent constraint forces perpendicular to the constraint manifold.

The key theorem (Leimkuhler-Skeel, 1994) states that RATTLE is symplectic on the constraint manifold—the restricted 2-form ω|M is preserved by the discrete flow. This combines exact constraint satisfaction with the long-term stability properties of symplectic integration, making RATTLE the gold standard for constrained molecular dynamics.

Computational cost involves iterative solution of constraint equations. Standard SHAKE converges linearly, typically requiring 5-10 iterations for tolerances around 10⁻⁷. Improved variants offer significant speedups: P-SHAKE (preconditioned) achieves 3× faster convergence than matrix-inversion SHAKE, while LINCS provides 3-4× speedup over SHAKE through a non-iterative approach.

Lie group integrators preserve rotational structure intrinsically

For systems evolving on Lie groups such as rigid body dynamics on SO(3), Lie group integrators exploit the mathematical structure to guarantee solutions remain on the group manifold while preserving associated invariants. Rather than integrating rotation matrices R ∈ SO(3) directly—which would accumulate numerical drift violating orthogonality—these methods integrate in the Lie algebra so(3) and map back via the exponential.

The fundamental update takes the form Rₙ₊₁ = Rₙ·exp(Ωₙ) where Ωₙ ∈ so(3) is computed by numerical integration in the linear Lie algebra space. For SO(3), the Rodrigues formula provides a closed-form exponential at cost of approximately 20-30 floating-point operations, making this approach computationally practical.

RKMK methods (Runge-Kutta-Munthe-Kaas) generalize classical Runge-Kutta to Lie groups. For an ODE dg/dt = g·ξ(g,t) on group G, the s-stage RKMK algorithm computes stages kᵢ = h·dexp⁻¹_{uᵢ}(ξ(gₙ·exp(uᵢ), tₙ + cᵢh)) using the inverse of the exponential derivative, then updates gₙ₊₁ = gₙ·exp(Σᵢbᵢkᵢ). The dexp⁻¹ operation for SO(3) expands as dexp⁻¹_Ω(ξ) = ξ - ½[Ω,ξ] + (1/12)[Ω,[Ω,ξ]] + O(‖Ω‖⁴).

Crouch-Grossman methods avoid the dexp⁻¹ computation entirely by using group operations directly, trading some achievable order for implementation simplicity. Magnus integrators expand the solution Y(t) = exp(Ω(t))Y(0) using the Magnus series Ω = Ω₁ + Ω₂ + ···, particularly effective for time-dependent systems.

For angular momentum preservation, Lie group variational integrators (Lee, Leok, McClamroch 2007; Bou-Rabee and Marsden 2009) combine the Lie group structure with discrete variational principles. The discrete Lagrangian on SO(3) takes form Lₐ(Rₖ,Rₖ₊₁) = (1/2h)tr((Rₖ₊₁ - Rₖ)ᵀJ(Rₖ₊₁ - Rₖ)), yielding discrete Euler-Poincaré equations that exactly preserve body angular momentum Π = IΩ via the discrete Noether theorem.

Energy-momentum methods sacrifice symplecticity for exact conservation

The Simo-Tarnow-Wong energy-momentum methods, developed for computational mechanics in 1991-1995, provide an alternative philosophy: preserve both energy and momentum exactly at the cost of symplecticity. These methods emerged from the observation that widely used implicit schemes (the Newmark family) fail to conserve angular momentum for nonlinear Hamiltonian systems.

The mathematical foundation uses discrete derivative operations defining algorithmic internal forces that guarantee conservation. For nonlinear elastodynamics with hyperelastic materials (Simo and Gonzalez), stress tensors are modified at equilibrium points as solutions to local convex projections, maintaining second-order accuracy through a simple two-step solution scheme.

What is sacrificed: the symplectic form. Phase space volume may not be exactly conserved, long-term qualitative dynamics differ from symplectic methods, and KAM-theory stability guarantees do not apply. However, for engineering applications where energy balance must be monitored and verified, this trade-off is often acceptable.

Average vector field (AVF) methods replace the vector field f(y) with its path average: yₙ₊₁ = yₙ + h∫₀¹f(yₙ + τ(yₙ₊₁ - yₙ))dτ. This construction achieves exact energy preservation for Hamiltonian systems, second-order accuracy, and time-symmetry—but not symplecticity. The improved AVF (IAVF) extends this to fourth order.

Discrete gradient methods construct custom update rules preserving first integrals. The discrete gradient ∇̄H satisfies H(yₙ₊₁) - H(yₙ) = ∇̄H(yₙ,yₙ₊₁)·(yₙ₊₁ - yₙ) and ∇̄H(y,y) = ∇H(y), guaranteeing exact energy conservation regardless of step size. Angular momentum follows when spatial discretization respects rotational symmetry.

Hamiltonian Boundary Value Methods (HBVMs) use discrete line integrals to achieve arbitrary-order energy conservation. For polynomial Hamiltonians of degree ν, HBVM(k,s) with k Gauss-Legendre nodes conserves energy exactly when 2k ≥ νs, providing high-order schemes that can simultaneously preserve multiple independent invariants.

Computational trade-offs shape method selection

Performance comparisons reveal significant differences. For the outer solar system integrated over 1 billion years, standard Wisdom-Holman symplectic integrators require over one year of computation to achieve 10⁻¹⁰ relative energy error, while SABA(10,6,4) higher-order methods achieve 10⁻¹⁴ error in approximately 4 hours—a dramatic efficiency gain from order elevation.

Force evaluations dominate computational cost. Second-order leapfrog requires 1 evaluation per step; fourth-order composition methods need 3-6; eighth-order methods require 15-17. For moderate precision, the Blanes-Moan 6-stage fourth-order method offers an excellent error-to-cost ratio. Higher-order methods win decisively for high-precision requirements.

| Method | Order | Energy | Momentum | Symplectic | Relative Cost | |--------|-------|--------|----------|------------|---------------| | Leapfrog/Verlet | 2 | Oscillating | Exact | Yes | Low | | RATTLE | 2 | Oscillating | Exact (constraint) | On manifold | Medium | | Variational (midpoint) | 2 | Oscillating | Exact | Yes | Medium | | Gauss-Legendre RK | 2s | Quadratic only | Quadratic only | Yes | High | | Simo-Tarnow | 2 | Exact | Exact | No | Medium | | Discrete gradient | 2+ | Exact | Exact* | No | High | | HBVM(k,s) | 2k | Exact | Partial | No | High |

Constraint solver convergence presents implementation challenges. SHAKE converges linearly; angle constraints converge twice as slowly as bond constraints. Ring topologies can cause non-convergence requiring matrix-inversion approaches. Maximum iteration limits of 100-1000 are typical before failure flags.

For Lie group methods, the cost of exponential maps varies significantly. SO(3) exponentials via Rodrigues formula cost ~20-30 flops (closed form). General n×n matrices require O(n³) via Padé approximation. Commutator-free methods minimize the number of exponentials per stage—optimal configurations use 5 commutators for RKF45-type methods and 21 for seventh-order Butcher methods.

Velocity-dependent forces require modified approaches

Standard symplectic integrators assume separable Hamiltonians H(q,p) = T(p) + V(q). Velocity-dependent forces like the Lorentz force F = q(E + v × B) violate this assumption, creating the non-separable electromagnetic Hamiltonian H = (1/2)(p - A(x))² + φ(x).

The Boris algorithm (1970) remains the standard for charged particle integration despite not being symplectic. Its effectiveness stems from being volume-preserving (proved by Qin et al. 2013) and time-reversible, which bounds global energy error over long times. The algorithm applies half-step electric kicks around a magnetic rotation: v⁻ = vⁿ + (qΔt/2m)E, then v⁺ via rotation, then vⁿ⁺¹ = v⁺ + (qΔt/2m)E.

Variants address different requirements. Boris-SDC uses spectral deferred corrections for high-order accuracy. The Vay integrator handles relativistic particles with large E×B drift better. Hyper Boris integrators extend to arbitrary high orders using Chebyshev polynomials.

K-symplectic methods exploit a non-canonical symplectic structure natural to charged particles: Ω = d(v + A) ∧ dx with H = (1/2)v² + φ. This enables explicit symplectic methods via 4-way Hamiltonian splitting, each part producing an analytically solvable flow.

For Weber electrodynamics, where F = (q₁q₂/4πε₀r²)[1 - ṙ²/2c² + rr̈/c²]r̂ includes acceleration dependence, numerical challenges multiply. The force depends on r̈, creating implicit equations requiring iterative solution. However, Weber's force is always central (directed along r̂), automatically conserving particle angular momentum—a significant theoretical advantage over Lorentz electrodynamics where field angular momentum must be included.

Integration strategies for Weber forces include implicit ODE methods solving F(t,y,y',y'') = 0 at each step, extrapolation methods estimating r̈ from previous positions, and energy-based weak formulations exploiting Weber's Lagrangian structure.

Software implementations provide practical starting points

GeometricIntegrators.jl (Julia) offers the most comprehensive collection of structure-preserving methods including symplectic Runge-Kutta, Lie group integrators, and variational integrators with active development and thorough documentation. REBOUND (C with Python bindings) specializes in N-body dynamics, implementing WHFast, IAS15, and high-order SABA-family integrators with SIMD optimization achieving up to 4.7× speedup.

For molecular dynamics, GROMACS provides LINCS (3-4× faster than SHAKE) and P-LINCS for parallel systems. LAMMPS offers configurable SHAKE/RATTLE with flexible tolerances. DifferentialEquations.jl includes ManifoldProjection callbacks for constraint enforcement.

Implementation best practices include timestep selection starting at τ = T_min/20 (T_min being the shortest timescale), monitoring bounded energy oscillations to verify stability, and running convergence studies to verify expected order. Shadow Hamiltonian monitoring—available to 24th order via Skeel and Hardy formulas—provides 100,000× more sensitivity to integration quality than direct energy measurement.

Conclusions and method selection guidance

The fundamental tension between symplecticity and exact conservation, formalized by the Ge-Marsden theorem, means method selection depends critically on application requirements. Symplectic-momentum methods (variational integrators) excel for long-time astronomical simulations where bounded energy error and preserved phase space structure matter more than exact conservation. Energy-momentum methods (Simo-Tarnow, discrete gradient) suit engineering applications requiring exact energy balance verification.

For angular momentum specifically: symplectic methods preserve it exactly when the discrete Lagrangian inherits rotational symmetry—no augmentation needed. RATTLE/SHAKE enforce it as a constraint while maintaining symplecticity on the constraint manifold. Lie group integrators preserve rotational structure intrinsically for systems on SO(3).

The most significant recent advances come from recognizing that different force types require different symplectic structures. K-symplectic methods for electromagnetic systems, extended phase space methods for general non-separable Hamiltonians, and specialized integrators for shearing-sheet dynamics in astrophysical disks each exploit problem-specific structure for dramatic efficiency gains. The choice between volume-preserving methods (explicit, robust) and strictly symplectic methods (stronger theoretical guarantees) depends on whether phase space structure preservation or computational efficiency takes priority.