In physics, specifically classical mechanics, the three-body problem is to take the initial positions and velocities (or momenta) of three point masses orbiting each other in space and then to calculate their subsequent trajectories using Newton's laws of motion and Newton's law of universal gravitation. Unlike the two-body problem, the three-body problem has no general closed-form analytic solution. The differential equations that govern the motions of three gravitating bodies are not integrable and cannot be solved to give explicit formulas for the positions of the bodies as a function of time. For most initial conditions, the dynamical system for three orbiting bodies is chaotic, and the only way to predict their motions is to estimate them using numerical methods. The three-body problem is a special case of the n-body problem. Historically, the first specific three-body problem to receive extended study was the one involving the Earth, the Moon, and the Sun. In an extended modern sense, a three-body problem is any problem in classical mechanics or quantum mechanics that models the motion of three particles.
Mathematical description The mathematical statement of the three-body problem can be given in terms of the Newtonian equations of motion for vector positions r i = ( x i , y i , z i ) {\displaystyle \ \mathbf {r} _{i}=(x_{i},y_{i},z_{i})\ } of three gravitationally interacting bodies with masses m i {\displaystyle m_{i}} :
r ¨ 1 = − G m 2 ( r 1 − r 2 ) | r 1 − r 2 | 3 − G m 3 ( r 1 − r 3 ) | r 1 − r 3 | 3 , r ¨ 2 = − G m 3 ( r 2 − r 3 ) | r 2 − r 3 | 3 − G m 1 ( r 2 − r 1 ) | r 2 − r 1 | 3 , r ¨ 3 = − G m 1 ( r 3 − r 1 ) | r 3 − r 1 | 3 − G m 2 ( r 3 − r 2 ) | r 3 − r 2 | 3 . {\displaystyle {\begin{aligned}{\ddot {\mathbf {r} }}_{1}&=-Gm_{2}{\frac {\left(\mathbf {r} _{1}-\mathbf {r} _{2}\right)}{\ \left|\mathbf {r} _{1}-\mathbf {r} _{2}\right|^{3}}}-Gm_{3}{\frac {\left(\mathbf {r} _{1}-\mathbf {r} _{3}\right)}{\ \left|\mathbf {r} _{1}-\mathbf {r} _{3}\right|^{3}}}\ ,\\{\ddot {\mathbf {r} }}_{2}&=-Gm_{3}{\frac {\left(\mathbf {r} _{2}-\mathbf {r} _{3}\right)}{\ \left|\mathbf {r} _{2}-\mathbf {r} _{3}\right|^{3}}}-Gm_{1}{\frac {\left(\mathbf {r} _{2}-\mathbf {r} _{1}\right)}{\ \left|\mathbf {r} _{2}-\mathbf {r} _{1}\right|^{3}}}\ ,\\{\ddot {\mathbf {r} }}_{3}&=-Gm_{1}{\frac {\left(\mathbf {r} _{3}-\mathbf {r} _{1}\right)}{\ \left|\mathbf {r} _{3}-\mathbf {r} _{1}\right|^{3}}}-Gm_{2}{\frac {\left(\mathbf {r} _{3}-\mathbf {r} _{2}\right)}{\ \left|\mathbf {r} _{3}-\mathbf {r} _{2}\right|^{3}}}~.\end{aligned}}} where G {\displaystyle \ G\ } is the gravitational constant. As astronomer Juhan Frank describes, "These three second-order vector differential equations are equivalent to 18 first order scalar differential equations." As June Barrow-Green notes with regard to an alternative presentation, if P i {\displaystyle P_{i}} represent three particles with masses m i {\displaystyle m_{i}} , distances P i P j = r i j , {\displaystyle \ P_{i}P_{j}=r_{ij}\ ,} and coordinates q i j {\displaystyle \ q_{ij}\ } ( i , j = 1 , 2 , 3 ) {\displaystyle \ (i,j=1,2,3)\ } in an inertial coordinate system ... the problem is described by nine second-order differential equations. The problem can also be stated equivalently in the Hamiltonian formalism, in which case it is described by a set of 18 first-order differential equations, one for each component of the positions r i {\displaystyle \ \mathbf {r} _{i}\ } and momenta p i {\displaystyle \ \mathbf {p} _{i}\ } :
d r i d t = ∂ H ∂ p i , d p i d t = − ∂ H ∂ r i , {\displaystyle {\frac {\mathrm {d} \ \mathbf {r} _{i}}{\mathrm {d} \ t}}={\frac {\partial \ {\mathcal {H}}}{\partial \ \mathbf {p} _{i}}}\ ,\qquad {\frac {\mathrm {d} \ \mathbf {p} _{i}}{\mathrm {d} \ t}}=-{\frac {\partial \ {\mathcal {H}}}{\partial \ \mathbf {r} _{i}}}\ ,}
where H {\displaystyle {\mathcal {H}}} is the Hamiltonian:
H = − G m 1 m 2 | r 1 − r 2 | − G m 2 m 3 | r 3 − r 2 | − G m 3 m 1 | r 3 − r 1 | + | p 1 | 2 2 m 1 + | p 2 | 2 2 m 2 + | p 3 | 2 2 m 3 . {\displaystyle {\mathcal {H}}\ =\ -{\frac {Gm_{1}m_{2}}{\left|\mathbf {r} _{1}-\mathbf {r} _{2}\right|}}\ -\ {\frac {Gm_{2}m_{3}}{\left|\mathbf {r} _{3}-\mathbf {r} _{2}\right|}}\ -\ {\frac {Gm_{3}m_{1}}{\left|\mathbf {r} _{3}-\mathbf {r} _{1}\right|}}\ +\ {\frac {\left|\mathbf {p} _{1}\right|^{2}}{2m_{1}}}\ +\ {\frac {\left|\mathbf {p} _{2}\right|^{2}}{2m_{2}}}\ +\ {\frac {\left|\mathbf {p} _{3}\right|^{2}}{2m_{3}}}~.}
In this case, H {\displaystyle {\mathcal {H}}} is simply the total energy of the system, gravitational plus kinetic.
Restricted three-body problem
In the restricted three-body problem formulation, in the description of Barrow-Green,two... bodies revolve around their centre of mass in circular orbits under the influence of their mutual gravitational attraction, and... form a two-body system... [whose] motion is known. A third body (generally known as a planetoid), assumed massless with respect to the other two, moves in the plane defined by the two revolving bodies and, while being gravitationally influenced by them, exerts no influence of its own. Per Barrow-Green, "[t]he problem is then to ascertain the motion of the third body." The restricted three-body problem is easier to analyze theoretically than the full problem. It is of practical interest as well since it accurately describes many real-world problems, the most important example being the Earth–Moon–Sun system. For these reasons, it has occupied an important role in the historical development of the three-body problem. The restricted 3-body problem has a 4-dimensional phase space, but only one conserved quantity, the Jacobi integral. It was shown by Heinrich Bruns that there are no more algebraic conserved quantities, and by Henri Poincaré in 1889 that there are no more analytic conserved quantities. Therefore, since the dimension of the phase space is larger than the number of constants of motion, the system is not exactly solvable; in fact, it is chaotic. Depending on the value of the Jacobi integral, a body initially orbiting the larger mass may be able to be captured by the secondary mass or be ejected via Lagrange points L2 or L3. A variant of this problem, where the two large bodies both exert radiation pressure, results in the addition of four additional equilibrium points in addition to the five classical Lagrange points. More recent work on the circular restricted problem has drawn on symplectic geometry and contact geometry. A 2025 survey by Agustín Moreno describes the use of global surfaces of section and related geometric methods, and their connection with the design of spacecraft trajectories. Such methods also yield existence results: for energies slightly above the first critical value, Jungsoo Kang and Kevin Ruck used Rabinowitz Floer homology to show that the planar problem admits either a periodic symmetric collision orbit or infinitely many symmetric consecutive-collision orbits, and that for generic mass ratios and energies at least two geometrically distinct symmetric consecutive-collision orbits exist for each primary.
Solutions
General solution
There is no general closed-form solution to the three-body problem. In other words, it does not have a general solution that can be expressed in terms of a finite number of standard mathematical operations. Moreover, the motion of three bodies is generally non-repeating, except in special cases. However, in 1912 the Finnish mathematician Karl Fritiof Sundman proved that there exists an analytic solution to the three-body problem in the form of a Puiseux series, specifically a power series in terms of powers of t1/3. This series converges for all real t, except for initial conditions corresponding to zero angular momentum. In practice, the latter restriction is insignificant since initial conditions with zero angular momentum are rare, having Lebesgue measure zero. An important issue in proving this result is the fact that the radius of convergence for this series is determined by the distance to the nearest singularity. Therefore, it is necessary to study the possible singularities of the three-body problem. As is briefly discussed below, the only singularities in the three-body problem are binary collisions (collisions between two particles at an instant) and triple collisions (collisions between three particles at an instant). Collisions of any number are somewhat improbable, since it has been shown that they correspond to a set of initial conditions of measure zero. But there is no criterion known to be put on the initial state to avoid collisions for the corresponding solution. So Sundman's strategy consisted of the following steps:
Using an appropriate change of variables to continue analyzing the solution beyond the binary collision, in a process known as regularization. Proving that triple collisions only occur when the angular momentum L vanishes. By restricting the initial data to L ≠ 0, he removed all real singularities from the transformed equations for the three-body problem. Showing that if L ≠ 0, then not only can there be no triple collision, but the system is strictly bounded away from a triple collision. This implies, by Cauchy's existence theorem for differential equations, that there are no complex singularities in a strip (depending on the value of L) in the complex plane centered around the real axis (related to the Cauchy–Kovalevskaya theorem). Find a conformal transformation that maps this strip into the unit disc. For example, if s = t1/3 (the new variable after the regularization) and if |ln s| ≤ β, then this map is given by σ = e π s 2 β − 1 e π s 2 β + 1 . {\displaystyle \sigma ={\frac {e^{\frac {\pi s}{2\beta }}-1}{e^{\frac {\pi s}{2\beta }}+1}}.}
This finishes the proof of Sundman's theorem. The corresponding series converges extremely slowly. That is, obtaining a value of meaningful precision requires so many terms that this solution is of little practical use. Indeed, in 1930, David Beloriszky calculated that if Sundman's series were to be used for astronomical observations, then the computations would involve at least 108000000 terms.
Special-case solutions
In 1767, Leonhard Euler found three families of periodic solutions in which the three masses are collinear at each instant. In 1772, Lagrange found a family of solutions in which the three masses form an equilateral triangle at each instant. Together with Euler's collinear solutions, these solutions form the central configurations for the three-body problem. These solutions are valid for any mass ratios, and the masses move on Keplerian ellipses. These four families are the only known solutions for which there are explicit analytic formulas. In the special case of the circular restricted three-body problem, these solutions, viewed in a frame rotating with the primaries, become points called Lagrangian points and labeled L1, L2, L3, L4, and L5, with L4 and L5 being symmetric instances of Lagrange's solution. In work summarized in 1892–1899, Henri Poincaré established the existence of an infinite number of periodic solutions to the restricted three-body problem, together with techniques for continuing these solutions into the general three-body problem. In 1893, Meissel stated what is now called the Pythagorean three-body problem: three masses in the ratio 3:4:5 are placed at rest at the vertices of a 3:4:5 right triangle, with the heaviest body at the right angle and the lightest at the smaller acute angle. Burrau further investigated this problem in 1913. In 1967 Victor Szebehely and C. Frederick Peters established eventual escape of the lightest body for this problem using numerical integration, while at the same time finding a nearby periodic solution. In the 1970s, Michel Hénon and Roger A. Broucke each found a set of solutions that form part of the same family of solutions: the Broucke–Hénon–Hadjidemetriou family. In this family, the three objects all have the same mass and can exhibit both retrograde and direct forms. In some of Broucke's solutions, two of the bodies follow the same path. In 1993, physicist Cris Moore at the Santa Fe Institute found a zero angular momentum solution with three equal masses moving around a figure-eight shape. In 2000, mathematicians Alain Chenciner and Richard Montgomery proved its formal existence. The solution has been shown numerically to be stable for small perturbations of the mass and orbital parameters, which makes it possible for such orbits to be observed in the physical universe. But it has been argued that this is unlikely since the domain of stability is small. For instance, the probability of a binary–binary scattering event resulting in a figure-8 orbit has been estimated to be a small fraction of a percent. In 2013, physicists Milovan Šuvakov and Veljko Dmitrašinović at the Institute of Physics in Belgrade discovered 13 new families of solutions for the equal-mass zero-angular-momentum three-body problem. In 2015, physicist Ana Hudomal discovered 14 new families of solutions for the equal-mass zero-angular-momentum three-body problem. In 2017, researchers Xiaoming Li and Shijun Liao found 669 new periodic orbits of the equal-mass zero-angular-momentum three-body problem. This was followed in 2018 by an additional 1,223 new solutions for a zero-angular-momentum system of unequal masses. In 2018, Li and Liao reported 234 solutions to the unequal-mass "free-fall" three-body problem. The free-fall formulation starts with all three bodies at rest. Because of this, the masses in a free-fall configuration do not orbit in a closed "loop", but travel forward and backward along an open "track". In 2021, Xiaoming Li, Xiaochen Li and Shijun Liao reported a single family containing 135,445 periodic orbits of non-hierarchical unequal-mass triple systems, of which 13,315 were found to be linearly stable. In 2023, Ivan Hristov, Radoslava Hristova, Dmitrašinović, and Kiyotaka Tanikawa published a search for "periodic free-fall orbits" in the three-body problem, limited to the equal-mass case, and found 12,409 distinct solutions. Because a free-fall periodic orbit generally corresponds to two distinct sets of initial conditions, the 24,582 sets of initial conditions located by that search reduce to 12,409 distinct solutions, 236 of which are self-dual. A further search in 2025, restricted to free-fall orbits possessing central symmetry, enlarged this class; all of the orbits found in it were linearly unstable. In 2026, Hristov, Hristova and Tanikawa established four regions of stability for the equal-mass zero-angular-momentum problem and reported 971 verified sets of initial conditions for linearly stable collisionless periodic orbits. They describe these as candidates for orbits that are stable in the sense of the Kolmogorov–Arnold–Moser (KAM) theorem, rather than as orbits whose nonlinear stability has been established. Most of these large searches were restricted to planar motion. In a 2025 preprint, Li and Liao reported 10,059 new three-dimensional periodic orbits of the general three-body problem, computed for two unit masses together with a third mass taking the twenty values m3 = 0.1n for integers 1 ≤ n ≤ 20; 1,996 of them were found to be linearly stable. The set includes 21 spatial choreographies of three equal masses, in which all three bodies follow a single closed curve, and 273 orbits the authors call "piano trios", in which two equal masses share one closed curve while the third, of a different mass, moves along another. Whether such periodic orbits arise in nature has also been examined numerically. In simulations of binary–binary and triple–single encounters published in 2026, Simon Portegies Zwart, Arjen Doelman and Jelmer Sein found that about 9% of their calculations produced periodic three-body systems, which generally occurred as transients surviving for tens to hundreds of orbital periods rather than as permanent configurations.
Numerical approaches Using a computer, the problem may be solved to arbitrarily high precision using numerical integration. There have been attempts at creating computer programs that numerically solve the three-body problem (and by extension, the n-body problem) involving both electromagnetic and gravitational interactions, and incorporating modern theories of physics such as special relativity. In addition, using the theory of random walks, an approximate probability of different outcomes may be computed. Machine-learning methods have also been applied, mainly to the restricted problem and to the search for periodic orbits rather than as general-purpose integrators. A convolutional variational autoencoder has been used to represent families of periodic orbits of the circular restricted problem in a low-dimensional latent space, from which candidate trajectories are generated and then converged onto periodic orbits by numerical continuation. Physics-informed neural networks in which the boundary conditions are imposed as hard constraints have been used to compute Lyapunov orbits about the L1 and L2 points, distant retrograde orbits, and three-dimensional halo orbits.
Statistical approaches Because individual trajectories in the chaotic regime cannot be predicted over long times, a complementary approach seeks the statistical distribution of outcomes of a three-body interaction—in particular which body escapes, and the binding energy and angular momentum of the binary that is left behind. Building on earlier statistical escape theories, Barak Kol and collaborators developed a flux-based formulation in which outcome distributions are derived from a regularized phase-space volume. In 2024, Viraj Manwadkar, Alessandro Trani and Kol measured the associated "chaotic absorptivity" from millions of scattering events, following each only until it could be classified as regular or chaotic, and used it to predict the distribution of chaotic outcomes over both binary binding energy and angular momentum. The regularized phase-volume was extended from a function of energy alone to a joint distribution over energy and angular momentum by Yogesh Dandekar and Kol in 2025. Such treatments generally assume that the accessible phase space is explored ergodically. In 2024, Trani, Nathan Leigh, Tjarda Boekholt and Portegies Zwart reported that regular, non-chaotic trajectories occupy a substantial fraction of the phase space of the gravitational three-body problem—between 28% and 84% of the sampled initial conditions, depending on the configuration—and that regular and chaotic regions coexist at all scales with a multifractal structure. Statistical escape theories that exclude the outcomes of these regular trajectories may therefore give biased results.
History The gravitational problem of three bodies in its traditional sense dates in substance from 1687, when Isaac Newton published his Philosophiæ Naturalis Principia Mathematica. Newton, having solved the two-body problem, tried to discover whether any long-term stability is possible for a system such as the Earth, the Moon, and the Sun. Guided by major Renaissance astronomers Nicolaus Copernicus, Tycho Brahe and Johannes Kepler, Newton introduced later generations to the beginning of the gravitational three-body problem. In Proposition 66 of Book 1 of the Principia, and its 22 Corollaries, Newton took the first steps in the definition and study of the problem of the movements of three massive bodies subject to their mutually perturbing gravitational attractions. In Propositions 25 to 35 of Book 3, Newton also took the first steps in applying his results of Proposition 66 to the lunar theory, the motion of the Moon under the gravitational influence of Earth and the Sun. Later, this problem was also applied to other planets' interactions with the Earth and the Sun. The physical problem was first addressed by Amerigo Vespucci and subsequently by Galileo Galilei, as well as Simon Stevin, but they did not realize what they contributed. Though Galileo determined that the speed of fall of all bodies changes uniformly and in the same way, he did not apply it to planetary motions. Whereas in 1499, Vespucci used knowledge of the position of the Moon to determine his position in Brazil. It became of technical importance in the 1720s, as an accurate solution would apply to navigation, specifically for the determination of longitude at sea, solved in practice by John Harrison's invention of the marine chronometer. However, the accuracy of the lunar theory was low, due to the perturbing effect of the Sun and planets on the motion of the Moon around Earth. Jean le Rond d'Alembert and Alexis Clairaut, who developed a longstanding rivalry, both attempted to anal