Theoretical and Observational Constraints on Lunar Orbital Evolution in the Three-Body Earth-Moon-Sun System

Extremely slow recession of the Moon from the Earth has been recently proposed and attributed to conversion of Earth’s axial spin to lunar orbital momentum. This hypothesis is inconsistent with long-standing recognition that the Moon’s orbit involves three-body interactions. This and other short-comings, such as Earth’s spin loss being internal, are summarized here. Considering point-masses is justified by theory and observational data on other moons. We deduce that torque in the Earth-Moon-Sun system increases eccentricity of the lunar orbit but decreases its inclination over time. Consequently, the average lunar orbital radius is decreasing. We also show that lunar drift is too small to be constrained through lunar laser ranging measurements, mainly because atmospheric refraction corrections are comparatively large and variations in lunar cycles are under-sampled. Our findings support co-accretion and explain how orbits evolve in many-body point-mass systems.


Introduction
Kepler's elliptical orbits (the reduced two-body problem) reasonably approximates much of the solar system, but neither describes its evolutionary behavior, nor accurately portrays the orbital dynamics of the Moon, which only appears to orbit the Earth. Rather, both bodies co-operatively orbit the Sun (e.g., [1]). Gravitational interactions of three unequal masses, considered by Euler and Lagrange, should pertain. However, the Moon does not occupy any of the Lagrangian points, which are stable positions when two minor bodies move in a plane fixed in space about a very large body ( Figure 1). Hence, instability of the lunar orbit is of long standing and current interest (e.g., [2]).
Recent studies consider the lunar orbital radius to be very slowly growing. The tiny change, alleged to be +38 mm y −1 from modelling recent lunar laser ranging (LLR) measurements (e.g., [3][4][5][6]), is attributed to conversion of Earth's spin, which is decreasing, into lunar orbital angular momentum, which is perceived as increasing (e.g., [7]). This transfer hypothesis approximates the Earth and Moon as being an isolated two-body system [8], a case where no torque can exist [9], so angular momentum of the system should be conserved.
The hypothesis is poorly reasoned, as follows: • Neither body is isolated from the Sun [1], as illustrated in Figure 2.
• Torque is required to modify angular momentum [9], yet torque (τ = → r × → F ) on the Moon originating from the Earth is zero, because their vectors for separation distance (r) and their attractive gravitational force (F) are collinear.

•
No mechanism has been proposed for transforming axial spin of a body to the orbital angular momentum of a distant body. • Equations for the gravitational origin of orbits do not depend on the spin of the central body: indeed, a body need not even exist at the system barycenter, as exemplified by Pluto and Charon. , which reaches a maximum of~405,000 km from Earth (e.g., [10]). The actual excursion is slightly less than the maximum shown, due to the orbit being elliptical. The lunar monthly period is idealized as 1/12th of a year; (c) Orientation of Earth's spin and various orbital elements. The barycenter and ecliptic planes are distinct. Pink indicates components affected by the time-varying lunar orbit. The directions of the barycenter axis and the spin axis vary from~0 to 28.5 • over the lunar cycling, and precess independently. The radial position of the barycenter (B, pink double dots) varies by~500 km over a month inside the lower mantle. Its longitude varies by 360 • every day. Earth-Moon distance is not to scale, but body sizes and angles are. Part c was modified after [11], which has a Creative Commons 4 license.
Hence, revisiting stability and evolution of the lunar orbit is warranted. The purpose of this report is to provide a physically plausible cause of the current changes, which could provide constraints on the origin of the moon.
The paper is based on Newtonian physics. Section 2 describes the complex lunar orbit, its multiple cycles, and how orbital parameters, namely the semi-major axis (a), eccentricity (e), and inclination (i), are ascertained. Section 3 explains why the lunar orbit can be treated as gravitational interactions between point-masses. Section 4 discusses reduction of body spin through internal dissipation. Section 5 sets limits on changes in lunar a, e, and i using conservation laws, and accesses the possible changes based on force imbalances and torques experienced by the Moon during the three-body point-mass interactions. Section 6 covers the problems in LLR models which prevent quantifying lunar orbital evolution. Section 7 discusses lunar evolution and origin, and compares our findings to behavior of the planets, which are also perturbed by three-body interactions [12]. Section 8 concludes.

Observations and Parameterization of the Unique Lunar Orbit and its Cycles
In this work, "lunar orbit" refers to the Moon encircling the Earth, whereas "sinuous lunar path" is applied to the actual motions of the Moon about the Sun in the heliocentric frame of reference. The lunar orbit is around the barycenter of the Earth-Moon system. The enormous Solar mass allows us to treat the Sun as stationary in our theoretical assessment, and so "barycenter" as used in this report refers to the barycenter of the Earth-Moon system, which orbits the Sun in a nearly circular path.

Long-Standing Techniques for Measuring Motions in the Solar System
Motions of the Moon, Sun, and planets have been reliably recorded from the Earth against the stars for~200 years. Only motions perpendicular to the line-of-sight are wellresolved. Consequently, mathematical modeling and multicomponent fitting are used to describe the apparent motions of astronomical bodies against the background of stars and to parameterize the orbital motions [13][14][15]. Although the Earth-Moon barycenter moves in an elliptical orbit around the Sun, the position of Earth's center in space is controlled by the varying Earth-Moon distance: Thus, the geoenter does not move around the Sun in a simple elliptical orbit ( Figure 1). Neither are geocenter cycles annual because the Julian year is not an integer multiple of any lunar monthly cycle (see below). Due to this complexity, observations are based on Earth crossing the ecliptic in the spring (the vernal equinox) which defines the ecliptic plane for that particular year, and also an imaginary grid of longitude and latitude against the stars.
Projecting motions forward in time requires revising model parameters. Newcomb's [13] meticulous calculations were used for many decades. Such models describe observations of the Moon from the Earth as polynomial expansions with five to nine terms [14,15], and similarly depict motions of the planets. From these models, which are affected by the Earth's complex motions in a time-varying plane, orbital parameters are extracted. Trade-offs exist between providing accuracy for the present century and a large time span of applicability [15].

