Compare trajectories without confusing three different errors

Define the same initial-value problem twice

The primary and companion each have unit mass. The inner planet is massless, so it does not change their circular relative orbit. At the initial time both relative position vectors point along the positive first axis. The inner orbit lies in the first two coordinate directions; the outer orbit is tilted about the first axis.

The companion-primary relative angular speed uses their total mass. Its prescribed displacement is therefore the expression below. In the primary-centered frame, the planet feels the primary, the companion and the indirect acceleration of the primary. Every term is three-dimensional.

The saved direct reference was actually integrated with REBOUND 5.2.2 and IAS15 at an accuracy parameter of one hundred-billionth. It stores Cartesian states every 0.005 secular times. The browser uses downsampled states for the named presets and an adaptive Dormand–Prince method for new parameter choices. The averaged prediction starts from the angular-momentum and eccentricity vectors derived from those same initial Kepler elements.

\boldsymbol R(t)=a_2(\cos n_2t,\ \sin n_2t\cos i_0,\ \sin n_2t\sin i_0),\qquad n_2=\sqrt{\frac{2}{a_2^3}}

The adopted units set both stellar masses and the gravitational constant to one.

\ddot{\boldsymbol r}=-\frac{\boldsymbol r}{r^3}+\frac{\boldsymbol R-\boldsymbol r}{|\boldsymbol R-\boldsymbol r|^3}-\frac{\boldsymbol R}{R^3}

The last term is required because the origin accelerates. The axes themselves do not rotate.

\boldsymbol r_0=(1-e_0,0,0),\quad \boldsymbol v_0=\left(0,\sqrt{\frac{1+e_0}{1-e_0}},0\right)

These initial states have unit inner semimajor axis and periapsis on the positive first axis.

Reading: REBOUND 5.2.2 official Kozai example and IAS15 integrator · Tremaine & Yavetz, Why do Earth satellites stay up?, arXiv:1309.5244v1

Extract orbit elements from the saved states

At each saved time, take the position-velocity cross product and form the eccentricity vector below. Its length is the osculating eccentricity; the angle of angular momentum to the outer normal is the mutual inclination. No sine curve or separately drawn eccentricity history enters this extraction.

The inner two-body energy varies because the companion exchanges energy with the planet. For a circular companion there is instead a conserved rotating Hamiltonian. Its final term subtracts the outer angular speed times the axial angular momentum. The indirect potential remains present.

This distinction makes a useful error test. A nearly constant inner Kepler energy is not the right criterion for a correct restricted three-body calculation. The circular rotating Hamiltonian should remain constant to numerical accuracy while the osculating elements oscillate. The provided scientific tests also differentiate that Hamiltonian along the full force and check that its derivative vanishes.

\boldsymbol h=\boldsymbol r\times\boldsymbol v,\quad \boldsymbol e=\boldsymbol v\times\boldsymbol h-\frac{\boldsymbol r}{r},\quad \cos i=\frac{\boldsymbol h\cdot\hat{\boldsymbol n}_2}{h}

The formula uses the unit primary gravitational parameter. Inclination is extracted from the same saved Cartesian state.

K=\frac{v^2}{2}-\frac1r-\frac1{|\boldsymbol R-\boldsymbol r|}+\frac{\boldsymbol R\cdot\boldsymbol r}{R^3}-n_2\boldsymbol h\cdot\hat{\boldsymbol n}_2

This invariant applies to the circular, Newtonian, massless restricted problem used by the direct reference.

Reading: REBOUND 5.2.2 official Kozai example and IAS15 integrator · Tremaine & Yavetz, Why do Earth satellites stay up?, arXiv:1309.5244v1

Check the integrator before judging the approximation

For the default first three secular times, the browser’s adaptive integrator was compared with the independent saved IAS15 states. The maximum eccentricity difference was about four ten-millionths. The largest position difference was about eight ten-thousandths of the inner semimajor axis. Accumulated fast orbital phase error contributes to that residual. The rotating-Hamiltonian drift was about four hundred-millionths in the adopted energy units.

