The Lunar Fossil Figure in a Cassini State

The present ellipsoidal figure of the Moon requires a deformation that is significantly larger than the hydrostatic deformation in response to the present rotational and tidal potentials. This has long been explained as due to a fossil rotational and tidal deformation from a time when the Moon was closer to Earth. Previous studies constraining the orbital parameters at the time the fossil deformation was established find that high orbit eccentricities (e ≳ 0.2) are required at this ancient time, which is difficult to reconcile with the freezing of a fossil figure owing to the expected large tidal heating. We extend previous fossil deformation studies in several ways. First, we consider the effect of removing South Pole−Aitken (SPA) contributions from the present observed deformation using a nonaxially symmetric SPA model. Second, we use the assumption of an equilibrium Cassini state as an additional constraint, which allows us to consider the fossil deformation due to nonzero obliquity self-consistently. A fossil deformation established during Cassini state 1, 2, or 4 is consistent with the SPA-corrected present deformation. However, a fossil deformation established during Cassini state 2 or 4 requires large obliquity and orbit eccentricity (ϵ ∼ 68° and e ∼ 0.65), which are difficult to reconcile with the corresponding strong tidal heating. The most likely explanation is a fossil deformation established during Cassini state 1, with a small obliquity (ϵ ∼ −0.2°) and an orbit eccentricity range that includes zero eccentricity (0 ≤ e ≲ 0.3).


Introduction
Planetary bodies are expected to be in a state near hydrostatic equilibrium, with larger bodies having smaller nonhydrostatic perturbations. However, the Moon is unique in that its present ellipsoidal figure requires a deformation that is significantly larger than the hydrostatic deformation in response to the present rotational and tidal potentials. This has long been suggested as due to a fossil figure preserving the Moon's past rotational and orbital state (Sedgwick 1898;Jeffreys 1915;Lambeck & Pullan 1980).
More recent studies consider the orbital parameters that can explain the fossil figure in more detail. Garrick-Bethell et al. (2006) found that the present lunar figure, as constrained by its degree-2 gravity, can be explained by a fossil figure assembled at semimajor axis a = 22.9 R ⊕ , where R ⊕ is Earth's radius, and orbit eccentricity e = 0.49 assuming synchronous rotation. They also found two solutions at a = 24.8 R ⊕ , e = 0.17 and a = 26.7 R ⊕ , e = 0.61, assuming 3:2 spin-orbit resonance. The model used to derive these solutions assumes that once the fossil figure is established, the Moon does not deform as the orbit expands and the rotational and tidal potentials change, which requires infinite rigidity after the fossil figure forms. Matsuyama (2013) developed an extended model that accounts for the effects of finite rigidity and the corresponding deformation and found solutions with similar eccentricities but smaller semimajor axes. The latter is expected because larger rotational and tidal potentials at the time the fossil figure forms are needed if the subsequent deformation is taken into account. Keane & Matsuyama (2014) considered the effects of removing mass anomalies from the present gravity field and found a significant contribution to the degree-2 gravity field associated with the South Pole−Aitken (SPA) basin. They found mass-anomaly-corrected fossil figure solutions with an initial eccentricity e ∼ 0.2 assuming synchronous rotation, and no self-consistent solutions assuming a fossil figure established during 3:2 spin-orbit resonance. Similar to the solutions of Matsuyama (2013), the required initial semimajor axis in Keane & Matsuyama (2014) depends on the elastic lithosphere thickness and rigidity, which determines its effective strength and ability to preserve a fossil figure (e.g., a ∼ 15 R ⊕ assuming a 25 km thick elastic lithosphere with a rigidity of 40 GPa).
In this paper, we extend previous studies by considering the effects of nonaxially symmetric SPA contributions and Cassini states. We derive SPA-corrected gravity constraints on the fossil figure in Section 2 using a nonaxially symmetric SPA model derived from the crustal thickness, and we consider the possible orbital parameters when the fossil figure forms in Section 3. An additional constraint can be used assuming that the Moon remains in a Cassini state as the orbital and rotational parameters, and the corresponding lunar shape, evolve. We develop a Cassini state model that incorporates a fossil figure and consider the effects of a fossil figure and obliquity feedback on Cassini states in Section 4. We combine the SPA-corrected gravity and Cassini state constraints to infer the possible orbital and rotational parameters when the fossil figure is established and discuss the evolution of the obliquity and gravity coefficients in Section 5. We end with conclusions and a discussion of the path forward in Section 6.

