Democratic Heliocentric Coordinates

• Democratic heliocentric coordinates (DHC) are a canonical system in symplectic N-body integrations that treat each planet equivalently in a heliocentric reference frame.
• DHC splits the Hamiltonian into Keplerian, interplanetary, and solar jump components, allowing operator-splitting methods but introducing artificial precession under typical timestep settings.
• Numerical experiments reveal that reducing the integration timestep restores physical instability rates, while Jacobi-coordinate splits provide enhanced stability for high-eccentricity planetary orbits.

Democratic heliocentric coordinates (DHC) are a canonical coordinate system employed in symplectic N-body integrations, most notably within Wisdom–Holman (WH) integrators, to model the long-term evolution of planetary systems. DHC are constructed to treat all planets equivalently with respect to the heliocentric frame, and enable operator splitting of the Hamiltonian into analytically tractable Keplerian, interplanetary, and so-called “jump” (solar) terms. Recent numerical experiments demonstrate that DHC introduce an eccentricity-dependent artificial precession in planetary orbits, significantly affecting the rates of secular instabilities, such as those involving Mercury in the Solar System, unless the integration timestep is reduced to values much smaller than are typically used (Rein et al., 12 Jan 2026).

1. Canonical Variables and DHC Transformations

The DHC construction proceeds from standard barycentric N-body variables—positions $r_i$, momenta $p_i$, masses $m_0$ (Sun), $m_1,\ldots,m_N$ (planets), and total mass $M = \sum_{i=0}^N m_i$. The canonical change of variables consists of:

• Center-of-mass coordinate and momentum:
  • $Q_0 = \frac{1}{M} \sum_{i=0}^N m_i r_i$
  • $P_0 = \sum_{i=0}^N p_i$
• Heliocentric “relative” coordinates and their conjugate momenta for $i=1,\ldots,N$:
  • $Q_i = r_i - r_0$
  • $P_i = p_i - \frac{m_i}{M} \sum_{j=0}^N p_j$

These transformed variables are canonical, i.e., $p_i$0, and their inverses express barycentric positions and momenta in terms of $p_i$1, and the center of mass variables. The total momentum $p_i$2 is conserved, so the barycentric motion trivially decouples from the internal planetary evolution (Rein et al., 12 Jan 2026).

2. Wisdom–Holman Hamiltonian Splitting in DHC

In DHC, the full N-body Hamiltonian is split exactly into three terms:

• $p_i$3: Kepler Hamiltonian for each planet’s heliocentric motion,

$p_i$4

• $p_i$5: Interplanetary potential,

$p_i$6

• $p_i$7: Solar “jump” term,

$p_i$8

The standard second-order WH map advances the system over timestep $p_i$9 via a sequence ("kick–kick–drift–kick–kick"):

1. Kick: update $m_0$0 by $m_0$1,
2. Kick: update $m_0$2 by $m_0$3,
3. Drift: exactly solve each heliocentric two-body problem for $m_0$4,
4. Kick: repeat step 2,
5. Kick: repeat step 1.

Each component of the Hamiltonian is solved exactly within its substep; only their noncommutativity induces integration errors (Rein et al., 12 Jan 2026).

3. Eccentricity-Dependent Artificial Precession in DHC

Due to the Baker–Campbell–Hausdorff expansion, second-order WH splitting in DHC generates leading errors that couple $m_0$5, $m_0$6, and notably $m_0$7. This manifests as an $m_0$8 secular term in each heliocentric Kepler subsystem, generating an artificial apsidal precession. Numerical results show:

• The artificial precession rate $m_0$9 is proportional to $m_1,\ldots,m_N$0 and increases rapidly with orbital eccentricity $m_1,\ldots,m_N$1.
• This precession exceeds the physical general-relativistic precession $m_1,\ldots,m_N$2 cy$m_1,\ldots,m_N$3 when $m_1,\ldots,m_N$4, with

$m_1,\ldots,m_N$5

where $m_1,\ldots,m_N$6 is the planet’s mean motion.

• For $m_1,\ldots,m_N$7–$m_1,\ldots,m_N$8, the DHC-induced precession grows rapidly at $m_1,\ldots,m_N$9, undermining accuracy at typical symplectic timesteps; Jacobi-coordinate splits remain accurate for $M = \sum_{i=0}^N m_i$0 up to $M = \sum_{i=0}^N m_i$1 days or greater, even at high $M = \sum_{i=0}^N m_i$2 (Rein et al., 12 Jan 2026).

