Excitation of Tumbling in Phobos and Deimos

TL;DR

This study uses mass-spring simulations to explore how Phobos and Deimos could have entered tumbling states, revealing that such tumbling can persist and influence their orbital eccentricities, with impacts possibly playing a key role.

Contribution

It demonstrates that crossing spin-orbit or mean motion resonances does not induce tumbling, and highlights the prolonged tumbling states and their effects on orbital eccentricity in Phobos and Deimos.

Findings

Tumbling states can persist for extended periods in Phobos and Deimos.

Tumbling can significantly reduce orbital eccentricity.

Impacts may be responsible for spin excitation rather than tidal resonance.

Abstract

Mass-spring model simulations are used to investigate past spin states of a viscoelastic Phobos and Deimos. From an initially tidally locked state, we find crossing of a spin-orbit resonance with Mars or a mean motion resonance with each other does not excite tumbling in Phobos or Deimos. However, once tumbling our simulations show that these moons can remain so for an extended period and during this time their orbital eccentricity can be substantially reduced. We attribute the tendency for simulations of an initially tumbling viscoelastic body to drop into spin-synchronous state at very low eccentricity to the insensitivity of the tumbling chaotic zone volume to eccentricity. After a tumbling body enters the spin synchronous resonance, it can exhibit long lived non-principal axis rotation and this too can prolong the period of time with enhanced tidally generated energy dissipation. The low orbital eccentricities of Phobos and Deimos could in part be due to spin excitation by nearly catastrophic impacts rather than tidal evolution following orbital resonance excitation.

Tables

Table 1. Table 1: Quantities for Phobos and Deimos

		Phobos	Deimos
Dimensions	$2a_b,2b_b,2c_b$	$27\times 22\times 18$ km	$15\times 12.2\times 11$ km
Body axis ratios	$(b_b/a_b,c_b/a_b)$	(0.82,0.67)	(0.82,0.73)
Mean Radius	$R$	11.2667 km	6.2 km
Mass	$M$	$1.0659\times {10}^{16}$ kg	$1.4762\times {10}^{15}$ kg
Mass ratio	$M_{\mathrm{Mars}}/M$	$6.0\times {10}^7$	$4.3\times {10}^8$
Eccentricity	$e$	0.0151	0.00033
Semi-major axis	$a$	$2.76R_{\mathrm{Mars}}$	$6.91R_{\mathrm{Mars}}$
Orbital period	$P=\frac{2\pi }{\omega }$	0.3189 day	1.263 day
Gravitational timescale	$t_g$	1555 s	1418 s
Spin	$\omega t_g$	0.32	0.09
Surface rotational velocity	$\omega R$	2.5 m/s	0.36 m/s
Semi-Major axis/Radius	$a/R$	830.3	3777.6
Mars spin	${\omega }_{\mathrm{Mars}}{t}_g$	0.10	0.11
Radius ratio	$R_{\mathrm{Mars}}/R$	300.8	546.7
Despin time	$t_{\mathrm{despin}}$	3000 yr	$4\times {10}^5$ yr
Eccentricity damping time	$t_{\mathrm{edamp}}$	$0.4\times {10}^9$ yr	$2\times {10}^{12}$ yr
Asphericity parameter	$\alpha$	0.82	0.65
Dissipation parameter	$K_d$	${10}^{-7}$	${10}^{-9}$
		Mars
Mars Radius	$R_{\mathrm{Mars}}$	3389.5 km
Mars Mass	$M_{\mathrm{Mars}}$	$6.4171\times {10}^{23}$ kg
Mars Rotation period	$P_{\mathrm{Mars}}$	1.025957 day
Synchron. semi-maj. axis	$a_{\mathrm{sync}}$	$6R_{\mathrm{Mars}}$
Moon Mass ratio	$M_{\mathrm{Phobos}}/M_{\mathrm{Deimos}}$	7.2

Table 2. Table 2: Base Simulation Parameters

Minimum interparticle distance	$d_I$	0.5
Max spring length	$d_S$	2.3 $d_I$
Number of nodes	$N$	$40$
Number of springs per node	$NS/N$	$8$
Spring constant	$k_s$	$0.22$
Young’s modulus	$E$	$\sim 2$
Time step	$dt$	0.001
Damping time	$t_{\mathrm{damp}}$	100

Table 3. Table 3: Simulations for Resonance Crossings

Simulation	PhM2:1	DeP2:1
Resolved body	Phobos	Deimos
Approximate Body axis ratios	0.8,0.7	0.8,0.7
$M_{\mathrm{Mars}}/M$	$6.0\times {10}^7$	$4.3\times {10}^8$
$R_{\mathrm{Mars}}/R$	301	547
Spin ${\omega }_{\mathrm{Mars}}{t}_g$	0.10	0.11
Quadrupole $c_{22}$	${10}^{-4}$	${10}^{-4}$
Initial $a_{\mathrm{Phobos}}/R$	1200	2390
Initial $a_{\mathrm{Deimos}}/R$		3778
Initial $a_{\mathrm{Phobos}}/R_{\mathrm{Mars}}$	4.0	4.7
Initial $a_{\mathrm{Deimos}}/R_{\mathrm{Mars}}$	-	6.9
$M_{\mathrm{Phobos}}/M_{\mathrm{Deimos}}$		700
Inward drift rate ${\left|{\tau }_{a,\mathrm{Phobos}}\right|}^{-1}{t}_g$	$2\times {10}^{-6}$	${10}^{-7}$
Initial eccentricity Phobos	0.001	0.0
Initial eccentricity Deimos	-	0.0
Spring damping parameter ${\gamma }_s$	3	3
Figures	2	3

Table 4. Table 4: Simulations of Tumbling

Mass ratio $M_{\ast }/M$			${10}^4$
Initial semi major axis $a/R$			48
Initial mean motion $n{t}_g$			0.30
Orbital period $P=2\pi /\left(n{t}_g\right)$			20.9
Number of simulations in each series			10
Resolved body	Deimos		Phobos
Simulation series	DeF	DeS	PhF	PhS
Asphericity $\alpha$	0.7	0.7	0.85	0.85
Initial eccentricity	0.2	0.05	0.2	0.05
Spring damping ${\gamma }_s$	0.1	0.01	0.1	0.01
Figures	4a,6a	4b,6b	5a,7a	5b,7b,8

Equations

$$\frac{dE}{dt}$$

$$T=\frac{3}{2}\frac{G{M}_{\ast }^2}{a}{\left(\frac{R}{a}\right)}^5{k}_2\sin {ϵ}_2.$$

$$t_{\mathrm{despin}}\approx P\frac{1}{15\pi }{\left(\frac{M}{M_{\ast }}\right)}^{\frac{3}{2}}{\left(\frac{a}{R}\right)}^{\frac{9}{2}}\frac{Q}{k_2}$$

$$e_g\equiv \frac{G{M}^2}{R^4}$$

$$t_{\mathrm{wobble}}\approx t_{\mathrm{despin}}{\left(\frac{n}{\omega }\right)}^4$$

$$t_{\mathrm{edamp}}\approx \frac{2}{21}\frac{M}{M_{\mathrm{Mars}}}{\left(\frac{a}{R}\right)}^5\frac{Q}{k_2}{n}^{-1}$$

$$\ddot{x}+\frac{{\alpha }^2}{2}{\left(\frac{a}{r}\right)}^3\sin (2x-2f)=0$$

$$K_d\equiv 3\frac{k_2}{\xi Q}{\left(\frac{R}{a}\right)}^3\frac{M_{\ast }}{M}$$

$$t_g$$

$$U(r,\varphi ,\psi )$$

$$\displaystyle\qquad\left(\tilde{C}_{lm}\cos(m\phi)+\tilde{S}_{lm}\sin(m\phi)\right)P_{l}^{m}(\sin\psi)\bigg{)}.$$

$$\displaystyle\begin{array}[]{lrll}\tilde{C}_{20,M}=&-\texttt{8.75021132354E-04}\pm\texttt{1.252311E-11}\\ \tilde{C}_{22,M}=&-\texttt{8.46359038694E-05}\pm\texttt{2.411019E-12}\\ \tilde{S}_{22,M}=&\texttt{4.89346258602E-05}\pm\texttt{2.415168E-12}\end{array}.$$

$$\frac{de}{dt}\approx n\frac{3}{2}{\left(\frac{R}{a}\right)}^2{c}_{22}\sin \varphi$$

$$\varphi \approx (n-2{\omega }_{\mathrm{Mars}})t+\dot{n}\frac{t^2}{2}+\mathrm{constant},$$

$$\Delta e\approx \sqrt{\frac{4\pi }{3}{\tau }_{a}n}\frac{3}{2}{\left(\frac{R}{a}\right)}^2{c}_{22}.$$

$$\ddot{x}+\frac{{\alpha }^2}{2}{V}_x(x,t)=0$$

$$V_x(x,t)\equiv \sum _{m\ne 0,m=-\infty }^{\infty }W\left(\frac{m}{2},e\right)\sin (2x-mt).$$

$$\alpha >\frac{1}{2+\sqrt{14e}}$$

$$\alpha >\frac{1}{1+e\sqrt{17/2}}.$$

$${\dot{E}}_{\mathrm{ee}}(e)$$

$$C_e$$

$${\dot{E}}_{\mathrm{despin}}$$

$${\dot{E}}_m(e)\approx -C_{e}n\left|\frac{m}{2}-1\right|.$$

$${\dot{E}}_{\mathrm{wob}}(e)=A_w{C}_{e}ne$$

