Tidal Evolution of the Earth-Moon System with a High Initial Obliquity

A giant impact origin for the Moon is generally accepted, but many aspects of lunar formation remain poorly understood and debated. \'Cuk et al. (2016) proposed that an impact that left the Earth-Moon system with high obliquity and angular momentum could explain the Moon's orbital inclination and isotopic similarity to Earth. In this scenario, instability during the Laplace Plane transition, when the Moon's orbit transitions from the gravitational influence of Earth's figure to that of the Sun, would both lower the system's angular momentum to its present-day value and generate the Moon's orbital inclination. Recently, Tian and Wisdom (2020) discovered new dynamical constraints on the Laplace Plane transition and concluded that the Earth-Moon system could not have evolved from an initial state with high obliquity. Here we demonstrate that the Earth-Moon system with an initially high obliquity can evolve into the present state, and we identify a spin-orbit secular resonance as a key dynamical mechanism in the later stages of the Laplace Plane transition. Some of the simulations by Tian and Wisdom (2020) did not encounter this late secular resonance, as their model suppressed obliquity tides and the resulting inclination damping. Our results demonstrate that a giant impact that left Earth with high angular momentum and high obliquity ($\theta>61^{\circ}$) is a promising scenario for explaining many properties of the Earth-Moon system, including its angular momentum and obliquity, the geochemistry of Earth and the Moon, and the lunar inclination.


