Gravitational N-Body Problem

The gravitational N-body problem is the foundational dynamical problem of astrophysics: given $N$ point masses $m_i$ at positions $\mathbf{r}_i$ , interacting only through their mutual Newtonian attraction, find their motion for all time. Each mass obeys $\ddot{\mathbf{r}}_i = -G\sum_{j\neq i} m_j (\mathbf{r}_i-\mathbf{r}_j)/|\mathbf{r}_i-\mathbf{r}_j|^3$ — a set of $3N$ coupled, nonlinear, second-order ordinary differential equations. Only the $N=2$ case admits a general closed-form solution, the elliptical, parabolic, and hyperbolic Kepler orbits. For $N=3$ Poincaré proved in the 1890s that no further single-valued analytic integrals beyond the classical ten exist, and that the phase space is threaded by homoclinic tangles — founding modern chaos theory. What decides the physics is not the force law but the number of bodies. When $N$ is small, every close encounter matters and the graininess of the gravitational field drives the evolution. When $N$ is enormous — a galaxy of $\sim10^{11}$ stars — individual encounters are negligible over the age of the Universe, and the lumpy potential can be replaced by a single smooth mean field $\Phi(\mathbf{x},t)$ . That is the collisionless limit, and it legitimizes the whole smooth-potential machinery of galactic dynamics: the collisionless Boltzmann equation, Jeans' theorem, the Jeans equations. Between the extremes sit globular clusters and galactic nuclei ( $N\sim10^5$ – $10^7$ ), where two-body relaxation drives core collapse, mass segregation, and evaporation on cosmologically interesting timescales. This article is the hinge between celestial mechanics and stellar dynamics.

Key takeaways

The gravitational N-body problem seeks the trajectories of $N$ point masses coupled only by Newtonian gravity, governed by $3N$ second-order ODEs — one line to write, unsolvable in closed form for $N\geq3$ .
For $N=2$ it collapses to the exactly integrable Kepler problem; for $N\geq3$ Poincaré showed no general analytic solution exists and the motion is generically chaotic.
Exactly ten classical integrals — energy, total momentum, center-of-mass drift, angular momentum — reduce the dimensionality but fall short of the $6N-1$ needed to integrate the flow once $N\geq3$ .
Two timescales decide everything: the crossing time $t_{\mathrm{cross}}\sim R/v$ and the two-body relaxation time $t_{\mathrm{relax}}\sim(0.1N/\ln N)\,t_{\mathrm{cross}}$ , dominated by many distant weak encounters through the Coulomb logarithm $\ln\Lambda\approx\ln N$ .
Galaxies ( $N\sim10^{11}$ ) have $t_{\mathrm{relax}}\sim10^{16}$ yr, far longer than their age, so they behave as a smooth collisionless fluid; star clusters are genuinely collisional. The practical tool is numerical integration, with cost scaling from $O(N^2)$ down to $O(N\log N)$ .
Key idea: the character of an N-body system is fixed by the ratio of its two-body relaxation time to its age — when $t_{\mathrm{relax}}\gg t_{\mathrm{age}}$ (galaxies) the discrete masses act as a smooth collisionless fluid in phase space; when $t_{\mathrm{relax}}\lesssim t_{\mathrm{age}}$ (star clusters) individual encounters reshape the system.

The core insight

The character of an N-body system is set by the ratio of its two-body relaxation time to its age: when $t_{\mathrm{relax}}\gg t_{\mathrm{age}}$ (galaxies) the discrete masses behave as a smooth collisionless fluid; when $t_{\mathrm{relax}}\lesssim t_{\mathrm{age}}$ (star clusters) individual encounters reshape the system.

Where this lives today