$$H(p,x,t)=\frac{p^2}{2}+\frac{{\alpha }^2}{2}{V}_x(x,t)$$

$$\ddot{x}+\frac{{\alpha }^2}{2}{V}_x(x,t)=-K_d(L(e)\dot{x}-N(e))$$

$$K_e\equiv K_d\xi {\left(\frac{R}{a}\right)}^2$$

$$L\approx I\omega +\beta R{m}_p{v}_{\mathrm{impact}},$$

$$\frac{m_p}{M}\sim \frac{1}{10}\frac{R\omega }{\beta v_{\mathrm{impact}}}.$$

Full text

1 Introduction

Phobos orbits within Mars’ geosynchronous orbit, and it is currently drifting inward due to tides excited by Phobos that are dissipated in Mars (Burns, 1992; Bills et al., 2005). Tides excited by Mars that are dissipated in Phobos (satellite tides) increase Phobos’ inward drift rate in semi-major axis and reduce Phobos’ orbital eccentricity. If Phobos’ eccentricity is only affected by tidal evolution, then Phobos must have had a significantly higher orbital eccentricity in the past (Singer, 1966; Lambeck, 1979; Cazenave et al., 1980; Yoder, 1982). Because the semi-major axis drift rate is increased by satellite tides, Phobos’ orbit would have crossed Deimos’ orbit only 1 to 2 Gyrs ago (Cazenave et al., 1980; Szeto, 1983). This would have caused orbital instability and collisions between the two moons (Szeto, 1983). An epoch of high eccentricity and swift migration is not ruled out if Phobos and Deimos are remnants of earlier moons that collided (Efroimsky, private communication). The tidal drift rates assumed by previous work (e.g., Cazenave et al. 1980) assumed a tidally locked state for Phobos throughout tidal evolution. However, the inferred previous high orbital eccentricities, $e>0.2$, would have caused Phobos to tumble prior to entering the spin synchronous state (Wisdom, 1987). Tumbling would have elevated the energy dissipation rate in Phobos and the associated eccentricity damping rate by orders of magnitude (Wisdom, 1987), exacerbating the problem caused by Phobos’ short tidal evolution timescale.

Alternatively, Phobos could have remained at low eccentricity and its current eccentricity value is due to excitation by other dynamical processes such as crossing spin-orbit resonances (Yoder, 1982). Phobos’ current semi-major axis is $a=2.76R_{\rm Mars}$, where $R_{\rm Mars}$ is the mean equatorial radius of Mars. When Phobos was at a semi-major axis of $a\approx 3.8R_{\rm Mars}$, its orbital period would have been half that of Mars’ spin rotation period. This is a 2:1 spin-orbit resonance between Phobos’ orbit and Mars’ rotation, with strength dependent on the quadrupole gravitational moment of Mars (Yoder, 1982), (also see appendix D by Frouard and Efroimsky 2017). If initially low, Phobos’ orbital eccentricity after passing through this spin-orbit resonance is estimated to be 0.032 and above its current value of 0.015 (Yoder, 1982). Because the satellite spin angular momentum is negligible compared to its orbital angular momentum, the rate of tidal energy dissipation $dE/dt$ is related to the rate of eccentricity change $\dot{e}$ with the constraint that angular momentum is conserved. At low eccentricity

[TABLE]

where $M$ is moon mass, and $n$ is its mean motion. Phobos’ current eccentricity value is consistent with an eccentricity excited previously while crossing this resonance and tidal eccentricity damping since then (Yoder, 1982). If Phobos remained at low eccentricity during the past few Gyr then its semi-major axis drift rate would have remained slow enough to ensure its orbital stability and Deimos’s orbit need never have been encountered.

Deimos’ anomalously low eccentricity, $e=0.00033$, presents a similar problem. Inward tidal migration of Phobos would have increased Deimos’ orbital eccentricity as they crossed the 2:1 orbital mean motion resonance where Deimos’ orbital period is twice that of Phobos (Yoder, 1982). As the time for its excited orbital eccentricity to damp (in the tidally locked state) is greater than the age of the solar system, Deimos would have retained a higher eccentricity than it currently has. This implies that there must be an additional mechanism damping or reducing Deimos’ eccentricity. Wisdom (1987) speculated “Perhaps the resonance passage occurred much earlier when Deimos was more likely to still be chaotically tumbling. ….. The small orbital eccentricity of Deimos could also be a consequence of an episode of chaotic tumbling.” Thus the low eccentricities of both Phobos and Deimos may be a result of prior episodes of tumbling.

The low orbital inclinations of Phobos and Deimos with respect to Mars’ equatorial plane suggest that they were born in or from a circumplanetary disk (e.g., Goldreich 1965; Burns 1992; Rosenblatt 2011; Craddock 2011; Rosenblatt et al. 2016; Hesselbrock and Minton 2017; Canup and Salmon 2018), whereas an inferred previous higher orbital eccentricity for Phobos has been used to argue for capture of Phobos from heliocentric orbit (Burns, 1992; Lambeck, 1979). Phobos and Deimos should remain within $1^{\circ}$ of their Laplace plane during tidal evolution (Cazenave et al., 1980), and this would be independent of Mars’ obliquity variations as they take place adiabatically compared to Phobos’ and Deimos’ nodal precessional motion. A low eccentricity history for Phobos and Deimos is concordant with their formation in and from a circumplanetary disk.

Both Phobos and Deimos are currently in a spin synchronous or tidally locked rotation states where their angular rotation rates match their mean motions. However due to their elongation, Phobos and Deimos can chaotically tumble for long periods of time even at their quite low orbital eccentricities (Wisdom, 1987). Wisdom (1987) speculated that prior to ending in a tidally locked or state, “Phobos must have spent many hundreds of thousands, probably millions, of years in the chaotic zone”. Wisdom (1987) estimated that the timescale for the damping of orbital eccentricity can be a few orders of magnitude shorter for a chaotically tumbling satellite than for a synchronously rotating satellite. However, due to resonant substructure in the chaotic region and its dependence on attitude, it’s difficult to estimate when or how an initially tumbling and tidally dissipating body would exit a tumbling state. Eventually, the rotation could remain close to one of the attitude-stable islands bordering the chaotic zone long enough that tidal dissipation would let the satellite exit the chaotic tumbling dynamical state (Wisdom, 1987).

Impacts significantly affect the spins distributions of asteroids (Harris, 1979) and their rotation states (Henych and Pravec, 2013), excepting the very largest ones, or those with diameters above 100 km, (Farinella et al., 1992). Mars’ moons lie well beneath this size limit, so their spins are likely to have been significantly perturbed by rare energetic impacts from asteroids crossing Mars’ orbit or secondary impacts from crater ejecta leaving Mars. A high fraction, 50% of rocky asteroids and satellites exhibit craters with diameters of 1 times their mean radius ($R_{m}$), and the great majority have more than two craters with diameters greater than $R_{m}/2$ (Thomas, 1999). The impact that formed the Stickney crater on Phobos was probably energetic enough to knock Phobos out of its tidally locked state (Weidenschilling, 1979; Ramsley and Head, 2019). Phobos would have tumbled prior to reentering the tidally locked state in which it currently resides (Wisdom, 1987).

In this paper we reconsider processes for spin excitation and evolution of non-spherical moons Phobos and Deimos. With mass-spring model viscoelastic dynamical simulations, we simulate spin-down, spin-orbit resonance and tumbling in orbit about host planet Mars. Prior studies (e.g., Wisdom 1987) have integrated the equations of motion for dissipationless rigid bodies. Our numerical study work differs from prior work as we use a viscoelastic model for the spinning body. In our simulations, tidal deformation, torque and dissipation are generated within the interior of the spinning body, directly coupling orbital and spin evolution. Our simulations are described in section 2. In section 3 we test the possibility that orbital resonance crossing could have excited the spins of Mars’ moons. In section 4 we explore spin evolution from an initial tumbling state. We look at eccentricities reached from the tumbling state during the transition to the tidally locked state and tidal dissipation rates while tumbling. In section 5 we discuss collisional excitation as a mechanism for tumbling excitation. A summary and discussion follow.

1.1 Physical Quantities

In Table 1 we list orbital and physical values for Phobos and Deimos. Ratios and quantities relevant for our numerical simulations are also listed.

Phobos and Deimos are both in spin synchronous or tidally locked spin states. We estimate the time it could have taken spherical bodies in similar orbits to reach this state. The spin-down time $t_{\rm despin}$ of an orbiting spherical moon can be estimated from its moment of inertia, $I$ and initial spin $\omega$, giving $t_{\rm despin}\sim I\omega/T$ (Peale, 1977) (see Eq. (9) by Gladman et al. 1996). The tidal torque is $T$. The secular part of the semi-diurnal ($l=2$) term in the Fourier expansion of the perturbing potential from point mass $M_{*}$, (here Mars) gives a tidally induced torque on a spherical body of mass $M$, radius $R$ and semi-major axis $a$ (e.g., Kaula 1964; Goldreich 1963; Efroimsky and Makarov 2013).

[TABLE]

Here $k_{2}\sin\epsilon_{2}$ is frequency dependent and is known as a quality function. It is often approximated as $k_{2}/Q$ with energy dissipation factor $Q$ and Love number $k_{2}$ describing the deformation and dissipation in mass $M$. See Efroimsky (2015) for discussion on the accuracy of this approximation. This torque gives a time for a spherical body to spin down

[TABLE]

where $P$ is the spin period. For an incompressible homogeneous spherical elastic body, the Love number $k_{2}\sim 0.038e_{g}/\mu$ where $\mu$ is the moon’s rigidity and