INTRODUCTION
The Moon is widely thought to have formed in the aftermath of a giant impact between the proto-Earth and another planetary body close to the end of terrestrial planet formation (Hartmann & Davis 1975  Incorrect handling of the axial precession of Earth was the main source of error in the results ofĆuk et al. (2016). Their integrator simply gave the Earth's axis a single "kick" per timestep (calculated using the positions of the Sun and the Moon), in a leapfrog scheme where perturbations were alternating with the Keplerian motions. While this approach did not lead to any secular errors when integrating the current Earth-Moon system or during the lunar Cassini state transition, significant drift in obliquity occurs when the Moon is in one of the secular resonances observed during the LPT and so usual symmetries do not apply. Since our integrator conserved the total spin of Earth, but kept artificially lowering the obliquity, this led to a secular increase in the ecliptic angular momentum. Figure 1 shows that about a third of the total change in AM inĆuk et al. (2016) Figure 1 was due to this "obliquity leakage". We corrected this problem by replacing this simple mapping by a fourth-order Runge-Kutta step (Press et al. 1992) which approximates the precession of Earth under lunar and solar torques. Unlike our treatment of orbital motions and lunar rotation, this approach is not symplectic, so absolute conservation of ecliptic AM is not expected. However, ecliptic AM is now approximately conserved and its variation in our simulations of the Laplace plane transition (LPT) is approximately 1% (Fig. 2). 4Ćuk et al. Cuk et al. (2016) assumed that Earth's oblateness moment depends on its spin rate ω R as J 2 ∝ ω 2 R . While this is accurate at low rotation rates, early Earth is significantly deformed in some high-AM scenarios (Ćuk & Stewart 2012;) and so we cannot assume oblateness is a small correction to the planet's figure. Based on the results of , we now use the relationship J 2 ∝ L 2 (1 − 0.64(L/L C )) where L is the Earth's rotational AM, and L C = α √ GM 3 R is the critical AM at which the equator of a (non-deformable) Earth would reach orbital velocity (M , R and α are Earth's present-day mass, radius and dimensionless moment of inertia, and G is the gravitational constant). Note that here we are using the Earth's AM rather than its spin rate, as our precession routine only needs to "know" about J 2 and L, and we do not need to track the spin rate explicitly. The fact that we encounter the Laplace Plane instability at the same distances as Tian & Wisdom (2020) indicates that our approaches yield approximately similar results. Cuk et al. (2016) designed their integrator to model the Cassini State transition of the lunar spin axis, which occurs later in lunar tidal evolution than the dynamical events studied in this work. Therefore, the earlier version of rsistem assumed a constant figure for the Moon. This is unlikely to be true during the earliest stages of lunar tidal evolution, when the Moon may have still had a magma ocean, and the lunar figure was likely in (or close to) hydrostatic equilibrium. WhileĆuk et al. (2016) periodically manually adjusted the lunar figure in their simulations of the LPT, here we implement automatic adjustment of the Moon to the triaxial figure calculated by Keane & Matsuyama (2014) at the relevant semimajor axis (assuming zero eccentricity and obliquity). As the Lie-Poisson approach we adopted (Touma & Wisdom 1994a) assumes a rigid-body Moon, we tried to minimize the number of figure adjustments, so we included it into a routine activated about every 1000 yr. Therefore, we adjusted the lunar figure about 10 5 times in the simulations shown in Figs. 1 and 3, which is only a small fraction of about 2 × 10 11 steps taken by the integrator. To make sure that our reshaping events do not lead to leakage of AM or energy, we re-ran a part of the simulation shown in Fig. 3 with a 10,000 year interval (rather than 1000 yr interval) of adjusting lunar shape and found no noticeable change in outcome.
We use essentially the same tidal model asĆuk et al. (2016). The similarity of our results to those of Tian & Wisdom (2020) suggests that our fundamentally different approaches produce consistent outcomes. However, in their supplementary discussion, Tian & Wisdom (2020) make several incorrect statements about our tidal model that we need to clarify. Our tidal model is based on instantaneous motions, both angular and radial, of Earth, the Moon and the Sun relative to each other. We do not use orbital elements in calculating tidal accelerations. This choice was made so we would not have to make any assumptions about lunar rotation, which can become asynchronous during the Cassini State transition (Ćuk et al. 2016, 2019). Therefore, when modeling lunar librational eccentricity tides, we use the instantaneous rotational motion of the Moon relative to the Earth-Moon line as the main input parameter. When calculating tidal despinning of a non-synchronous satellite, the "constant Q" approach means that the tidal bulge is assumed to be at a fixed angular distance ǫ from the sub-planet point, with synodic rotation rate of the planet determining only the sign of this angle, but not its magnitude ǫ = (2Q) −1 (Murray & Dermott 1999). Applying the same approach to lunar librations would make the amplitude of librations (proportional to eccentricity in absence of other dynamical excitation) irrelevant for the resultant force on Earth and the opposing torque on lunar figure. In this situation, the magnitude of eccentricity damping would be independent of eccentricity, which is not physically correct. Because of this, our integrator assumes linear frequency dependence of the tidal response for rotational rates comparable to the mean motion, as described inĆuk et al. (2016). The fact that we cannot use a constant-Q approach to calculate the lunar librational motion is not in contradiction with the more traditional approach of Tian & Wisdom (2020), who use a constant Q model of tides which assumes synchronous rotation and explicitly involves lunar orbital eccentricity. Figure 3 shows a simulation of the evolution of an Earth-Moon system with an initial obliquity of θ = 65 • through the LPT using our improved numerical integrator. The initial AM of the Earth-Moon system was 2.37 times the current value, so that the ecliptic AM is the same as it is now. In this example we used constant tidal parameters for Earth and the Moon, with Q E /k 2E = Q M /k 2M = 100, where k 2 and Q are respectively the tidal Love number and the tidal quality factor (Murray & Dermott 1999), and subscripts E and M refer to Earth and the Moon. At the end of the simulation, Earth has a low obliquity θ = 5 • , the Moon has an inclination of just over 30 • , and the AM is close to that of the present-day system. The small excess in total AM is associated with the lunar inclination which can be damped during the lunar Cassini State transition (Chen & Nimmo 2016;Ćuk et al. 2016)  Total AM (now=1) Figure 3. Evolution of the Earth-Moon system starting with a 65 • obliquity and 2.37 times its present-day total AM (chosen to match the present ecliptic AM, L = 0.35LC ; see supplementary text). Tidal parameters were set to QM /k2M = QE/k2E = 100, and the figures of Earth and the Moon were assumed to be in hydrostatic equilibrium. The left-hand panels plot (top to bottom) the lunar semimajor axis, eccentricity, and inclination with respect to the ecliptic. Also top to bottom, the right hand panels plot the system's AM (solid curve is total AM, nearly-horizontal line ecliptic AM), spin AM of Earth (roughly equivalent to number of rotations per day), and Earth's obliquity relative to the ecliptic. The solid red lines plot the averaged tilt to the ecliptic of the lunar Laplace plane (Tremaine et al. 2009, Eq. 22). The beginning and the end of each secular resonance (see description in text) is marked by vertical lines, and the resonances are numbered in the middle panels: 1. the inner 3/2 secular resonance, 2. Quasi-Kozai-Lidov resonance, 3. the outer 3/2 secular resonance.

DYNAMICAL STAGES OF THE LAPLACE PLANE INSTABILITY
There are two distinct phases to the LPT in Fig. 3. Before about 40 Myr into the simulation, the node of the lunar orbit with respect to the ecliptic is librating around the line of the equinoxes, rather than circulating (Fig. 4,middle panel). At 40 Myr, when the lunar inclination with respect to the ecliptic reaches zero, the lunar line of nodes stops being tied to the line of equinoxes and starts circulating. This transition was identified as critical by Tian & Wisdom (2020) and they referred to this event as the termination of the LPT. While an important dynamical event, we must point out that the timing of the transition from libration to circulation of the angle between the lines of lunar nodes and Earth's equinoxes (which we designate h = Ω − Ω spin after Tian & Wisdom (2020)) depends not only on the usual determinants of the Laplace Plane (Earth's mass, obliquity, oblateness etc.) but also on the Moon's free inclination. The circulation in h simply indicates that the tilt of the lunar Laplace Plane to the ecliptic (i.e. lunar forced inclination) is now smaller than the Moon's free inclination. However, the Laplace Plane Transition is certainly not over, as the Moon's instantaneous Laplace plane at this point is still inclined to the ecliptic by about 30 • (red line, Fig. 3).
The transition between libration and circulation in h leads to a temporary drop in lunar eccentricity (Fig. 3, middle left panel). Excitation of eccentricity is generally due to secular resonances present during the LPT, and this event represents the end of the second of three major secular resonances that dictate the dynamics of the orbital evolution shown in Fig. 3. We discuss the other secular resonances below. This middle resonance lasts from about 13 Myr to 40 Myr, and the top panel in Fig. 4 shows that during this period (labeled "2" in Figs argument of perigee with respect to the ecliptic, ω, is librating, indicating a form of Kozai-Lidov resonance (Lidov 1962;Kozai 1962). Unlike in the usual Kozai-Lidov problem, where the outer perturber dominates, here the precession is mostly driven by the Earth's oblateness. While h librates, the argument of pericenter referenced to the ecliptic is ω = ̟ − Ω spin , where ̟ and Ω spin are the longitudes of the lunar perigee and Earth's vernal equinox, respectively. SinceΩ spin <<Ω eq (where Ω eq is the lunar node with respect to Earth's equator), this "Quasi-Kozai-Lidov" (QKL) resonance requires thatω eq ≈ −Ω eq , a kind of secular resonance previously studied in the context of artificial satellite dynamics (Rosengren et al. 2015). The Sun, through the Kozai-Lidov mechanism keeps exciting the lunar eccentricity, effectively removing AM from the system, while preserving the ecliptic AM. Given the criteria for this resonance, it will necessarily break once the nodal parameter h starts circulating (middle panel in Fig. 4). The QKL resonance found in our simulations differs from other variations of the Kozai-Lidov resonance that are observed when there is also an inner perturber present in the system (the original Kozai-Lidov resonance relies exclusively on the outer perturber). When there is little or no obliquity of the central body (or the inner perturber is coplanar with the outer), the definition of the argument of pericenter are the same in equatorial and ecliptic reference frames. In that case, a single threshold for instability (i.e., excitation of eccentricity for initially circular orbits) connects Molniya orbits at 63.4 • (close to the planet) and the Kozai-Lidov threshold at 39.2 • far from the planet (Tremaine & Yavetz 2014;Petrovich et al. 2020). Since our inner perturber (Earth) has a large obliquity, our secular dynamics is more complex. QKL resonance in the present paper relies on the longitude of the node being almost stationary, as the orbital plane is dominated by the tilt of the planet's equatorial plane. In this configuration, the argument of pericenter (with respect to the Sun) precesses at rates similar to the longitude of pericenter, and the stationary ω with respect to the Sun can be reached at lower free inclinations than in the low obliquity case. Once h begins to circulate due to decreasing obliquity, this form of instability disappears without free inclination crossing any critical values. Similarly, while our QKL resonance is in the regime where precession is dominated by Earth's oblateness, the resonant interaction is with the Sun (as our Earth is azimuthally symmetric). This sets the QKL apart from the outer Kozai-Lidov resonance (Naoz et al. 2017) in which the resonance arises with interactions with an eccentric inner perturber.
After a brief period (∼ 5 Myr in the example in Fig. 3) of low inclination, the system enters another secular resonance. This final resonance is absent in the results of Tian & Wisdom (2020) and makes our simulated evolution diverge fundamentally from theirs. Lasting from about 44 Myr to 114 Myr (labeled "3" in Figs. 3 and 4), this secular resonance excites lunar eccentricity and therefore arrests further outward tidal evolution. In this dynamical state, the lunar inclination remains approximately constant (i ≈ 33 • ), while Earth's obliquity decreases, while conserving ecliptic AM (Fig. 3). We find that this resonance is caused by the periods of variation of lunar inclination (due to perturbations from Earth's oblateness) and lunar eccentricity (due to the Kozai-Lidov mechanism) being exact multiples of each other. For each complete cycle in lunar inclination (driven by the circulation of the nodal angle h) there are three cycles in eccentricity (dependent on the angle 2ω). Therefore, the resonant argument is k = 3h + 2ω or k = 3Ω − 3Ω spin + 2ω. This secular resonance requires thatω ≈ −(3/2)Ω, which is consistent with the lunar inclination in our simulations (see Discussion). The bottom panel in Fig. 4 plots the resonant argument k during the simulation shown in Fig. 3, with libration of k confirming the presence of the resonance. This resonance was found bý Cuk et al. (2016), but its effects were unfortunately entangled with the erroneous obliquity drift in their simulations. Our corrected simulations here confirm that this secular resonance takes non-ecliptic AM from Earth's rotation and, through the Kozai-Lidov mechanism, passes it to the heliocentric orbit of Earth. Fig. 3 shows that this resonance can last until Earth's obliquity is almost completely removed, at which point the resonance breaks.
The earliest, short-lived secular resonance seen within the first 10 Myr in Fig. 3 (labeled "1") is related to the late secular resonance we just described, likewise requiring the precession rate of the argument of perigee to be about -150% of that of the node. However, these orbital elements need to be measured relative to Earth's equator, rather than the ecliptic, as the perturbations from Earth's oblatenesss still dominates the lunar orbit at this early time. The relevant resonant argument is j = 3Ω eq + 2ω eq − 2Ω spin . Note that j does not obey D'Alembert's rules (Murray & Dermott 1999) because it is written using variables from two different coordinate systems 1 . Figure 5 shows a short snapshot Earth's obliquity (°) Figure 5. A short simulation illustrating the state of the system 10 Myr into the simulation shown in Fig. 3, illustrating the inner 3/2 secular resonance ( Table 1). The top left panel plots the resonant argument j = 3Ωeq + 2ωeq − 2Ωspin, where Ωeq and ωeq are the lunar longitude of the ascending node and argument of perigee in the equatorial reference frame, while Ωspin is the longitude of vernal equinox. Left-hand panels show Earth's obliquity and lunar inclination relative to the ecliptic.
of the first secular resonance 10 Myr into the simulation, illustrating the libration of j and the coupling between different dynamical quantities. To avoid confusion, from here on we will refer to the first and third secular resonances encountered during our example lunar orbital evolution as the inner and outer 3/2 secular resonances, respectively, while we will continue to refer to the second resonance as the quasi-Kozai-Lidov (QKL) resonance (see Table 1).

HIGH ECCENTRICITY EPISODES
The sequence of secular resonances encountered in the simulation shown in Fig. 3, ending with a low-obliquity Earth, is typical of most of our integrations that have explored starting with a compact Earth-Moon system that has the present-day ecliptic AM and obliquities of 62 • or 65 • , with various tidal properties of Earth and the Moon. For some integrations restricted to initial obqliuities of 61 • − 62 • , relatively early in the LPT the lunar orbit is captured in a secular resonance that forces lunar eccentricity to become exceptionally high (e > 0.5), as seen in Fig. 6. In such scenarios, extremely strong satellite tides quickly collapse the lunar semimajor axis, leading to simulations crashing as our fixed timestep cannot handle much smaller a. When repeated with a smaller timestep, simulations do not suffer crashes, but the Moon effectively re-starts its tidal evolution at a smaller a with completely damped eccentricity and inclination. As some of the total AM has already been lost and ecliptic AM is conserved, Earth's obliquity is now lower than it was initially, and the system can evolve without encountering the Laplace Plane instability. The end result would be a high AM, high obliquity system with a non-inclined Moon, inconsistent with the present-day Earth-Moon system.
We do not expect such an evolutionary track is a realistic possibility for the Earth-Moon system. An increase in the tidal dissipation within the Moon avoids the secular resonance found to produce very high eccentricities, resulting in an orbital evolution like the one shown in Fig. 7 and leading to a low-obliquity Earth. High eccentricity leads to intense tidal heating of the Moon and, as shown by Tian et al. (2017), tidal heating affects the Moon's tidal properties, typically leading to the increase in the parameter A, which measures the relative strength of tides within the Moon and Earth (Canup et al. 1999): where M and R are mass and radius, and subscripts E and M refer to Earth and the Moon. This dependence of tidal properties on tidal dissipation would act to moderate or eliminate high-e episodes of lunar tidal evolution (Tian et al.  Figure 6. Part of an integration that experienced rapid collapse of semimajor axis. Gray line plots the initial simulation with initial obliquity of θ = 62 • , the correct ecliptic AM, and QM /k2M = QE/kE2 = 100. At the point in orbital evolution when the quasi-Lidov-Kozai resonance is encountered by other simulations, this integration appears to be captured into a very high-e secular resonance, leading to rapid collapse of semimajor axis. Changing the lunar tidal dissipation to QM /k2M = 50 (black line) avoids this instability. We argue that this instability is probably not relevant for the real Earth-Moon system, as it requires the Moon to maintain the same tidal parameters despite having very high eccentricity and so experiencing intense tidal heating.
2017). Therefore, we believe that cases of lunar semimajor axis collapse due to high eccentricities are an artifact of our use of constant tidal properties for the Moon, rather than a plausible orbital evolution pathway. However, we think that the use of constant tidal properties is generally justified in other cases, such as evolutions shown in Figs. 3 and 7, as the Moon does not reach very high eccentricities during those simulations. The high-eccentricity evolution shown here is restricted to cases with Earth's initial obliquities in the 61 • −62 • range. For θ initial = 61 • we find instability using both A = 10 and A = 20, while for θ initial = 62 • high-e crash happens only when A = 10. Apart from the above discussion of how variable tidal properties may affect this instability, we note that Tian & Wisdom (2020) integrated systems with θ initial = 61 • without suffering a high-e crash. Clearly this phenomenon is at least somewhat dependent on the tidal model used, and therefore more work is needed to evaluate its intrinsic relevance to the evolution of a high-obliquity Earth-Moon system.

OBLIQUITY OF EARTH AFTER THE LAPLACE PLANE TRANSITION
The evolution through the three secular resonances plotted in Fig. 3 produces a system with an AM similar to that of the present-day Earth-Moon system but with a very low obliquity of Earth. Further tidal evolution would have made Earth's obliquity increase as the Moon absorbed more and more ecliptic AM from Earth (Touma & Wisdom 1994b). Based on the constraints of the current lunar orbit and terrestrial rotation rate, the post-LPT obliquity of Earth would need to have been in the range 10 − 15 • (Touma & Wisdom 1994b;Rubincam 2016), somewhat higher than that found at the end of the simulation in Fig. 3. Simple extrapolation from work of Touma & Wisdom (1994b) would imply that this simulations would lead to present-day Earth obliquity of ≤ 10 • . We find that the dynamics of the end of the LPT are quite complex and dependent on the tidal properties of Earth and the Moon and that different choices of tidal properties can lead to systems with higher final obliquities. Figure 7 shows a tidal evolution similar to that in Fig. 3, except that the Earth begins with somewhat lower obliquity and AM (although with the same ecliptic AM), and that the lunar tidal response is twice as strong (A = 20 vs. A = 10). The initial conditions do not seem to affect the final outcome, as simulations with the current ecliptic 10Ćuk et al. Total AM (now=1) Figure 7. Evolution of the Earth-Moon system starting with a 62 • obliquity and 2.13 times its current total AM. We used constant tidal parameters QE/k2E = 100 and QM /k2M = 50 for Earth and the Moon, respectively. The panels plot the same quantities as in Fig. 3. The main difference from the simulation in Fig. 3 is that the Moon's orbit is never librating within the outer 3/2 secular resonance, but is circulating outside the resonance (Fig. 8), resulting in a higher obliquity of Earth (θ ≈ 20 • ) at the end of the simulation.
AM tend to pass through a "keyhole" in phase space at the end of the QKL resonance (labeled "2" in Fig. 3 and 4), with a total AM about 50% greater than the present-day and obliquity θ ≈ 50 • . This "keyhole" is equivalent to the end of LPT proposed by Tian & Wisdom (2020), except that in our simulations it is followed by the outer 3/2 secular resonance. Fig. 7 shows that this resonance is significantly affected by the lunar tidal response, with lunar eccentricity periodically going through zero, indicating a lack of libration. Figure 8 shows a snapshot of this near-resonance at 50 Myr into the simulation shown in Fig. 7. Fig. 8 indicates that the same resonant angle k as in Fig. 4 is responsible for the excitation of lunar eccentricity, despite not being in libration. Given that the system is not as deep in the outer 3/2 secular resonance in Fig. 7 as in Fig. 3, the Moon exits the resonance when Earth's obliquity is about θ ≈ 20 • . This value is probably too high to match the present Earth-Moon system and would result in 30 • obliquity for Earth. However, a full numerical study of the tidal evolution out to 60R E (updating Touma & Wisdom 1994b) is probably necessary in order to quantify the correspondence between post-LPT obliquities and the current day values for Earth-Moon system histories with large-scale lunar inclination damping. The lowest initial obliquity for which the Laplace Plane instability occurs is about θ = 61 • , while the upper limit of obliquities is determined by the plausible upper limit on the total AM of the post-giant-impact Earth. Obliquities larger than θ ≈ 72 • would require initial AM that are typically not found in the final states of numerical giant impact simulations (Lock & Stewart 2017;Lock et al. 2018;Ćuk & Stewart 2012;Canup 2012;. Our integrations indicate that the final obliquity of Earth depends on the tidal parameters more than the initial obliquity and AM. In particular, early indications are that the tidal response of Earth and the Moon during and after capture of the Moon into the outer 3/2 secular resonance determines the final obliquity of Earth. The dynamics at resonance capture determine whether further evolution will involve libration (Fig. 4)  Earth's obliquity (°) Figure 8. A short integration illustrating the dynamical state of the system 50 Myr into the simulation shown in Fig. 7. The top left panel plots the resonant argument k = 3Ω − 3Ωspin + 2ω, where Ω and ω are lunar longitude of the ascending node and argument of perigee, and Ωspin is the longitude of the vernal equinox (in an ecliptic reference frame). Other three panels plot Earth's obliquity to ecliptic (top right), and the Moon's eccentricity (bottom left) and inclination relative to the ecliptic (bottom right). While the system spends most of the time at k ≈ 0 • , the argument k is circulating, rather than librating as in Fig. 4. 8) of the resonant argument k. We will explore the sensitivity of the evolution to tidal parameters further in future work. The fact that the reconstructed obliquity of Earth (Touma & Wisdom 1994b) is within the range of our early results (Figs. 3, 7) gives us confidence that a high-obliquity pathway to today's exact configuration will be found.