The N-body problem is now a computational science, spanning three bodies to trillions. At the small- $N$ end, mission designers exploit the chaotic manifolds of the restricted three-body problem — the same tangles Poincaré found — to route probes on almost no fuel, as the Genesis and GRAIL missions did. At intermediate $N$ , direct integrators such as Aarseth's NBODY6 follow globular clusters star by star through core collapse, binary hardening, and black-hole ejection — the physics of worked example 1. At the largest scale, cosmological simulations — Millennium, IllustrisTNG, trillion-particle runs on machines like Frontier — evolve dark matter as a collisionless N-body fluid using tree and particle-mesh solvers that cut the cost from $O(N^2)$ to $O(N\log N)$ , and their halos are weighed with the same virial theorem that weighed Coma. The through-line is unbroken: because no pen can write $\mathbf{r}_i(t)$ for $N\geq3$ , the field advances by numerical integration and by statistical descriptions wherever $N$ permits.

The physical picture

Picture $N$ stars as points threading a shared potential well. To an excellent approximation each star's orbit is set by the smooth mass distribution of all the others — a mean field felt as a continuous force — and only occasionally does a near neighbour deliver a discrete kick. The subject is the competition between the smooth field, which conserves each star's orbital integrals, and these granular kicks, which slowly randomize them. The first figure separates the two: a rosette orbit in the smooth well, and a single deflecting encounter whose transverse impulse is $\delta v_\perp = 2Gm/(bv)$ .

Figure

Two faces of the N-body problem. Left: a test star threads the smooth mean-field potential on a rosette orbit that conserves energy and angular momentum. Right: a single passage at impact parameter b delivers a transverse kick δv⊥ = 2Gm/(bv) — the granular deflection that drives two-body relaxation

The two-body problem is integrable because, in the center-of-mass frame, it has one relative coordinate and just enough conserved quantities — energy plus the angular-momentum vector — to lock the orbit into a closed conic. Add a third body and that closure is gone: too few isolating integrals remain, so trajectories wander over the energetically allowed region and nearby initial conditions diverge exponentially. The escape is statistical. Relaxation is a random walk in velocity: each distant star deflects the test star by a tiny angle, the deflections add incoherently, and the variance of the velocity change grows linearly with the number of encounters — it takes of order $N/(8\ln N)$ crossings for the accumulated change to rival the star's own speed. The second figure plots this ratio and shows why the same physics makes a globular cluster collisional and a galaxy collisionless.

Figure

The number of crossings a star survives before encounters randomize its velocity, t_relax ∕ t_cross = 0.1 N ∕ ln N, rising steeply with N. A globular cluster (N ≈ 10⁵) relaxes in ≈ 870 crossings ≈ 0.85 Gyr, less than its age, so it is collisional; a galaxy (N ≈ 10¹¹) needs ≈ 4×10⁸ crossings ≈ 4×10¹⁶ yr, far longer than cosmic age, so it is collisionless

The core equations

Because the potential is generated by the very bodies moving in it, the N-body problem is intrinsically nonlinear and self-gravitating — the source of the field is also its response. The equations of motion are Newton's second law for each body under the summed attraction of the rest.

\ddot{\mathbf{r}}_i = -G\sum_{j\neq i}^{N} m_j\,\frac{\mathbf{r}_i-\mathbf{r}_j}{|\mathbf{r}_i-\mathbf{r}_j|^{3}},\qquad i=1,\dots,N \qquad (1)

where: $\mathbf{r}_i$ — position vector of body $i$ ; $m_j$ — mass of body $j$ ; $G$ — Newton's gravitational constant, $6.674\times10^{-11}\ \mathrm{m^3\,kg^{-1}\,s^{-2}}$ ; $N$ — number of bodies. These are $3N$ coupled second-order ODEs; the $j\neq i$ excision removes the divergent self-force. Bundling the sum into a scalar potential gives the equivalent form $\ddot{\mathbf{r}}_i=-\nabla\Phi(\mathbf{r}_i)$ with $\Phi(\mathbf{r})=-G\sum_j m_j/|\mathbf{r}-\mathbf{r}_j|$ , and replacing that discrete $\Phi$ by a smooth $\Phi(\mathbf{x},t)$ sourced through Poisson's equation is precisely the step that defines the collisionless limit. The system is conservative, so the total mechanical energy is a constant of the motion — the most robust classical integral and a stringent check on any numerical integrator.