[TABLE]

is a measure of the moon’s gravitational energy density (following Equation 6 by Quillen et al. 2017). Wobble damping timescales estimated from asteroid light curves give an estimate for asteroid viscoelastic material properties $\mu Q\sim 10^{11}$ Pa (Pravec et al., 2014). For energy dissipation factor $Q=100$ this gives a material with rigidity $\mu=1$ GPa, roughly consistent with ice or rocky rubble. Spin down times with $\mu Q=10^{11}$ Pa and $M_{*}=M_{\rm Mars}$ are listed in Table 1 and are computed using semi-major axes, masses and mean radii for Phobos and Deimos that also listed in this Table. The despin times are shorter than those used by Peale (1977) as he assumed a larger value for $\mu Q$.

A time for wobble or non-principal axis rotation to decay (following Peale 1977; Gladman et al. 1996)

[TABLE]

where $n$ is the mean motion and $\omega$ the spin. For a low eccentricity tidally locked body $n=\omega$ and the time for non-principal axis rotation to decay is the same as the spin down time. If tumbling is continuously excited in an elongated body then it may remain tumbling for far longer than the spin down time of Equation 3.

We check estimates for eccentricity damping due to satellite tides in Phobos and Deimos using a timescale

[TABLE]

based on the formula for eccentricity damping for a tidally locked body in eccentric orbit (following Peale and Cassen 1978; Yoder 1982 and Murray and Dermott 1999 Equation 4.198). The resulting eccentricity damping timescales for Phobos and Deimos are listed in Table 1. We confirm that Phobos would have tidally damped its eccentricity but Deimos would not have, confirming previous estimates (Lambeck, 1979; Yoder, 1982; Szeto, 1983).

The classic spin-orbit dynamics problem is that of a rigid non-spherical body undergoing principal axis rotation that is orbiting a massive central point mass at zero obliquity. The equation of motion

[TABLE]

(Wisdom, 1987; Celletti, 2010) with $a,r,f$ respectively, the semimajor axis, the orbital radius and the true anomaly. Here time $t$ is in units of the mean motion and $x$ is the angle between the body’s long axis and the pericenter line. The asphericity parameter $\alpha\equiv\sqrt{3(B-A)/C}$ and $A<B<C$ are the body’s principal axis moments of inertia. The ratio $(B-A)/C\approx(b_{b}/a_{b})^{2}-(c_{b}/a_{b})^{2}$ in terms of body axis ratios where $a_{b},b_{b},c_{b}$ are body semi-major axis lengths. Using the body axis ratios for Phobos and Deimos we compute asphericity parameters for each moon and these are listed in Table 1. The rotational dynamics will be discussed further in section 4. The dissipative spin-orbit problem depends on the unitless parameter

[TABLE]

(Correia and Laskar, 2004; Celletti, 2010; Celletti and Lhotka, 2014). The parameter $\xi$ depends on the moment of inertia, $C=\xi MR^{2}$. With $\mu Q=10^{11}$ Pa and $\xi\approx 0.4$ we include in Table 1 estimates for $K_{d}$ for Phobos and Deimos.

2 Numerical Simulations

To explore tidal and spin evolution in migrating multiple planet and moon systems we use our mass-spring and N-body code (Quillen et al., 2016a; Frouard et al., 2016). Viscoelastic response arises naturally through dissipation in the springs and this makes it possible to explore scenarios where orbital migration and tidal evolution occur (Quillen et al., 2017). The mass-spring model is a viscoelastic and N-body code of low computational complexity that can accurately measure small deformations and conserves angular momentum (Quillen et al., 2016a; Frouard et al., 2016; Quillen et al., 2016b, 2017, 2019c, 2019a, 2019b). Excited normal modes, excited non-principal axis rotation or tumbling, libration, attitude instability, spin-orbit, spin-binary, spin-secular resonance and spin-precession mean-motion capture, all occur naturally within our simulations. The code accuracy (in its ability to measure small strains), simplicity and speed, compared to more computationally intensive grid-based or finite element methods, are strengths of this type of code. The code has been used to measure tidal spin down (Frouard et al., 2016; Quillen et al., 2016b), directly connecting numerically measured torques to the underlying viscoelastic model and without any the underlying analytic assumptions commonly used to compute tidal drift rates. The mismatch between predicted and numerically measured spin down rates (Frouard et al., 2016) is not due to inaccuracy of the code, but due to neglect of bulk viscosity in the analytic computation. We have used the model to explore resonant obliquity evolution in long term integrations of rapidly spinning and non-round satellites (Quillen et al., 2017) and with shorter integrations we resolve the internal tidal heating distribution in a spin synchronous Moon in eccentric orbit (Quillen et al., 2019a). Comparison between numerically measured and analytically predicted wobble (or non-principal axis) damping rates illustrate that the code is accurate and that we understand how the simulation dissipation rates depend on the simulated rheology (Quillen et al., 2019b).

2.1 Sizes and units

Our simulations work with mass $M$ in units of the the mass of a single resolved spinning body. Distances in units of volumetric radius, $R=R_{\rm vol}$, the radius of a spherical body with the same volume. For our mass-spring simulations we work with time in units of a gravitational timescale

[TABLE]

Gravitational timescales for Phobos and Deimos are listed in Table 1. Their low mean density causes these times to be longer than for denser objects.

2.2 Simulations of Phobos and Deimos

As done previously (Quillen et al., 2017), we consider two or three masses in orbit under the influence of gravity. One mass is the non-spherical spinning body, here Phobos or Deimos, and it is resolved with masses and springs. A second mass is the central host planet, here Mars. In some of our simulations a third point mass is included that represents the other moon.

Due to the proximity of Mars to its moons, the rotation of its non-spherical figure can cause significant perturbations on its moon’s orbits (see Yoder 1982). We model the static gravitational field of Mars by taking into account the quadrupole terms in its gravitational potential,

[TABLE]

The normalized (or dimensionless) coefficients from recent measurements of Mars’ gravity field (Genova et al., 2016) are

[TABLE]

Here $\phi$ is longitude and $\psi$ is latitude on Mars, not the colatitude $\theta$ that is often used for spherical coordinates. The associated Legendre functions $P_{2}^{2}(x)=3(1-x^{2})$ and $P_{2}^{0}(x)=(3x^{2}-1)/2$. The $\tilde{C}_{22,M}$ and $\tilde{S}_{22,M}$ coefficients of Mars sum in quadrature to an amplitude $\tilde{c}_{22}\approx 10^{-4}$.

To take into account the rotation of Mars’ gravitational field we rotate the quadrupole components with an angular rotation rate $\omega_{\rm Mars}$. The azimuthal angle with respect to a fixed direction (a true longitude) $\phi_{in}=\phi+\omega_{\rm Mars}t$. The gradient of the quadrupolar gravitational potential is added as an additional force onto each point mass and mass node in the simulation, and equally and oppositely to the mass representing Mars.

As we only resolve Phobos or Deimos with the mass-spring model, tidal evolution due to tidal dissipation in Mars is neglected. To mimic inward or outward drift in orbital semi-major axis (migration) due to tidal dissipation in Mars, we apply small velocity kicks to the resolved body using the recipe for migration given in Equations 8–11 by Beauge et al. (2006) and as we did previously for the migrating Pluto system (Quillen et al. 2017). The migration rate, $\dot{a}$, depends on an exponential timescale $\tau_{a}$ (the parameter $A$ in Equation 9 by Beauge et al. 2006), giving $\dot{a}\sim a\tau_{a}^{-1}$. We adopt a convention $\tau_{a}<0$ corresponding to inward migration for Phobos as it lies within Mars’ synchronous orbit.

For our random spring model, particle or node positions for the resolved body are drawn from an isotropic uniform distribution but only accepted into the spring network if they are within the surface bounding a triaxial ellipsoid, $x^{2}/a_{b}^{2}+y^{2}/b_{b}^{2}+z^{2}/c_{b}^{2}=1$, and if they are more distant than $d_{I}$ from every other previously generated particle. Particles nodes that are not inside the ellipsoid are deleted. For more discussion on the different particle models see Quillen et al. (2016b). Here $a_{b},b_{b},c_{b}$ are the body’s semi-major axes. We use a subscript to differentiate between body semi-major axis and orbital semi-major axis. Once the particle positions have been generated, every pair of particles within $d_{S}$ of each other are connected with a single spring of spring constant $k_{s}$ and damping parameter $\gamma_{s}$. Springs are initiated at their rest lengths.

Because we need to simulate the spin and orbit dynamics for hundreds of orbital periods, we use only a few mass nodes to resolve Phobos or Deimos, similar to our study of the obliquity evolution of Pluto and Charon’s minor satellites (Quillen et al., 2017). The Young’s modulus and viscoelastic relaxation timescale are estimated from the spring parameters, (described previously by Frouard et al. 2016) and here these quantities are approximate because of the sparse network. Sensitivity to the number of springs and nodes in the network was discussed by Quillen et al. (2016b). A simulation snap shot is shown in Figure 1 to illustrate the spring network.

Because of self-gravity the body is not exactly in equilibrium at the beginning of the simulation and the body initially vibrates. We apply additional damping in the springs at the beginning of the simulation for a time $t_{\rm damp}$ to remove these vibrations. Only after $t_{\rm damp}$ do we process simulation outputs for measurements of dissipation. After this time spring parameters are not varied. The numbers of springs, nodes and masses remain fixed during the simulation. When initializing orbital motion, the mean anomalies, longitudes of the ascending node, and arguments of periapse are chosen using uniform random distributions. Simulation outputs are separated by the time interval $\Delta t=10$ in our N-body units.

