Early Dynamics of the Lunar Core

The Moon has a much smaller bulk density than Earth, indicating that the lunar iron core is much smaller than that of Earth. Lunar laser ranging (LLR) and GRAIL mission results indicate that the Moon’s core has a radius of about 340 km (with the 200-380 km diameter range allowed by data) and is (at least partially) liquid (Williams \BOthers., \APACyear2006, \APACyear2014). The presence of the core is further suggested by magnetic induction (Hood \BOthers., \APACyear1999; Shimizu \BOthers., \APACyear2013) and seismic modeling (Weber \BOthers., \APACyear2011; Garcia \BOthers., \APACyear2011). While the Moon does not have a global magnetic field, remanent magnetization observed from orbit and in lunar samples indicates that there was a lunar magnetic field in the distant past. Since the magnetic fields of terrestrial planets are thought to be generated in their liquid iron cores, the rotational dynamics of the Moon’s core throughout lunar history is directly relevant to studies of ancient lunar magnetism.

Lunar rotational dynamics is very different from that of Earth. The Moon is in synchronous rotation due to Earth’s tidal forces and the resulting tidal dissipation, and this state was likely established very soon after the Moon’s formation. Tidal dissipation on Earth over billions of years made the Moon’s semimajor axis expand from several Earth radii ($R_{E}$) at formation to the present value of 60.3 $R_{E}$. As the lunar orbit and orbital period grew, so did the Moon’s rotation period. While the Moon’s rotational period is equal to its orbital period, the Moon’s obliquity is determined by the so-called Cassini States. Cassini States are minimum-energy solutions for the forced obliquity of a rotating body that experiences tidal dissipation, and most synchronous satellites have evolved into stable Cassini States (Colombo, \APACyear1966; Peale, \APACyear1969). The spin axis of the satellite in a Cassini State is in the same plane as the normals to the orbital plane and the Laplace plane (the plane around which the orbit precesses). When the period of the axial precession is longer than that of orbital precession (as is the case for the present-day Moon), the satellite is in Cassini State 2, in which the spin axis is closer to the Laplace plane normal than it is to the orbit normal (in the present-day Earth-Moon system, the Laplace plane is the ecliptic). Ward (\APACyear1975) first established that the Moon underwent a transition between the Cassini States 1 and 2 at a semimajor axis of about 34$R_{E}$, assuming that the Moon’s current shape was already established at this time. The necessity of the Cassini State transition is dictated by the dependence of the lunar spin and orbital precession on Earth-Moon separation. The lunar spin precession is driven by Earth’s tidal torque and therefore slowed down with the Moon’s tidal recession. However, beyond the Earth-Moon distance of about $15R_{E}$ (within which Earth’s oblateness was an important perturbation) the tidal recession led to speeding up of lunar apsidal precession, as it is driven by solar perturbations. Therefore, to reach the present state where the natural periods of nodal and spin precession are about 19 and 80 yr, the Moon had to experience equality of these two periods at some point in its evolution, which Ward (\APACyear1975) placed at the Earth-Moon distance of about $34~R_{E}$. Ward (\APACyear1975) also found that the Moon would experience very large forced obliquities during the Cassini State transition.

Meyer \BBA Wisdom (\APACyear2011) examined how the Moon’s liquid core would behave during the obliquity evolution established by Ward (\APACyear1975). Meyer \BBA Wisdom (\APACyear2011) show that the Moon’s free core nutation (FCN) period should be several centuries, much longer than the current nodal precession period of the lunar orbit (18.6 yr). Since the Moon’s spin axis is in a Cassini State and its precession follows that of the orbit, the natural frequency of nutation of the core around the mantle is too slow to allow the core to follow the mantle’s motion. Instead, the lunar core spin axis stays close to the ecliptic normal, in an arrangement somewhat equivalent to Cassini State 2 (the core spin axis is in the plane defined by the ecliptic normal and the mantle spin axis, and the two spin axes are on the opposite sides of the ecliptic normal; Fig. 1). This is in contrast to the case of Earth, where FCN takes only about 400 days, coupling the core to the mantle, which itself is precessing much more slowly, with a 26 kyr period. Meyer \BBA Wisdom (\APACyear2011) find that this decoupling of the lunar core’s rotation axis from that of the mantle should have happened relatively early in lunar orbital evolution, when the Moon had semimajor axis $a=26-29R_{E}$, before the mantle’s Cassini State transition. Note that the core being “de-coupled” from the mantle means that the two are spinning around different axes with different precessional motion; this terminology does not imply that there is no contact or no interaction between the core and the mantle. Dwyer \BOthers. (\APACyear2011) recognized that the results of Meyer \BBA Wisdom (\APACyear2011) imply a large past mutual obliquities between the core and the mantle, and Dwyer \BOthers. (\APACyear2011) argued this differential precession may have been the ultimate driving force behind the lunar dynamo.