SPA Correction
The long-term lunar gravitational potential at an exterior point with spherical coordinates (r, θ, f) can be expanded in spherical harmonics as where G is the gravitational constant, M = 7.346 × 10 22 kg is the lunar mass, R 0 = 1738 km is a reference radius, Φ ℓm are dimensionless expansion coefficients, and we adopt the sign convention of geodesy and astronomy in which the gravitational potential is positive. The spherical harmonics Y ℓm are given by where P ℓm is the unnormalized associated Legendre function (e.g., Arfken & Weber 1995;Wieczorek 2015). In addition to the spherically symmetric contribution corresponding to ℓ = 0, the potential contains nonspherical contributions associated with rotational and tidal deformation, as well as mass anomalies. Below, we adopt the common notations Φ ℓm ≡ C ℓm if m 0, Φ ℓm ≡ S ℓm if m < 0, and C 20 ≡ − J 2 , and we rescale the gravity coefficients using the mean radius R = 1737 km instead of the reference radius R 0 . Matsuyama et al. (2021) found best-fit SPA ejecta thickness models by assuming an elastic lithosphere thickness and rescaling the inferred crustal thickness associated with the SPA ejecta using the observed gravity harmonics as constraints. In this paper, we take a different approach and assume that the observed SPA crustal thickness is the result of ejecta loading alone, with the elastic lithosphere thickness at the time of SPA loading as the unknown parameter. The former approach allows for nonloading contributions to the observed SPA crustal thickness; however, it requires prescribing the elastic lithosphere thickness. The latter approach does not allow for nonloading contributions; however, it allows us to infer the elastic lithosphere thickness at the time of SPA loading. Both approaches capture the complex geometry of the SPA ejecta (Matsuyama et al. 2021, Figure 2) and yield nearly identical solutions for the SPA-corrected gravity field, indicating that nonloading contributions are small.
We expand the SPA crustal thickness in spherical harmonics as where R = 1737 km is the lunar mean radius and H ℓm are dimensionless expansion coefficients. The gravity coefficients corresponding to the direct SPA ejecta contribution and the deformation associated with loading are where k ℓ L is the degree-ℓ load Love number, ρ L is the density of the load, and ρ is the mean density of the Moon. The load Love numbers describe the deformation in response to mass loading and depend on the interior structure. Since we are interested in the long-term deformation associated with SPA ejecta loading, the relevant Love numbers in Equation (4) are the long-term ones. We compute these Love numbers using the classical propagator matrix method (e.g., Sabadini et al. 2016) and assuming the interior structure parameters summarized in Table 1. The long-term load Love numbers are sensitive to the elastic lithosphere strength, which depends on its thickness and rigidity, and we consider a range of elastic lithosphere thicknesses.
We find best-fit SPA ejecta loading models, which corresponds to finding a best-fit elastic lithosphere thickness, using the long-wavelength gravity field as a constraint. That is, we minimize the misfit are the observed gravity coefficients, ℓm OBS dF are their uncertainty (Konopliv et al. 2013), and ℓ max is the expansion truncation degree. We exclude degree-2 gravity because it is dominated by rotational and tidal deformation. Figure 1 shows the normalized probability distribution (proportional to e SPA 2 c with SPA 2 c given by Equation (5)) for a range of ℓ max using the SPA crustal thickness expansion coefficients H ℓm of Matsuyama et al. (2021). When higher degrees are included, the probability distribution becomes wider owing to contributions of smaller mass anomalies not associated with SPA. We adopt the 1σ confidence region for the ℓ 3 max = solution, T e = 12.8 ± 3.4 km, as a constraint for the elastic lithosphere thickness at the time of SPA loading. Assuming the interior parameters summarized in Table 1 and T e = 12.8 km, k 0.83 L 2 = -, and the corresponding SPA gravity coefficients are given by Equation (4).
We are ignoring the effect of lithosphere thickness variations due to tidal heating, which can result in perturbations at degrees 2 and 4 (Beuthe 2013;Garrick-Bethell et al. 2014). However, the gravitational potential perturbations associated with these perturbations are small owing to the large isostatic compensation expected during a time of strong tidal heating. Regardless of their magnitude, degree-2 and degree-4 perturbations do not affect the best-fit SPA model because we use the best-fit solution limited to degree 3 alone.