Base parameters for most of the simulations discussed here are listed in Table 2.

3 Simulations of Resonance Crossings

In this section we test the hypothesis that crossing of an orbital resonance could excite the spins of Phobos or Deimos. Two orbital resonances were discussed by Yoder (1982); Wisdom (1987), a 2:1 resonance between Phobos mean motion and Mars’ spin and a 2:1 mean motion resonance between Phobos and Deimos. Phobos’ eccentricity would have been increased by the 2:1 resonance with Mars’ rotation where Phobos’ orbital period is half that of Mars’ spin period. Phobos would have crossed this resonance as it drifted inwards and when it was near an orbital semi-major axis of $3.8R_{\rm Mars}$. The other resonance is the 2:1 mean motion resonance between Deimos and Phobos which would have been crossed when Phobos was near an orbital semi-major axis of $4.3R_{\rm Mars}$. This resonance would have significantly increased the eccentricity of Deimos.

We show two simulations, one with a resolved Phobos, labelled PhM2:1 and crossing the 2:1 resonance with Mars’ rotation and the other with a resolved Deimos and labelled DeP2:1 crossing the 2:1 mean motion resonance with Phobos. Simulation parameters for these two simulations are listed in Table 3. Initial orbital inclinations are zero. For both of these simulations Phobos is forced to migrate inward with a migration timescale $\tau_{a}$, and mimicking drift due to tidal dissipation in Mars. The mass of Mars in units of the spinning body is used for the central body. The radius of Mars in units of the volumetric radius of the spinning body, and the spin rate of Mars in units of $t_{g}$ (for Phobos or for Deimos) are used to compute the quadrupole force terms from Mars.

Figure 2 shows the simulation of Phobos crossing the 2:1 resonance with Mars’ figure and with parameters listed in Tables 2 and 3. In Figure 2 we plot quantities as a function of time in units of the initial orbital period of the spinning moon (here Phobos).

In Figure 2 orbital elements for the spinning moon are computed using its center of mass position and velocity and the position, velocity and mass of Mars and neglecting its oblateness. The angular momentum of the spinning body, ${\bf L}_{s}$, is computed at each simulation output by summing the angular momentum of each particle node, using node positions and velocities measured with respect to the center of mass of the spinning body. The moment of inertia matrix of the spinning body, $\bf I$, in the fixed reference frame, is similarly computed from node positions and masses. The moment of inertia matrix is diagonalized and the eigenvector corresponding to the largest eigenvalue is identified as the body’s principal axis.

The instantaneous spin vector $\boldsymbol{\omega}$ is computed by multiplying the spin angular momentum vector by the inverse of the moment of inertia matrix, ${\boldsymbol{\omega}}={\bf I}^{-1}{\bf L}_{s}$. The spin $\omega=|{\boldsymbol{\omega}}|$ is the magnitude of this spin vector and shown in the top panel in Figure 2 in units of $t_{g}^{-1}$. At each simulation output, the obliquity is computed from the angle between the body’s spin vector and its orbit normal. One of the Andoyer-Deprit variables, the non-principal angle, $J$, is computed from the angle between angular momentum vector and the body’s principal axis with the largest moment of inertia. The energy dissipation rate is measured from the strain rates in the springs (see Quillen et al. 2019a, b). The energy dissipation rate in the figure is in N-body units or $e_{g}/t_{g}$ (as defined in Equations 4 and 9). Here, we only compare relative variations in the energy dissipation rates. In studies with a larger number of nodes and springs we have quantitatively measured energy dissipation directly from the simulations (Frouard et al., 2016; Quillen et al., 2019b).

From top to bottom, the panels in Figure 2 show the simulated Phobos’ spin, obliquity in degrees, semi-major axis in units of Mars’ radius, energy dissipation rate (in Phobos), non-principal angle in degrees and orbital eccentricity. In the top panel, the green points show spin and the brown dotted line shows the location of the 2:1 orbital resonance with Mars’ figure. The resonance is crossed where the green points cross the brown dotted line. The 2:1 resonance with Mars’ rotation is crossed at $t\approx 0.75\times 10^{3}$ orbital periods and increases the orbital eccentricity of Phobos from 0.001 to $\approx 0.008$. The energy dissipation rate (in Phobos) is approximately 64 times higher after crossing the resonance ($\log_{10}(64)=1.8$) and this ratio is consistent with a tidal heating rate in the tidally locked state that is approximately proportional to $e^{2}$, as expected (Kaula, 1964). The libration caused by the increased orbital eccentricity is evident as a thickening in the breadth of the green points, showing spin, after the resonance is crossed, seen in the top panel of Figure 2.

The eccentricity jump (0.008) caused by the resonance in the PhM2:1 simulation is smaller than the 0.03 estimated by Yoder (1982). We attribute our lower final eccentricity values to the drift rate, not to errors in the calculation by Yoder (1982).111The exponents should be 2 (not 12) for the coefficients in the top two lines in the left column of page 336 by Yoder (1982) and his coefficient $A_{22}$ depends on the ratio $C_{22}/a$ not $C_{22}/A$ as written. Using modern measurements of Mars, we agree with the sizes of the resonant coefficient $C_{22}$ and predicted eccentricity jump estimated by Yoder (1982). By varying Phobos’ drift rate $\tau_{a}$, we find that the eccentricity jump is smaller if Phobos is migrating more quickly. The eccentricity jump is lower than predicted by Yoder (1982) because our simulations are not drifting sufficiently slowly to be well inside the adiabatic regime (e.g., Quillen 2006). The size of the eccentricity jump crossing this resonance can be estimated using Fresnel integrals, following appendix A by Laskar et al. (2004). Following equation 30 by Yoder (1982)

[TABLE]

where the resonant angle $\phi=\lambda-2\theta_{\rm Mars}-\varpi$, the angle $\lambda$ is Phobos’ mean longitude, $\varpi$ is Phobos’ longitude of pericenter, and $\theta_{\rm Mars}$ is a rotation orientation angle for Mars. We approximate the time variation of the resonant angle as

[TABLE]

neglecting the resonant perturbation during resonant crossing. The drift rate in Phobos’ mean motion $\dot{n}=\frac{3}{2}\frac{n}{\tau_{a}}$ in terms of our drift parameter. Neglecting the sensitivity to initial phase, Equation A.8 by Laskar et al. (2004) gives an eccentricity jump

[TABLE]

Using the drift rate and mean motion from our simulation, we computed $\tau_{a}n=9\times 10^{4}$. This, the quadrupolar moment and the above equation gives $\Delta e\sim 0.006$ which is close to that observed in our simulation (0.008).

After an initial spin down period (at $t<10^{4}$), Phobos remains tidally locked, at low obliquity and rotating about a principal body axis. Passage through the resonance does increase the energy dissipation rate as the orbital eccentricity increases, but not as much as it would have if the body had started tumbling. Because the body spins down into a tidally locked state prior to encountering the orbital resonance, the final spin state is insensitive to initial conditions. However there is some variation between different simulations in the size of the eccentricity jump that takes place when the resonance is crossed, as it depends on the phase of the resonant angle when the resonance is encountered (Yoder, 1979; Quillen, 2006). A similar simulation but with a three times larger $c_{22}$ coefficient gives a higher eccentricity jump and also fails to push the simulated Phobos out of the tidally locked state. We conclude that crossing the 2:1 orbital resonance with Mars’ figure rotation was unlikely to push Phobos out of its tidally locked state or cause it to tumble.

In Figure 3 we show a simulation, denoted DeP2:1, where the spinning body represents Deimos. This simulation is designed to see how the 2:1 mean motion resonance between Phobos and Deimos would affect the spin state of Deimos. As Deimos is near Mars’ synchronous radius (its mean motion is similar to the rotation rate of Mars), Deimos’ semi-major axis drifts slowly due to tides excited in Mars. It is Phobos’ inward migration that lets the two moons cross this mean motion resonance. Due to the low mass ratio $m_{\rm Phobos}/m_{\rm Mars}\sim 10^{-8}$, the resonance is weak. To see the effect of the resonance, Phobos must drift slowly compared to the resonance libration time (e.g., Quillen 2006). To reduce the time required for the simulation, we increased the mass of Phobos by a factor of 100 relative to its actual mass ratio with Mars.

Figure 3 of a spinning Deimos is similar to Figure 2 for Phobos, and its parameters are also listed in Table 3. In this simulation, Phobos is included as a third point mass and is again forced to migration inward. In the third panel from top, the semi-major axis times $2^{-\frac{2}{3}}$ of Phobos is shown in brown and the semi-major axis of Deimos is shown in black. The 2:1 mean motion resonance occurs where these two lines cross. In the bottom panel the eccentricities for Deimos and Phobos are shown in black and dark red, respectively. The ratio of Phobos’ to Deimos’ mass is 100 times larger than their actual mass ratio, to reduce the drift rate needed to see the effect of this resonance. The resonance does increase the eccentricity of Deimos, but slowly. Even with a much higher Phobos mass, crossing the resonance does not excite Deimos’ spin.