More recently, Chen \BBA Nimmo (\APACyear2016) have shown that lunar obliquity tides should have significantly affected the past lunar inclination. Ćuk \BOthers. (\APACyear2016) proposed a model in which the Earth has a very fast spin ($<$2.5 hour period) and a high obliquity ( 70^∘) after lunar formation. Independent of the constraints from lunar inclination, this very fast initial spin of Earth is a requirement by some of the current lunar formation models in order to explain Earth-Moon isotopic similarity (Lock \BOthers., \APACyear2018). As the lunar orbit grows in the Ćuk \BOthers. (\APACyear2016) model, it encounters an unstable transition between equatorial and ecliptic Laplace planes (Tremaine \BOthers., \APACyear2009; Tamayo \BOthers., \APACyear2013). The lunar orbit is temporarily trapped in this zone, where solar perturbations induce large eccentricities, loss of the system’s angular momentum (which is transferred to the heliocentric orbit), a reduction of Earth’s obliquity, and a large lunar inclination ($\simeq 30^{\circ}$). In this picture, the Moon has a high inclination as it enters the Cassini State transition, with strong obliquity tides during the transition reducing the inclination to close to the observed values. In this paper we will consider the dynamics of the lunar core in the context of the Ćuk \BOthers. (\APACyear2016) model, and determine the consequences for the core’s orientation and the Moon’s tidal evolution.

In this work we use a custom made numerical integrator cr-sistem, which combines mixed-variable symplectic integrator for orbital motion and a Lie-Poisson mapping for the lunar rotation. In all aspects except the treatment of lunar rotation (including the motions of the core), cr-sistem is identical to r-sistem used by Ćuk \BOthers. (\APACyear2016). Therefore we direct the reader to the Methods section of Ćuk \BOthers. (\APACyear2016) for all details regarding the treatment of mutual perturbation between the bodies, tides within Earth and the Moon, and the shape and precession of Earth. In this section we will describe our approach for integrating the rotations of the lunar mantle and core, and their mutual interaction.

Ćuk \BOthers. (\APACyear2016) used the approach of Touma \BBA Wisdom (\APACyear1994), who advance the rotation of the Moon in the body-fixed reference frame implicitly, while the changes to this motion arising from oblateness and triaxiality are integrated explicitly (this algorithm is ultimately based on the Poincaré-Hough method Hough, \APACyear1895; Poincaré, \APACyear1910). The fact that these perturbations to the uniform rotation are relatively small enabled Ćuk \BOthers. (\APACyear2016) to efficiently integrate the rotational dynamics with only dozens of time steps required for each orbital/rotational period. However, inclusion of a triaxial core within a triaxial mantle makes the Hamiltonian much more complex.

Touma \BBA Wisdom (\APACyear2001) discuss in some detail the optimal algorithms for numerically integrating the rotational dynamics of a planet with a solid mantle and a liquid core. While Touma \BBA Wisdom (\APACyear2001) are able to derive a reasonably efficient algorithm for the dynamics of an oblate (i.e. rotationally symmetric) planet with a core, there is no such option for a triaxial core-mantle system. As suggested by Touma \BBA Wisdom (\APACyear2001), we divide our Hamiltonian into three parts, which correspond to rotations around the three Cartesian axes:

where $A,B,C$ are the principal moments of inertia and $P,Q,R$ are the components of angular momentum measured with respect to the mantle’s principal axes (values without subscript refer to the whole body, and those with the subscript ”C” to the core. As Touma \BBA Wisdom (\APACyear2001) show, $P_{C}$, $Q_{C}$ and $R_{C}$ are not really the angular momentum components of the core, but function equivalently. Additionally, the core moments of inertia are $F_{C}={2\over 5}m_{C}bc$, $G_{C}={2\over 5}m_{C}ac$, $H_{C}={2\over 5}m_{C}ab$, and we define $\alpha=AA_{C}-F_{C}^{2}$, $\beta=BB_{C}-G_{C}^{2}$ and $\gamma=CC_{C}-H_{C}^{2}$, where $m_{C}$ is the core’s mass and $a$, $b$ and $c$ are the core’s principal axes (our approach and notation follows that of Touma \BBA Wisdom, \APACyear2001). These three Hamiltonians each generate separate rotations of total and core angular momentum vectors ${\bf M}$ and ${{\bf M}_{C}}$ around the corresponding mantle axes. For example, the angular velocities generated by $\mathcal{H}_{z}$ for ${\bf M}$ and ${\bf M}_{C}$ are $w_{z}=(C_{C}R-H_{C}R_{C})/\gamma$ and $w_{Cz}=(CR_{C}-H_{C}R)/\gamma$ (with the sense of rotations reversed for $w$ and $w_{C}$ due to our mantle-fixed reference frame). Using these angular velocities, we can define rotation matrices that convert vectors between the mantle-fixed and inertial frames, and calculate periodic “kicks” on the rotational motion due to gravitational torques and tidal dissipation.

It is clear from our definition of $w$ that the rotation frequencies of the mantle around all three axis include terms independent of the core’s properties (at least to the lowest order) and therefore $w$ does not vanish even in the limit of the core’s size and mass being zero. The core-mantle cross-terms do become negligible, but the main components of $\dot{\bf M}$ remain, and correspond to the rotation of ${\bf M}$ around the mantle (or, in an inertial frame, rotation of the mantle around the total angular momentum vector). This means that the integrator needs to resolve this basic rotational motion which is handled explicitly, in contrast to the treatment in Touma \BBA Wisdom (\APACyear1994) (and subsequently Ćuk \BOthers., \APACyear2016) where only the slower nutational motion (separate from principal axis rotation) had to be resolved. Our numerical model which includes the lunar core therefore requires significantly larger number of time-steps per rotation period. There is a direct parallel here to orbital integrations, with our integrator being analogous to a basic $T+U$ symplectic integrator (Forest \BBA Ruth, \APACyear1990), which requires enough time-steps to resolve basic Keplerian motion, while the Touma \BBA Wisdom (\APACyear1994) scheme is analogous to mixed-variable symplectic integrators (Wisdom \BBA Holman, \APACyear1991) which assume Keplerian orbits, so fewer time-steps are needed per orbit. Therefore, we are not in a good position to replicate the long-term integrations of Ćuk \BOthers. (\APACyear2016), as inclusion of the lunar core would increase integration time by at least an order of magnitude. Therefore we will restrict ourselves to examining how the presence of the lunar core affects the Cassini States of the Moon at different Earth-Moon distances, and determining the mutual motions of the core and the mantle for the purposes of understanding the past lunar dynamo.

One remaining choice we had to make was the introduction of core-mantle friction. We calculate the vector difference between the rotational velocities of the mantle and the core $\Delta{\vec{\omega}}$, and then apply a torque to $\bf{M_{C}}$ every time-step:

where $dt$ is the time-step and $K_{CM}=8.43\times 10^{-3}$ year^-1 is the magnitude of spin-orbit friction. We based this form of core-mantle torque that is linear in the core-mantle relative motion on Eq. 19 in Pavlov \BOthers. (\APACyear2016). Our value for $K_{CM}$ corresponds to a damping timescale of $(1000/2\pi)$ sidereal months, or 120 years, and is equivalent to Pavlov \BOthers. (\APACyear2016) $k_{v}/C_{T}=16\times 10^{-9}$ day^-1. Our use of linear core-mantle boundary (CMB) friction follows the usage in LLR community, but the flow at the CMB boundary is likely to be turbulent and non-linear (Williams \BOthers., \APACyear2001), so our model should be taken as a first approximation. As we ignored any tidal forces from Earth or the Sun on the core, the core-mantle friction term (Eq. 4) is the only dissipative influence on the core’s rotation, and it acts to match it to the rotation of the mantle. However, non-dissipative core-mantle torques inherent in our integrator make the core precess at a different angle and phase from the mantle, which can lead to equilibrium points very different from core-mantle co-rotation, as we show in the following sections.

After completing our integrator, the first situation we wanted to explore is the present-day dynamics of the lunar core. We used the same lunar shape parameters as Ćuk \BOthers. (\APACyear2016), and we decided to set the mass and radius of the lunar core exactly to 0.01 $M_{M}$ and 0.2$R_{M}$, and we will use these values throughout the paper. Note that the ellipticity of the core is by far the most important parameter for its dynamics, and the mass and density have only second-order effects on the core’s nutation. Here we used shape parameters $(C_{C}-A_{C})/C_{C}=2.5\times 10^{-4}$ and $(C_{C}-B_{C})/C_{C}=2\times 10^{-4}$ for the core (ellipticity was based on Pavlov \BOthers., \APACyear2016). We used long-term average lunar tides (quantified by tidal quality factor $Q_{M}=60$ and Love number $k_{2M}=0.024$).

Figure 2 shows a $10^{6}$ yr simulation in which the mantle was put in synchronous rotation and Cassini State 2 with obliquity to ecliptic $\epsilon_{M}=1.5^{\circ}$, while the core is made to approximately co-rotate with the mantle. For comparison, we also conducted a separate experiment using an entirely solid Moon with no fluid core. In about 300 years the core decouples from the mantle’s spin direction and settles in an apparently stable state with obliquity to ecliptic of $\epsilon_{C}=0.1^{\circ}$. The spin rate of the core on the other hand settles close to that of the mantle, firmly within the range of mantle’s rotational librations (bottom). The mantle is meanwhile damping librations around its Cassini state close to $1.58^{\circ}$, with the damping timescale being much faster than for a solid moon, indicating that the damping is dominated by CMB friction.

Figure 3 shows the situation at the end of the same simulation as shown in Fig. 2. Here we plot the $y$-component the core and mantle spin vectors in the polar representation $\epsilon_{i}\cos(\Omega_{i})$, where $\epsilon$ is the obliquity with respect to ecliptic and $\Omega$ the longitude of the node (i.e. intersection between the plane normal to the spin axis and the ecliptic) for either the mantle and the core. This plot enables us to see how the obliquity vectors are oriented in space. From Fig. 3, it is clear that the mantle and the core spins are both precessing around the ecliptic normal with the 18.6 yr period dictated by the (orbital) nodal precession. While the core’s obliquity is much smaller, it is precessing exactly out-of-phase (i.e. $180^{\circ}$ away) from the mantle, as expected from a ”quasi-Cassini State 2” (see Fig. 1) in which the natural precession frequency of the core (around the mantle) is much slower than the mantle’s precession around the reference direction (in this case, the ecliptic normal). As explained in Meyer \BBA Wisdom (\APACyear2011), the core should precess around the mantle with the frequency:

where $\omega$ is the lunar spin rate, and $C_{M}=C-C_{C}$ the mantle moment of inertia (for our purposes, $C/C_{M}\simeq 1$). Since we use $(C_{C}-A_{C})/(C_{C})=2.5\times 10^{-4}$, and the lunar spin period is 27.3 days, then the period for the core nutation should be $\simeq 300$ yr. Since $\epsilon_{C}<<\epsilon_{M}$, we can estimate $\epsilon_{C}/\epsilon_{M}\simeq\omega_{FCN}/\dot{\Omega}$, where $\dot{\Omega}$ is the nodal precession rate for the Moon (based with analogy with Cassini State 2 Ward, \APACyear1975). Therefore, the expected forced obliquity of the lunar core should be $\epsilon_{C}=1.6^{\circ}(18.6{\rm yr}/300{\rm yr})\simeq 0.1^{\circ}$, fully consistent with the results of our simulation. These results are also very similar to the findings of Dumberry \BBA Wieczorek (\APACyear2016) and Stys \BBA Dumberry (\APACyear2018) who additionally took into account the existence of the inner solid core. So we can conclude that for the present Earth-Moon system, our numerical model returns an answer consistent with the analytical expectations, with the lunar core precessing separately from the mantle with a much smaller obliquity to the ecliptic.

We now turn our attention to the transition between the core being coupled and uncoupled from the mantle, which is in many ways analogous to the Cassini State transition that the mantle experienced somewhat later, with the difference that the former involved the core and the mantle, and the latter the mantle and the lunar orbit. Meyer \BBA Wisdom (\APACyear2011) calculated that the core should have uncoupled at around $26R_{E}$, assuming that the mantle was not in hydrostatic equilibrium but had a figure identical to the present one. As explained in Ćuk \BOthers. (\APACyear2016), this is probably the smallest Earth-Moon distance at which a non-hydrostatic mantle would be a reasonable assumption. However, in our simulations of the Moon at this distance from Earth, we did not exactly reproduce Meyer \BBA Wisdom (\APACyear2011) result, and we kept getting an uncoupled core in quasi-Cassini State 2 (note that the lunar mantle is at this point in Cassini State 1 relative to the orbit). We tried simulations at even smaller Earth-Moon distances, but even for $a=18R_{E}$ we were getting a decoupled core (Figs. 4 and 5), using $(C_{C}-A_{C})/C_{C}=2.5\times 10^{-4}$ (as elsewhere in the paper). However, when we reduced the semimajor axis to $a=15R_{E}$, the core in the simulation is now coupled to the mantle (Figs. 6 and 7). Note that these simulations only intended to demonstrate that the core-mantle coupling is possible for small Earth-Moon distances, and are not in any way an accurate representation of the Moon’s dynamics at the time, as the mantle was likely not rigid at this time, and we did not try to accurately match Earth’s rotation and oblateness which affect the Moon’s orbital precession.

We conclude that the criterion used by Meyer \BBA Wisdom (\APACyear2011) for the core-mantle coupling, $\omega_{C}>\dot{\Omega}$, is approximate and cannot be used to accurately determine the distance at which the lunar core dynamically decoupled from the mantle. Numerical integrations suggest that the decoupling probably happened at lunar semimajor axis somewhat smaller than $26R_{E}$ estimated by Meyer \BBA Wisdom (\APACyear2011). This is a relatively minor discrepancy, as the nutation period of the core is about 50% and 33% of the mantle precession period in our Figs. 4 and 6, respectively. Therefore, our simulations suggest that dynamical decoupling happens when the core nutation rate is about twice the nodal precession rate. A more accurate determination would require many more simulations and is of limited interest to this paper. Dumberry \BBA Wieczorek (\APACyear2016) have conducted an analysis of core-mantle coupling (but without dissipation) and found that the exact identity of the mantle precession rate and the FCN is in fact a resonance which forces a large core-mantle tilt; they find that FCN must be significantly faster than mantle precession to enable dynamical coupling. Therefore, the main conclusion of this section is that, assuming a rigid lunar mantle with the present shape, the core should be dynamically decoupled from the mantle for lunar semimajor axes $a>25R_{E}$ that are relevant for the Cassini State Transition and lunar paleomagnetism.

In order to have a self-consistent model of past lunar core dynamics, we need to re-examine the dynamical history of the lunar mantle with the core’s influence included. Lunar inclination damping due to obliquity tides, as modeled by Ćuk \BOthers. (\APACyear2016), directly depends on the mantle’s rotational state, which is in principle affected by the presence of the core. It is important to revisit these calculations and determine whether the results obtained by modeling the Moon as a solid body are still valid when core-mantle interaction is included. However, due to the limitations of our current integrator relative to that used by Ćuk \BOthers. (\APACyear2016), we are not able to directly integrate the lunar tidal evolution, as complete simulations would take unreasonably long time.

In this section we will instead present a series of simulations of relaxation to equilibrium of the lunar core and mantle. These simulations are designed to follow the evolutionary track in semimajor axis and inclination found by Ćuk \BOthers. (\APACyear2016). We started these simulations with a sub-synchronous Moon, let it damp to a stable state, and recorded the final rotational state of the core and the mantle. While this is not a model of lunar orbital history independent of Ćuk \BOthers. (\APACyear2016), is an important test of their model. If the rotational state of the mantle, especially the obliquity, is fundamentally changed by the presence of the core, the results of Ćuk \BOthers. (\APACyear2016) would not be valid. On the other hand, if the simulations with the core-mantle interactions included produce lunar obliquities close to those in the core-less model, the overall nature of the lunar tidal evolution is likely unchanged by core-mantle interactions.

In Fig. 8, we plot the results of our equilibrium solutions along the track of the lunar tidal evolution found by Ćuk \BOthers. (\APACyear2016). We generated these solutions by placing the Moon at a sequence of semimajor axes, with the spacing of 1 $R_{E}$. At each semimajor axis, the Moon was given orbital inclination based on the orbital evolution track taken from Figure 4 in Ćuk \BOthers. (\APACyear2016) (the same data are plotted by small red dots in Fig. 8), while the eccentricity was set to $e=0.01$. The solid blue circles show the end-states of the simulations after 1 Myr (extended to 2 Myr for the two simulations at 30-31 $R_{E}$ which showed slowest convergence) averaged over the final 1000 years. Due to short-period variations in orbital elements, our initial conditions and averaged final states do not perfectly match in either inclination or semimajor axis, causing the blue circles to be slightly misaligned with the red dots they track.

In panel B, we compare our solutions for the mantle obliquity for the specific $a$ and $i$ plotted in panel A (solid blue circles) to the core-less Moon’s obliquity found by Ćuk \BOthers. (\APACyear2016) (small red dots). The solid black line plots the Cassini State obliquities corresponding to the inclinations in panel A (unlike the numerical simulations, the Cassini State calculation, by definition, assumes synchronous rotation). In panel C we plot the corresponding rotational rates from our mantle equilibrium solutions (solid blue circles) and the Ćuk \BOthers. (\APACyear2016) core-less Moon simulation (small red dots), with the thin black line plotting the synchronous state.

Our solutions clearly show similarity to the previously found evolution, including the divergence from a Cassini State. However, while the divergence from the Cassini state for a core-less moon was limited to the immediate aftermath of the Cassini State Transition, our simulations suggest that the Moon was in non-synchronous rotation for all semimajor axes in the $a=25-40R_{E}$ range. Particularly interesting is the divergence from Ćuk \BOthers. (\APACyear2016) results in semimajor axis range 26-29 $R_{E}$, where our mantle is in a non-synchronous state, while Ćuk \BOthers. (\APACyear2016) find the Moon to be in synchronous rotation. There are two distinct reasons for this divergence. One is from the different ways the two solutions were calculated: Ćuk \BOthers. (\APACyear2016) integrated a single continuous simulation of lunar tidal evolution, while we performed a number of short simulations that started with the Moon in non-synchronous rotation. Therefore, in this segment of Ćuk \BOthers. (\APACyear2016) simulations the Moon is likely in a quasi-stable Cassini State 1, which is temporarily stable against Earth’s torques, but would not be re-established if it were to be broken by an outside factor (such as a large impact). In the Ćuk \BOthers. (\APACyear2016) simulation, the synchronous rotation is broken by an annual libration resonance at 29.7 $R_{E}$, after which the Moon settles in a sub-synchronous state adjacent to Cassini State 2 (their Extended Data Figs. 8 and 9). In our core-mantle simulations the Moon at $a=27-29R_{E}$ is in a sub-synchronous state adjacent to Cassini State 1 (Fig. 8), which apparently does not exist in the Ćuk \BOthers. (\APACyear2016) model (their Extended Data Figs. 7 and 8).

Apart from the divergence prior to the Cassini State transition, we also observe that our simulations at semimajor axes 35-40 $R_{E}$ have the Moon in a non-synchronous state, while the Ćuk \BOthers. (\APACyear2016) find the Moon in the synchronous Cassini State 2 at this distance (beyond 40 $R_{E}$ we also find that the Moon is mostly in synchronous rotation). The main reason for the prevalence of non-synchronous rotation in our simulations is that the core, which lags behind the mantle in its spin rate, exerts a relatively substantial spin-changing torque on the mantle. This can move the equilibrium spin value outside of the synchronous spin-orbit resonance, producing a stable sub-synchronous state. While for some semimajor axes $a=30-35R_{E}$ the synchronous state itself may be unstable (due to too-high Cassini State obliquity; Beletskii, \APACyear1972), we find that in the $a=35-40R_{E}$ range the synchronous rotation is sometimes reachable by spin-down but not by spin-up (i.e. initial conditions may determine the end point rotation state reached by the Moon).

In general, the presence of the core drives the sub-synchronous states further away from synchronicity, and as a result the corresponding Cassini States have somewhat higher obliquities. This can be understood by recognizing that the pole of a slower-rotating Moon will precess more quickly, so to match the precession rate of the orbit, axial precession must be slowed down by high obliquity. On the other hand, we find that the Cassini State obliquities for synchronous cases (such as our $a=40R_{E}$ simulation) are lower when the core is included in the simulation. The reason for this change is that the core, which is significantly misaligned with the mantle, induces additional librations of the mantle in synchronous rotation. Librations make the precession slower (as the Moon presents a different average figure to Earth’s torques), requiring lower obliquities to reach the same precession rate as in the non-librating case.

Despite the modification of synchronous rotation caused by the presence of the fluid core, we find that the lunar orbital evolution history found by Ćuk \BOthers. (\APACyear2016) is unlikely to be fundamentally changed. Using the lunar inclination model published by Ćuk \BOthers. (\APACyear2016), we find that the lunar mantle has had a large obliquity for the significant part of the Moon’s history. Obliquities and associated tides in our model are sometimes lower, sometimes higher than those found by Ćuk \BOthers. (\APACyear2016), but there is every reason to think that large damping of lunar eccentricity took place. We also note that core-mantle friction also damps obliquity relative to the Laplace plane (Rochester, \APACyear1976), which leads to damping of inclination (as the obliquity is forced). Further more extensive numerical simulations of the full evolution of lunar spin and orbit in the presence of the core are needed to determine the exact contributions of different inclination-damping mechanisms.

Our results are dependent on our choice for the functional form of CMB friction (Eq. 4), in which the torque is directly proportional to the difference between core and mantle spin rates, with a constant damping timescale of about 120 yr. Other functional forms are possible, including those that are quadratic in relative core-mantle motion (Yoder, \APACyear1981). Given that the current core-mantle tilt is small, these alternative forms would produce much shorter timescales for damping of core-mantle motion when the lunar obliquity to ecliptic was high. We tried such a functional form of CMB friction in some simulations, in order to determine whether such strong friction could couple core to the mantle. Surprisingly, we find that beyond about $30R_{E}$ the CMB torque on the mantle is too strong for mantle to stay in the Cassini State, and the Moon evolves into a very slow-rotating, high-obliquity state with the core and the mantle still decoupled. We expect that this state would persist until the lunar inclination is damped, ultimately allowing core-mantle coupling in the absence of mantle’s forced obliquity. As a substantial lunar inclination has survived the Cassini State transition, this puts an upper limit on the intensity of past CMB friction, and may indicate that past quadratic CMB friction is not consistent with the present lunar orbit. However, our model includes a number of approximations, and we only explored relatively short-term dynamics of the core and the mantle at fixed Earth-Moon distances; full-scale evolutionary integrations are necessary to fully constrain long-term inclination damping due to quadratic CMB friction.

For the purposes of the generation of the lunar magnetic field, we are chiefly interested in the rotational state of the liquid core. In Fig. 9 we plot the final rotational parameters of the core in the simulations shown in Fig. 8. The solid squares and circles plot the results of our simulations featuring the core, while small dots plot the Ćuk \BOthers. (\APACyear2016) simulations of a core-less Moon. In the top panel we plot the core and mantle obliquities relative to the ecliptic, a convenient system of reference because the core’s spin axis is close to the ecliptic for the duration of the simulations. Unlike the core, the mantle has a much larger tilt to the ecliptic, both due to large orbital inclination (most important before the Cassini State transition) and the large forced obliquity (dominant after the Cassini State transition). As the core is in a non-synchronous rotation, we also plot the predicted core rotation rate based on that of the mantle using the solid magenta line. This rate is calculated as $\omega_{C}=\omega_{M}\cos(\psi_{CM})$, where $\omega_{C}$ and $\omega_{M}$ are the core and mantle spin rates, and $\psi_{CM}$ is the core-mantle obliquity. This calculated spin rate assumes that the core spin matches the component of the mantle’s angular velocity parallel to the core’s spin axis. Comparison between these predicted and numerically simulated spin rates of the core (Fig. 9, bottom panel) shows this to be a good approximation for the core’s final spin rate in the simulations.

Figure 10 re-plots some of our results for the core in terms of quantities measured relative to the mantle, and extends the calculation to the present day using a semi-analytical model. Since the Moon was most likely in Cassini State 2 for $a>40R_{E}$, and direct simulations take much longer for large semimajor axes, it is more efficient to calculate core and mantle rotational properties using a simplified model. We find the core obliquity as the solution to the equation:

where $\dot{\Omega}$ is the precession rate of the longitude of the Moon’s ascending node, and $\epsilon_{M}$ is the obliquity of the mantle with respect to the ecliptic. We obtain Eq. 6 from Eq. 5 (in this paper) and Eq. 1 in Ward (\APACyear1975), and assumed $\omega_{C}=\omega_{M}\cos(\psi_{CM})$ for the core’s rotation rate. Eq. 6 can alternatively be derived from the model of Stys \BBA Dumberry (\APACyear2018), by setting the size of the inner core to zero. Eq. 6 is solved numerically by simply advancing through obliquities and finding the value for which the sum of terms goes through zero. The mantle obliquity and lunar orbital precession are obtained from the semi analytical model of lunar tidal evolution described by Eqs. M1-M8 in the Methods section of Ćuk \BOthers. (\APACyear2016). We used $i=10^{\circ}$ and $e=0.01$ at $a=40R_{E}$ as initial conditions, and advanced the model forward, with Earth tidal $Q$ set at lifetime average of $Q_{E}=35$, and the lunar tidal $Q$, as well as Earth and Moon tidal Love numbers, set at present values $Q_{M}=38$, $k_{2}^{E}=0.3$ and $k_{2}^{M}=0.024$ (Williams \BBA Boggs, \APACyear2015). This model evolution gives us an inclination just over $5^{\circ}$ at $a=60R_{E}$, consistent with the present lunar orbit. Thick lines in Fig. 10 plot the core obliquity and spin rate derived from this model.

Fig. 10 shows that, while the numerical and semi-analytical obliquities mostly agree, there is a dramatic divergence for points at $a=43R_{E}$ and $a=50R_{E}$, for which the numerical model does not predict synchronous rotation. This is related to the shifting of the equilibrium spin rate from exact synchronicity due to the CMB friction torque on the mantle. It appears that the equilibrium rotation is just at the outer limit of the synchronous spin-orbit resonance, making reestablishment of synchronous rotation sensitive to initial conditions and possible secondary resonances. Here we used a low lunar eccentricity ($e=0.01$), but a more substantial eccentricity would shift the equilibrium spin value to higher values, counteracting the effects of CMB friction and shifting the equilibrium back into the synchronous rotation. Lunar eccentricity was likely affected by planetary orbital resonances at $a=46R_{E}$ and $a=53R_{E}$ (Ćuk, \APACyear2007), with the latter resonance likely bringing the lunar eccentricity close to the present value. The actual history of the lunar rotation depends on the timing of synchronous-lock-breaking impacts, the relative strength of satellite tides and spin-orbit friction, the time-evolution of lunar eccentricity and the (poorly understood but likely) influence of secondary resonances

It is evident from Fig. 10 that, in our model, when the Moon had semimajor axis $a<45R_{E}$ there was a large mutual obliquity between the lunar core and the mantle, and a significant difference between the rotation rates of the two. In order to compare the models of orbital and rotational dynamics to the record of lunar magnetism, we would need to know the absolute timing of the lunar tidal evolution. Unfortunately, there is currently no direct way to constrain the speed of early lunar tidal evolution. While the age of the Earth-Moon system implies Earth’s long-term average tidal $Q_{E}\simeq 34$, the current tidal $Q$ of Earth is about 12, meaning that there must have been at least some intervals of higher tidal $Q$ in the past. Bills \BBA Ray (\APACyear1999) show that often contradictory claims of tidal periods captured in the geological record only concern the last billion years or so, and the timeline of the early history of the Earth-Moon system is completely unconstrained. Some models of ocean dissipation (Webb, \APACyear1982) indicate general reduction of $Q$ over geological time, with tidal $Q$ being a factor of several higher than now in the first billion years of the Solar System, but so far there is no widely accepted correspondence between lunar distance and time.

In Table 1, we present two possible timelines for lunar tidal evolution as a function of Earth’s tidal $Q$ during the system’s early history. These numbers do not establish a preferred time-obliquity correspondence, but rather illustrate that a wide range of reasonable parameters would predict a large lunar core-mantle obliquity at the epoch when the lunar magnetic field is thought to have been strongest, 4.25-3.6 Gyr ago (Tikoo \BOthers., \APACyear2014; Weiss \BBA Tikoo, \APACyear2014). In Table 1 we used $\dot{a}\simeq a^{-11/2}Q_{E}^{-1}$ (Murray \BBA Dermott, \APACyear1999), i.e. we assumed a constant $Q$ for early Earth, and ignored lunar eccentricity and obliquity tides, as well as temporary halts of tidal evolution due to the Laplace Plane transition (Ćuk \BOthers., \APACyear2016).

The obliquities shown in Figs. 8-10 are model-dependent, and lunar tidal evolution models that produced lower past lunar inclinations (Touma \BBA Wisdom, \APACyear1994) would also result in smaller obliquities (although high obliquities close to the Cassini State transition are inevitable). However, lower past lunar inclinations in models predating work by Chen \BBA Nimmo (\APACyear2013) are primarily an artifact of the obliquity tides being ignored. While Peale \BBA Cassen (\APACyear1978) considered tidal heating from obliquity tides, they did not investigate associated inclination damping; Williams \BOthers. (\APACyear2001) proposed that tidal heating during the Cassini State transition powered the lunar dynamo, but did not address orbital consequences. Therefore, pre-2013 models of lunar tidal evolution cannot be directly used to derive past lunar mantle obliquities. More recently, the model of Pahlevan \BBA Morbidelli (\APACyear2015) gets around the inclination damping at the Cassini State transition by exciting lunar inclination through planetesimal flybys after a rapid early lunar orbital expansion. Since Pahlevan \BBA Morbidelli (\APACyear2015) propose that the Moon underwent a rapid orbital expansion in the first 10 Myr, and had a close-to-present inclination at $a\simeq 45R_{E}$ after the encounters were over ($\simeq 100$ Myr), it appears unlikely the Moon would have had large enough obliquity to power a dynamo at the “magnetic epoch” 4.25-3.6 Gya. While the effects of late planetesimal encounters need more exploration (we ignored them in this work as they happen before the Cassini State transition for tidal parameters we considered), it is likely that the Pahlevan \BBA Morbidelli (\APACyear2015) model is incompatible with a precession-driven dynamo, and that a different mechanism would be needed to power the ancient lunar dynamo.

In our calculations, we ignored the presence of the solid inner core. While the size of the inner core is poorly constrained (Williams \BOthers., \APACyear2014), it should have significant effects on the core dynamics by introducing new nutation frequencies (Dumberry \BBA Wieczorek, \APACyear2016) and it may in some cases significantly change the tilt of the outer core (Stys \BBA Dumberry, \APACyear2018). It is thought that the solid inner core grew over time (Scheinberg \BOthers., \APACyear2015), so it would have been much smaller (or not present at all) during the early history of the Earth-Moon system that we study here. However, assuming a later formation for the inner core may not be viable, as some models of the magnetic dynamo generation require the presence of the solid inner core (Stanley \BOthers., \APACyear2017), so improved dynamical models will have to include a two-component core. Based on our results for the outer core, the inner core is unlikely to be synchronously locked as Stanley \BOthers. (\APACyear2017) have proposed, as friction with the sub-synchronous outer core would likely be a more important factor for determining the inner core’s spin than Earth’s tidal torque. Additionally, Stys \BBA Dumberry (\APACyear2018) have found that the free inner core nutation period is currently in the 10-40 yr range, making a large forced tilt due to resonant interaction with the 18.6 yr mantle precession period likely.

Finally, we note that a magnetic field that is locked to the core, or changes its geometry relatively slowly relative to the core, would not be stationary relative to the lunar mantle. Earth’s core is practically co-rotating with the mantle and its spin axis is not significantly offset from that of the mantle (as the free core nutation is much, much faster than Earth’s axial precession). Therefore, experience from Earth is not directly applicable to the Moon, and a lunar dipole field generated in the core would drift across the lunar surface due to different rotations of the core and the mantle. For example, at $a=40R_{E}$, when we find the core-mantle obliquity to be 25^∘, a magnetic pole that is aligned with the core spin axis would sweep on the surface a circle at $65^{\circ}$ latitude in a course of a lunar orbital period ($\simeq 15$ days at that time). The pole of a dipole field that is significantly tilted relative to the core spin axis would trace a much more complex track on the timescales determined by the differential rotation of the core and the mantle. Therefore, directional magnetization of lunar rocks does not record a long-term orientation of the lunar magnetic field, but possibly only a snapshot taken at the critical points of the material’s thermal evolution. Some of the recent work on the locations of proposed paleopoles does indicate a great variability in the inferred orientation of the lunar dynamo (Oliveira \BBA Wieczorek, \APACyear2017; Nayak \BOthers., \APACyear2017), and relative motions of the core and the mantle may at least partially explain these findings.

In this work we numerically modeled the rotational dynamics of the Moon’s core and the mantle during the early history of the Earth-Moon system, and we report the following conclusions:

1. The lunar core was precessing independently of the mantle for practically all of lunar history.

2. Our numerical model also suggests that in the distant past the lunar core rotated significantly slower than the mantle, assuming that the main dissipative torque on the core comes from friction at the core-mantle boundary.

3. The presence of the small lunar core makes the non-synchronous rotation before and after the Cassini State transition even more likely, and the CMB friction should have caused additional damping of lunar inclination, separate from obliquity tides.

4. For much of the early history of the Earth-Moon system there was a large mutual obliquity between the core and the mantle. This has direct implications for the generation of the ancient lunar magnetic field, and we discuss plausible absolute timings of events in lunar tidal evolution.

Finally, we note that this is only the first step in numerically modeling the dynamics of the lunar core. Due to numerical limitations, we only explored stationary solutions for the lunar rotation, and in the future we hope to produce a self-consistent lunar tidal evolution model that includes the effects of core-mantle interactions. We also hope that additional factors such as the presence of the solid inner core and the possible non-rigid behavior of the core-mantle boundary will be addressed in future work.