Fossil Figure Solutions
Any interior region with long-term elastic strength can preserve a fossil figure. The theory outlined below is general enough to be applicable to arbitrary interior structures, including arbitrary regions with long-term rigidity, through the use of Love numbers. We assume a fossil figure preserved by an elastic lithosphere because it is the first region to develop Note. Parameters are based on mass, moment of inertia, gravity, and topography constraints (Matsuyama et al. 2016). We consider a range of elastic lithosphere thicknesses, T e . The shear modulus of all layers is set to zero except for the elastic lithosphere because we are interested in the long-term Love numbers.
long-term elastic strength as the Moon cools. We investigate the effect of increasing the elastic lithosphere thickness, which has a similar effect to assuming long-term elastic strength in other interior regions in addition to the elastic lithosphere, in Section B.1. For a given elastic lithosphere thickness, we remove the SPA contribution and the contribution associated with the deformation in response to the present rotational and tidal potentials (hereafter referred to as the equilibrium contribution) from the gravity field to obtain a gravity constraint on the fossil figure.
The SPA contribution is given by Equation (4), and the equilibrium contribution can be written as where k 2 is the long-term degree-2 Love number describing the long-term deformation due to gravitational forcing and ℓm á ñ  are time-averaged gravity coefficients of the sum of the rotational and tidal potentials (Appendix A). The time average is taken because we are interested in the long-term rotational and tidal response.
The emplacement of SPA ejecta alters the inertia tensor and induces a reorientation of the rotation and tidal axes relative to the surface, or true polar wander. We use the method of Matsuyama et al. (2014) to partition the inertia tensor into SPA, equilibrium, and fossil figure contributions, where the latter is associated with past rotational and tidal potentials. Removing the SPA and equilibrium contributions from the inertia tensor ( Table 2)  10 5 -. The difference arises because Keane & Matsuyama (2014) impose an axially symmetric geometry for the SPA ejecta while we capture the complex geometry of the ejecta using a model based on the SPA crustal thickness, as described above. Keane & Matsuyama (2014) removed SPA and other mass anomalies from the observed gravity field, while we only remove SPA contributions; however, this contributes little to the discrepancy because Keane & Matsuyama (2014) found that SPA contains 96% of the total degree-2 power from mass anomalies. E) prior to the emplacement of SPA ejecta, where the uncertainties correspond to the 1σ uncertainty for the elastic lithosphere thickness (T e = 12.8 ± 3.4 km). For comparison, the true polar wander solutions of Keane & Matsuyama (2014) for the rotation pole and tidal axis are (71.7 ± 2.0°N, 221.6 ± 2.1°E) and (14.3 ± 1.4°N, 1.9 ± 0.9°E), respectively. The rotation pole solutions are consistent with each other given the uncertainties, while there are large differences for the tidal axis solutions. As discussed above, the discrepancy arises owing to differences in the SPA ejecta model.
Assuming a fossil figure established when the rotational and tidal potentials were larger and a reference frame aligned with the rotation and tidal axes at that time, the fossil figure gravity coefficient can be written as where quantities with an asterisk are evaluated at the time the fossil figure is established. We are implicitly assuming that a fossil figure is established as an elastic lithosphere forms over a timescale much shorter than the orbital evolution timescale. That is, * k 2 is computed assuming an interior without an elastic lithosphere, and k 2 is computed assuming an interior with an elastic lithosphere. If an elastic lithosphere does not form, * k k 2 2 = , and the fossil figure contribution to the gravity field vanishes, as expected.
It is worth noting that the long-term Love numbers used here differ from the Love numbers describing the periodic deformation at the tidal forcing period. As discussed above for the load Love numbers, we compute the long-term tidal Love numbers using the classical propagator matrix method (e.g., Sabadini et al. 2016) and assuming the interior structure parameters summarized in Table 1. Assuming the best-fit solution described above with T e = 12.8 km, the values of the Love numbers are * k 1.44 2 = and k 2 = 1.23. The rotational potential contribution in Equations (6) and (7) is where w is the rotation rate (Appendix A). Assuming synchronous rotation and averaging over the orbital and  , Table S1). c SPA (Equation (4)) and equilibrium (Equation (6)) contributions assuming the interior structure parameters in Table 1, the best-fit elastic lithosphere thickness, and the corresponding Love numbers (T e = 12.8 km, k 0.86 L 2 = -, and k 2 = 1.23). Figure 1. Normalized probability distribution of the elastic lithosphere thickness at the time of SPA ejecta loading assuming different truncation degrees ℓ max (Equation (5)). The light-red area indicates the 1σ confidence region for the ℓ 3 max = solution.
precession periods, the tidal potential contributions are á ñ=´- to eighth order in eccentricity and arbitrary obliquity (Correia & Rodríguez 2013, Appendix A). For the largest eccentricity e = 0.7 considered below, the typical precision of Equations (9) and (10) is 0.7 9 ∼ 5%. The total time-averaged gravity coefficients in Equation (7)  The obliquity is usually positive (between 0°and 180°), but since the nodes of the Laplace, orbital, and equatorial planes on one another coincide in Cassini states, we allow negative obliquities for mathematical convenience. A positive obliquity means that the ascending node of the orbital plane on the Laplace plane coincides with the ascending node of the equatorial plane on the orbital plane. A negative obliquity (between -180°and 0°) means that the ascending node of the orbital plane on the Laplace plane coincides with the descending node of the equatorial plane on the orbital plane (or, equivalently, with the ascending node of the orbital plane on the equatorial plane). Equations (7)-(10) extend the fossil figure solutions of Matsuyama (2013) by taking the effects of nonzero obliquity into account, which is critical for considering Cassini states as described below. The effects of a finite rigidity interior are captured by the long-term Love numbers as described above. Combining Equations (7)-(9) and assuming synchronous rotation, zero eccentricity, and a homogeneous fluid interior at all times ( * k k 3 2 2 2 = = ), we can recover Equations (7)-(8) of Jankowski et al. (1989) for the hydrostatic moment of inertia differences under these assumptions. Figure 2 shows the orbital and rotational parameters at the time the fossil figure is established, as constrained by the present SPA-corrected degree-2 gravity coefficients. Here we focus on zero-eccentricity or zero-obliquity solutions, and solutions where both the eccentricity and obliquity do not vanish are considered in Section 5.
Assuming synchronous rotation and zero obliquity, a fossil figure formed at a semimajor axis a * ∼ 13 R ⊕ and eccentricity e * ∼ 0 can explain the present gravity coefficients (Figure 2(a)). Allowing for only requiring a small eccentricity is significant because the freezing of a fossil figure may not be selfconsistent with the large energy dissipation driven by high eccentricity. For comparison, assuming a fossil figure established during synchronous rotation, the solution of Keane & Matsuyama (2014) requires e * ∼ 0.2, and the solutions of Garrick-Bethell et al. (2006) and Matsuyama (2013) require e * ∼ 0.5.
If the fossil figure is established during 3:2 spin-orbit resonance and zero obliquity, only high-eccentricity solutions (e * ∼ 0.25 − 0.45) are possible (Figure 2(c)). For comparison, the 3:2 spin-orbit resonance solutions of Garrick-Bethell et al. (2006) and Matsuyama (2013) require an eccentricity e * ∼ 0.2 or ∼0.6, and Keane & Matsuyama (2014) find no selfconsistent solutions in this case. There are no solutions for 3:2 spin-orbit resonance and zero eccentricity because the fossil figure contribution to C 22 vanishes in this case (Equation (10)).
For the solutions assuming zero obliquity, the difference between our solutions and those of Garrick-Bethell et al. (2006) and Matsuyama (2013) arises as a result of the effects of removing the SPA contribution from the degree-2 gravity coefficients. The discrepancy between our solutions and those of Keane & Matsuyama (2014) is due to differences in the adopted SPA model ( Table 2). The low-obliquity, zero-eccentricity solutions (Figure 2(b)) and the low-eccentricity, zero-obliquity solutions (Figure 2(a)) are equivalent, both constraining the semimajor axis at the fossil figure formation time to ∼13 R ⊕ . An interesting case that has not been previously considered is a fossil figure established during high obliquity (ò * ∼ 100°), zero eccentricity, and synchronous rotation (Figure 2(b)). This case is intriguing because the required obliquity is similar to the Cassini state 2 obliquity (Ward 1975;Siegler et al. 2011), a minimum energy state expected owing to tidal energy dissipation (Peale 1969;Ward 1975). However, previous Cassini state studies do not take into account the effects of a fossil figure preserved by an elastic lithosphere and the feedback between obliquity and the gravity coefficients controlling the Cassini states. We develop a model that includes these effects in Section 4, and we consider possible fossil figure solutions assuming that the fossil figure is established during a Cassini state in Section 5.