Using scaling for first order resonances (Quillen, 2006; Mustill and Wyatt, 2011), we can check to see if Phobos’ semi-major axis drift rate would have been sufficiently slow that the 2:1 mean motion resonance between Phobos and Deimos resonance was crossed adiabatically. The resonance is crossed adiabatically if the drift rate in Phobos’ mean motion $\dot{n}/n^{2}$ is significantly smaller than the mass ratio $(m_{\rm Phobos}/m_{\rm Mars})^{\frac{4}{3}}=4.2\times 10^{-11}$ (Quillen, 2006; Mustill and Wyatt, 2011). Using Phobos’ current semi-major axis drift rate $\dot{n}/n^{2}=5\times 10^{-13}$ (Bills et al., 2005) and scaling to the semi-major axis of the resonance at $a=4.3R_{\rm Mars}$, we estimate that Phobos’ drift rate would have been about 1000 times slower than the critical value defining the adiabatic limit for this resonance. Thus Deimos’ low eccentricity is not because the 2:1 mean motion resonance was crossed too quickly to increase Deimos’ eccentricity.

In both simulations, crossing orbital resonances due to an inward drifting Phobos did not excite the spin of a tidally locked Deimos or Phobos. We rule out the possibility that crossing the strongest orbital resonances also excited Phobos’ or Deimos’ spin.

4 Evolution of tumbling bodies

In the previous section (section 3) orbital migration was forced to mimic drift caused by tides dissipated in Mars. Dissipation within the spinning body primarily served to damp the body into the tidally locked state and keep it there. However, in this section orbital migration and eccentricity damping arise solely due to satellite tides and internal dissipation coming from damping in the springs.

We have run two series of simulations of a spinning resolved body, with axis ratios similar to Deimos and two series with axis ratios similar to Phobos. The spinning body is alone in orbit about a point mass. The DeF and PhF series have initial eccentricity $e=0.2$ whereas the DeS and PhS series start with $e=0.05$. The DeF and PhF series have 10 times higher damping in the springs (the $F$ represents fast) than the DeS and PhS series (the $S$ represents slow). The DeF and DeS simulations of Deimos have the same node and spring network and the PhF and PhS simulations of Phobos have the same network. The parameters for these simulations are listed in Tables 2 and 4.

By initially setting the body’s spin above the mean motion and at moderate obliquity, we ensure that the body begins to tumble soon afterwards. Individual simulations in each series are begun at different spin, obliquities and mean anomaly. In these simulations there is no forced migration $\tau_{a}^{-1}=0$ and the quadrupole moment of the central body is zero. As the energy is tidally dissipated, the eccentricity decreases and eventually the body drops into a spin synchronous state. We measure the energy dissipation rates during the simulations and the eccentricities reached when the simulated elongated moon finally enters into the tidally locked state.

As the tidal dissipation depends strongly on the ratio of semi-major axis to radius, tidal dissipation is stronger if we put the spinning body closer to the central mass. To reduce simulation times, we reduced the central mass $M_{*}$ to lower than that of Mars in units of the spinning body. This is equivalent to inflating the radius of the spinning body and it should not significantly affect the dynamics of the spin-orbit problem which is primarily sensitive to the body axis ratios (Quillen et al., 2017). The central mass $M_{*}$ and spring damping parameter $\gamma_{s}$ were adjusted to allow significant orbital eccentricity evolution during the simulations.

The DeF and PhF simulations were run to a maximum time of $t=10^{6}$ and the DeS and PhS simulations to $2\times 10^{6}$ in our gravitational units. For the DeF and PhF simulations, this maximum length of time is equivalent to about $5\times 10^{4}$ orbital periods. As our mass-spring network is course and randomly generated, the moments of inertia are not the same as the triaxial boundary used to confine the randomly generated node positions. We have included in Table 4 the asphericity parameter computed from the nodes in the simulations. For the DeF and DeS simulations the asphericity parameter is similar to that of Deimos and for the PhF and PhS simulations, similar to that of Phobos.

Eccentricity and spin evolution for these series of simulations are shown as a function of time in orbital periods in Figures 4 and 5. We find that tumbling can last thousands of orbital periods, confirming the work by Wisdom (1987). Simulations with large semi-major axes or lower levels of spring damping (the DeS simulations compared to the DeF simulations) or stiffer springs spend even longer time in tumbling states, as expected. Prior to tidal lock, the simulations primarily exhibited tumbling and only rarely experience short intervals of spin-orbit resonance capture.

Only the DeF simulations show intervals captured into the 2:1 spin-orbit resonance (see Figure 4a). The DeF simulations were begun at higher eccentricity than the DeS simulations and the spin reaches higher values at the beginning of the simulation. The DeF simulations also have a higher spring damping parameter than the DeS simulations. We ran a series of simulations at the lower damping rate but starting at $e=0.2$ and found that capture into the 2:1 resonance can occur. We suspect that the reason the DeS simulations don’t capture into the 2:1 spin-orbit resonance is not because these simulations have a lower dissipation rate, but because the mean spin values while tumbling are too low at lower eccentricity. At lower eccentricity, the tumbling body is usually spinning slower than the 2:1 spin-orbit resonance. For more elongated bodies, as seen in the PhF and PhS simulations, capture into the 2:1 spin-orbit resonance was also unlikely, though two of the PhF simulations had intervals in the 3:1 spin-orbit resonance (see Figure 5a).

The equation of motion for the classic conservative spin-orbit dynamical problem from Equation 7 can be written (following Celletti 2010; Celletti and Lhotka 2014)

[TABLE]

with asphericity parameter $\alpha$ and

[TABLE]

For $m=2,3,4$, the coefficients are $W(1,e)=1$, $W(3/2,e)=\frac{7e}{2}$ and $W(2,e)=\frac{17}{2}e^{2}$ (Cayley, 1859; Goldreich and Peale, 1966) and describe the strengths of the 1:1, 3:2 and 2:1 spin-orbit resonances, respectively. Conditions for resonance overlap are computed using a distance between spin-orbit resonances in units of the mean motion and resonant widths. The resonant widths are $\alpha$ for the 1:1 (spin synchronous) resonance, $\alpha\sqrt{7e/2}$ for the 3:2 spin-orbit resonance and $\alpha e\sqrt{17/2}$ for the 2:1 spin-orbit resonance. The condition for overlap of the 1:1 and 3:2 spin-orbit resonances at low obliquity and eccentricity is

[TABLE]

and predicts when chaotic tumbling takes place (Wisdom et al., 1984).

Equation 20 implies that for asphericity parameter $\alpha>1/2$ even a value of zero eccentricity would give chaotic behavior. Asphericity parameters for Phobos and Deimos are both above 1/2. For $\alpha>1/2$, the width of the spin synchronous or 1:1 spin-orbit resonance exceeds the distance between the 1:1 and 3:2 spin-orbit resonances. This means that the 3:2 spin-orbit resonance is diminished or swallowed up by the chaotic zone associated with the separatrix of the 1:1 resonance. This would explain why we do not see capture into the 3:2 spin-orbit resonance in these 4 series of simulations. We have run similar simulations of rounder bodies, and those do spend time trapped in the 3:2 spin-orbit resonance.

A similarly derived condition for resonance overlap between 1:1 and 2:1 spin-orbit resonances gives

[TABLE]

For the asphericity of Phobos the overlap condition is satisfied at $e>0.08$ and for Deimos at $e>0.19$. Our DeF and DeS simulations have asphericity parameter similar to that of Deimos. Contrary to expectation, the 2:1 resonance island might increase with decreasing eccentricity rather than shrink as it emerges from the tumbling chaotic region. This may explain why we tended to see capture into the 2:1 resonance after the eccentricity had time to drop in the DeF simulations and why similar simulations of more elongated bodies (the PhF and PhS simulations) are unlikely to spend time intervals in the 2:1 spin-orbit resonance.

At the beginning of the simulations and at high eccentricities the PhF simulations (see Figure 5a) have high spin values, sometimes even above 1. The spin in these simulations in gravitational units is 0.3 which is approximately equal to that of Phobos in its gravitational units. A spin value of 1 is spinning so fast that material could leave the surface. If Phobos were ever at large eccentricity $e\gtrsim 0.15$ it would have tumbled and during tumbling it would have experienced episodes where material was spun off its surface.

4.1 Dissipation rates

In Figure 6 we show normalized energy dissipation rates versus orbital eccentricity for the DeF and DeS series of simulations and in Figure 7 for the PhF and PhS series. We box averaged the energy dissipation over intervals $10^{2}$ long and these points are shown as solid points in the Figure. Open points show measurements from the simulation output without average, though we only plotted 1 out of 10 points to reduce confusion. Each color in these figures shows a different simulation and points show simulation measurements at different times. Labels in the key on these figures refer to their order in Figures 4 and 5.

In Figures 6 an 7, as eccentricity decreases due to tidal dissipation, a simulation begins on the upper right side in one of and drifts to the lower left. The dissipation rate drops abruptly when the system falls into a spin synchronous state. After reaching the spin synchronous resonance, obliquity, libration and non-principal axis rotation are damped and the energy dissipation rate can continue to drop. Afterwards the dissipation rate is extremely low and so is the rate of eccentricity change. After the tidally locked state is reached, tidal evolution slows. All the simulations end in the spin synchronous state, though some still exhibit non-principal axis rotation, non-zero obliquity and free libration amplitudes which take longer to decay.

The energy tidally dissipated for a body that is spinning about a principal body axis and in the spin synchronous state with little free libration is

[TABLE]

(Kaula, 1964; Cassen et al., 1980), where $I$ is the obliquity. Damping in eccentricity is accompanied by inward drift in orbital semi-major axis and tidal heating rates are extremely sensitive to semi-major axis. Equations 22, 23 and $n\propto a^{-3/2}$ implies that $\dot{E}_{\rm ee}\propto a^{-\frac{15}{2}}$. In Figures 6 and 7 we have corrected the energy dissipation rates by this factor of semi-major axis to take into account the slight decrease in semi-major axis occurring during the simulations. This correction slightly lifts the low eccentricity points on Figures 6 and 7.