4. Numerical Integration Outcomes in the Solar System

Long-term ensemble integrations serve to illustrate the practical consequences of DHC’s artificial precession. Key findings from $M = \sum_{i=0}^N m_i$3 Gyr Solar System runs are:

Integration Setting	Δt	Instability Rate (%)	Notes
DHC	≈ 6 days	0	Mercury instabilities artificially suppressed
DHC	≈ 0.6 days	≈ 1	Agrees with Jacobi results
Jacobi (with GR)	6 days	100	Instabilities present (converged)
DHC (restart, high $M = \sum_{i=0}^N m_i$4)	6 days	0	Artificially stabilized
DHC (restart, high $M = \sum_{i=0}^N m_i$5)	0.6 days	100	Converged (instabilities appear)

These results demonstrate that with standard DHC timesteps (several days), secular increases in eccentricity (e.g., Mercury’s $M = \sum_{i=0}^N m_i$6 surge via $M = \sum_{i=0}^N m_i$7–$M = \sum_{i=0}^N m_i$8 resonance) are suppressed by the numerical artifact, leading to an unrealistically low instability rate. Reducing the timestep by an order of magnitude restores physical instability rates and agreement with Jacobi splits (Rein et al., 12 Jan 2026).

5. Comparison with Jacobi Coordinate Wisdom–Holman Integrators

Jacobi coordinates employ a hierarchical splitting wherein planet $M = \sum_{i=0}^N m_i$9 orbits the cumulative mass $Q_0 = \frac{1}{M} \sum_{i=0}^N m_i r_i$0. The Kepler and interaction Hamiltonians are:

• $Q_0 = \frac{1}{M} \sum_{i=0}^N m_i r_i$1, with reduced mass $Q_0 = \frac{1}{M} \sum_{i=0}^N m_i r_i$2,
• $Q_0 = \frac{1}{M} \sum_{i=0}^N m_i r_i$3.

No separate “jump” term appears; all barycentric effects for inner bodies are handled exactly by the Kepler solver. As a result, the Jacobi split does not introduce an $Q_0 = \frac{1}{M} \sum_{i=0}^N m_i r_i$4-dependent precessional error, and integration remains accurate for large $Q_0 = \frac{1}{M} \sum_{i=0}^N m_i r_i$5, even at high eccentricity. There is no requirement to reduce $Q_0 = \frac{1}{M} \sum_{i=0}^N m_i r_i$6 as $Q_0 = \frac{1}{M} \sum_{i=0}^N m_i r_i$7 increases, which is reflected in the convergence of Mercury’s instability rate at much larger timesteps (Rein et al., 12 Jan 2026).

6. Practical Implementation Guidelines

For reliable long-term N-body integrations in planetary systems:

• Prefer Jacobi-coordinate WH integrators whenever the goal is to capture large-$Q_0 = \frac{1}{M} \sum_{i=0}^N m_i r_i$8 secular instabilities (e.g., Mercury’s $Q_0 = \frac{1}{M} \sum_{i=0}^N m_i r_i$9–$P_0 = \sum_{i=0}^N p_i$0 resonance). These allow stable integration with timesteps of several days without introducing artificial precession.
• Using DHC: If required (e.g., for hybrid or parallel-in-time algorithms), ensure $P_0 = \sum_{i=0}^N p_i$1 or, more conservatively, as recommended by Wisdom (2015), $P_0 = \sum_{i=0}^N p_i$2, with $P_0 = \sum_{i=0}^N p_i$3 as defined above.
• High-e switching: Once eccentricities reach $P_0 = \sum_{i=0}^N p_i$4–$P_0 = \sum_{i=0}^N p_i$5, integrations should switch to either hybrid symplectic algorithm (including a Bulirsch–Stoer pericenter solver) or a high-accuracy non-symplectic scheme to preserve physical fidelity.

The strong $P_0 = \sum_{i=0}^N p_i$6-dependence and slow convergence of DHC’s “jump” term error impose stringent constraints on timestep selection and justify the continued preference for Jacobi splits in high-fidelity Solar System chaos studies (Rein et al., 12 Jan 2026).