---
title: "Derive the quadrupole eccentricity cycle"
description: "Average the companion’s tidal potential, derive the vector equations and recover the critical inclination and maximum eccentricity."
canonical_url: "https://sandboxphysics.com/stories/secular-orbits/derive-the-quadrupole-cycle/"
language: "en"
page_type: "story"
alternate_languages:
  zh-CN: "https://sandboxphysics.com/zh/stories/secular-orbits/derive-the-quadrupole-cycle/"
last_modified: "2026-10-28"
---

[Kozai–Lidov cycles from torque to planetary orbits](https://sandboxphysics.com/stories/secular-orbits/)

AC09 · L2

# Derive the quadrupole eccentricity cycle

Sandbox Physics · 2026-10-28

The quadrupole cycle follows from averaging the companion’s differential gravitational potential over both orbits. Its conserved energy and axial angular momentum predict the eccentricity maximum. The vector equations evolve the same quantities without singular orbital angles.

## Before you begin

Chapter L1 supplies the normalized angular momentum, eccentricity and secular time. Some derivatives and averages over a Kepler ellipse appear below.

## Remove the acceleration shared by star and planet

Let the companion’s displacement from the primary be the outer vector, and the planet’s displacement be the inner vector. The primary is accelerating toward its companion. In a nonrotating frame centered on the primary, an indirect potential must remove that common acceleration.

Expand the companion potential for an inner distance much smaller than the outer distance. The term linear in the inner position cancels the indirect potential. A constant changes no force. The first surviving term is quadratic in the inner position, the quadrupole tidal potential.

Now average over a circular outer orbit. Its direction sweeps uniformly in a plane. The average of its squared projection onto any inner position is half the squared part of that position lying in the outer plane. This gives an axisymmetric potential about the outer normal. At this step we have averaged the outer clock but have not yet averaged the inner orbit.

$$
U=-\frac{Gm_2}{|\boldsymbol R-\boldsymbol r|}+\frac{Gm_2\boldsymbol R\cdot\boldsymbol r}{R^3}
$$

The gradient of the second term subtracts the primary’s acceleration.

$$
U_Q=-\frac{Gm_2}{2R^3}\left[3(\boldsymbol r\cdot\hat{\boldsymbol R})^2-r^2\right],\qquad \langle U_Q\rangle_2=-\frac{Gm_2}{4a_2^3}\left[r^2-3(\boldsymbol r\cdot\hat{\boldsymbol n}_2)^2\right]
$$

The first expression keeps the outer phase; the second replaces its circular motion by an orbit average.

Reading: Tremaine & Yavetz, Why do Earth satellites stay up?, arXiv:1309.5244v1 · Liu, Muñoz & Lai, Short-range forces in triple systems, arXiv:1409.6717v2

## Average around the ellipse with the right clock

Choose a unit vector along periapsis and a second unit vector in the ellipse plane perpendicular to it. Eccentric anomaly parameterizes the position using these two vectors. Equal steps of eccentric anomaly do not represent equal amounts of time: the averaging measure contains the extra factor below.

For example, integrate the squared distance over this measure. Odd cosine terms vanish, while the average of cosine squared is one-half. The result contains an eccentricity-squared correction. Applying the same procedure to the squared projection onto the outer axis yields the second identity.

Substitute both integrals into the outer-averaged potential, divide by the potential scale defined in L1, and collect terms. The resulting dimensionless Hamiltonian depends only on the eccentricity vector and the angular-momentum vector. The eccentricity vector points toward periapsis and has eccentricity as its length. Its direction therefore carries the orientation of the ellipse that a scalar eccentricity alone would lose.

$$
\boldsymbol r=a_1\left[(\cos E-e)\hat{\boldsymbol p}+\sqrt{1-e^2}\sin E\,\hat{\boldsymbol q}\right],\quad dM=(1-e\cos E)\, dE
$$

A time average is an integral over mean anomaly, divided by one full revolution.

$$
\langle r^2\rangle_1=a_1^2\left(1+\frac32e^2\right),\qquad \left\langle(\boldsymbol r\cdot\hat{\boldsymbol n}_2)^2\right\rangle_1=\frac{a_1^2}{2}\left[1-e^2-j_\parallel^2+5e_\parallel^2\right]
$$

Here the parallel eccentricity is the eccentricity vector’s projection along the outer normal.

$$
\phi_Q=\frac{\langle U_Q\rangle_{1,2}}{\Phi_0}=\frac18\left[1-6e^2-3j_\parallel^2+15e_\parallel^2\right]
$$

This is potential per unit test-particle mass. A convention using the negative disturbing function reverses its sign.

Reading: Tremaine & Yavetz, Why do Earth satellites stay up?, arXiv:1309.5244v1

## Turn one Hamiltonian into two vector equations

The angular-momentum and eccentricity vectors describe the same ellipse, so their components are constrained. They are perpendicular, and their squared lengths add to one. Their orbital Poisson brackets produce the cross-product evolution written below; using ordinary Cartesian position-momentum Hamilton equations on these constrained vectors would give the wrong dynamics.

Differentiate the polynomial Hamiltonian with respect to both vectors, insert the derivatives into the cross-product rule and normalize time by the secular time. The explicit equations below are exactly the averaged equations integrated in the Lab.

Taking scalar products verifies the constraints directly: a vector has zero scalar product with its own cross product. The same calculation preserves the axial angular momentum, while the autonomous Hamiltonian preserves its value. The numerical solver retains the raw vectors so these conditions can test the integration; it does not project an inaccurate state back onto the constraints.

$$
\frac{d\boldsymbol j}{d\tau}=-\boldsymbol j\times\nabla_j\phi-\boldsymbol e\times\nabla_e\phi,\qquad \frac{d\boldsymbol e}{d\tau}=-\boldsymbol j\times\nabla_e\phi-\boldsymbol e\times\nabla_j\phi
$$

The vector brackets and potential convention fix both signs.

$$
\boldsymbol j\cdot\boldsymbol e=0,\quad j^2+e^2=1,\quad \boldsymbol j\cdot\hat{\boldsymbol n}_2=\mathrm{constant},\quad \phi_Q=\mathrm{constant}
$$

These are independent diagnostics of the raw averaged trajectory.

$$
\boldsymbol j^{\prime}=\frac34\left[j_\parallel\boldsymbol j\times\hat{\boldsymbol n}_2-5e_\parallel\boldsymbol e\times\hat{\boldsymbol n}_2\right]
$$

A prime denotes a derivative with respect to the dimensionless secular time.

$$
\boldsymbol e^{\prime}=\frac34\left[2\boldsymbol j\times\boldsymbol e-5e_\parallel\boldsymbol j\times\hat{\boldsymbol n}_2+j_\parallel\boldsymbol e\times\hat{\boldsymbol n}_2\right]
$$

The changing eccentricity direction describes periapsis rotation as well as stretching.

Reading: Tremaine & Yavetz, Why do Earth satellites stay up?, arXiv:1309.5244v1 · Liu, Muñoz & Lai, Short-range forces in triple systems, arXiv:1409.6717v2

## Recover the critical angle and an eccentricity maximum

Take the limit of an initially almost circular orbit. Its axial angular momentum approaches the initial inclination cosine. On the growing branch, the eccentricity maximum occurs with periapsis ninety degrees from the line of nodes. The eccentricity vector’s axial projection then has the magnitude of eccentricity times the inclination sine.

Use the conserved axial momentum to eliminate inclination from the Hamiltonian. Equating its value at the initial state and at the maximum gives the factored expression below. The nonzero solution exists above the critical mutual inclination of about 39.23 degrees, on the prograde branch.

For an initial inclination of sixty degrees the predicted maximum is about 0.763763, giving a periapsis of about 0.236237 times the semimajor axis. This is an analytic limiting example. The default finite initial eccentricity of 0.05 has a slightly different conserved Hamiltonian and reaches about 0.838909 at sixty-five degrees. The zero-seed limit gives 0.838047.

An exactly circular vector state remains exactly circular under these equations. The unstable solution grows from an arbitrarily small perturbation, and its growth time depends on that seed. A solver that silently adds eccentricity would change the initial-value problem. The Lab displays and saves its finite initial seed.

$$
0=e_{\max}^2\left[9-\frac{15\cos^2 i_0}{1-e_{\max}^2}\right],\qquad e_{\max}=\sqrt{1-\frac53\cos^2i_0}
$$

The expression uses the near-circular limit and the quadrupole test-particle Hamiltonian.

$$
i_{\mathrm{crit}}=\arccos\sqrt{\frac35}=39.2315^\circ
$$

A larger inclination opens the growing branch for this particular circular-outer approximation.

Reading: Tremaine & Yavetz, Why do Earth satellites stay up?, arXiv:1309.5244v1

## Which slow angle can resonate?

The argument of periapsis measures the direction of the eccentricity vector from the ascending node within the inner orbit plane. The ascending node is where the planet crosses the outer plane toward its positive normal. With that convention, the eccentricity vector’s axial component is eccentricity times the inclination sine times the periapsis-angle sine.

Substitution puts the Hamiltonian into an eccentricity-and-angle form. Contours of fixed Hamiltonian describe the allowed slow motion. Some trajectories circulate through every periapsis direction; others librate around ninety or two hundred seventy degrees. Libration means the slow angle swings around a preferred value instead of winding through full turns, the resonance behavior usually associated with the Kozai–Lidov mechanism.

A large eccentricity cycle alone does not distinguish these possibilities. With the finite seed and periapsis direction used here, the default averaged record circulates: its unwrapped periapsis angle advances about 408 degrees over twelve secular times. The exported eccentricity and angular-momentum vectors contain the information needed to reconstruct that angle. Close to a circle, the vector remains useful even while the angle itself becomes poorly defined.

$$
\phi_Q=\frac18\left[1-6e^2-3j_\parallel^2+15e^2\left(1-\frac{j_\parallel^2}{1-e^2}\right)\sin^2\omega\right]
$$

The argument of periapsis is measured from the ascending node in the inner plane; the conserved axial projection eliminates inclination.

Reading: Tremaine & Yavetz, Why do Earth satellites stay up?, arXiv:1309.5244v1 · Liu, Muñoz & Lai, Short-range forces in triple systems, arXiv:1409.6717v2

## Add conservative relativistic periapsis precession

An additional potential inversely proportional to angular-momentum length generates the leading orbit-averaged relativistic periapsis precession. It contributes no angular-momentum torque and rotates the eccentricity vector within the orbit plane. At high eccentricity its rate becomes rapid, interrupting the coherent Newtonian torque in this doubly averaged model.

The default orbit with the dimensionless relativistic strength set to 0.8 reaches a maximum eccentricity about 0.629278. This value follows from the modified conserved Hamiltonian and agrees with a separately solved energy root in the tests. No relativistic force is added to the saved direct reference; the amber curve remains the explicitly labelled Newtonian comparison.

The precession term assumes a weak gravitational field and slow motion. The inner orbit must remain far outside the gravitational radius, and the integration stops before the stated near-primary boundary or before the gravitational radius reaches one percent of periapsis. The latter is a conservative chosen validity cutoff, not a sharp physical transition. Dissipative tides and gravitational radiation are omitted. Beyond double averaging, the outer orbital phase can interact with precession and change the peak from one cycle to the next. The final paper guide follows this behavior.

$$
\phi=\phi_Q-\frac{\epsilon_{\mathrm{GR}}}{j},\quad \left.\boldsymbol e^{\prime}\right|_{\mathrm{GR}}=\frac{\epsilon_{\mathrm{GR}}}{j^3}\boldsymbol j\times\boldsymbol e,\quad \dot\omega_{\mathrm{GR}}=\frac{\epsilon_{\mathrm{GR}}}{t_Kj^2}
$$

The coefficient is dimensionless; its relation to the physical speed of light is fixed below.

$$
\epsilon_{\mathrm{GR}}=\frac{3GM}{a_1c^2}\frac{M}{m_2}\left(\frac{a_2}{a_1}\right)^3,\qquad \frac{GM}{c^2q}\ll1
$$

Changing the independent teaching coefficient at fixed masses and lengths changes the implied relativistic compactness.

Reading: Liu, Muñoz & Lai, Short-range forces in triple systems, arXiv:1409.6717v2 · Klein & Hamilton, Short-range forces beyond double averaging, arXiv:2606.19456v1

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.

| Finite initial eccentricity, initial inclination | $e_0=0.0500,\quad i_0=65.0^\circ$ |
| --- | --- |
| Default maximum averaged eccentricity | $e_{\max}=0.838909$ |
| Near-circular analytic limiting maximum | $e_{\max}(e_0\to0)=0.838047$ |
| Maximum with conservative averaged GR | $\epsilon_{\mathrm{GR}}=0.800,\quad e_{\max}=0.629278$ |
| Representative averaged Hamiltonian | $\phi_Q=0.056315120$ |

The representative frame is taken at 1.8 secular times. The finite initial eccentricity and the near-circular analytic limit are kept distinct. The downloadable analysis reproduces the maxima and supplies actual direct-state records.

## [Try the orbit exchange](https://sandboxphysics.com/labs/kozai-lidov-cycles/)

1. Compare the exchange and relativistic-precession presets. The averaged eccentricity maximum falls while the Newtonian direct reference retains its stated force model.
2. Increase the initial tilt toward the near-polar preset and watch the record stop at its exclusion boundary. A partial high-eccentricity excursion has no measured full cycle period.

## Exercises with explanations

1. Why would an average uniformly sampled in eccentric anomaly give the wrong tidal potential?

   Read the explanation

   It overweights the quickly traversed part near periapsis and underweights the time spent near apoapsis. The mean-anomaly measure supplies the missing time factor; that factor changes the eccentricity-squared coefficients.

2. Use the near-circular result for an initial inclination of thirty degrees. What does the negative square-root argument mean?

   Read the explanation

   It means the nonzero growing-branch solution is unavailable. It is not an imaginary physical eccentricity. The orbit lies on the stable low-eccentricity branch of this approximation.

3. A numerical run keeps eccentricity bounded but its vector orthogonality drifts. Is that enough to trust its maximum eccentricity?

   Read the explanation

   No. The vectors no longer describe one physical Kepler ellipse. Check both geometric constraints, the conserved projection and Hamiltonian. A smaller timestep or corrected integrator must recover their accuracy; clipping a scalar leaves the vector geometry broken.

[Previous: Angular momentum connects orbital shape and tilt](https://sandboxphysics.com/stories/secular-orbits/angular-momentum-and-slow-torques/) [Next: Compare trajectories without confusing three different errors](https://sandboxphysics.com/stories/secular-orbits/compare-trajectories-and-averaging/)

## Sources and calculation scope

1. [Tremaine & Yavetz, Why do Earth satellites stay up?, arXiv:1309.5244v1](https://arxiv.org/abs/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.

2. [Liu, Muñoz & Lai, Short-range forces in triple systems, arXiv:1409.6717v2](https://arxiv.org/abs/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.

3. [Klein & Hamilton, Short-range forces beyond double averaging, arXiv:2606.19456v1](https://arxiv.org/html/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.