Cassini States
The Moon is currently in a Cassini state, a configuration where the spin and orbit angular momenta are misaligned and precess synchronously (Colombo 1966;Peale 1969;Ward 1975). The Cassini state condition can be written as where i is the orbit inclination, p is the ratio of the orbital mean motion n to the orbit plane precession rate dΩ orb /dt, and c ≡ C/(MR 2 ) is a normalized principal moment of inertia (Peale 1969;Ward 1975;Bills & Nimmo 2008). We follow the approach of Chen & Nimmo (2016) to compute the orbit precession rate, which has contributions associated with torques from the solar tides and the rotational flattening of Earth. Equation (11) provides a constraint on the possible Cassini states for a rigid body with a constant inertia tensor and corresponding constant gravity coefficients. We are assuming that inertia tensor perturbations due to deformation in response to changes in the tidal potential occur over timescales that are significantly longer than the tidal forcing period, which is a reasonable assumption for the Moon. In this regard, the hydrostatic solutions considered below also involve deformation timescales that are significantly longer than the tidal forcing period because we are assuming hydrostatic deformations over geologic timescales.
In the calculations below we consider the evolution of the Moon to its current state using the observed degree-2 gravity coefficients as a constraint. Therefore, SPA formation and its effect on the degree-2 gravity field and Cassini states must be taken into account. The lunar age is constrained to 4.48-4.54 Gyr by crustal age constraints (Elkins-Tanton et al. 2011), and the SPA age is constrained to 4.25 ± 0.01 Gyr assuming that one of the Apollo Mg-suite rocks was excavated by SPA (Garrick-Bethell et al. 2020), which corresponds to SPA formation ∼250 Myr after the formation of the Moon. If the Apollo Mg-suite rock was excavated from the Procellarum KREEP Terrane (Elardo et al. 2020), its age would be associated with one of the nearside basins, making SPA older. However, this does not affect the fossil figure constraints considered below because SPA formation likely occurred after the fossil figure formed, as discussed in Section B.3.
The SPA formation semimajor axis is weakly constrained owing to uncertainties in the energy dissipation rate that controls the orbit recession rate. The lunar semimajor axis as a function of time can be written as (e.g., Murray & Dermott 1999, Equations (4.214) and (4.216) where R ⊕ , k 2⊕ , and Q á ñ Å are Earth's mean radius, degree-2 tidal Love number at the tidal forcing period, and average tidal quality factor, respectively. Assuming that the present k 2⊕ ∼ 0.3 (Munk & MacDonald 1960, Equation (10.3.2)) and Q ⊕ ∼ 12 (Murray & Dermott 1999, Equation (4.225)) have always prevailed, Q 12 á ñÅ , and the lunar age would only be ∼1.6 Gyr. This indicates that the present Q ⊕ is higher than the average value, which can be explained by lower tidal energy dissipation in Earth's ocean in the past owing to Earth's faster rotation rate (Webb 1982). Assuming a constant k 2⊕ = 0.3, a lunar age of 4.5 Gyr corresponds to an average tidal quality factor Q 33.7 á ñ = Å and an SPA formation semimajor axis a SPA = 38.7 R ⊕ if SPA formed 250 Myr after lunar formation. We adopt a SPA = 38.7 R ⊕ as a fiducial value and investigate the dependence of our results on a SPA in Section B.3.
Before SPA formation and the corresponding true polar wander, the degree-2 gravity coefficients in the pre-SPA principal axes reference frame can be written as where the first terms correspond to the fossil figure contributions (Equation (7)) and the second terms correspond to the equilibrium contribution (Equation (6)).
After SPA formation and the corresponding true polar wander, the degree-2 gravity coefficients in the post-SPA principal axis reference frame can be written as where the first term corresponds to the fossil figure contribution, the second term corresponds to the SPA contribution, and the third term corresponds to the equilibrium contribution (Equation (6)). Given the observed gravity coefficients C ℓm OBS at the present semimajor axis a OBS , the gravity coefficients after SPA formation can be written as Previous studies considering the spin axis evolution as the orbit expands and the Cassini states evolve have made some simplifying assumptions. Ward (1975) assumes a permanent fossil figure with gravity coefficients that remain constant as the orbit expands, which requires an interior with infinite rigidity. Chen & Nimmo (2016) assume a hydrostatic response to the rotational and tidal potentials, which is not consistent with the presence of a fossil figure. Siegler et al. (2011) assume a combination of a hydrostatic and a frozen-in component; however, they do not describe how this partitioning is determined. Our model provides a natural partitioning into frozen-in and hydrostatic components that depends on the elastic lithosphere strength, which is sensitive to the elastic lithosphere thickness and rigidity.
The deformation in response to the tidal potential depends on the obliquity, and the Cassini states determining the obliquity depend on this deformation. Therefore, an obliquity feedback between tidal deformation and the Cassini states arises. This feedback has been considered for the Cassini states of Triton assuming a hydrostatic deformation and a homogeneous interior structure . Our model takes this feedback into account by including the obliquity dependence of the gravity coefficients (Equations (9) and (10)) and allows for the presence of a fossil figure and arbitrary interior structures through the use of Love numbers. Figure 3 illustrates the effects of a fossil figure and obliquity feedback on Cassini states. Taking obliquity feedback into account reduces the obliquity amplitude for Cassini states 2 and 4 significantly. This effect is more pronounced for solutions assuming a hydrostatic deformation at all times. The presence of a fossil figure reduces the effects of obliquity feedback because it preserves a memory of the orbital and rotational state, including the obliquity, when the fossil figure was established. Cassini state solutions with a fossil figure differ significantly from solutions assuming a hydrostatic deformation at all times owing to the large differences in the corresponding degree-2 gravity coefficients. In particular, the annihilation of Cassini states 1 and 4, as well as the corresponding expected transition to Cassini state 2, occurs at a larger semimajor axis when the fossil figure is taken into account. The formation of a fossil figure increases the Cassini state transition semimajor axis from a ∼ 27 R ⊕ (Figure 3, orange curve) to a ∼ 34 R ⊕ (Figure 3, red curve).
As discussed above, the observed lunar deformation is not consistent with hydrostatic deformation given the present rotational and tidal potentials, requiring a fossil figure contribution. Nevertheless, hydrostatic deformation is relevant for the Cassini state evolution before the fossil figure is established. Therefore, we consider the hydrostatic Cassini state evolution in more detail in Figure 4. An interesting effect of obliquity feedback is the appearance of new Cassini states at small semimajor axes and high eccentricities. Without obliquity feedback, the obliquity magnitude of Cassini states 2 and 4 decreases as the eccentricity increases (Figure 4(a)). The effects of eccentricity on Cassini states are significantly more pronounced when obliquity feedback is taken into account: Cassini states 2 and 4 change significantly when e  0.5 (Figure 4(b)), and four new Cassini states appear when e  0.6 ( Figure 4(c)).
The evection resonance between the lunar orbit precession and the yearly orbital motion of Earth can increase the eccentricity to e ∼ 0.5 at a ∼ 4.6 R ⊕ , when the lunar orbit precession period is near 1 yr, and the eccentricity continues to increase after this resonance (Touma & Wisdom 1998). Therefore, the post-evection orbital parameters can be consistent with one of the new Cassini states; however, as discussed below, the new Cassini states are not compatible with the fossil figure constraint.

Fossil Figure Solutions Consistent with Cassini States
In Section 3 we find possible fossil figure solutions assuming zero obliquity (Figures 2(a), (c)) or zero eccentricity (Figure 2(b)). As discussed above, under the assumption of an equilibrium Cassini state, the obliquity is no longer a free parameter. In this section, we use the Cassini state model described above to find fossil figure solutions with an obliquity that is consistent with a Cassini state. We explore all the possible Cassini state solutions, including the new Cassini states described in Section 4. Figure 5 illustrates that a fossil figure formed during Cassini state 1, 2, or 4 is consistent with the SPA-corrected degree-2 gravity coefficients. On the other hand, a fossil figure established during Cassini state 3 or one of the new Cassini states described above is incompatible with the SPA-corrected degree-2 gravity coefficients. The Cassini state 1 solution requires a * ∼ 13 R ⊕ , e *  0.3 (Figure 5(a)), and ò * ∼ − 0.2°F igure 3. Cassini state solutions assuming hydrostatic deformation at all times or a combination of hydrostatic and frozen-in fossil figure components. Cassini states 1-4 are labeled as S1-S4. Hydrostatic deformation solutions are found by assuming an interior without an elastic lithosphere ( * k k 2 2 = in Equation (7)) at all times. Solutions without obliquity feedback are found by assuming zero obliquity in the gravity coefficients at all times (ò * = ò = 0 in Equation (9)). For the fossil figure solutions, we assume the best-fit elastic lithosphere thickness at the time of SPA loading, T e = 12.8 km (Section 2), and best-fit orbital and rotational parameters constrained by the SPA-corrected gravity coefficients, a * = 13.1 R ⊕ , e * = 0, and ò * = −0.16°(Section 5). In all cases, we assume synchronous rotation, the present orbit inclination (i = 5.1°), and zero eccentricity at all times.
( Figure 4(b)). The Cassini state 2 and 4 solutions require a * ∼ 3.8 R ⊕ , e * ∼ 0.65 (Figures 5(b) and (c)), and ò * ∼ 68° (  Figure 4(c)), with a significantly smaller range of possible parameter values. A fossil figure formed during Cassini state 2 or 4 is problematic for two reasons. First, it requires a very high eccentricity before it can get excited by the evection resonance at a * ∼ 4.6 R ⊕ (Touma & Wisdom 1998). Second, the combination of a large eccentricity, a large obliquity, and a small semimajor axis produces strong tidal heating that is difficult to reconcile with the formation of a fossil figure. Therefore, hereafter we focus on the Cassini state 1 solution.
Given the orbital parameters at the time the fossil figure is established, we can compute the subsequent evolution of the obliquity and gravity coefficients. Figure 6 shows a representative evolution assuming a fossil figure established during Cassini state 1 and zero eccentricity. As discussed above, assuming a fossil figure formed during Cassini state 1, an eccentricity e *  0.3 is consistent with the present SPAcorrected degree-2 gravity coefficients ( Figure 5(a)). Before the fossil figure is established, the spin axis and gravity coefficients follow the hydrostatic solution shown as blue curves ( * k k 2 2 = in Equation (9)). A fossil figure is established at a * = 13.1 R ⊕ , which corresponds to a Cassini state 1 obliquity ò * = −0.16°. The subsequent evolution follows the fossil figure solution shown as orange curves ( * k k 2 2 ¹ in Equation (9)). The obliquity transitions from Cassini state 1 to state 2 when the former disappears at a ∼ 34 R ⊕ . The formation of SPA perturbs the degree-2 gravity field (Figures 6(b) and (c)), which perturbs the Cassini state 2 (Figure 6(a)) and includes true polar wander.

Conclusions
We have developed a new lunar fossil figure model that extends previous models in several ways. First, we remove SPA contributions from the observed degree-2 gravity field using a nonaxially symmetric SPA model, which provides a new gravity constraint on the possible orbital and rotational parameters at the time the fossil figure is established. Second, we use the assumption of an equilibrium Cassini state as an additional constraint. To accomplish this, we construct a Cassini state model that takes the effects of a fossil figure and = in Equation (7)). Solutions without obliquity feedback are found by assuming zero obliquity in the gravity coefficients at all times (ò * = ò = 0 in Equation (9)). In all cases, we assume the present orbit inclination (i = 5.1°) and synchronous rotation at all times. Figure 5. Constraints on the orbital parameters at the time the fossil figure is established from the SPA-corrected gravity coefficients (J 2 = (1.32 ± 0.17) × 10 −4 and C 22 = (3.71 ± 0.35) × 10 −5 ). Solutions in panels (a)−(c) are found assuming the obliquity of Cassini state 1, 2, or 4, as labeled. In all cases, we assume the best-fit elastic lithosphere thickness at the time of SPA loading (T e = 12.8 km, Section 2), the present orbit inclination (i = 5.1°), and synchronous rotation. obliquity feedback into account, and we compute the Cassini state obliquity self-consistently given the SPA-corrected degree-2 gravity coefficients.
There are four previously unidentified high-obliquity Cassini states when the semimajor axis a  12 R ⊕ and the orbit eccentricity e  0.6 (Figure 4(c)); however, a fossil figure formed during one of these states is incompatible with the SPAcorrected gravity constraints. The most likely solution is a fossil figure established during Cassini state 1, when the semimajor axis a ∼ 13 R ⊕ , the eccentricity e  0.3, and the obliquity ò ∼ −0.16°. The low eccentricity and low obliquity in this case allow for the formation of a fossil figure during a time of low tidal energy dissipation. Assuming zero eccentricity, high-obliquity Cassini states (|ò|  68°) are unstable when an energy dissipation term is included in the equations of motion describing Cassini states (Gladman et al. 1996, Appendix B). The obliquity evolution of the most likely solution described above satisfies this stability condition (Figure 6(a)).
We assume the best-fit elastic lithosphere thickness at the time of SPA loading, T e = 12.8 km (Figure 1), and the present orbit inclination (i = 5.1°) in all the calculations above. Assuming a thinner/thicker elastic lithosphere thickness decreases/increases the required semimajor axis at the time the fossil figure forms (Appendix B.1, Figure 8), as expected owing to the ability of a thicker lithosphere to preserve a larger fraction of the initial fossil figure. The formation of the fossil figure at different semimajor axes when considering different elastic lithosphere thickness does not affect the semimajor axis and amplitude of the Cassini state transition (Appendix B.1, Figure 9) because this transition is determined by the degree-2 gravity coefficients and the fossil figure solutions are found by imposing the same SPAcorrected degree-2 gravity constraint. Varying the orbit inclination has a small effect on the possible orbital parameters at the time the fossil figure forms (Appendix B.2, Figure 10) and affects the amplitude of the Cassini state transition significantly, with larger transition amplitudes for larger orbit inclinations (Appendix B.2, Figure 11).
The elastic lithosphere thickness is expected to change as the orbital parameters and the corresponding tidal energy dissipation evolve and a magma ocean cools. Therefore, considering these effects requires a model that couples the thermal and orbital evolution. Chen & Nimmo (2016) developed a similar coupled thermalorbital model but including the effects of tidal dissipation in the magma ocean due to obliquity tides and using the present orbit inclination instead of the present fossil figure as a constraint. They found a rapid inclination damping due to large obliquity tidal heating associated with a Cassini state transition at a ∼ 30 R ⊕ , which places strong constraints on the possible orbital evolution. The Cassini state model of Chen & Nimmo (2016) does not take the effects of a fossil figure and obliquity feedback into account, which affects the Cassini state transition significantly (Figure 3). Our model can be used to extend previous coupled thermal-orbital models by taking into account the effects of removing SPA contributions from the present gravity field, the effects of finite rigidity on preserving a fossil figure, and the effects of a fossil figure and obliquity feedback on Cassini states. The large obliquity associated with the Cassini state transition produces strong tidal heating, damping the orbit inclination (Chen & Nimmo 2016;Hamilton et al. 2016). The feasibility of the fossil figure solutions discussed in this paper, as constrained by the present orbit inclination, is the subject of ongoing work. This work was supported by the National Aeronautics and Space Administration (NASA) under grant No. 80NSSC20K0334 issued through the NASA Emerging Worlds program. A portion of this research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (80NM0018D0004).