DISCUSSION
In the previous sections, we have demonstrated that an initially high-obliquity, high-AM Earth-Moon system (θ > 61 • ) naturally evolves into a state with much lower AM, a low obliquity of Earth and a lunar inclination of i ≈ 30 • . We confirm the finding of Tian & Wisdom (2020) that the ecliptic AM is conserved during the LPT and agree with them that the original simulations ofĆuk et al. (2016) suffered from significant numerical errors. However, our findings of low final obliquity and AM, consistent with the present Earth-Moon system, directly contradict the conclusions of Tian & Wisdom (2020) that the Laplace Plane instability cannot lead to the system's present configuration.
Comparing our simulations with those of Tian & Wisdom (2020), it is evident that the crucial difference in our simulations is the presence or absence of the dynamical feature we label the outer 3/2 secular resonance, in which the resonant argument k either librates (Fig. 4) or is at the boundary between circulation and libration (Fig. 8). In order to understand when this resonance is encountered, we explored where it is located in orbital element space. While the dynamics of a satellite in the LPT can be very complex (Tremaine et al. 2009;Tamayo et al. 2013), Figures 3 and  7 show that the outer 3/2 secular resonance is encountered at a point where the Sun is the principal perturber on the lunar orbit. This is the well-known Kozai-Lidov regime (Lidov 1962;Kozai 1962) and we can use the following expressions for the orbital precession rates (Innanen et al. 1997;Carruba et al. 2002): where τ is dimensionless time, and all orbital elements are in an ecliptic coordinate system. Since e ≃ 0.1 during much of the duration of the outer 3/2 secular resonance, while sin i ≃ 0.5, we will ignore terms involving e 2 but keep those with sin 2 i. We will also assume that ω is circulating faster than the resonant argument k (Figs. 4 and 8), so that < sin 2 ω >= 0.5, which is not true during Kozai-Lidov resonance but is acceptable here, as ω has to circulate when h is circulating and k is librating. Therefore, if we assume that Earth's axial precession is much slower than the Moon's orbital precession, the condition for the libration of k becomesω = −(3/2)Ω. Combining the above equations, we find that 5 cos 2 i − 3 cos i − 1 = 0. The solution to this equation is cos i = (3 + √ 29)/10, or i ≈ 33 • . This result matches the average inclination of the Moon while in the outer 3/2 secular resonance in Figs. 3 and 7 rather well.
The lunar inclination at the end of the simulations shown in Tian & Wisdom (2020) Fig. 2 is only slightly larger than that expected during the outer 3/2 secular resonance, i ≈ 34 • . While this difference in inclination is apparently sufficient to avoid the secular resonance, the Moon has a similarly high inclination in our integrations at the end of the QKL resonance and therefore this cannot be the cause of the different results of our studies. All our simulations that feature the QKL resonance have the Moon subsequently experience either resonant or near-resonant evolution along the 3/2 outer secular resonance, evolving the Earth to low obliquities. However, in the simulations of Tian & Wisdom (2020), the Moon avoids the 3/2 outer secular resonance because its inclination remains above 33 • despite the expectation of inclination damping (Chen & Nimmo 2016). The limited inclination damping is the result of Tian & Wisdom (2020) using a constant figure of the Moon, based on the hydrostatic shape at the distance of 6 R E . This assumption makes the Moon unrealistically non-spherical at the end of QKL resonance, and makes the model of Tian & Wisdom (2020) underestimate by orders of magnitude the lunar obliquity tides and the extent of the resulting inclination damping.
The lunar inclination in our Figs. 3 and 7 is constant only while the Moon is in the outer 3/2 secular resonance. Once the resonance is broken (at the very end of both simulations), lunar inclination decreases noticeably over time. This decrease is due to obliquity tides and is rather dramatic as we assumed strong lunar dissipation in these simulations (A = 10 and A = 20 similar to that used by Tian & Wisdom (2020)). Inclination damping is related to the tidal evolution rate as (Chyba et al. 1989): where θ M is the lunar obliquity, and we assumed thatȧ is dominated by tides raised on Earth by the Moon. In our simulations, the Moon is in Cassini state 1 around a = 23R E , and (θ M /i) ≈ 0.25. Using analytical estimates for the lunar figure (Garrick-Bethell et al. 2006), axial precession (Ward 1975) and both solar-induced (Innanen et al. 1997) and oblateness-induced (Danby 1992) nodal precession, we find axial and nodal precession periods of about 20 yr and 80 yr, respectively, consistent with this θ M /i ratio. Therefore, we would expect sin i −1 (di/dt) ≈ −0.3(ȧ/a), so when the Moon evolves from 18 to 24R E at the end of simulations in Fig. 2 of Tian & Wisdom (2020), we would expect lunar inclination to decrease noticeably, by several degrees, which is not seen in their simulations. Note that this calculation does not depend on assumptions about lunar shape, as the Moon's current oblateness and the oblateness expected from hydrostatic equilibrium are equal at about 23R E . To illustrate the importance of taking into account the correct lunar figure, we ran a control simulation using the same fixed lunar shape as Tian & Wisdom (2020). For initial conditions we used the state of the system 42 Myr into the simulation shown in Fig. 3, which is between residence in the QKL and the outer 3/2 resonances. Results are plotted in Fig. 9; the Moon maintains a constant orbital inclination and does not encounter outer 3/2 resonance as it moves out beyond 25 R E .
We tentatively conclude that the main difference between our model and that of Tian & Wisdom (2020) is their apparent lack of significant inclination damping due to obliquity tides upon leaving the QKL resonance, which makes the encounter with the outer 3/2 secular resonance much less likely. In our integrations shown in Figs 3 and 7, the lunar inclination is above i = 33 • when the QKL resonance is broken, as the QKL resonance requires slowerω (and therefore higher i) than the 3/2 outer secular resonance. Through the action of obliquity tides, the lunar inclination steadily decreases and enables capture into (or just outside) the outer 3/2 secular resonance. We note that this resonance is present in some of Tian & Wisdom (2020) simulations included in their Supplementary material. Since the outer 3/2 secular resonance is driven primarily by solar perturbations, it is not finely sensitive to Earth's spin and is therefore also present in cases with a lower Earth-Moon system AM. In Tian & Wisdom (2020) Fig. S1, the eccentricity of the Moon is weakly excited after h starts circulating, most likely because of proximity to the outer 3/2 secular resonance. While the effect on the total AM is small, there is a noticeable decrease in Earth's obliquity while the resonance is active. Additionally, some of the final obliquities of Earth in Fig. S3 of Tian & Wisdom (2020) are notably smaller than the corresponding lunar inclinations, indicating that some obliquity-reducing mechanism must have been acting after h started circulating.  Fig. 3 (gray). The panels plot the same quantities as in Fig. 3. As in the simulations of Tian & Wisdom (2020), the Moon does not encounter the outer 3/2 secular resonance as it evolves outward on a circular orbit.
The evolution of the Earth-Moon system's angular momentum, Earth's obliquity, lunar inclination, lunar tidal response and lunar shape are all closely connected when the Moon is at 20 − 25R E . A deformable, and therefore highly dissipative, Moon is required for AM loss and a reduction in Earth's obliquity through the outer 3/2 secular resonance. However, a continuing strong tidal response of the Moon as it approaches the Cassini state transition (at 30-34R E ; Ward 1975;Ćuk et al. 2016, 2019 would completely remove all of the lunar inclination, which is in conflict with the present day lunar orbital tilt (Chen & Nimmo 2016) 2 .Ćuk et al. (2016 showed that a Moon with close to present-day tidal properties can damp the lunar inclination from values around i = 30 • at 25R E to the present day value, i = 5 • . Therefore, the Moon must have transitioned from a very strong tidal response (expected from a fluid body) to a state close to its current solid-body response around 25R E , or soon after the end of the outer 3/2 secular resonance, according to our results. We also note that the Earth-Moon distance at which the Moon's current shape was frozen in is likely in the a = 23 − 30R E range (Garrick-Bethell et al. 2006), but the question of lunar shape is notoriously complex, and any connection to LPT will need to be studied separately. In any case, our model makes definite predictions about the Moon's early orbital, geophysical and thermal evolution, and we hope that further work on lunar geological history and paleomagnetism will be able to further constrain early lunar tidal evolution.
Finally, the requirement of initial obliquity in the 61 − 71 • range is actually more likely by a factor of ≈ 5 (assuming random distribution on a sphere) than the initial obliquity of < 10 • . The canonical model (Canup & Asphaug 2001) requires θ 0 = 10 • , but even lower obliquities are needed to match AM loss through the evection resonance (Ćuk & Stewart 2012). The constraint of conservation of ecliptic AM requires that systems with obliquities of 61 − 71 • have total initial AM of 2.1-3.1 times the present-day AM of the Earth-Moon system. It has been shown that many such high-AM impacts can produce post-impact structures with significant amounts of mass and AM beyond the Roche limit and with silicate vapor pressures at the Roche-limit of tens of bars, consistent with the moderately-volatile depleted composition of the Moon (see Figures 14 and 17 in Lock et al. 2018).