The tidal energy dissipation rate for an object undergoing principal axis rotation and at zero obliquity that is spinning down (and not in a spin synchronous state)

[TABLE]

(e.g., Murray and Dermott 1999 Equation 4.151 and 4.160) and using the coefficient $C_{e}$ defined in Equation 23. This rate is expected to characterize the energy dissipation rate for tumbling states (Wisdom, 1987). By setting $\omega$ equal to the value expected in a spin-orbit resonance, Equation 24 gives an estimate for energy dissipation rates in a spin-orbit resonance. For example in the 3:2 spin-orbit resonance set $\omega-n=3/2n-n=n/2$. The mean tidal dissipation rate in spin-orbit resonance at low obliquity is approximately characterized by spin $\omega/n=m/2$ for $m\neq 2$, and dissipation rate

[TABLE]

The lower left hand side of Figure 6a shows that simulations of a low obliquity, tidally locked body undergoing principal axis rotation show $dE/dt\propto e^{2}$, as given in Equation 22. Because it is well studied, we normalize the dissipation rates using the tidally locked states. Checking to make sure the body is rotating about a principal axis (at low angle $J$), we measure the energy dissipation rate, obliquity and eccentricity in a simulation after it has reached a tidally locked state. We use that state and Equation 22 to estimate $C_{e}n$. The energy dissipation rates are then divided by $C_{e}n$. In the 2:1 spin-orbit resonance $|\omega-n|/n=1$ so our normalization gives $|\dot{E}|/(C_{e}n)\approx 1$ in this resonance as long as the obliquity is low. The grey dotted lines in Figures 6 and 7 show $|\dot{E}|/(C_{e}n)=7e^{2}$, consistent with Equation 22 for a tidally locked state at $I=0$. The minimum energy dissipation rates in the tidally locked states all lie near and parallel to this line, as expected. Light green dot-dashed lines in Figure 6 and 7 show estimates for the 3:2, 2:1 and 5:2 spin-orbit resonances. Four simulations in the DeF series have intervals in the 2:1 spin-orbit resonance and while they are in that resonance, their energy dissipation rates lie near the green dot-dashed line for that resonance in Figure 6a.

Figures 6 and 7 show that the energy dissipation rate during tumbling is 1 to 6 orders of magnitude above that in the tidally locked states, quantitatively confirming the rough estimate by Wisdom (1987). There is no analytical formula for the energy dissipation rate in a tumbling state, though an approximate value was estimated by Wisdom (1987) using Equation 24 and replacing $|\omega-n|$ with $n$. Figure 6a shows that this is a good approximation as long as the eccentricity is above about 0.05. Figures 6 and 7 show that the energy dissipation rate is weakly dependent on eccentricity with somewhat lower rates at lower eccentricities. To roughly characterize the mean energy dissipation rate seen in the wobbling states we plotted as a brown lines in Figures 6 and 7

[TABLE]

using coefficient $A_{w}=14$ in Figure 6 and $A_{w}=7$ in Figure 7.

We lack a proposed theoretical explanation for the dependence of dissipation rate while tumbling on eccentricity, but the brown lines are a decent match to the energy dissipation rates exhibited by the tumbling bodies. In the classic spin-orbit problem (Equation 7), the chaotic tumbling zone seen in Poincaré maps of the planar system becomes increasingly narrow with decreasing eccentricity (Wisdom, 1987). At low eccentricity, we would not expect the average dissipation rate to be sensitive to eccentricity. Equation 7 for the spin-orbit problem only allows a single angle for body rotation and tumbling bodies rotate in three-dimensions. For the three-dimensional problem chaotic regions are present even at zero eccentricity (Wisdom, 1987). The dependence on eccentricity is related to the size of the configuration space that the tumbling body can explore in the three-dimensional problem.

With principal axis rotation, our previous work has estimated how the tidal spin down rate depends on body shape (at most a factor of 2 difference from the volume equivalent sphere; Quillen et al. 2016b). However the energy dissipation rate in wobbling or tidally locked states probably does not depend on body shape in the same way. The difference between the dissipation rate in the 2:1 spin orbit resonance and the green dashed line in Figure 6a may be due to this discrepancy or other degrees of freedom that have not damped away, such as free libration, obliquity and non-principal axis rotation. Similarly the differences in the estimated coefficients $A_{w}$ for the brown lines might be sensitive to body shape. The lines giving dissipation rates in spin-orbit resonances, tidally locked and tumbling states on Figures 6 and 7 could in future be corrected for a dependence on body shape.

Our simulations, shown in Figures 6 and 7, show that a tumbling body can remain tumbling long enough that it reaches a very low orbital eccentricity, confirming the prediction by Wisdom (1987). If Deimos or Phobos were excited into a tumbling state, they could end up at quite low eccentricity, so excitation of tumbling could account for their current low eccentricities. The eccentricity damping time for Deimos is currently of order $\sim 10^{12}$ years (as we listed in Table 1; see section 1.1). From comparing the height of the brown line and the grey dotted line at an eccentricity of 0.003, we see that the dissipation rate in the tumbling state is about 5 orders of magnitude higher than that in the tidally locked state. An increase in the dissipation rate by only 3 orders of magnitude could lower the eccentricity damping timescale for Deimos enough to damp its eccentricity within the age of the Solar system (assuming $\mu Q=10^{11}$ Pa). With dissipation 5 orders of magnitude higher than in the tidally locked state, Deimos would have had to be tumbling only for tens of millions of years before falling into its current tidally locked rotation state in order for its eccentricity to have dropped significantly due to tidal dissipation.

We examine the tidal heat power that could be present in a tumbling state. Using Equation 25, with $m=2$, representative of the tidal power in the tumbling state, and using values for $k_{2},Q$ discussed previously in section 1.1 and listed in Table 1, we estimate the tidal heating power for Phobos and Deimos. We find $\dot{E}\approx 4.6\times 10^{5}$ W for Phobos and $\dot{E}\approx 5$ W for Deimos. Using a thermal conductivity of $K_{T}=1$ W/(m K) which is within a factor of a few similar to that of ice, rock or porous rock, we estimate a core to surface temperature difference of $\Delta T=\dot{E}/(K_{T}4\pi R)=3.2$ K for Phobos. These low powers and small core to surface temperature difference imply that tidal heating is not likely to affect the composition or interior structure of Phobos or Deimos, even if they spent extended periods of time tumbling.

We consider the sensitivity of our estimates to rheology. The tidal frequency $2(n-\omega)$ can be normalized by the viscoelastic relaxation time $\tau_{\rm relax}\equiv\eta/\mu$ with $\eta$ the viscosity. This gives a normalized tidal frequency ${\bar{\chi}}=2(n-\omega)\tau_{\rm relax}$. We run our simulations in the regime with frequency $\bar{\chi}\ll 1$ where the quality function is proportional to $\bar{\chi}$ for both Maxwell rheology and that exhibited by our numerical simulations (Efroimsky and Makarov, 2013; Frouard et al., 2016). However, the short rotation periods of Phobos and Deimos may put them in the opposite limit ($\bar{\chi}\gg 1$) where the quality function is instead inversely proportional to ${\bar{\chi}}$ (again for both simulation and Maxwell rheologies). In the tumbling state, the spin contains multiple frequency components, with average near the mean motion. In the spin-synchronous state, the tidal perturbations are at multiples of the mean motion with strongest term set by the mean motion. We expect that the ratio of dissipation rates, between spin-synchronous and tumbling, would be similar for the high frequency limit as for the low frequency limit. As a consequence, the relative dissipation rates for the spin-synchronous and tumbling states for the low frequency limit, which is shown in our figures, can be used to estimate the relative dissipation rates for bodies that are in the high frequency limit. However, the higher spin rate spin-orbit resonances would have larger normalized tidal frequency $\bar{\chi}$ than the 2:1 spin-orbit resonance and would have lower dissipation rates for a Maxwell model, rather than higher ones as we estimated in Equation 25. For the spin-orbit resonances, a Maxwell model would show different relative dissipation rates.

4.2 Long-lived non-principal axis rotation

The grey dotted lines in Figures 6 and 7 representing the dissipation rates in the tidally locked states are corrected for non-zero obliquity but not for long lasting non-principal axis rotation or free libration. The reason that points lie above the grey dotted line at $e\sim 0.01$ in Figure 7a showing the PhF series is that the simulations have obliquity of a few degrees. Frouard and Efroimsky (2017) have estimated the sensitivity of the dissipation rate to free libration. Less straightforward is how the dissipation rate depends on non-principal axis rotation.

In some simulations we saw long lived non-principal axis rotation in the spin-synchronous state. This is evident from looking at their oscillations in spin, eccentricity and from the non-zero value of the non-principal angle $J$. An example of long lived non-principal axis rotation is shown in Figure 8 which shows the 10-th simulation from the PhS series. The panels in Figure 8 are similar to those in Figure 2. In Figure 8, the second panel from the bottom shows the non-principal angle $J$. After $3\times 10^{4}$ orbital periods, the body has entered the spin synchronous state, but $J$ continues to be above zero. After $7\times 10^{4}$ orbital periods, the mean $J$ stops decreasing and the body remains in a long lived non-principal axis rotation state. Many of the PhS simulations show long lived non-principal axis rotation but with slowly decaying amplitudes.