E = \sum_i \tfrac{1}{2}m_i|\dot{\mathbf{r}}_i|^2 - G\sum_{i<j}\frac{m_i m_j}{|\mathbf{r}_i-\mathbf{r}_j|} = \text{const} \qquad (2)

where: $E$ — total mechanical energy (kinetic plus potential); the $i<j$ sum counts each interacting pair once. Energy conservation constrains the flow but cannot integrate it. The symmetries of the Newtonian action supply nine more conserved quantities alongside energy — the ten classical integrals, whose count is exactly why $N=2$ closes and $N\geq3$ does not.

\mathbf{P}=\sum_i m_i\dot{\mathbf{r}}_i,\quad \mathbf{L}=\sum_i m_i\,\mathbf{r}_i\times\dot{\mathbf{r}}_i,\quad \mathbf{R}_{\rm cm}=\frac{\sum_i m_i\mathbf{r}_i}{\sum_i m_i} \qquad (3)

where: $\mathbf{P}$ — total linear momentum (3 integrals); $\mathbf{L}$ — total angular momentum (3 integrals); $\mathbf{R}_{\rm cm}$ — center of mass, giving 6 integrals jointly with $\mathbf{P}$ . With $E$ these total ten — enough for $N=2$ but not $N\geq3$ . The dynamical clock of any stellar system is its crossing time — how long a typical star takes to traverse it once, the unit against which every slower process is measured.

t_{\mathrm{cross}} \sim \frac{R}{v} \sim \left(\frac{R^3}{GM}\right)^{1/2} \qquad (4)

where: $R$ — characteristic system size; $v$ — typical stellar speed; $M$ — total mass. The second form uses $v^2\sim GM/R$ , so $t_{\mathrm{cross}}$ is essentially the free-fall time. The decisive quantity is the two-body relaxation time — how many crossings a star survives before the cumulative graininess of the potential randomizes its velocity; its ratio to the crossing time grows almost linearly with $N$ .

t_{\mathrm{relax}} \approx \frac{0.1\,N}{\ln N}\, t_{\mathrm{cross}} \qquad (5)

where: $t_{\mathrm{relax}}$ — time for accumulated encounters to deflect a star's velocity by order itself; $\ln N$ — the Coulomb logarithm $\ln\Lambda$ . The prefactor $0.1$ absorbs the $1/8$ and geometric factors. For a bound steady state, the time-averaged energies obey the scalar virial theorem — the most useful equilibrium relation in extragalactic astronomy, turning a velocity dispersion and a size into a mass.

2K + W = 0,\qquad K=\sum_i\tfrac12 m_i v_i^2,\quad W=-G\sum_{i<j}\frac{m_i m_j}{r_{ij}} \qquad (6)

where: $K$ — total kinetic energy in a steady bound state; $W$ — total gravitational potential energy. It follows that $E=K+W=-K<0$ : a self-gravitating system in equilibrium is always bound and heats up when it loses energy (negative heat capacity). In the large- $N$ limit the discrete swarm becomes a phase-space fluid whose density $f$ is conserved along orbits in the smooth potential — the collisionless Boltzmann equation, gateway to galactic dynamics.

\frac{\partial f}{\partial t} + \mathbf{v}\cdot\nabla f - \nabla\Phi\cdot\frac{\partial f}{\partial\mathbf{v}} = 0 \qquad (7)

where: $f(\mathbf{x},\mathbf{v},t)$ — phase-space distribution function; $\Phi$ — smooth self-consistent potential from $\nabla^2\Phi=4\pi G\rho$ . The right-hand side is zero precisely because two-body encounters are negligible when $t_{\mathrm{relax}}\gg t_{\mathrm{age}}$ .

Derivation AWhy N equals two integrates and N equals three does not

The goal is to make concrete the claim that the same ten integrals which fully solve the two-body problem leave the three-body problem open. We separate the center-of-mass and relative motions and watch the integrals peel away the degrees of freedom until an integrable one-dimensional problem remains for $N=2$ — then repeat the count for $N=3$ and show the reduction stalls with too many degrees of freedom still coupled.