Those errors are smaller than the percent-level direct-versus-averaged eccentricity differences in the default comparison. This does not certify every possible long integration. It establishes a specific reference check over the stated interval, with both an orbital-element diagnostic and a more phase-sensitive Cartesian diagnostic.

The averaged solver is checked differently. Finite numerical integrals of the instantaneous tidal potential reproduce its orbit average; the raw vectors preserve their geometry, Hamiltonian and axial projection. These tests probe the force and normalization as well as timestep convergence.

Reading: REBOUND 5.2.2 official Kozai example and IAS15 integrator · Tremaine & Yavetz, Why do Earth satellites stay up?, arXiv:1309.5244v1

Use a third force model to isolate averaging

A full Newtonian direct calculation retains all orders of the companion potential. The doubly averaged curve also drops terms beyond quadrupole order. Their difference therefore combines two physical approximations.

The public materials add a controlled intermediate calculation: IAS15 integrates the instantaneous quadrupole tidal acceleration, retaining both fast orbital clocks. Comparing it with the doubly averaged quadrupole changes averaging alone. Comparing it with the full Newtonian force changes multipole truncation at the same initial state.

At a mutual inclination of seventy-five degrees and an outer separation six times the inner semimajor axis, the measured periods are 5.005, 4.810 and 4.670 secular times for the full force, instantaneous quadrupole and double average. The outer-orbit smoothing window changes the quadrupole estimate between 4.745 and 4.810. At a separation ratio of twenty, the corresponding periods are 4.675, 4.665 and 4.670. These are finite-record measurements, not universal correction factors.

As a worked comparison, subtract the averaged compact period from the instantaneous-quadrupole period and divide by the averaged period. The relative shift is about three percent. Repeating with the full force gives about seven percent, so calling the whole shift an averaging error would assign the omitted spatial terms to the wrong cause.

The compact case visibly accumulates phase differences. Because nonlinear evolution changes which part of a cycle is sampled, the two period shifts are not guaranteed to add linearly. Use the three-model experiment to identify which approximations matter, and retain the raw trajectory for another choice of summary statistic.

\boldsymbol a_{Q,\mathrm{tidal}}=\frac{1}{R^3}\left[3(\boldsymbol r\cdot\hat{\boldsymbol R})\hat{\boldsymbol R}-\boldsymbol r\right]

This intermediate direct reference truncates the spatial force expansion but keeps the outer phase. Its central primary attraction is still included.

Reading: Liu, Muñoz & Lai, Short-range forces in triple systems, arXiv:1409.6717v2 · REBOUND 5.2.2 official Kozai example and IAS15 integrator · Klein & Hamilton, Short-range forces beyond double averaging, arXiv:2606.19456v1

Stop where the declared physical problem ends

The near-polar preset reaches the near-primary exclusion before completing a cycle. Its saved full-force record ends at 2.180 secular times; the averaged record ends near 2.160. The last successful state is retained and no later curve is invented. A full period cannot be measured from this partial excursion.

The browser also stops a custom direct record for a close companion approach, loss of a bound inner osculating ellipse, nonfinite state or exhausted integration budget. These are reported boundaries, not collision or tidal-disruption models. Body radii in the scene are schematic.

The main scene displays the averaged ellipse and its changing plane. Its companion marker is a grip on the representative outer plane; the fast orbital phases have been removed. Inner and outer paths have separately stated display magnifications; the numerical records retain their actual units. Use the data row for the direct reference and its residual.

A convincing comparison states which model evolved each curve, what initial state it used, which time interval survived, and which observable was measured. Those details allow another reader to ask a sharper physical question with the same saved experiment.

Reading: REBOUND 5.2.2 official Kozai example and IAS15 integrator · Liu, Muñoz & Lai, Short-range forces in triple systems, arXiv:1409.6717v2