CONCLUSIONS
We confirm the finding of Tian & Wisdom (2020) that the ecliptic angular momentum of the Earth-Moon system is approximately conserved during the Laplace-plane transition. However, contrary to their claims, we find that Earth-Moon system with a high initial obliquity (θ > 61 • ) and AM (combined to produce the correct ecliptic AM) naturally evolves into a state with the current total AM and a low obliquity of Earth. The crucial dynamical mechanism that allows for the evolution to low obliquities is a secular resonance (the "outer 3/2 secular resonance") that is invariably encountered by the (then highly inclined) lunar orbit during the later stages of the Laplace Plane transition. The capture into this secular resonance requires some inclination damping, which is expected for realistic obliquity tides, but is absent in the results of Tian & Wisdom (2020) due to their treatment of the lunar figure.
We conclude that the hypothesis of an initially highly tilted Earth with a high AM is viable and offers much promise in explaining the implied common source for terrestrial and lunar materials (Ćuk & Stewart 2012;Canup 2012;Lock et al. 2018), the moderately volatile depleted composition of the Moon (Lock et al. 2018), and subsequent AM loss. It also naturally explains the current lunar inclination while allowing for large scale core-mantle friction during the lunar Cassini State transition, potentially powering the ancient lunar dynamo (Ćuk et al. 2019).
MĆ is supported by NASA Emerging Worlds Program award 80NSSC19K0512. SJL gratefully acknowledges funding from NSF through award EAR-1947614. Comments by Daniel Tamayo and two anonymous reviewers greatly improved the manuscript. The authors would like to thank Jihad Touma and the faculty and staff of the American University in Beirut for hosting a workshop on lunar formation in October 2019, which greatly advanced our understanding of early lunar evolution.