Appendix A Rotational and Tidal Potentials
In a reference frame with the rotation vector directed along the z-axis, the rotational potential at the lunar surface (r = R) is given by where w is the magnitude of the rotation vector. In terms of a spherical harmonics expansion similar to the one used in Equation (1), á ñ =   . Using Kaula's expansion (Kaula 1964;Peale & Cassen 1978), the gravitational potential at a surface point with spherical coordinates (R, θ, f) due to tides raised by Earth can be written as M ⊕ is Earth's mass, a is the lunar orbit semimajor axis, ω is the argument of pericenter, Ω is the longitude of the ascending node, M is mean anomaly, and Ψ is the angle between the axis of minimum moment of inertia and a reference inertial axis at the equator. As noted by Chyba et al. (1989), the argument of the F ℓmp functions is the lunar obliquity ò when describing the tidal potential acting on the Moon. Setting t = 0 as the time of pericenter passage, wt, w Y = W + ¢ + where tan tan cos w w ¢ =  , and where M and R are the lunar mass and mean radius, respectively. Rewriting this potential using the commonly adopted notation for the expansion coefficients, Assuming ℓ = 2 and synchronous rotation and averaging over the orbital and precession periods yields Equations (9) and (10) in the main text, which we cross-validated using the symbolic computation package TenGSHui for Mathematica (Trinh 2019).