\mu\ddot{\mathbf{r}} = -\frac{G m_1 m_2}{r^2}\hat{\mathbf{r}},\qquad \mu=\frac{m_1 m_2}{m_1+m_2}\qquad\text{(Kepler reduction)}

Step 1. Introduce the center of mass $\mathbf{R}_{\rm cm}=(m_1\mathbf{r}_1+m_2\mathbf{r}_2)/M$ ; since $\mathbf{P}$ is conserved it drifts uniformly and we work in its rest frame. Why: this uses 6 of the 10 integrals ( $\mathbf{P}$ and $\mathbf{R}_0$ ) to strip 6 of the 12 phase-space dimensions at a stroke.

Step 2. Form the relative coordinate $\mathbf{r}=\mathbf{r}_1-\mathbf{r}_2$ ; the two equations combine into one effective one-body problem with reduced mass $\mu$ , the boxed equation. Why: the CM and relative motions decouple exactly for an inverse-square pair force, so the hard part is a single particle in a fixed central field.

Step 3. Conservation of $\mathbf{L}=\mu\,\mathbf{r}\times\dot{\mathbf{r}}$ confines the motion to a plane, and its magnitude eliminates $\theta$ via $\dot\theta=L/(\mu r^2)$ . Why: angular momentum flattens the problem to 2D and reduces the radial motion to a single coordinate in an effective potential.

Step 4. Energy conservation $E=\tfrac12\mu\dot r^2 + L^2/(2\mu r^2) - GM\mu/r$ then gives $\dot r(r)$ , so $r(t)$ follows from one quadrature and the orbit is a conic. Why: with $E$ and $L$ fixed the last degree of freedom is integrable by separation of variables — the problem is solved in closed form.

Step 5. For $N=3$ the same $\mathbf{P},\mathbf{R}_0,\mathbf{L},E$ reductions (plus the ascending node) bring the phase space down to Poincaré's reduced problem, but no further single-valued analytic integral exists to close it. Why: too few isolating integrals leave the degrees of freedom coupled, the flow develops chaotic regions, and no closed-form $\mathbf{r}_i(t)$ can be written — Poincaré's non-integrability theorem.

Go deeper · Derivation B · 1 minFrom a single flyby to the relaxation time

Derivation BFrom a single flyby to the relaxation time

The goal is to build the relaxation time from the ground up: compute the transverse kick one distant star gives a passing test star, add such kicks as an incoherent random walk over all impact parameters, extract the Coulomb logarithm, and set the accumulated velocity change equal to the star's own speed to read off how many crossings relaxation takes — the calculation that ultimately tells a galaxy from a globular cluster.

|\delta v_\perp| \approx \frac{2Gm}{b\,v}\qquad\text{(impulse deflection)}

Step 1. For a nearly straight passage only the perpendicular force $F_\perp=Gm^2 b/(b^2+v^2t^2)^{3/2}$ accumulates; integrating, $\delta v_\perp=(1/m)\int_{-\infty}^{\infty}F_\perp\,\mathrm{d}t=(Gmb/v)\int_{-\infty}^{\infty}(b^2+x^2)^{-3/2}\mathrm{d}x=(Gmb/v)(2/b^2)=2Gm/(bv)$ . Why: the parallel force reverses across closest approach and cancels, so the net effect is the transverse impulse, and the integral evaluates to exactly $2/b^2$ .

Step 2. A single crossing brings the star past $\delta n = n\,(2\pi b\,\mathrm{d}b)\,(v t_{\mathrm{cross}})$ field stars near impact parameter $b$ ; since the kicks are random in direction they add in quadrature, so $\sum(\delta v_\perp)^2 \approx \int (2Gm/bv)^2\,\delta n$ . Why: with $\langle\delta v_\perp\rangle=0$ but $\langle\delta v_\perp^2\rangle\neq0$ , velocity diffuses as a random walk — the variance grows, not the mean.