Lunar Orbital Parameters and Velocities
Available data show that the Moon's orbit around the barycenter, on average, resembles a Keplerian ellipse. Table 1 provides orbital and body parameters. Tangential velocities are~1 km s −1 for the Moon's apparent path around the Earth, or 30 km s −1 for the barycenter around the Sun (e.g., [10]). To "keep pace with the Earth" the Moon moves faster when it is outboard of the barycenter, than inboard, but because the barycenter orbit is curved, the moon spends more time outboard (Figure 2).

Lunar Cycles: "Year" vs. "Month"
To observers on Earth, phases of the Moon and eclipses define key lunar cycles. Phases of the Moon repeat on average every 29.53 days, providing the well-known synodic month. Full moons occur at the extreme outboard position with respect to the Sun and the barycenter orbit, whereas new moons occur at the extreme inboard position (Figure 2). Quarter moons roughly occur at medial times along the barycenter path, due to its curvature. Tilt of the lunar orbital plane with respect to the barycenter plane ( Figure 1) and its varying orientation make eclipses infrequent, rather than monthly [10, [20][21][22].
Other definitions exist for a lunar month [10]. The important cycle as regards orbital evolution is the anomalistic period of 27.55 days, which represents the interval from perigee to perigee (or apogee to apogee).
Due to inclination, eccentricity, the variable orientation of the lunar orbital plane with respect to the fixed reference frame of the barycenter plane, and the anomalistic cycle not being some simple fraction of a year (i.e., the barycenter orbital period of 365.25 days), the orbit of the Moon around the Sun involves several long cycles. Many are on human time-scales [10, [20][21][22]. Apsidal precession takes 8.85 years ( Figure 2). The~19 year Metonic cycle links the synodic month to the tropical year, but the Callippic cycle of 76 years is more exact. The Saros cycle of 18.029 years concerns eclipses, and thus a special repeated orientation of the lunar and barycenter planes.

Motivation to Directly Measure the Moon's Orbital Radius
Various scenarios for formation suggest slow growth of the average radius of the lunar orbit. For example, if the Moon was ejected from the Earth or fission occurred [23], its current average orbital radius suggests growth at~90 mm y −1 assuming zero radius at solar system beginnings. This upper limit for recession is 12 or 13 orders-of-magnitude smaller than orbital velocities (Section 2.2), and is far too small to quantify based on observations against the stars. Alternatively, as proposed by Laplace, the Moon could have formed from a ring of dust around the Earth (see, e.g., [23,24] and references therein). Co-accretion points to an originally circular orbit at some finite radius, permitting negative values for lunar drift. Section 5 provides details.
Crucially, annual resetting of the geocentric reference frame at the vernal equinox greatly complicates, if not prohibits, quantification of multi-year changes. Thus, to measure possible tiny perturbations of the Moon's orbit, laser light is reflected from retro-reflectors (denoted "mirrors" for brevity) placed on the Moon by astronauts (e.g., [3]). See Section 6 for details.

Gravitational Attraction to Oblate Planets and to the Pluto-Charon Binary
Two facts pertain to the motion of a small mass in the gravitational potential well of large body. First, at some distance from any mass distribution, the gravitational potential must reduce to the inverse radial dependence associated with a point mass. Second, the radial force of gravity, combined with forces during axial spin, renders large self-gravitating bodies oblate: this reduction in symmetry makes the force field non-central [25]. This section discusses how gravitational attraction to the oblate and lower symmetry distribution affect natural satellites in our solar system.

Gravitational Attraction to Non-Spherical Mass Distributions
The gravitational potential exterior to an oblate body, first derived by Gauss and Dirichlet, is complex even for a body of homogenous density. Hofmeister et al. [25] provide simpler, but exact, equations along the special axes of an axisymmetric body for the potential (ϕ) and force (F) as functions of the radial (r) and axial (z) directions: where m is the mass of an exterior test body, and M is the mass of the central body with equatorial radius α and ellipticity ε. MacMillan [26] provides a formula equivalent to Equation (1). Figure 3a shows that the axial, equatorial and spherical potentials converge rapidly with increasing distance even if the oblate body is very flat. Thus, beyond a few body radii, the attraction is close to that of a spherical distribution. Convergence becomes more rapid as ellipticity decreases. For the general point exterior to the body: where: The term in square brackets in Equation (5) is dimensionless and reveals the noncentral nature of gravitational attraction to an oblate body [25]. Specifically, at each and every external point, the potential is not inversely proportional to the spherical radial distance s = (r 2 + z 2 ) 1 2 . Moreover, at a general point the force is not even directed towards the object's center (Figure 3b: see [25] for additional examples).
Regarding the special directions, symmetry requires that force vectors along the polar axis and above the equatorial plane do point towards the center, but again force magnitudes neither vary as 1/r 2 nor as 1/s 2 per Equations (1)-(6).

Comparison with the Generalized Potential Used in Fitting
The "generalized" potential describes attraction of a test mass to a large body as a perturbation on results for a sphere [27]: α ave s n P k n (sin θ)(C nk cos kλ + S nk sin kλ) where α ave is mean equatorial body radius; s, θ, and λ are the spherical orbital radius, latitude, and longitude of the test mass; P n k are associated Legendre functions of the first kind; and C nk and S nk are numerical coefficients obtained from fitting. As a potential for a planet must necessarily include their predominantly oblate shape, limitations of the generalized potential are understood by comparing Equation (7) to the exact results for an oblate body, Equations (1)- (6).
For an axially symmetric body, C nk and S nk terms involve longitude and k = 0, where the zonal harmonic coefficients (J n , defined as −C n0 ) are used in fitting, e.g., of motions of satellites around a gas giant. Equation (7) can only be reconciled with results for the special axes, Equations (1) and (2), if and only if the relationship for each term: is used [25]. Furthermore, this equivalence also presumes that many terms in each series are used to achieve convergence, which is atypical of fitting procedures using Equation (7), see, e.g., [28]. Furthermore, the series describing points off the special axes of the oblate body is incompatible with the generalized potential. The meaning of the coefficients is lost in fitting, because attraction to an oblate body is not central as in Equation (7), as illustrated in Figure 4. In short, the properties of Equation (7) are incompatible with those of oblate bodies. The theoretical and observed shape of spinning, self-gravitating, astronomical entities is the oblate spheroid: However, the generalized potential cannot accurately describe attraction to bodies such as gas giants and stars.
Here, we note that the generalized potential with many terms can be used to fit orbital data for satellites proximal to Earth. As discussed by Transtrum et al. [29], multiplied parameters permit fitting data with incorrect physics. The fitting approach necessitates that the empirical fits are adjusted as orbits evolve, and is useful for satellites that typically orbit at heights of only 1.02 to 1.3 Earth radii, and experience additional short-term accelerations from regional asymmetries such as continents with their massive mountain ranges. However, the physical meaning of the most important terms, J 2 in particular, is lost in applying Equation (7) to the giant planets [25] which are distinctly oblate, rather than approximately spherical with bumps.

Orbits involving A Non-Central Potential
Far from a spinning astronomical body, where the point mass approximation is valid, orbits can be inclined, eccentric, and follow Keplerian dynamics. Close in, where forces are not central, stable orbits are either equatorial and circular or polar and elliptical. Bodies close-in with inclined, but non-polar, orbits are unstable, and will be pulled toward the equatorial plane.
Orbital parameters are controlled by distance, plus central body size and ellipticity [25]. For bodies with close-in, low inclination orbits, equating the centripetal and gravitational forces yields: where a is the radius and v the tangential speed of a circular orbit. Additionally, µ = 0 for a constant density body, increasing to 3 /5 for the limiting case of a point or spherical mass where ε = 0 [25]. For the above reasons, Section 3.2 analyzes data on natural satellites in terms of the equatorial radius (α) of the central planet, and discusses satellite systems individually, to account for ε varying among the gas giants, and to address unique Pluto's system.

Examples of Orbits around Uniaxial Central Mass Distributions
The giant planets are significantly flattened, with ellipticities of 0.432 for Saturn, 0.354 for Jupiter, 0.213 for Uranus, and 0.184 for Neptune [16]. The smallest axial ratio (for Saturn) is 0.902, which is far more spherical than the example in Figure 3. In contrast, rocky Earth and Mars have ε = 0.00335 and 0.00589, respectively, whereas tiny Pluto, which is spin-orbit coupled with nearly as large Charon, has ε = 0.0000.
Orbital characteristics of the moons of the outer planets and Pluto occupy three distinct trends with a few outliers (Figure 4a). Our Moon falls on a trend with 20 others. The selfgravitating close moons almost all have nearly circular orbits (e.g., eccentricity e < 0.018 for the 8 close Jovian moons) that lie in the equatorial planes of these planets (e.g., inclination i < 0.47 • for the Jovian moons). More distant satellites have large i and e, with both covering a wide range, and include many (e.g.,~60 around Jupiter) in retrograde orbits ( Figure 4a). Overall, inclination correlates with eccentricity for the moons of all planets, and those with many moons show fits with similar slopes (Figure 4b). Figure 5 shows how e and i depend on the semi-major axis, normalized to central body radius. Orbits of distant moons have high eccentricities and widely ranging inclinations ( Figure 5a). Inclinations are larger for the close moons of Saturn, reaching 1.86 • , but the same bimodal pattern exists for e and i. Moons of Uranus follow similar trends to the gas giant systems, whereas the few moons of Neptune are distributed bimodally. Mars is quite round, like Earth, and its tiny moons have orbital characteristics similar to the close moons of the giants (Figures 4 and 5). The Plutonian system is discussed below. Pluto's tiny moons [30], which actually orbit the Pluto-Charon barycenter are not shown: see text for discussion of this unique satellite system.
Iapetus of Saturn resembles our Moon in all three orbital characteristics (Figures 4 and 5a). Its average radius (~1420 km) is also similar, and Iapetus is fairly round, but its composition provides a lower density, in accord with most (icy) moons in the outer solar system. The inclination of Iapetus is also variable like our Moon (not shown).

A Triaxially Shaped and Time-Varying Central Mass Distribution
Pluto and Charon essentially compose a binary system, where the remaining 4 tiny satellites are spin-locked, resonating, and co-planar [30], even at large a/α of 60 ( Figure 5a). These 4 tiny distant moons are forced into equatorial orbits because Pluto and Charon together compose the central mass: Their combined gravitational potential is highly flattened in their barycenter plane and time-dependent with a repeat near 6 days. As occurs for a single oblate body, the direction of pull is towards their barycenter only in their orbital plane and along the special, perpendicular, axis above their barycenter. Nowhere in this fascinating system is the gravitational attraction proportional to the distance to the Pluto-Charon barycenter. On average and roughly, Pluto and Charon together act as a highly flattened oblate body over their~6-day orbital period.
Hence, the relevant parameter is not Pluto's α, used in Figure 5, but Charon's orbital radius. For this extreme case, the uniaxial force operates out to a factor of 3.3 × larger than the apparent limit of a/α = 27 seen for the oblate giants. This observation, in view of Earth's low ellipticity, supports modelling the Moon-Earth binary as point mass interactions. The remainder of Section 3 concerns the slightly oblate planets.

Ellipticity and Spin Rate
Todhunter [31] devised the formula for the dependence of angular velocity (ω) of a spinning, homogeneous body of homogeneous density (ρ) on ellipticity: based on the geometrical constructions of MacLaurin. This classical result for homogeneous uniaxial bodies simplifies to a linear relationship between ω and ε for ε < 0.5: [32], which therefore holds to significantly larger ε than the planets of our solar system. The small deviation between the exact relationship and the approximation (Figure 6a) is insignificant, given that real bodies have radially varying density. As body spin declines, the object rounds up (Figure 6a). Distant and close in orbits are affected differently from theoretical considerations, which are consistent with the observations of moons in the solar system, as follows:

Effect of Slowing Planet Spin on nearby Satellites
For moons with close-in orbits, which occur near a planet's equatorial plane (Figures 4 and 5), spin down of the planet will effectively weaken the exterior attractive force per Equations (2), (4), and (9). As the attractive force changes not only in magnitude, but also in direction, in response to spin down, orbital evolution is complicated and should depends both on parameters of the planet and the satellite's initial orbital parameters.

Effect of Slowing Planet Spin on Distant Satellites
Moons in distant orbits (a > 27α) about individual planets will not be affected by spin decline because the central potential acts as that of a point mass at such distances (Figures 4 and 5). Considering Figure 3a and results for round oblates [25], convergence occurs at a small, roughly~0.05%, difference in the actual potential from that of a sphere. Our Moon is described by this situation because its orbital radius is~60 Earth radii. As Earth's tidal bulge is miniscule compared to the difference between its polar and equatorial radii, this bulge negligibly affects the Moon's motions. Thus, the lunar orbit is a point mass problem.

Internal Dissipation of Body Spin and Implications for Conservation Laws
Loss of spin makes planetary bodies less oblate, per Equation (10ab). The rocky Earth is a differentially spinning, mostly solid, layered body [34,35]. Implications of friction existing between the layers are as follows:

Oversimplified Case: Conservation of Spin Angular Momentum Holds
Earth's~60 km thick outermost layer (the solid lithosphere), moves westward at about 60 mm y −1 relative to the mantle and below (~6300 km radius, sketched in Figure 6b). From a classical physics standpoint, friction exists between these layers and degrades the spin while heating their interface, but (angular) momentum is alleged to be conserved nevertheless. Given that the lithosphere is only 0.6 wt% of the Earth [36], recent slowdown of the surface at 2 × 10 −6 s y −1 [37,38] means that if spin angular momentum were indeed conserved, then the interior would accelerate by a mere~10 −9 s y −1 .
The lithosphere lies on a low velocity zone, detected in seismology, and inferred to be weak (e.g., [36]). The base of the lithosphere is partially melted, which reduces resistance to the shearing motion of differential rotation. However, Earth's subducting tectonic plates extend to depths of 300 to 660 km [34,35] and thereby provide resistance: Some plates break while others bend, depending on whether the subduction direction shares the same sense as spin. This geometry shows that differential rotation exists between the outermost layer and the interior [34,35].
Another prominent weak zone exists: Earth's outer core does not propagate shear waves and is therefore at least partially liquid, whereas the inner core is solid due to compression [36]. Therefore, differential rotation could exist between the strong mantle and the inner core. Super-rotation of the core has been considered, but its small size (<1 • y −1 ) makes quantification difficult (see [32] and the review of [39]). A tangential velocity similar to continental drift is well below the limit of detection. Notably, super-rotation of the whole core differs from super-rotation of the inner core. In particular, the molten outer core flows (shears) without resistance, so a faster rotating, more oblate solid inner core could be accommodated while the mantle rotates more slowly and is rounder. Criss [32] provides analytical formulae for spinning stratified bodies along with calculations for several endmember cases.
The roughly horizontal frictional force of viscous drag provides the torque needed to change relative angular momentum of the layers (Figure 6). Conservation of momentum would cause spin to decrease towards the surface, while increasing towards the interior, thereby preserving the original angular momentum of the whole Earth.
Thus, transfer of spin angular momentum to the Moon is unnecessary because internal processes redistribute spin during conservative behavior.

Frictional Forces Decrease Spin
Previous discussions of friction and momentum conservation implicitly assume that the force applied is external. Regarding the textbook example of a block sliding down a hill, the two objects are considered to be distinct and the force is a constant, because it is visualized as stemming from the nearly constant gravitational acceleration on Earth's surface. Earth acts as the external object (surroundings) in the sliding block problem. In contrast, the internal layers of a spinning planet compose the object itself.
Spin energy declines due to friction between internal planetary layers, which produces heat. Loss of energy depends on material properties such as the coefficient of friction and thermal conductivity. The amount may be large, as the model of Na and Lee [40] suggests that 97% of energy loss during spin down becomes heat. Thus, Earth's average spin rate is declining.
Earth's spin angular momentum decreases during heat production because friction provides torque ( Figure 6). This decline is at a slower rate than spin energy decreases (Iω vs. 1 /2 Iω 2 , where I is the moment of inertia). For additional discussion of large-scale processes occur inside the Earth, see [33]. Microscopic processes that dissipate momentum must be intimately tied with those converting mechanical energy to heat, but a detailed discussion is beyond the scope of this report.
Thus, dissipation of spin is also an internal process, and does not affect the Moon.

Permissible Changes in Orbits around a Central Point Mass
The above discussion of Earth's declining spin shows that changes in the interior rotational motions depend on whether or not heat is produced, i.e., whether mechanical energy is identical to the total energy, which is the quantity actually conserved. Similarities and differences of spin with orbital changes are covered here. This analysis extends our earlier analysis of the weak perturbations of planetary orbits by other planets [12].

Conservation Laws for Orbits Are Stringent
The Virial theorem of Clausius, predicated on conservative forces, results from linear momentum cancelling in all directions over the restricted space that defines some bound state [41]. This balance is readily demonstrated for elliptical orbits. Balance involves the canonical conjugate pairs of position and linear momentum, and is independent of conservation of energy and angular momentum. The extra stipulation is attributable to energy being consumed to form the bound state, which by definition is restricted in space. Orbits are cyclical bound states. The Virial theorem involving Newtonian gravitational potential energy (P.E.) tightly restricts kinetic energy (K.E.) of the orbit: where the angle brackets indicate averaging over the orbit. Conservation of energy provides an additional constraint: K.E. + P.E. = E tot = a constant (12) which holds at any given time. Multiple constraints make bound states resistant to change. Equations (11) and (12) hold for repeating paths, such as aperiodic orbits [42]. Thus, the Virial constrains hold for the Moon's path which is elliptical over a long term average.

Energy Conservation Is Key
Constant total energy provides a stringent restriction. Regarding orbits, friction does not exist for an object moving in the rarefied medium of space. Although collisions could alter the energy of an orbiting object, these have been few and small in the solar system over the past several billion years [43], so are not considered here.

Changes in the Lunar Orbit Permitted under Energy Conservation
Constant E tot of an orbiting object does not permit changes in its semi-major axis, but allows eccentricity to vary [42]. Importantly, for a 2-body system, any inclination is permissible, as this factor does not enter into the description of the orbit around a point mass or sphere. Thus, orbital inclination may also change. Figure 7a compares the current, quite circular, average lunar orbit to hypothetical orbits in the same plane that have the same energy, with the Earth at one focus. The radial distance as a function of orbital position is [42] (p. 131): where the angle θ is defined by a point on the orbit and the focus, with respect to a the line connecting the perigee and apogee. The angular average of r in Figure 7b was computed by integrating Equation (13): Equation (14) provides the geometrical average of an ellipse, and rests on infinitesimally small angles. An angular average is appropriate because a is constant under energy conservation and so evolution involves changes in shape (i.e., in e or b) over long periods of time.
Elongation drives the Moon towards Earth at perigee, but away at apogee, while shrinking the semi-minor axis (b). Consequently, the average radius decreases during elongation. If the endmember condition of a perfectly circular orbit describes formation 4.52 billion years ago, then the rate of lunar drift is −0.13 mm y −1 . A circular orbit must describe the starting point because torque is elongating the orbit, as follows:

Solar Torques Change Eccentriciy and Inclination of the Lunar Orbit
If the Earth and Moon were truly an isolated 2-body system, the lunar orbit could not change. Changes result from interactions with the Sun (the actual 3-body system), which provides a force on the Moon that is 2.2× greater than does Earth:  (15), lists their numerical values, and describes our notation for subscripts.
The Sun holds the Earth-Moon pair in a heliocentric orbit. The barycenter orbit is nearly circular, indicating long-term stability (e.g., [12]). With postulated lunar drift being~12 ordersof-magnitude smaller than barycenter or lunar orbital velocities, force imbalances must also be small relative to those largely controlling the orbit, i.e., those in Equation (15). Thus, the Sun's influence on the Moon is deduced by considering solar forces as perturbing the lunar orbit around the Earth. The planets slightly perturb the lunar orbit, but since the same conservation principles hold and the largest force on the Moon (from Jupiter at its closest pass) is only 0.00017 times solar [12], this minor effect can be neglected.
Essential behaviors, such as the shape of spinning planets and orbital characteristics, have been inferred by comparing an attractive gravitational pull to opposing centrifugal forces [42]. Following this approach does not idealize the orbits of these 3-bodies as coplanar, but for simplicity we consider both the barycenter and lunar orbits to be circular. Hence, gravitational pull (Equation (15)) and radial centripetal acceleration of the barycenter are balanced: F G,SB = F R,SB , and similarly for a nearly stable lunar orbit: F G,Em ∼ = F R,Em (Figure 8). In contrast, F G,Sm and F R,Sm are imbalanced even in this idealization.  (Figures 2 and 9). Rarely is the force applied by the Sun to the Moon collinear with the Sun-Moon direction. The geometry of the 3-body system thus provides torque which necessarily changes angular momentum, yet lunar orbital energy must remain constant. Thus, the imbalanced solar force and consequent torque can affect only e and i of the lunar orbit. Figure 9. Characteristics of the lunar orbit over the perigee-to-apogee half-cycle: (a) Distances associated with the lunar and terrestrial Keplerian orbits around the barycenter. Red dots = lunar radii spaced at equal times of 0.54 days, which were calculated from Equation (13) using 0.2 • increments, and a and e in Table 1. These points are reproduced by a 4th order polynomial fit (not shown). The average radius (black dashes) is slightly lower than the semi-major axis (grey heavy line). Blue X and right axis = distance to the geocenter, calculated from the fit; (b) Time derivatives of the calculated distances. Red solid curve = lunar variations over the half-cycle. Blue dots near x axis = the barycenter variation on the same scale. Long dashes and right axis = lunar drift, as modelled from LLR delay times [3], scaled to a daily basis.

Permissible Changes in Eccentricity
In the plane of the lunar orbit, the Sun can accelerate the Moon both radially (perpendicular) and along (parallel) to the orbit (Figure 2b). Centripetal acceleration of the Moon by the Sun is evident in the variation in velocity along the month. From the NASA factsheets [16], the maximum orbital velocity of the moon around the Earth is 1.082 km s −1 while the minimum is 0.970 km s −1 . Parallel acceleration should not change eccentricity. Orbital precession, i.e., rotation of the orbital ellipse as the barycenter moves around the Sun (Figure 2d), probably links to parallel acceleration.
As the Sun pulls the Moon towards the Sun in an inboard position (Figure 2ab), this elongates the orbit, and is denoted as grab. In an outboard position, the Sun pulls the Moon towards the Earth (and Sun) which brings the Moon closer to the lunar orbit foci, which also elongates the orbit (Figure 2). Thus, perpendicular acceleration associated with the Sun increases e of the lunar orbit, while decreasing the average radius (−0.13 mm y −1 from Equation (14) and Figure 7b).

Permissible Changes in Inclination
The Moon lies~5 • out of the barycenter orbital plane. Hence, the direction of pull on the Moon towards the Sun is infrequently co-linear with the centrifugal force associated with the Moon's path around the Sun (Figure 8a,b). Stability of the lunar orbital plane is impossible because pull towards the barycenter plane provides torque at all times except during the instant of the~bimonthly crossings (Figure 2b). Hence, inclination of the Moon's orbital plane changes unilaterally, because it is slowly but continuously pulled into the barycenter plane.
Inclination of the lunar plane to the ecliptic is currently 5.14 • , while varying from 18 to 28 o with respect of Earth's equatorial plane (Figure 1c). The latter is relevant (Section 3.2). Assuming that an equatorial configuration existed at formation suggests a progression of~5 × 10 −9 • y −1 . Jupiter also applies torque, but to both Earth and its Moon, thereby moving the barycenter orbital plane (Section 7).

Uncertainties in Modelling Lunar Drift from LLR
Difficulties in establishing slow evolutionary changes in the lunar orbit are linked to large and rapid variations in the orbit of the Moon around the barycenter, and the many correction terms needed to convert surface-to-surface measurements to orbital parameters. This section begins by describing known lunar motions, then discusses the database [3] provided by laser lunar ranging (LLR), and how drift is ascertained by combining several independent models.

Daily Variations
In a perfect elliptical orbit, the orbital radius changes throughout each cycle. Given the average values in Table 1, the idealized lunar orbital radius around the barycenter increases on average by~3000 km d −1 from perigee to apogee, and then shrinks in the second half of the month. Simultaneously, the distance from the barycenter to the geocenter changes proportionately: From Table 1, the associated change of ±500 km inside the Earth is a factor of~83 smaller than the lunar variation over a month (Figure 9).
Lunar drift modelled from LLR [3-6,44] as 38 mm y −1 is nearly 10-orders of magnitude smaller than the average daily change in the lunar radius over a half-month ( Figure 9). Accuracy in radius better than a parts-per-billion level is required to established drift. This level of precision is not achievable, as follows:

Variations in Tangential Velocity during the Lunar Orbit
Because tangential velocity varies continuously with position, intermittent measurements of the orbital radius are insufficient to establish lunar drift. Evenly spaced time intervals of 0.54 days were used to construct Figure 9 mainly because LLR measurements can, at best, be made over the 1 2 day period when a given LLR station spins across the side of Earth facing the Moon. Additionally, most years in the LLR data base have~730 data points or fewer, so intervals between measurements are, on average,~1/2 day or longer.
Fixed angular intervals provide additional information. Histograms of radius for the average lunar ellipse constructed over an anomalous month at narrow (0.5 • ) intervals provide a smooth distribution (Figure 10a). A wider increment of 5 • returns the expected mean radius, but the distribution has some irregularities. With a 13 • increment, which is similar to~1 day intervals common in the LLR database, the mean is not reproduced, which shows that large increments are insufficient to provide the mean radial distance. Laser ranging measurements are made at irregular intervals. Our results for wide intervals show that this approach is unlikely to constrain the average radius of an elliptical lunar orbit.

Lunar Laser Ranging Data and How Drift Is Ascertained
Muller et al. [6] review Earth-based LLR measurements, upon which the 38 mm y −1 recession value is based. Satellite experiments began afterwards and so are not covered here. Models used to extract drift are discussed after we describe LLR raw data and compare these to distances ascertained from ephemeris tables [19,45,46].

Available Station-to-Mirror Travel Times
For details on LLR measurements, see Samain et al. [47]. Battat et al. [48] discuss the Apache point station in particular. Murphy [49] provide additional discussion and include a map of retroreflector sites on the moon, established by astronauts. For brevity, "mirror" is used here.
Data collected from various terrestrial stations from late 1969 to early 2019 are compiled on the LLR website [3]. Data from the two longest running stations resemble the aggregate data, which is presented below. The Apollo 15 site was used in 70% of the measurements [44]. About 40% of the data were collected near quarter moons [44]. Since 1984, 58% of data since 1984 are from Grasse station [6], which contributes to~1/2 day gaps frequently existing in the database.
Round-trip travel times were converted to distances by dividing by 2 and multiplying by lightspeed. The figures below show these raw data. Travel times are reported at the level of 10 −13 s, which corresponds to an uncertainty of ±0.03 mm for each uncorrected station-to-mirror distance.

LLR Data Collection Is Too Infrequent to Describe Monthly and Longer Lunar Cycles
Most years are under-sampled: Fewer than 2 measurements per day, on average, have been made over the program lifetime ( Figure 11). Sampling may be sufficient to provide reasonable data for 1985, 1995 through 2001, and 2015 through 2018. However, this presumes reasonable coverage of the monthly cycles throughout any given year, which is not the case (Figure 12). An additional problem exists because the anomalous monthly cycle is not related to the Julian year by an integer, see Figure 12. Considering an apsidal precession cycle (8.85 years) might improve statistics. However, dense sampling only occurred from~1990 to 2001 and again after 2011, but the post 2011 interval is less than 8.85 years ( Figure 11). Additionally, the complete apsidal cycle under-samples the closer distances (cf. average in Figures 12b and 13). This option for ascertaining changes in radius with time is not viable. This problem is compounded by the anomalous monthly cycle not being related to the Julian year by an integer, see Figure 12.   [45,46] less the combined body radii of the Moon and Earth (fine black line): (a) Example of data collection at an average interval of~1 day. Purple vertical bar shows the 8109 km body radii contribution, which is substantial; (b) Example of data collected at an average interval near 1 2 day. Blue curve shows the monthly average distance. For both panels,~90% of the data were collected at Grasse [7,44]. Note that distances longer than expected occur when the station is not collinear with the line containing the geocenter, barycenter, and Moon's center. Distance varying greatly over any given anomalous cycle (Section 6.1) exacerbates the sampling problem. As shown in Figure 11, the monthly average distance also cycles over circa one year. Capturing this interesting variation, which may be connected with 3-body dynamics, requires adequate sampling of both the perigee and apogee, which is not available (Figure 13), perhaps due to preferential data collection at the quarter moons. Distances are neither constrained over the anomalistic month, nor over the anomalistic year, nor over the apsidal cycle, even in the most densely sampled years of LLR data acquisition.

Lunar Drift Is a Model Value
Lunar drift ascertained from LLR involves combining several different models and applying many correction terms, each of which dwarfs the magnitude of the quantity being sought ( Table 2). The models compute the barycenter orbit about the Sun, the lunar orbit around the geocenter, lunar orientation (mirror position), station positioning, and atmospheric refraction. The~1 cm weighted residual for 1987 to the present, mentioned in several papers (e.g., [6,50]), refers to differences between the combined multi-parameter models and LLR measurements. Given the findings of Transtrum et al. [29], who showed that multiplied parameters can reproduce data with incorrect physics, the small residual demonstrates neither accuracy, nor precision, nor reliability.

5.
Non-uniform distributions in the data set are one contributor to correlations between solution parameters, as stated by [44] in summarizing William et al. [52]. Basically, LLR data are insufficient to constrain the models. 6.
Importantly, apparent agreement exists between LLR raw data and ephemeris tables roughly midway between perigee and apogee ( Figure 12). This section of the orbit is unimportant because an elliptical orbit is defined by its apogee and perigee. The apsides are under-sampled by LLR, and are most sensitive to corrections. This serious problem is not discussed in available reports. 7.
Atmospheric effects, producing refraction, are highly variable. This poorly constrained correction term is circa 2 m, independent of all other modelling efforts. 8.
Last, but not least, elliptical orbits are described by two parameters: the semi-major axis and the eccentricity. Measurements of radius alone are insufficient, as two measurements are needed to solve for two unknowns. This serious problem has been overlooked.

Discussion
Motions in the Earth-Moon-Sun system are complex. Short-term behavior is not captured by modeling LLR measurements (Section 6), so long-term evolutionary behavior is unconstrained. Lunar orbital evolution was evaluated theoretically in Sections 4 and 5. Implications of these findings are covered next.

The Lunar Orbital Radius Is Decreasing, Not Increasing
Changes in Earth's spin with time are dissipative, arising during differential rotation of its layers which has been detected in seismologic studies [34,35]. There is no need to invoke transfer of angular momentum to the Moon, because all the interesting action involving spin occurs inside Earth [33].
The difference between Earth's polar and equatorial radii is too small for this body to exert non-central forces on the distant Moon, based on analytical results for attraction to the oblate. It is far less possible for the significantly smaller tidal bulges to affect the lunar orbit. That the Moon's orbit is a point-mass problem is supported by the lunar orbital plane being substantially inclined to Earth's equatorial plane, and by trends displayed by the moons of the giant planets (Section 3). It follows that the lunar orbit is governed by point mass interactions.
Astronomical observations [19,45,46] show that the current, average lunar orbit, even 4.5 × 10 9 y after formation, remains nearly circular (Figure 7a), despite Solar forces acting on the Moon being double Earth's. This observation testifies to stability of the lunar orbit. Conservation of energy, and Solar forces providing torque in the 3-body geometry cause eccentricity to increase, which contracts the semi-minor axis while maintaining the semimajor axis (Figure 7a). Consequently, the average lunar radius has decreased over geologic time (−0.13 mm y −1 ; Figure 7b). Inclination of the lunar orbit from Earth's equatorial plane is consistent with Solar torque existing, unopposed, but the inclination is not part of the orbital energetics. The change in i is rather small and is consistent with stability and tiny changes in eccentricity.

Comparative Planetology and Evolutionary Behavior throughout Our Solar System
Orbits around a single central body are controlled by its mass, size, and ellipticity, plus the distance to the small mass. For separation distances exceeding~27 planetary body radii, forces are effectively central and orbits are governed by point mass interactions (Figures 3-5). Inward densification of planetary mass makes these objects more similar to point masses than the constant density oblate bodies considered in Section 3, and therefore promotes Keplerian behavior over much of the solar system.
Our recent study of interplanetary interactions [12] shows that planets act as point masses on each other in pairs, resulting in multiple periodic surges in forces and also torques. Previous models [53] treated the perturbing planet as a ring that exerts a central force, which does not occur.
The variable influences on any given planet of the 7 other planets produce many cycles for each of the 8 planets: these effects sum, yielding patterns more complex than perturbations of the Moon. Yet, from energy conservation and symmetry arguments, eccentricities are increasing whereas inclinations decrease as the planets are drawn to Jupiter's orbital plane. Solar system data provide confirmation [12]. Jupiter is highly influential because it is most massive and is near the rocky planets. In the outer solar system, the more distant and second most massive planet (Saturn) is also influential, as evidenced by the inclinations of Neptune and Uranus. Mercury's orbit is the most elongated and is changing the fastest, because it is the smallest and innermost planet. Rates of all planets (de/dt) vary from 1.5 × 10 −12 per year to 4.5 × 10 −11 which are similar to and bracket the Moon's rate, established here as 1.2 × 10 −11 .
Elongation appears to set the Moon on a collision course with the Earth (Figure 7a). Before this occurs, elongation of the lunar orbit may be sufficient that the reduced 2-body approximation no longer reasonably approximates the dynamics. For example, at apogee for e~0.98 (Figure 7a), the Solar force would be~10 times Earth's per Equation (14). This configuration seems unlikely to persist. Essentially, the Sun is capturing the Moon. The absence of moons for Venus and Mercury suggests that the Sun captured their moons long ago [33]. In contrast, planets in the distant outer solar system have retained their moons, as suggested by their large numbers.
For Uranus, di/dt is negative with a magnitude below 2 × 10 −8 • y −1 , assuming an initially polar orbit, which is consistent with orientation of its current spin axis and with its equatorial satellite orbits. Assuming that the Moon's orbit was in Earth's equatorial plane at formation suggests a similar progression of~0.5 × 10 −8 • y −1 . Both Earth and its Moon are being pulled to the orbital plane of Jupiter around the Sun, as are Mercury, Venus, and Mars. However, much stronger forces from the Sun has caused wobbling of the lunar plane about the barycenter plane ( Figure 1) which is still not parallel to Jupiter's orbital plane.
After 4.5 billion years, eccentricities and inclinations actually have changed very little, based on the rates deduced from spherical orbits. Observations today testify to the stability of the Solar System, which includes our Moon.

Implications for Formation in General and the Moon in Particular
Gravitational forces are the key driver of planetary formation, as discussed by many authors. Gravitational forces also govern their subsequent evolution, via perturbations, as discussed here and earlier [12]. Governing physical principles are summarized as follows: • Large perturbing forces limit the applicability of the reduced 2-body approximation. • Due to geometry, the fixed-plane three-body approximation is likewise insufficient to describe the Earth-Moon-Sun system.

•
The highly variable 3-dimensional geometry of the lunar orbit permits the Sun to apply torque, which changes angular momentum. • Forces on the Moon imperfectly balance, mainly due to the time-varying distance of the Moon from the Sun. • However, without dissipation of orbital energy, which requires a non-conservative force such as friction, orbital evolution is limited to changes in eccentricity and inclination, as long as the 2-body approximation is reasonably accurate.
The unique dynamics of the Earth-Moon-Sun system reveals information regarding its origin and evolution. The initial orbit was more circular, because the Sun elongates the lunar orbit in most positions (Figure 2). The Moon's orbit was probably initially in Earth's equatorial plane, but has been modified by the Sun consistently pulling the Moon toward the orbital plane of the barycenter while Jupiter simultaneously has pulled bodies in the inner solar system towards is orbital plane [12]. Jupiter's action on the Moon is very small compared to that of the Sun, but because this torque is unopposed, Jupiter has some influence.
Circular orbits describing formation is consistent with contraction of a spinning cloud of dust and gas, see, e.g., [23]. Notably, conservation of angular momentum causes the radius of the outer shell of a cloud to grow as the interior shrinks: see discussion of galaxies by Criss and Hofmeister [54]. For this reason, the Earth has a distant moon. The many distant moons in the outer solar system satisfy conservation of angular momentum during formation of the large planets. Moreover, division of the outer and inner solar systems is a consequence of contraction of the solar nebula to the center being balanced by expulsion of the outer shell of gas and dust. Conservation of angular momentum is not part of currently popular fractal accretion models (e.g., [55]).
Importantly, circular orbits are the only solution for a truly isolated 2-body point mass system. Close to an oblate central body, polar elliptical orbits are stable, but none are observed in the solar system. Tilted orbits close to an oblate could also be elliptical but these are not stable, and are very rare ( Figure 5). Multi-body interactions are largely, if not entirely, the cause of the currently elliptical orbits in our Solar system.
The above findings support, if not require, co-accretion of the Earth and the Moon. Formation takes finite time and so larger bodies necessarily form more slowly. Existence of miniature solar systems indicates that the nuclei of the planets were present during contraction of the overall nebula to the center, forming the Sun. Planets survive today due to the stability of their orbits, whereas the unconsolidated matter and objects in unstable orbits were sucked to the center. As the small dust-gas clouds of the planets contracted, their outer parts expanded, yielding satellite systems. Earth may have had more distant moons, which were robbed by the Sun long ago. Like the Earth, Saturn and Neptune likewise each have one moon that is much larger than the others. Unlike Earth, the outer planets are mostly composed of volatiles, and very light H and He gas, so these seem much larger than their moons, which are quite rocky. Thus, the size of the Moon compared to the Earth is not unusual, as the rock and metal cores of the giants are considerably smaller than their enshrouding atmospheres (see, e.g., [56]). Our Moon is simply not as unusual as currently perceived.

Relationship of Our Work to Previous Studies
Prior to LLR measurements which began in late 1969, lunar drift was unknown since radial motions from astronomical observations were insufficiently accurate to constrain small changes (Section 2). Hence, earlier estimates of lunar recession were maxima inferred from assumed formation conditions. Subsequently, the small positive drift suggested by LLR models were well received because they supported popular models for the Moon's origin.
Inconsistencies and assumptions in models (Sections 1 and 6) motivated us to closely examine the physical principles underlying evolution of the lunar orbit. Our theoretical approach contrasts with the LLR modeling, which stems from Newcomb's mathematical approach to deciphering orbital parameters, to which several additional correction terms are added.

Conclusions
Conservative conditions govern temporal changes in orbits because friction between distal bodies is null. In contrast, axial spin of bodies is dissipated non-conservatively and internally due to friction, torque, and heat production accompanying differential rotation of the internal layers (Section 4) and has no impact on orbital evolution. Thus, orbital evolution consists of increases in eccentricity and variations in inclination as a result of three-body interactions: These tiny perturbing forces and torques arise from distant objects. Even for the Moon, where the perturbing object, the Sun, exerts twice the force of the Earth, the changes are slow (Section 5) and similar to rates at which e and i of planets have changed since formation [12].
We have shown that long-term (~10 9 year) changes in orbits are not constrained by lunar laser ranging measurements (Section 6). Instead, we use conservation laws to deduce the following:

•
The lunar orbital radius is shrinking, so drift is about −0.13 mm y −1 .

•
The lunar orbit was originally circular and probably about Earth's equator.
• Co-accretion is favored, if not proven, since these findings rule out all other contenders for the lunar beginnings (a giant impact, fission, or capture).
Regarding broader implications of this study, the nearly circular paths remaining after 4.5 billion years on diverse scales in the solar system demonstrate the radial nature of gravitational pull. Large distances between small bodies make the two-body approximation valid over much of the solar system, which allowed Kepler to analyze the motions and to set the stage for development of the theory of gravity [56]. Our results point to a need for further exploration of Newtonian forces and torques that arise in three-body, pair-wise point-mass interactions.