A.1. Hydrostatic Figure
The gravitational field at the time the fossil figure was established can be expressed to the first order (in the polar and equatorial flattenings) in terms of the Love number * k 2 and is restricted to degree 2: * * * C k m m 2 2 2 = á ñ  . As the fossil figure freezes closer to Earth, second-order contributions * C m 2 D and * C m 4 to the degree-2 and degree-4 gravitational field, required to maintain the shape in hydrostatic equilibrium, become more important. These corrections can be computed using second-order theory of figures (see, e.g., Section 3.2 of Beuthe et al. (2016)). Figure 7 shows the first-order and second-order contributions to the gravitational potential coefficients as a function of the lunar semimajor axis assuming hydrostatic equilibrium and an interior without an elastic lithosphere. Second-order corrections do not exceed 1.2% and 3.9% for J 2 and C 22 , respectively, over the range of fossil figure semimajor axes (i.e., > 5 R ⊕ ) considered in the main text. For this reason, and since we use the best-fit SPA solution limited to degree 3 alone (ℓ 3 max = in Equation (5)), we choose to ignore second-order effects on the hydrostatic figure. This leaves the best-fit SPA model unaffected.

Appendix B Effects of Elastic Lithosphere Thickness, Orbit Inclination, and SPA Formation Semimajor Axis
In the main text, we assume the best-fit elastic lithosphere thickness at the time of SPA loading (T e = 12.8 km, Figure 1), the present orbit inclination (i = 5.1°), and SPA formation when the semimajor axis is a SPA = 38.7 R ⊕ . In this section we consider the effects of varying these parameters on the fossil figure solution, the corresponding constraints on the orbital parameters, and the subsequent evolution of the gravity coefficients and obliquity. Figure 8 illustrates how the orbital parameter constraints vary for different elastic lithosphere thicknesses and Cassini states given the SPA-corrected gravity coefficients. A thicker elastic lithosphere can retain a larger fossil figure. Therefore, smaller rotational and tidal potentials at the time the fossil figure forms are required as the elastic lithosphere thickness increases, which corresponds to larger semimajor axes. Varying the elastic lithosphere thickness does not modify the required orbit eccentricity at the time the fossil figure forms. Figure 9 shows the evolution of the obliquity and gravity coefficients for representative solutions assuming different elastic Figure 8. Constraints on orbital parameters at the time the fossil figure is established from the SPA-corrected gravity coefficients (J 2 = (1.32 ± 0.17) × 10 −4 and C 22 = (3.71 ± 0.35) × 10 −5 ) assuming different elastic lithosphere thicknesses (T e ) and Cassini states at this ancient time, as labeled. Additionally, we assume the present orbit inclination (i = 5.1°) and synchronous rotation in all cases.