Saved numerical example and method

Synthetic dynamics in units with gravitational constant, primary mass, companion mass and initial inner semimajor axis all equal to one. Raw-vector RK4 uses a maximum secular step of 0.001; the direct reference uses frozen REBOUND 5.2.2 IAS15 states at 0.005 secular-time cadence.

Default full-force peaks	$\tau_1=2.510,\quad\tau_2=7.565$
Default direct and averaged periods	$P_{\mathrm{full}}/t_K=5.055,\quad P_{\mathrm{DA}}/t_K=5.025$
Compact, same 75-degree initial tilt	$A=6:\quad P/t_K=5.005,\ 4.810,\ 4.670$
Wide, same 75-degree initial tilt	$A=20:\quad P/t_K=4.675,\ 4.665,\ 4.670$
Near-polar completed intervals	$\tau_{\mathrm{full}}=2.180,\quad\tau_{\mathrm{DA}}=2.160$

The compact and wide rows list full Newtonian, instantaneous quadrupole and doubly averaged quadrupole periods in that order. Periods come from grouped eccentricity peaks after one-outer-period smoothing. A stopped near-polar excursion supplies no full period.

Try the orbit exchange

Compare the wide-hierarchy and compact-hierarchy presets. Both begin at seventy-five degrees; their different direct-versus-averaged agreement is therefore easier to interpret.
Export a record after changing inclination. Check the direct force description, completed interval and averaged relativistic coefficient before plotting the two eccentricities together.

Exercises with explanations

The full direct curve differs from the double average, but an instantaneous quadrupole curve agrees with the double average. Which approximation is implicated?
Read the explanation
For that observable and interval, truncating the spatial force expansion is implicated. Averaging has not produced a comparable discrepancy in the intermediate control. Numerical convergence must still be checked before interpreting either result.
A position residual grows while the eccentricity residual remains tiny. What extra diagnostic would distinguish numerical failure from phase accumulation?
Read the explanation
Check the circular rotating Hamiltonian, repeat at tighter tolerances and compare orbital phase separately from osculating shape. Fast phase is more sensitive than the slowly varying eccentricity.
Transfer the comparison to an eccentric outer orbit. Can the same rotating Hamiltonian and circular outer clock be retained?
Read the explanation
No. The companion’s speed and separation now vary, the circular rotating invariant no longer applies in that form, and an octupole contribution may survive averaging. Start by specifying the new Kepler clock and choosing diagnostics appropriate to the nonautonomous force.

Sources and calculation scope

Tremaine & Yavetz, Why do Earth satellites stay up?, arXiv:1309.5244v1
Fixed 20 September 2013 version, DOI 10.1119/1.4874853. Potential 18, vector equations 23, constraints 24 and section IV give the quadrupole model and near-circular critical-angle limit.
Liu, Muñoz & Lai, Short-range forces in triple systems, arXiv:1409.6717v2
Fixed 7 January 2015 version, DOI 10.1093/mnras/stu2396. Sections 2–4, quadrupole potential 7, octupole 8, GR 32–33/42 and Figures 3–4/7. The browser retains the circular-outer, test-particle quadrupole; the guide reads the wider model.
REBOUND 5.2.2 official Kozai example and IAS15 integrator
The official example motivates an independent direct reference. Frozen trajectories use REBOUND 5.2.2, source commit c34134716cd484ed9cdb9f7d54ef8d859bdf547d, IAS15 accuracy parameter one hundred-billionth, with separately stated initial conditions and stopping rules. REBOUND software is GPL-3.0-only; no publisher code is copied into the Lab.
Klein & Hamilton, Short-range forces beyond double averaging, arXiv:2606.19456v1
Fixed 17 June 2026 preprint version; an August v2 also exists. Figure 1 separates averaging/force choices, Figures 2–3 show cycle changes, and Appendix B compares approximations with direct N-body examples. This is a reading guide; the browser does not reproduce its single-averaged GR/tidal calculations.