Step 3. The integral over $b$ contains $\int \mathrm{d}b/b=\ln(b_{\max}/b_{\min})\equiv\ln\Lambda$ ; with $b_{\max}\sim R$ and $b_{\min}\sim Gm/v^2\sim R/N$ this gives $\ln\Lambda\approx\ln N$ , and collecting factors $\sum(\delta v_\perp)^2\approx 8(Gm/Rv)^2 N\ln N$ per crossing. Why: each decade of impact parameter contributes equally, so distant weak encounters dominate and the result is only logarithmically sensitive to the cutoffs — the Coulomb logarithm.

Step 4. Relaxation is complete when the accumulated $\sum(\delta v_\perp)^2$ reaches $v^2$ ; dividing $v^2$ by the per-crossing value gives $n_{\mathrm{relax}}\sim v^2/[8(Gm/Rv)^2 N\ln N]\sim N/(8\ln N)$ crossings, hence $t_{\mathrm{relax}}\sim[N/(8\ln N)]\,t_{\mathrm{cross}}$ . Why: writing $Gm/Rv\sim v/N$ for a virialized system makes the $N$ -dependence explicit — the more bodies, the smoother the field, the longer relaxation takes. *The key results are: the integral $\int(b^2+x^2)^{-3/2}\mathrm{d}x=2/b^2$ and the scaling $t_{\mathrm{relax}}\sim[N/(8\ln N)]t_{\mathrm{cross}}$ .*

Go deeper · Derivation C · 1 minThe virial theorem from the moment of inertia

Derivation CThe virial theorem from the moment of inertia

The goal is to derive $2K+W=0$ for a bound self-gravitating system directly from the equations of motion, by taking two time derivatives of the moment of inertia and time-averaging over a bounded orbit. The key move is that the $-1$ homogeneity of the Newtonian potential lets the awkward $\sum m_i\,\mathbf{r}_i\cdot\ddot{\mathbf{r}}_i$ collapse cleanly into the total potential energy $W$ .

\langle 2K + W\rangle = 0\qquad\text{(scalar virial)}

Step 1. Define the moment of inertia $I=\sum_i m_i|\mathbf{r}_i|^2$ , a positive scalar measuring the system's spatial extent. Why: its second derivative packages the kinetic and potential energies together, and its boundedness in equilibrium is what makes the theorem work.

Step 2. Differentiate twice: $\tfrac12\ddot I=\sum_i m_i(\dot{\mathbf{r}}_i\cdot\dot{\mathbf{r}}_i+\mathbf{r}_i\cdot\ddot{\mathbf{r}}_i)=2K+\sum_i m_i\,\mathbf{r}_i\cdot\ddot{\mathbf{r}}_i$ . Why: the first term is twice the kinetic energy; the second holds the gravitational virial.