B.1. Elastic Lithosphere Thickness
lithosphere thicknesses, a fossil figure established during Cassini state 1, and zero eccentricity based on the orbital constraints in Figure 8. Before the fossil figure is established, the spin axis and gravity coefficients follow the hydrostatic solution (blue curves). Assuming T e = 12.8 km (top panels), 25 km (middle panels), and 50 km (bottom panels), the fossil figure forms at a * = 13.1, 15.6, and 18.4 R ⊕ , respectively. The subsequent evolution (orange curves) is not sensitive to the assumed elastic lithosphere thickness because the fossil figure solutions are found using the same SPA-corrected degree-2 gravity constraints.
A fossil figure can be preserved by any interior region with long-term elastic strength. An interior model containing an additional region with long-term elastic strength has a similar effect to that of a thicker elastic lithosphere. Figure 10 illustrates how the orbital parameter constraints vary for different orbit inclinations and Cassini states given the SPA-corrected gravity coefficients. Under the assumption of an equilibrium Cassini state, the obliquity and orbit inclination are coupled (Equation (11)), and therefore the degree-2 gravity coefficients become dependent on orbit inclination (Equation (9)). For the fossil figure solutions considered here requiring small semimajor axes at the time the fossil figure is established (a *  20 R ⊕ ), the Cassini state obliquities do not vary significantly with the semimajor axis (e.g., Figure 4). This results in a small effect of varying the orbit inclination on the semimajor axis at the time the fossil figure forms. Figure 11 shows the evolution of the obliquity and gravity coefficients for representative solutions assuming different orbit inclinations, a fossil figure established during Cassini state 1, and zero eccentricity based on the orbital constraints in Figure 10. As discussed above, varying the orbit inclination has a small effect on the semimajor axis and eccentricity constraints at the time the fossil figure forms. In all cases shown in Figure 11, the fossil figure forms at a * = 13.1 R ⊕ (indicated by the orange dashed lines) regardless of the assumed orbit inclination. The orbit inclination has a significant effect on the transition from Cassini state 1 to state 2. The transition semimajor axis slightly decreases and the obliquity amplitude increases with increasing orbit inclination.