The PhS simulation shown in Figure 8 is striking because at later times there is no evidence for decay in the non-principal axis rotation and this suggests that the body has been trapped in a resonance with body axis precession frequencies. Because the non-principal angle $J$ does not decay, the energy dissipation rate remains higher than it would be had the body damped into a state where it rotated about a principal body axis. These resonances might have previously been seen as stable non-chaotic islands in a surface of section for the out-of-plane motion of a prolate, axisymmetric body in a circular orbit that exist even at zero orbital eccentricity (see Figure 7 by Wisdom 1987). Such resonant islands probably exist for bodies that are not axisymmetric. Our PhS simulations suggest that tidally dissipating bodies can be captured into this type of resonance. This might be more likely for longer and more slowly tidally evolving bodies as we saw this behavior in the slower simulations of larger asphericity bodies.

The energy dissipation rate is higher in the long lived non-principal axis rotation states than in a tidally locked state where the body is undergoing rotation about a principal body axis. In the 4-th and 10-th PhS simulations the energy dissipation rate flattens out at $\sim 5\times 10^{-5}$ due to non-principal axis rotation. The energy dissipation rate is about an order of magnitude above that expected for the tidally locked states. The discrepancy is noticeable because of the low eccentricity. The other simulations in the PhF series also exhibit non-principal axis rotation but because the eccentricity is higher, the enhanced dissipation is less noticeable.

If the body stays captured in a resonance involving non-principal axis rotation for a long time then the classic formula for tidal evolution in the spin synchronous state would underestimate the energy dissipation rate and the associated eccentricity damping rate. The size of the underestimate would depend on the amplitude of non-principal axis rotation, the eccentricity and obliquity.

4.3 Chaotic Dynamics

Wisdom (1987) had speculated that once tumbling, eventually the body would exit the chaotic zone if it entered a sticky region near a stable resonant island. Our numerical study suggests that at low eccentricity, and for bodies with large asphericity such as Phobos and Deimos, stability is only achieved in the spin synchronous state. As the chaotic zone volume shrinks due to a slow decrease in eccentricity, a dissipating body can nevertheless remain for large periods of time in the tumbling state, even to eccentricities below 0.001. In Figures 6 and 7, the distributions of eccentricities reached when the simulations drop into the spin-synchronous states is fairly uniform and quite broad. We discuss possible qualitative explanations for this behavior.

The classic spin-orbit dynamics problem (shown in Equations 7 and 18) is that of a rigid non-spherical body undergoing principal axis rotation that is orbiting a central point mass at zero obliquity. The dynamics is conservative and described by a time dependent and periodic Hamiltonian that is a function of rotation angle and its time derivative (Celletti et al., 2007; Celletti, 2010),

[TABLE]

with canonical momentum $p=\dot{x}$. The Hamiltonian is a function of orbital eccentricity and the eccentricity is constant.

The dissipative spin-orbit problem includes in addition a small tidal torque that is due to tidally generated heat, giving equation of motion

[TABLE]

with unit-less dissipation factor $K_{d}$ defined in Equation 8. The functions $L(e),N(e)$ depend on eccentricity (Correia and Laskar, 2004; Celletti, 2010; Celletti and Lhotka, 2014). The dissipative spin-orbit problem also usually neglects variations in eccentricity.

For the DeS and PhS series of simulations we estimate dissipation parameter $K_{d}\sim 10^{-5}$ using estimates for the quality function, valid for the limit of low tidal frequency times viscoelastic relaxation time ($\bar{\chi}<1$), that is based on the tidal frequency, spring strengths and damping parameters ($k_{2}/Q\sim 0.038{\bar{\chi}}/\mu$; for details see Frouard et al. 2016). The DeF and PhF simulations have 10 times higher dissipation rate than the DeS and PhS simulations. For the actual Phobos and Deimos, $K_{d}\sim 10^{-7}$ and $10^{-9}$, respectively, (see Table 1). The numerical values are a few orders of magnitude larger than those of the actual moons. However, our numerical simulations have unit-less dissipation factor $K_{d}\ll 1$, as do Phobos and Deimos.

Because the tidal dissipation rate is low, the dynamics of the related conservative problem can help interpret the dynamics of the non-conservative system. The dissipative spin-orbit problem exhibits different kinds of attractors, including periodic, quasi-periodic and strange attractors (Celletti et al., 2007; Celletti and MacKay, 2007; Melnikov, 2014). Regions in phase space can be divided into neighborhoods of attraction. Some regions may not contain rotational invariant tori (this is referred to as non-existence of rotational invariant tori by Celletti and MacKay 2007). These regions can contain strange attractors so some orbits remain in them forever.

The problem we have studied here is more complex than the dissipative spin-orbit dynamical model as our body spins in three-dimensions and the tidal dissipation also causes slow variations in eccentricity. Nevertheless as the eccentricity damping rate is slow, the dynamics of the related low dimension dissipative and non-dissipative (classic) spin-orbit problems might help interpret our simulations. A difference between the classic spin-orbit problem and the full 3d but conservative spin dynamics problem is that the tumbling chaotic zone can persist even at zero eccentricity Wisdom (1987) when rotation out of the plane is allowed, but vanishes at zero eccentricity in the classic low dimensional spin-orbit problem.

Our simulations have faster eccentricity damping rates in units of the mean motion than Phobos and Deimos because we have inflated the radii of the spinning bodies. The eccentricity damping rate depends on dimensionless parameter

[TABLE]

(following Equations 22, 23, 25, 26) with $\frac{de^{2}}{dt}n^{-1}\sim K_{e}$ for the tumbling states. The dimensionless parameter $K_{e}$ characterizes the rate of eccentricity decrease. The two factors of $R/a$ imply that $K_{e}\ll K_{d}$ for our simulations of Phobos and Deimos (and as desired) but the PhS simulations have $K_{e}$ about $10^{3}$ times that of Phobos and the DeS simulations have $K_{e}$ about $10^{6}$ times that of Deimos.

The existence of attractors in the dissipative spin-orbit system implies that without eccentricity damping the dissipative spin-orbit system probably can remain in a tumbling state forever. The dissipation itself does not necessarily cause the system to exit the chaotic zone, even though it causes volume in phase space to shrink. In contrast, we can consider the evolution of the classic spin-orbit problem, without any dissipation but with an adiabatically decreasing eccentricity. This is a time dependent but conservative problem. As the eccentricity decreases, the width of the resonances decrease. A system initially in a chaotic zone could drop out of it, not because of dissipation, but because the volume in phase space that is covered by the chaotic region shrinks as the eccentricity drops.

A tumbling body continues to dissipate energy, giving an inward drift in orbital semi-major axis. At extremely low eccentricity, angular momentum conservation can no longer be maintained by an accompanying decrease in orbital eccentricity. Wisdom (1987) postulated that “if the chaotic tumbling did persist to the point of completely damping the eccentricity, the chaotic tumbling would then be damped at a somewhat faster rate through the secular damping of the rotational Jacobi integral”. Here the Jacobi integral is the rotational kinetic energy minus the rotation angular momentum times the mean motion.

Our simulations do not exhibit a transition to a faster mode of damping at low eccentricity. Rather the eccentricities reached when dropping into the spin synchronous state seem fairly evenly distributed. The eccentricity distribution seems to be related to a reduction in the accessible tumbling chaotic zone volume as the eccentricity drops. We use the chaotic zone width estimated from the classic low dimensional spin-orbit problem to estimate a dependence of chaotic zone volume on eccentricity, even though the low dimensional spin-orbit problem might poorly approximate the full 3D problem.

At low eccentricity, the 3:2 spin orbit resonance can be considered a weak perturbation primarily causing chaotic behavior on the separatrix of the 1:1 spin orbit resonance (Wisdom, 1987). The 3:2 resonance perturbation term in the classic spin-orbit Hamiltonian would be stronger than the other spin-orbit perturbations. The width of the chaotic zone in the classic spin-orbit problem depends on the strength of the 3:2 resonance or $W(3/2,e)\propto e$. The frequency of the perturbation angle (of the 3:2 resonant angle) compared to that of libration in the 1:1 resonance for order unity asphericity is $\lambda=\alpha/2$ which is of order unity. The perturbation is neither in the adiabatic or non-adiabatic limits, however Shevchenko (2012) estimates that the width of the chaotic zone near the separatrix has energy proportional to the perturbation strength in both limits. This implies that the separatrix width in energy is proportional to eccentricity. The width in units of $\dot{x}$ would be proportional to $\sqrt{e}$ giving a chaotic zone volume in 2-dimensional phase space ($\dot{x},x$) that is also proportional to $\sqrt{e}$. The weak dependence of chaotic zone width on eccentricity in the spin-orbit problem was noticed by Wisdom (1987) who commented in his Figure 3 caption regarding the chaotic zone for Deimos ‘Even with such a small eccentricity, the chaotic zone is not microscopic.’

We consider the volume in 2-dimensional phase space lost per decade of eccentricity decrease. We start with points evenly distributed in the chaotic zone at the separatrix with volume $\propto\sqrt{e}$. With eccentricity dropping adiabatically, the number of points leaving the chaotic zone in a decade in eccentricity would only be $\sqrt{10}\sim 3$ times more than in a consecutive lower eccentricity decade. For example the number of points leaving the chaotic zone at the separatrix between $e=0.01$ and $0.001$ should be about three times that leaving between $e=0.001$ and 0.0001. Such a weak dependence of chaotic zone volume on eccentricity would be consistent with the broad distributions of eccentricities reached when the simulations drop into the spin-synchronous states seen in Figures 6 and 7. We attribute the tendency for simulations of an initially tumbling body to drop into spin-synchronous state at very low eccentricity to the insensitivity of the chaotic zone volume to eccentricity.

5 Spin excitation by impacts