Step 3. Insert the equations of motion and use the homogeneity of the $1/r$ potential (Euler's theorem, degree $-1$ ): $\sum_i m_i\,\mathbf{r}_i\cdot\ddot{\mathbf{r}}_i=-G\sum_{i<j}m_i m_j/r_{ij}=W$ . Why: the pairwise force terms recombine so the double sum is exactly the total potential energy $W$ — special to the inverse-square law.

Step 4. Time-average: for a bound system $I$ stays bounded, so $\langle\ddot I\rangle\to0$ , leaving $2\langle K\rangle+\langle W\rangle=0$ . Why: a bounded quantity has no persistent nonzero average second derivative, so the equilibrium condition falls out — and $E=-\langle K\rangle<0$ confirms the system is bound.

Limiting cases

The behaviour is governed almost entirely by $N$ and the ratio $t_{\mathrm{relax}}/t_{\mathrm{age}}$ :

Case	Result	Meaning
Two bodies	$N=2$ , integrable	Relative orbit is a Kepler conic; the ten classical integrals fully determine the motion.
Three bodies	$N=3$ , chaotic	No general closed-form solution; sensitive dependence on initial conditions (Poincaré); only special periodic families are known.
Collisional	$t_{\mathrm{relax}}\lesssim t_{\mathrm{age}}$	Globular clusters, galactic nuclei ( $N\sim10^5$ – $10^7$ ): encounters drive core collapse, mass segregation, binary hardening, evaporation.
Collisionless	$t_{\mathrm{relax}}\gg t_{\mathrm{age}}$	Galaxies, dark-matter halos ( $N\gtrsim10^{11}$ ): mean-field limit; use the collisionless Boltzmann equation and Jeans machinery.
Cold / unstable	$\sigma^2 \ll GM/R$	Velocity dispersion too small to support self-gravity; Jeans instability drives collapse and structure formation.

Common trapClose encounters must dominate relaxation because gravity is strong up close

The intuition is natural and wrong. Because the mean-square deflection summed over impact parameter goes as $\int \mathrm{d}b/b=\ln\Lambda$ , every decade of impact parameter contributes equally — the many distant weak encounters collectively outweigh the rare strong ones, and relaxation is a diffusive random walk, not a sequence of dramatic collisions. The same logarithm explains why a galaxy is collisionless: not because encounters are absent, but because $t_{\mathrm{relax}}\sim(0.1N/\ln N)\,t_{\mathrm{cross}}$ scales almost linearly with $N$ , reaching $\sim10^{16}$ yr at $N\sim10^{11}$ . A related trap: larger $N$ at fixed size makes the potential smoother and relaxation slower, not faster — which is why simulations soften the force to suppress spurious two-body scattering.

Worked examples

1 · Is a globular cluster collisional? (fully worked)

A typical globular cluster has $N\approx10^5$ stars, a half-mass radius $R\approx10$ pc, and a velocity dispersion $\sigma\approx10\ \mathrm{km/s}$ . Estimate its crossing and relaxation times, and decide whether two-body encounters matter over its $\sim12$ -Gyr life.

Step 1. Crossing time. With $R=10\ \mathrm{pc}=3.09\times10^{17}\ \mathrm{m}$ and $v=\sigma=10^4\ \mathrm{m/s}$ , $t_{\mathrm{cross}}=R/v=3.09\times10^{13}\ \mathrm{s}$ ; converting ( $1\ \mathrm{yr}=3.16\times10^7\ \mathrm{s}$ ) gives $t_{\mathrm{cross}}\approx\mathbf{1.0\ Myr}$ .

Step 2. Relaxation factor. With $\ln N=\ln(10^5)=11.5$ , the number of crossings to relax is $0.1N/\ln N=0.1\times10^5/11.5=\mathbf{8.7\times10^2\ crossings}$ .

Step 3. Relaxation time. Multiply: $t_{\mathrm{relax}}=(0.1N/\ln N)\,t_{\mathrm{cross}}=870\times1.0\ \mathrm{Myr}\approx\mathbf{0.85\ Gyr}$ .

Step 4. Verdict. Since $t_{\mathrm{relax}}\approx0.85\ \mathrm{Gyr}\ll12\ \mathrm{Gyr}$ , the cluster has relaxed many times over, driving mass segregation, core collapse, and evaporation — it is genuinely collisional. By contrast, a galaxy's $N\approx10^{11}$ and $t_{\mathrm{cross}}\approx10^8$ yr give $t_{\mathrm{relax}}\approx(0.1\times10^{11}/25.3)\times10^8\ \mathrm{yr}\approx\mathbf{4\times10^{16}\ yr}$ , so the very same physics leaves a galaxy collisionless.

2 · Weighing the Coma cluster with the virial theorem (fully worked)

Galaxies in the Coma cluster show a line-of-sight velocity dispersion $\sigma\approx1000\ \mathrm{km/s}$ over a radius $R\approx1.5\ \mathrm{Mpc}$ . Use $2K+W=0$ to estimate the cluster's dynamical mass, and compare with its luminous mass of a few $\times10^{13}\,M_\odot$ .

Step 1. Turn the virial theorem into a mass. For a bound system $2K+W=0$ with $K=\tfrac12 M\langle v^2\rangle$ and $W\approx-GM^2/R$ gives $M\langle v^2\rangle=GM^2/R$ , hence $M\sim\langle v^2\rangle R/G$ ; for an order-of-magnitude estimate we use $M\sim\sigma^2 R/G$ .

Step 2. Insert numbers. With $\sigma=10^6\ \mathrm{m/s}$ , $R=1.5\ \mathrm{Mpc}=4.63\times10^{22}\ \mathrm{m}$ , and $G=6.674\times10^{-11}$ : $M\sim\sigma^2 R/G=(10^6)^2\times4.63\times10^{22}/6.674\times10^{-11}=\mathbf{6.9\times10^{44}\ kg}$ .

Step 3. Convert. Dividing by $M_\odot=1.99\times10^{30}\ \mathrm{kg}$ gives $M\approx\mathbf{3.5\times10^{14}\,M_\odot}$ ; the structure factor ( $\sim3$ ) pushes the true virial mass toward $\sim10^{15}\,M_\odot$ .

Step 4. Interpret. The dynamical mass exceeds the luminous mass of a few $\times10^{13}\,M_\odot$ by more than an order of magnitude. The galaxies move far too fast to be bound by the visible matter — exactly the discrepancy Zwicky found in 1933, the first dynamical evidence for dark matter, resting entirely on the N-body virial theorem.

Physics of the portrait

The N-body problem is where celestial mechanics turned into chaos theory. Its portrait belongs to the mathematician who proved that the three-body problem cannot be solved the way the two-body problem can — and who, in failing to solve it, discovered sensitive dependence on initial conditions.

Origins · 1 minHenri Poincaré and how this was discovered

Henri Poincaré

1854–1912 · celestial mechanics and the birth of chaos theory

Poincaré was the last universalist of mathematics, and the N-body problem was his proving ground. In 1887 King Oscar II of Sweden offered a prize for a solution to the $N$ -body problem — a proof of the stability of the Solar System. Poincaré's memoir won in 1889, but while the proofs were printing the editor Phragmén found an error. Correcting it, Poincaré discovered something far stranger than the stability he had sought: in the restricted three-body problem the stable and unstable manifolds of a periodic orbit intersect in a homoclinic tangle of such intricacy that "I shall not even try to draw it." Neighbouring trajectories diverge exponentially — deterministic yet unpredictable. He had found deterministic chaos, seventy years before Lorenz rediscovered it in a weather model. Along the way he proved that beyond the ten classical integrals the three-body problem admits no further single-valued analytic integrals, so the closed-form solution that works for $N=2$ cannot exist for $N\geq3$ . The modern N-body problem — attacked by numerical integration and statistical mechanics precisely because analytic solution is impossible — is the discipline Poincaré's negative result created.

History

1687

— Newton states universal gravitation and solves the two-body problem in the Principia

1772

— Euler and Lagrange find the collinear and triangular equilibrium points of the restricted three-body problem

1889

— Poincaré wins King Oscar II's prize and, correcting an error, discovers homoclinic chaos in the three-body problem

1912

— Sundman gives a convergent series solution of the three-body problem that converges far too slowly to be useful

1943

— Chandrasekhar formulates two-body relaxation and dynamical friction for stellar systems

Going further

Texts. Binney & Tremaine, Galactic Dynamics, 2nd ed., Ch. 1, 3, 4, 7 · Carroll & Ostlie, An Introduction to Modern Astrophysics, 2nd ed., Ch. 2, 25 · Aarseth, Gravitational N-Body Simulations (2003) · Heggie & Hut, The Gravitational Million-Body Problem (2003) · Mo, van den Bosch & White, Galaxy Formation and Evolution, Ch. 5 · Dodelson, Modern Cosmology (structure formation). The key checkable results in this article are: (1) the impulse integral $\int_{-\infty}^{\infty}(b^2+x^2)^{-3/2}\mathrm{d}x=2/b^2$ giving $\delta v_\perp=2Gm/(bv)$ ; (2) the scaling $t_{\mathrm{relax}}\sim[N/(8\ln N)]\,t_{\mathrm{cross}}$ ; (3) the virial identity $\tfrac12\ddot I=2K+W$ ; (4) the worked numbers $t_{\mathrm{cross}}\approx1.0$ Myr, $t_{\mathrm{relax}}\approx0.85$ Gyr for a globular cluster, and $M\approx3.5\times10^{14}\,M_\odot$ for Coma from $\sigma^2R/G$ .

Prerequisites. Gravitational Two-Body Problem · Distribution Function and Phase Space. Leads to. Collisionless Boltzmann Equation · Two-Body Relaxation Time · Virial Theorem in Stellar Systems · Jeans Equations.

Ask the tutor: Ask me to run the relaxation-time or virial-mass numbers for any real system and I will show you where it lands on the collisional-to-collisionless line.

Why can we solve the two-body problem exactly but not the three-body problem?

Both have the same ten classical integrals — energy, total momentum, center-of-mass motion, and angular momentum — but integrating a flow requires enough isolating integrals to reduce it to quadratures, roughly $6N-1$ of them. For $N=2$ the ten suffice: they peel the twelve-dimensional phase space down to a single integrable radial equation, and the orbit is a conic. For $N=3$ the same ten leave too many coupled degrees of freedom, and Poincaré proved that no further single-valued analytic integral exists to close the gap. Worse, the phase space contains chaotic regions where trajectories diverge exponentially, so no closed-form solution can exist — not for lack of cleverness, but as a theorem.

If there is no solution, how do astronomers predict orbits at all?

By numerical integration. The equations of motion (1) are well-posed initial-value problems, so given positions and velocities now, a computer can step them forward as accurately as you like — solar-system ephemerides are computed this way over centuries. What Poincaré ruled out is a formula $\mathbf{r}_i(t)$ valid for all time, not step-by-step prediction. The catch is chaos: in chaotic regions small errors grow exponentially, so predictions have a finite horizon (a few tens of millions of years for the inner Solar System). For large $N$ , astronomers switch to statistical or fluid descriptions instead of tracking every orbit.

What is the Coulomb logarithm and why does it appear in a gravity problem?

The Coulomb logarithm $\ln\Lambda=\ln(b_{\max}/b_{\min})$ is the factor you get when you add up the deflections a star suffers from encounters at all impact parameters $b$ . Because each encounter's squared kick scales as $1/b^2$ and the number of encounters at that $b$ scales as $b\,\mathrm{d}b$ , the contributions go as $\mathrm{d}b/b$ , whose integral is a logarithm. It is called "Coulomb" because the same mathematics governs a charged particle scattering in a plasma. Its value $\Lambda\approx N$ is only logarithmically sensitive to the cutoffs, which is why order-of-magnitude estimates are robust.

What does "collisionless" mean if the stars never actually touch?

"Collision" here is a term of art: a two-body gravitational encounter strong enough to noticeably deflect a star's velocity, not a physical impact. Stars are so tiny compared with their separations that literal collisions essentially never happen. A system is collisionless when even these gravitational deflections are negligible over its lifetime — that is, when $t_{\mathrm{relax}}\gg t_{\mathrm{age}}$ , so each star moves as if in a perfectly smooth potential. Galaxies are collisionless in exactly this sense; globular clusters are not.

How does the virial theorem let you weigh something you cannot put on a scale?

The relation $2K+W=0$ says the kinetic energy (from the velocity dispersion $\sigma$ ) balances the gravitational potential energy (from the mass $M$ and size $R$ ). Rearranged, $M\sim\sigma^2 R/G$ — so if you measure how fast the members move (via Doppler shifts) and how big the system is (via angular size and distance), gravity itself hands you the mass, including any invisible matter. This is how Zwicky discovered dark matter in Coma in 1933 and how cluster masses are still estimated today, up to an order-unity structure factor that depends on the density profile and must be calibrated for precision work.

Why do N-body simulations "soften" gravity, and doesn't that change the physics?

Softening replaces $1/r^2$ with $1/(r^2+\epsilon^2)$ for a small length $\epsilon$ , capping the force at tiny separations. It is done for two reasons: to prevent numerical blow-ups when particles pass extremely close, and — more subtly — because each simulation particle represents an enormous number of real stars, so it should not scatter its neighbours the way a single point mass would. Softening deliberately suppresses spurious two-body relaxation to keep a collisionless system collisionless. For genuinely collisional systems like star clusters, where you want the encounters, you use un-softened direct integrators (NBODY6) instead. Choosing $\epsilon$ is thus a physics decision, not just a numerical convenience.