B.3. SPA Formation Semimajor Axis
Given the weak constraints on Earth's energy dissipation controlling the Moon's recession rate, we consider the effect of varying the SPA formation semimajor axis. Figure 12 shows the evolution of the obliquity and gravity coefficients for representative solutions assuming different SPA formation semimajor axes and a fossil figure established during Cassini state 1. The fossil figure semimajor axis and eccentricity constraints are not affected by the formation of SPA because we are assuming an SPA formation semimajor axis larger than the fossil figure formation semimajor axis. This assumption is justified by the extremely large Earth's average tidal quality factor,  Q 4 10 4 á ñǺ , required to form SPA before the fossil figure forms at a * = 13.1 R ⊕ (Equation (12)). SPA formation and the corresponding true polar wander perturb the degree-2 gravity coefficients, increase/decrease the obliquity amplitude if it occurs after/before the CST, and increase the CST semimajor axis slightly if it occurs before the CST (e.g., ∼36 R ⊕ in Figure 12(a) compared with ∼34 R ⊕ in Figures 12(d) and (g)). Figure 10. Constraints on orbital parameters at the time the fossil figure is established from the SPA-corrected gravity coefficients (J 2 = (1.32 ± 0.17) × 10 −4 and C 22 = (3.71 ± 0.35) × 10 −5 ) assuming different orbit inclinations (i) and Cassini states at this ancient time, as labeled. Additionally, we assume the best-fit elastic lithosphere thickness at the time of SPA loading, T e = 12.8 km, and synchronous rotation in all cases.