Large but rare collisions can significantly excite the spins of asteroids with diameters less than 100 km (Harris, 1979; Farinella et al., 1998; Henych and Pravec, 2013). As a consequence rare but nearly catastrophic impacts could also have pushed Phobos and Deimos out of their spin synchronous states. Impacts below the catastrophic disruption threshold cannot significantly change orbital eccentricity. This follows because satellite spin angular momentum is negligible compared to satellite orbital angular momentum. However a moon that is put in a tumbling state would experience enhanced dissipation until it returns into the tidally locked state. During this time, the moon’s orbital eccentricity would be damped.

To push a moon out of tidal lock, the change in its spin $\delta\omega$ need only be $\omega/4$ as this is half of the distance to the 3:2 and 1:2 spin-orbit resonances. With a tangential impulse from a grazing impact, the spin angular momentum after impact is

[TABLE]

where $\omega$ is the moon’s initial spin, $I\sim\frac{2}{5}MR^{2}$ is the moon’s moment of inertia, $M,R$ are the moon’s mass and radius, $m_{p}$ is the projectile mass and $v_{\rm impact}$ is the relative velocity of projectile and target moon. The dimensionless parameter $\beta$ describes the total momentum transfer. Setting a change in angular momentum to $I\omega/4$ gives a minimum mass ratio

[TABLE]

The spin rotational velocity at the surface of the moon $R\omega$ is 2.5 m/s for Phobos and 0.4 m/s for Deimos. We adopt a velocity of $v_{\rm impact}=10$ km/s typical of collisions onto Mars (Neukum and Wise, 1976; Ivanov, 2001). This gives a minimum projectile to target moon mass ratio of order $2.5\times 10^{-5}$ corresponding to a projectile diameter of about 0.7 km onto Phobos and approximately consistent with the models by Ramsley and Head (2019) for the projectile that formed the Stickney crater on Phobos. To kick Deimos out of the tidally locked state, a projectile to target moon mass ratio mass ratio of $4\times 10^{-6}$ is required, corresponding to a 0.2 km diameter projectile. Deimos is spinning less quickly than Phobos, so a lower mass ratio projectile can knock it out of the spin synchronous spin state. These mass ratios lie in the regime of rare and nearly catastrophic collisions that can excite tumbling (see Figure 12 by Henych and Pravec 2013).

How rare are such impacts? The cross sectional areas of Phobos and Deimos are approximately 400 km2 and 120 km2. The mean time between impacts as a function of impactor diameter has been computed for Mars crossing asteroid 433 Eros (with area 1125 km2) by Richardson et al. (2005) based on the asteroid size distribution model by O’Brien and Greenberg (2005). We estimate the mean time between impactors by multiplying by the ratio of cross sections. This gives a few times $10^{9}$ years between 0.7 km diameter impactors on Phobos and $10^{9}$ years between 0.2 km diameter impacts on Deimos. These times neglect secondary impacts from Mars and that the bombardment rate was higher in the past. We conclude that as true for asteroids, rare and sub-catastrophic impacts capable of excitation of wobble would also have occurred on Phobos and Deimos.

6 Discussion and Summary

Using a mass-spring model for modeling spin and tidally generated dissipation within non-spherical elastic bodies, we have explored scenarios for spin excitation of Phobos and Deimos. We test the possibility that crossing the strongest orbital resonances (due to inwards drift of Phobos) could have excited the spins of Phobos and Deimos. Our simulations of an initially tidally locked Phobos crossing the 2:1 resonance with Mars’ figure rotation and an initially tidally locked Deimos crossing the 2:1 mean motion resonance with Phobos failed to show excitation of libration or wobbling in Phobos or Deimos.

We have run longer simulations of 2 drifting moons in orbit around Mars and also did not see tumbling excitation or large obliquity variations. However, weak resonances can only be seen in extremely long simulations with slow drift rates so longer simulations of more slowly drifting systems may yet reveal new resonant phenomena. We have not yet simulated secular phenomena such as the secular resonances that were explored by Yoder (1982). These would require longer integrations than run here and driving orbital precession frequencies or including the Sun and additional planets in the simulation. Due to their elongated shape, body axis precession and libration frequencies for Phobos and Deimos are not slow (as they are for nearly spherical bodies like the Moon and Mercury), so strong coupling via secular spin resonances would be unexpected. Our study rules out spin excitation by two of the strongest orbital resonances, however we do not necessarily rule out all types of spin resonances.

We have carried out simulations of initially tumbling bodies in eccentric orbit for bodies with axis ratios similar to Phobos and Deimos. Our simulations dissipate energy within the simulated elastic body, so tidal dissipation is directly tied to decreasing orbital semi-major axis and eccentricity. Our simulations show that even at quite low eccentricity, initially tumbling bodies remain so for thousands of orbital periods, confirming prior results by Wisdom (1987). When they finally fall into the spin synchronous state, the eccentricity can be substantially reduced from its initial value due to enhanced tidal dissipation while in the tumbling state. Depending upon the eccentricity, energy dissipation rates can be 1 to 6 orders of magnitude higher when tumbling than when tidally locked. Using our simulations we measure the mean dissipation rate in the tumbling state and find that it is weakly dependent on eccentricity.

We suspect that tidal dissipation does not necessarily cause the spinning body to leave the separatrix region or tumbling chaotic zone because this region contains attractors. From the conservative spin-orbit problem with one rotation angle at low eccentricity, we estimate that the volume of the chaotic zone is proportional to the square root of eccentricity. With eccentricity adiabatically dropping, the likelihood that a body leaves the separatrix region can be estimated from the rate of chaotic zone volume decrease. We attribute the tendency for simulations of an initially tumbling body to drop into spin-synchronous state at very low eccentricity to the insensitivity of the chaotic zone volume to eccentricity. These rough arguments were made using qualitative properties of the classic and dissipative spin-orbit dynamical problems (depending upon only a single rotation angle) and could be tested and extended to the full three dimensional spin dynamics problem.

Because of the higher tidal dissipation rates when tumbling, we support the low eccentricity scenario for Phobos proposed by Yoder (1982). Most likely Phobos has remained at low eccentricity during much of the age of the solar system. If Phobos had been at higher eccentricity, it would have spent time tumbling and it would have experienced spin out events. Its tidal evolution rate in semi-major axis would have been too fast to remain stable during the age of the Solar system. Impact followed by tumbling can decrease eccentricity, but would not increase it. Consequently orbital resonance crossing is still required to account for Phobos’ non-zero eccentricity (Yoder, 1982). As we support a low eccentricity history for Phobos, we prefer scenarios where Phobos and Deimos formed in or from a circumplanetary disk (Goldreich, 1965; Burns, 1992; Rosenblatt, 2011; Craddock, 2011; Rosenblatt et al., 2016), rather than were captured from heliocentric orbit (Burns, 1992; Lambeck, 1979). In the scenario explored by Rosenblatt et al. (2016) for Phobos’ and Deimos’ formation from a circumplanetary disk, larger moons migrated outward due to interactions with a circumplanetary disk. Interaction with the disk could also have damped orbital eccentricities, leaving them at low eccentricity when the circumplanetary disk dissipated and prior to long timescale tidal evolution.

Using the dissipation rates measured from our simulations, we find that an episode of tumbling could have increased the dissipation rate for Deimos sufficiently that its eccentricity could have been damped during the age of the Solar system. This is despite its current slow rate of tidal evolution. We estimate that both Deimos and Phobos likely experienced one or more sub-catastrophic impacts that would have caused them to tumble. As previously suggested by Wisdom (1987), Deimos’ low eccentricity may in part be due to an impact excited epoch of tumbling that then lowered its eccentricity. The formation of the Stickney crater might have excited tumbling on Phobos (Weidenschilling, 1979; Wisdom, 1987; Ramsley and Head, 2019). Because eccentricities can be substantially damped while in a long lived tumbling state, spin excitation via rare impacts is a mechanism that can reduce satellite orbital eccentricity.

We see some interesting phenomena in our simulations. We find that a tumbling body, after falling into the spin synchronous state, can spend time in a non-principal axis rotation state. These probably involve out of plane precession resonances. These states exhibit enhanced dissipation which is particularly noticeable at low orbital eccentricity. In some simulations we have also seen resonant excitation of free libration that can also be long lived. These types of excitation could cause enhanced tidal dissipation compared to that commonly predicted for a tidally locked low free libration, low obliquity body undergoing principal axis rotation.

In this paper we have focused on simulations of bodies with the axis ratios similar to Phobos and Deimos. With such extreme axis ratios, we primarily see tumbling or spin synchronous states. However, simulations of initially tumbling somewhat rounder bodies show that they can more easily be captured into and spend more time in spin-orbit resonances such as the 3:2 or 2:1 resonances. Except for a short discussion at the end of section 4.1, sensitivity to rheology was neglected in this study. This is because our simulations currently only exhibit a Kelvin-Voigt rheology and because we required long simulations, we compromised by poorly resolving the spinning bodies with our mass-spring model. Sensitivity to rheology could be explored in future studies using simulations of more accurately resolved bodies.

Acknowledgements.

This material is based upon work supported in part by NASA grant 80NSSC17K0771, and National Science Foundation Grant No. PHY-1757062. We acknowledge support from the NASA Summer Undergraduate Program for Planetary Research (SUPPR) internship program. We thank Juliana South, Tyler LaBree, John Siu, and Jonas Smucker for helpful discussions. We are grateful to Christoph Lhotka, Michael Efroimsky, Valeri Makarov, and James Kwiecinski for suggestions that inspired and significantly improved this manuscript.

Bibliography