Exoplanets Torqued by the Combined Tides of a Moon and Parent Star

It has been revealed in recent years that terrestrial extrasolar planets are basically ubiquitous in our Galaxy (Burke et al. 2015; Mulders et al. 2015), and, although extrasolar moons have yet to be discovered, certainly many of these exoplanets have moons just like the planets in our own solar system. Whether or not a planet has a moon is not just a minor curiosity but potentially fundamental to the question of whether these planets host life. In the case of the Earth, the Moon's relatively large size allows it to stabilize the Earth's obliquity. This may be crucial for stabilizing conditions on Earth for a sufficiently long time to allow life to develop (Ward & Brownlee 2000, although also see Lissauer et al. 2012; Li & Batygin 2014).

Another factor that may affect a planet's ability to host life is the presence of internal heating and volcanism. The associated tectonic activity due to the movement of plates sitting atop a fluid mantle can trap atmospheric gases such as carbon dioxide into rocks and help stabilize the climate (Walker et al. 1981; Sleep & Zahnle 2001; Foley & Driscoll 2016). Tectonic activity can also cycle fresh rock and minerals out of deeper regions of the planet, providing the building blocks and nutrients for life. Heating can provide the energy needed to drive biochemical reactions as is seen from hydrothermal vents on Earth. In the case of the Earth, this heating is driven by a combination of radioactivity, latent heat, and heat from formation; but in principle, a moon could also provide a source of heating through tidal interactions, as is seen for Io. The easiest way to drive such tidal heating would be for a moon via an eccentric orbit. In the absence of a large eccentricity though there could also be cases where a planet is being tidally torqued by both its moon and parent star. The competing torques ensure that the planet is never perfectly synchronized to either the moon or star, so that tidal heating can continue to persist.

Motivated by these possibilities, here I consider the evolution of a planet tidally torqued to a similar level by both its moon and parent star. Particular focus is on stars in the mass range of M_* ≈ 0.1–0.6 M_⊙, as they have two attractive properties: (1) there has been strong interest in studying exoplanets in the habitable zones of these stars, and (2) their habitable zones are closer-in, where tidal interactions with the star are more important. Of course, various aspects of star–planet–moon tidal interactions have been investigated previously. As early as in Counselman (1973), it was noted such systems would evolve toward one of three states, (1) the moon migrates inward until it reaches the planet, (2) the moon migrates outward until it escapes from the planet, and (3) the moon finds a stable state where its orbital frequency and the planetary spin frequency are at a mutual resonance. Since this work, the problem has continued to be studied in various ways over a wide variety of possible systems (e.g., Ward & Reid 1973; Touma & Wisdom 1994; Neron de Surgy & Laskar 1997; Barnes & O'Brien 2002; Sasaki et al. 2012; Sasaki & Barnes 2014; Adams & Bloch 2016). Here, I focus on particular aspects of this previous work by highlighting the possibility of tidal heating and the importance of the initial spin period of the planet in determining the resulting dynamics.

In Section 2, I motivate and summarize the parameter range for the star–planet–moon systems I will be considering. In Section 3, I present the set of equations I use to solve for the orbits and tides. In Section 4, I summarize the characteristic timescale of the problem and solve for the evolution of the exoplanet's spin in the limit of no orbital evolution. This provides some intuition that is useful for solving the more detailed evolution of the system. In Section 5, I explore the full evolution of the star–planet–moon system over a range of parameters, and I conclude in Section 6 with a summary.

Before diving into the details of tidal interactions, it is helpful to motivate the range of parameters that will be considered. In Figure 1, I plot the habitable zone for stars in the mass range of M_* = 0.1–0.6 M_⊙ (green shaded region). This is estimated as the range of distances at which an exoplanet would receive the same flux as at distances of 0.8–1.7 au from the Sun (this is the more optimistic range from the work of Kasting et al. 1993—a less optimistic range would be 0.95–1.4 au). In comparison to the Sun, the habitable zone must be fairly close to the parent star because the luminosity varies strongly with mass.

Zoom In Zoom Out Reset image size

If an exoplanet with an exomoon is within this range of distances from the star, it is possible that the star's gravity is sufficiently strong to unbind that moon. The critical distance within which a moon can remain is proportional to the so-called Hill sphere, R_H, so that

$$\begin{eqnarray}&&{a}_{\mathrm{crit},m}={{fR}}_{{\rm{H}}}={{fa}}_{* }{\left(\displaystyle \frac{{M}_{p}}{3{M}_{* }}\right)}^{1/3},\end{eqnarray} \tag{ 1 }$$

for which I use a constant factor f = 0.49 (Domingos et al. 2006, as appropriate for prograde orbits) for this work (although note that a value of f = 0.36 has also been argued for in the past by Holman & Wiegert 1999). This relation can be inverted to find a critical distance of the planet from the star, within which the moon would be stripped away. This is given by

$$\begin{eqnarray}{a}_{\mathrm{crit},* } & = & \displaystyle \frac{{a}_{m}}{f}{\left(\displaystyle \frac{3{M}_{* }}{{M}_{p}}\right)}^{1/3}\\ & = & 0.52\left(\displaystyle \frac{{a}_{m}}{{a}_{{\rm{M}}}}\right){\left(\displaystyle \frac{{M}_{p}}{{M}_{\oplus }}\right)}^{-1/3}{\left(\displaystyle \frac{{M}_{* }}{{M}_{\odot }}\right)}^{1/3}\mathrm{au},\end{eqnarray} \tag{ 2 }$$

where ${a}_{{\rm{M}}}=3.84\times {10}^{10}\,\mathrm{cm}$ is the distance between the Earth and Moon. In Figure 1, I plot ${a}_{\mathrm{crit},* }$ for different values of ${a}_{m}/{a}_{{\rm{M}}}$ (blue dashed lines). In all cases, I assume the mass of the exoplanet is M_p = M_⊕, where ${M}_{\oplus }=5.97\times {10}^{27}\,{\rm{g}}$ is the mass of the Earth, and the mass of the exomoon is ${M}_{m}={M}_{{\rm{M}}}$, where ${M}_{{\rm{M}}}=7.35\times {10}^{25}\,{\rm{g}}$ is the mass of the Moon. From comparing these critical distances to the habitable zone, one can see that the moon must be sufficiently close to the exoplanet to prevent being stripped. For example, for M_* = 0.5 M_⊙ and a_* = 0.3 au, then the moon must be at a distance ${a}_{m}\lesssim 0.7\,{a}_{{\rm{M}}}$ from the exoplanet to stay bound. These considerations roughly set the parameters I will be using for the initial conditions of the dynamical evolution calculations.

I next present the set of equations I will be using to solve for the dynamics of the planet's spin and the orbital separations. The basic strategy is to focus on the secular evolution of the star–planet–moon system rather than follow each orbit individually. In this sense, the conservation equations used represent averages over many orbits. This allows the evolution of these systems to be considered over much longer timescales of ∼10⁵–10¹⁰ years rather than being focused on the much shorter timescale of the orbits themselves. Further simplifications include assuming that the spin angular momentum of the planet is parallel to the orbital angular momentum of both the moon and planet.

Consider the configuration of a star, planet, and moon, with masses M_*, M_p, and M_m, respectively, as shown schematically in Figure 2. The semimajor axis of the orbit of the star and planet is a_* and the semimajor axis of the planet and moon is a_m. The orbits are all assumed to be circular. Tidal forces from the star and moon each generate two bulges on the planet on opposite sides. The frequency of these bulges in the frame of the planet are ${\hat{\omega }}_{* }=2({n}_{* }-{\rm{\Omega }})$ and ${\hat{\omega }}_{m}=2({n}_{m}-{\rm{\Omega }})$ due to the star and moon, respectively, where the factor of two represents the two bulges, Ω is the spin of the planet, and the orbital frequencies are ${n}_{* }^{2}={{GM}}_{* }/{a}_{* }^{3}$ and ${n}_{m}^{2}={{GM}}_{p}/{a}_{m}^{3}$.

Zoom In Zoom Out Reset image size

These bulges are not perfectly aligned with the position of the mass generating the bulge because the planet is not a perfect fluid. Instead, the bulge can proceed or lag behind. This is emphasized in Figure 2, where the bulge caused by the star proceeds the star by an angle α and the bulge caused by the moon lags behind by an angle β (note that these angle are exaggerated here and are much smaller in real systems). These angles depend on the "time lag" of the tidal forcing of the planet τ and tidal forcing frequencies such that $\alpha \propto \tau {\hat{\omega }}_{* }$ and $\beta \propto \tau {\hat{\omega }}_{m}$. In Figure 2, it is assumed that the planet is spinning faster than its orbit around the star, thus ${\hat{\omega }}_{* }\lt 0$ and the bulge from the star proceeds its position.

Due to the misalignment between the bulge and the object causing the tide, the planet's spin can be torqued up or down. These torques over secular timescales can be approximated as (Ogilvie 2014)

$$\begin{eqnarray}&&{N}_{* }=\displaystyle \frac{3}{2}{k}_{2}\tau {\hat{\omega }}_{* }\displaystyle \frac{{{GM}}_{* }^{2}{R}_{p}^{5}}{{a}_{* }^{6}},\end{eqnarray} \tag{ 3 }$$

and

$$\begin{eqnarray}&&{N}_{m}=\displaystyle \frac{3}{2}{k}_{2}\tau {\hat{\omega }}_{m}\displaystyle \frac{{{GM}}_{m}^{2}{R}_{p}^{5}}{{a}_{m}^{6}},\end{eqnarray} \tag{ 4 }$$

from the star and moon, respectively, where R_p is the radius of the planet and k₂ is the Love number which represents the planet's rigidity. In the specific case of Figure 2, ${N}_{* }\lt 0$ and ${N}_{m}\gt 0$ because of whether the bulge proceeds of lags. These torques can be rewritten in terms of the orbital frequencies as

$$\begin{eqnarray}&&{N}_{* }=\displaystyle \frac{3}{2}{k}_{2}\tau {\hat{\omega }}_{* }{M}_{* }{R}_{p}^{2}{\left(\displaystyle \frac{{R}_{p}}{{a}_{* }}\right)}^{3}{n}_{* }^{2},\end{eqnarray} \tag{ 5 }$$

and

$$\begin{eqnarray}&&{N}_{m}=\displaystyle \frac{3}{2}{k}_{2}\tau {\hat{\omega }}_{m}{M}_{m}{R}_{p}^{2}\left(\displaystyle \frac{{M}_{m}}{{M}_{p}}\right){\left(\displaystyle \frac{{R}_{p}}{{a}_{m}}\right)}^{3}{n}_{m}^{2}.\end{eqnarray} \tag{ 6 }$$

Note that because I am assuming that ${M}_{p}\gg {M}_{m}$, there is an additional factor of ${M}_{m}/{M}_{p}$ for N_m.

Note that for this work, I express the tidal torque using a constant time lag model (Mignard 1979, 1980, 1981; Hut 1981; Heller et al. 2011) rather than a constant geometric lag model (MacDonald 1964; Goldreich 1966; Murray & Dermott 1999), which would be parameterized with a quality factor Q. The relation between the two is ${k}_{2}\tau \hat{\omega }=\sigma {k}_{2}/Q$ with $\sigma =\mathrm{sgn}\,\hat{\omega }$ (see Efroimsky & Makarov 2013 for a discussion of the limitations and implications of different tidal implementations). I choose this formalism simply because it provides a more smooth transition from a positive torque (when the planet is spinning slowly) to zero torque (when the planet is tidal locked) to a negative torque (when the planet is spinning quickly). In contrast, if I leave Q fixed, the torque changes more abruptly as the planet's spin evolves, which causes numerical issues when evolving the orbits and spins.

Another issue is that I consider the torque by the star from the star's bulge and the torque by the moon from the moon's bulge but not vice versa (for example, the torque by the star from the moon's bulge). This is a reasonable approximation for the secular limit, because over long timescales, the moon's bulge, for example, will appear at all different positions with respect to the star and therefore will average to zero torque. If the orbits were followed individually, there would be more complicated behavior on shorter timescales, but this is beyond the secular focus of the present study.

The torques change the spin angular momentum of the planet as

$$\begin{eqnarray}&&\displaystyle \frac{d}{{dt}}({I}_{p}{\rm{\Omega }})={N}_{* }+{N}_{m},\end{eqnarray} \tag{ 7 }$$

where I_p and Ω are the moment of inertia and spin frequency of the planet, respectively. To conserve angular momentum, the orbital separations must also change. The differential equations that describe these are

$$\begin{eqnarray}&&\displaystyle \frac{d}{{dt}}[{M}_{p}{({{GM}}_{* }{a}_{* })}^{1/2}]=-{N}_{* },\end{eqnarray} \tag{ 8 }$$

and

$$\begin{eqnarray}&&\displaystyle \frac{d}{{dt}}[{M}_{m}{({{GM}}_{p}{a}_{m})}^{1/2}+{I}_{m}{n}_{m}]=-{N}_{m},\end{eqnarray} \tag{ 9 }$$

where I_m is the moment of inertia for the moon. For simplicity, I assume that the moon is tidally locked to the planet, which gives rise to the term ${I}_{m}{n}_{m}$ in Equation (9). As relatively little angular momentum is in the moon's spin, whether or not this tidal locking occurs does not impact my main conclusions. For the star, I assume that its spin is not changing appreciably from the small tidal torque of the planet and that it can be ignored.

From the above differential equations, a number of key timescales can be identified. Rewriting Equation (7), one finds

$$\begin{eqnarray}&&\displaystyle \frac{d{\rm{\Omega }}}{{dt}}=\displaystyle \frac{{n}_{* }}{{\tau }_{\mathrm{syn},* }}\left(1-\displaystyle \frac{{\rm{\Omega }}}{{n}_{* }}\right)+\displaystyle \frac{{n}_{m}}{{\tau }_{\mathrm{syn},m}}\left(1-\displaystyle \frac{{\rm{\Omega }}}{{n}_{m}}\right),\end{eqnarray} \tag{ 10 }$$

where the associated synchronization timescales are

$$\begin{eqnarray}&&{\tau }_{\mathrm{syn},* }\equiv \displaystyle \frac{\lambda }{3{k}_{2}\tau {n}_{* }}\left(\displaystyle \frac{{M}_{p}}{{M}_{* }}\right){\left(\displaystyle \frac{{a}_{* }}{{R}_{p}}\right)}^{3}{n}_{* }^{-1},\end{eqnarray} \tag{ 11 }$$

and

$$\begin{eqnarray}&&{\tau }_{\mathrm{syn},m}\equiv \displaystyle \frac{\lambda }{3{k}_{2}\tau {n}_{m}}{\left(\displaystyle \frac{{M}_{p}}{{M}_{m}}\right)}^{2}{\left(\displaystyle \frac{{a}_{m}}{{R}_{p}}\right)}^{3}{n}_{m}^{-1},\end{eqnarray} \tag{ 12 }$$

and $\lambda ={I}_{p}/{M}_{p}{R}_{p}^{2}$ is the radius of gyration for the planet.

The orbital separations also change as described by the differential equations

$$\begin{eqnarray}&&\displaystyle \frac{{{da}}_{* }}{{dt}}=-\displaystyle \frac{{a}_{* }}{{\tau }_{\mathrm{mig},* }}\left(1-\displaystyle \frac{{\rm{\Omega }}}{{n}_{* }}\right),\end{eqnarray} \tag{ 13 }$$

and

$$\begin{eqnarray}&&\displaystyle \frac{{{da}}_{m}}{{dt}}=-\displaystyle \frac{{a}_{m}}{{\tau }_{\mathrm{mig},m}}\left(1-\displaystyle \frac{{\rm{\Omega }}}{{n}_{m}}\right){\left[1-\displaystyle \frac{6}{5}{\left(\displaystyle \frac{{R}_{m}}{{a}_{m}}\right)}^{2}\right]}^{-1},\end{eqnarray} \tag{ 14 }$$

where the last term is due to synchronous spin of the moon and I assume ${I}_{m}=(2/5){M}_{m}{R}_{m}^{2}$. For the change in the orbital orbital separations,

$$\begin{eqnarray}&&{\tau }_{\mathrm{mig},* }\equiv {(6{k}_{2}\tau {n}_{* })}^{-1}\left(\displaystyle \frac{{M}_{p}}{{M}_{* }}\right){\left(\displaystyle \frac{{a}_{* }}{{R}_{p}}\right)}^{5}{n}_{* }^{-1},\end{eqnarray} \tag{ 15 }$$

and

$$\begin{eqnarray}&&{\tau }_{\mathrm{mig},m}\equiv {(6{k}_{2}\tau {n}_{m})}^{-1}\left(\displaystyle \frac{{M}_{p}}{{M}_{m}}\right){\left(\displaystyle \frac{{a}_{m}}{{R}_{p}}\right)}^{5}{n}_{m}^{-1}.\end{eqnarray} \tag{ 16 }$$

for the migration timescales.

Besides the planet's spin and orbital separations, tides will also change the eccentricity of the orbits, which is governed by the equation (Ogilvie 2014)

$$\begin{eqnarray}&&\displaystyle \frac{1}{e}\displaystyle \frac{{de}}{{dt}}=-\displaystyle \frac{3}{2}{k}_{2}\tau (18{n}_{m}-11{\rm{\Omega }})\displaystyle \frac{{M}_{m}}{{M}_{p}}{\left(\displaystyle \frac{{R}_{p}}{{a}_{m}}\right)}^{5}{n}_{m}.\end{eqnarray} \tag{ 17 }$$

This happens on a roughly similar timescale to ${\tau }_{\mathrm{mig},m}$. Nevertheless, I assume circular orbits in this work because one of my main goals is to highlight the fact that tidal heating is possible even without eccentricity.

Tidal heating can be present for circular orbits because when both the moon and star exert their tides on the planet, the planet is always forced to spin asynchronously with respect to each of them. To estimate the strength of this heating, first consider the total energy of the combined spinning planet plus gravitational interaction with a moon is

$$\begin{eqnarray}&&{E}_{m}=\displaystyle \frac{1}{2}{I}_{p}{{\rm{\Omega }}}^{2}-\displaystyle \frac{{{GM}}_{p}{M}_{m}}{2{a}_{m}}.\end{eqnarray} \tag{ 18 }$$

Taking the derivative of this expression results in

$$\begin{eqnarray}&&{\dot{E}}_{m}={I}_{p}{\rm{\Omega }}\displaystyle \frac{d{\rm{\Omega }}}{{dt}}+\displaystyle \frac{{{GM}}_{p}{M}_{m}}{2{a}_{m}^{2}}\displaystyle \frac{{{da}}_{m}}{{dt}}.\end{eqnarray} \tag{ 19 }$$

Substituting the part of the torque due to the moon's tide on $d{\rm{\Omega }}/{dt}$ and doing some algebra one finds

$$\begin{eqnarray}{\dot{E}}_{m} & = & -\displaystyle \frac{{\hat{\omega }}_{m}}{2}{N}_{m}=\displaystyle \frac{{I}_{p}{\hat{\omega }}_{m}^{2}}{4{\tau }_{\mathrm{syn},m}}\,= & \displaystyle \frac{3}{4}{k}_{2}\tau {\hat{\omega }}_{m}^{2}{M}_{m}{R}_{p}^{2}\left(\displaystyle \frac{{M}_{m}}{{M}_{p}}\right){\left(\displaystyle \frac{{R}_{p}}{{a}_{m}}\right)}^{3}{n}_{m}^{2}.\end{eqnarray} \tag{ 20 }$$

This change of energy represents the maximum amount of heating possible on the planet due to asynchronous rotation with respect to the moon. Performing a similar set of arguments on the tidal forcing from the star results in

$$\begin{eqnarray}&&{\dot{E}}_{* }=-\displaystyle \frac{{\hat{\omega }}_{* }}{2}{N}_{* }=\displaystyle \frac{{I}_{p}{\hat{\omega }}_{* }^{2}}{4{\tau }_{\mathrm{syn},* }}=\displaystyle \frac{3}{4}{k}_{2}\tau {\hat{\omega }}_{* }^{2}{M}_{* }{R}_{p}^{2}{\left(\displaystyle \frac{{R}_{p}}{{a}_{* }}\right)}^{3}{n}_{* }^{2}.\end{eqnarray} \tag{ 21 }$$

for the heating of the planet due to the star. I take the total heating rate on the exoplanet to be the sum of the two tidal heating contributions,

$$\begin{eqnarray}&&\dot{E}={\dot{E}}_{* }+{\dot{E}}_{m}.\end{eqnarray} \tag{ 22 }$$

In detail, the tidal heating can occur at different locations depending on how the tides are damped, but this requires a more sophisticated treatment of the tides that is outside the scope of this work.

The detailed evolution of the star–planet–moon system can be complicated because of the many parameters that can be varied. It is therefore helpful to consider some of the key timescales that will govern the evolution as well as solving a simpler set of evolution equations to provide some intuition on how the evolution will proceed.

4.1. Orbital and Tidal Timescales

First, we consider the four key timescales for the synchronization and migration of the orbits that were identified in Section 3. Substituting physical values, these are

$$\begin{eqnarray}{\tau }_{\mathrm{syn},* } & = & 1.8\times {10}^{10}{\left(\displaystyle \frac{{k}_{2}\tau /\lambda }{580{\rm{s}}}\right)}^{-1}\left(\displaystyle \frac{{M}_{p}}{{M}_{\oplus }}\right){\left(\displaystyle \frac{{M}_{* }}{{M}_{\odot }}\right)}^{-2}\\ & & \times {\left(\displaystyle \frac{{R}_{p}}{{R}_{\oplus }}\right)}^{-3}{\left(\displaystyle \frac{{a}_{* }}{1\mathrm{au}}\right)}^{6}\,\mathrm{years},\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}{\tau }_{\mathrm{syn},m} & = & 3.7\times {10}^{9}{\left(\displaystyle \frac{{k}_{2}\tau /\lambda }{580{\rm{s}}}\right)}^{-1}\left(\displaystyle \frac{{M}_{p}}{{M}_{\oplus }}\right){\left(\displaystyle \frac{{M}_{m}}{{M}_{{\rm{M}}}}\right)}^{-2}\\ & & \times {\left(\displaystyle \frac{{R}_{p}}{{R}_{\oplus }}\right)}^{-3}{\left(\displaystyle \frac{{a}_{m}}{{a}_{{\rm{M}}}}\right)}^{6}\,\mathrm{years},\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}{\tau }_{\mathrm{mig},* } & = & 2.5\times {10}^{17}{\left(\displaystyle \frac{{k}_{2}\tau }{191{\rm{s}}}\right)}^{-1}\left(\displaystyle \frac{{M}_{p}}{{M}_{\oplus }}\right){\left(\displaystyle \frac{{M}_{* }}{{M}_{\odot }}\right)}^{-2}\\ & & \times {\left(\displaystyle \frac{{R}_{p}}{{R}_{\oplus }}\right)}^{-5}{\left(\displaystyle \frac{{a}_{* }}{1\mathrm{au}}\right)}^{8}\,\mathrm{years},\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}{\tau }_{\mathrm{mig},m} & = & 2.5\times {10}^{11}{\left(\displaystyle \frac{{k}_{2}\tau }{191{\rm{s}}}\right)}^{-1}{\left(\displaystyle \frac{{M}_{m}}{{M}_{{\rm{M}}}}\right)}^{-1}\\ & & \times {\left(\displaystyle \frac{{R}_{p}}{{R}_{\oplus }}\right)}^{-5}{\left(\displaystyle \frac{{a}_{m}}{{a}_{{\rm{M}}}}\right)}^{8}\,\mathrm{years},\end{eqnarray} \tag{ 26 }$$

where ${R}_{\oplus }=6.37\times {10}^{8}\,\mathrm{cm}$ is the radius of the Earth, and I have used λ = 0.33 and k₂ = 0.3, which are motivated by empirical measurements and fits to the shear modulus and stiffness of the present day Earth (Williams 1994; Henning et al. 2009; Ray & Egbert 2012; Heller & Barnes 2013; Driscoll & Barnes 2015). The time lag is set to τ = 638 s (Lambeck 1977; Neron de Surgy & Laskar 1997), which for the current values of the Earth–Moon system gives a migration rate of the moon of ${{da}}_{m}/{dt}=3.8\,\mathrm{cm}\,{\mathrm{yr}}^{-1}$ and spindown rate of the planet of ${dP}/{dt}=-2.2\,\mathrm{ms}\,{\mathrm{century}}^{-1}$, where P = 2π/Ω is the planet's spin period. As a check on the general framework used here, these rates are both in agreement with the currently measured values for the Earth and Moon.

At least for the specific values used in Equations (23)–(26), the ordering of the timescales is ${\tau }_{\mathrm{syn},m}\lesssim {\tau }_{\mathrm{syn},* }\lesssim {\tau }_{\mathrm{mig},m}\ll {\tau }_{\mathrm{mig},* }$. Thus, for the Sun–Earth–Moon system, the synchronization of the Earth with the Moon's orbit is occurring the fastest of any of these processes. As emphasized in the discussion in Section 1, where things potentially get interesting is when ${\tau }_{\mathrm{syn},* }\approx {\tau }_{\mathrm{syn},m}$, so that both the tides acting from the star and moon are impacting the planet on similar timescales. Setting ${\tau }_{\mathrm{syn},* }\approx {\tau }_{\mathrm{syn},m}$ using Equations (23) and (24), one can estimate that the typical separation between the planet and moon where this occurs is

$$\begin{eqnarray}&&{a}_{m}\approx 0.5\left(\displaystyle \frac{{a}_{* }}{0.3\,\mathrm{au}}\right){\left(\displaystyle \frac{{M}_{m}}{{M}_{{\rm{M}}}}\right)}^{1/3}{\left(\displaystyle \frac{{M}_{* }}{0.5{M}_{\odot }}\right)}^{-1/3}{a}_{{\rm{M}}}.\end{eqnarray} \tag{ 27 }$$

Comparing this estimate for a_m with Figure 1, one can see that ${\tau }_{\mathrm{syn},* }\approx {\tau }_{\mathrm{syn},m}$ will naturally occur for low-mass stars over a wide range of the parameter space where moons can remain bound to planets.

4.2. Analytic Solutions without Orbital Migration

To get more intuition for what will happen to the system when ${\tau }_{\mathrm{syn},* }\approx {\tau }_{\mathrm{syn},m}$, I consider the simplified case when the orbital separations do not evolve. In other words, I make the approximation that ${\tau }_{\mathrm{mig},* },{\tau }_{\mathrm{mig},m}\gg {\tau }_{\mathrm{syn},* },{\tau }_{\mathrm{syn},m}$. It is found for this case that a particularly simple analytic solution is possible for the planet's spin evolution.

Assuming that a_m and a_* are constant, Equation (10) becomes a differential equation simply in Ω(t). Integrating this then results in

$$\begin{eqnarray}&&{\rm{\Omega }}(t)={{\rm{\Omega }}}_{\mathrm{eq}}+({{\rm{\Omega }}}_{0}-{{\rm{\Omega }}}_{\mathrm{eq}})\exp (-t/{\tau }_{\mathrm{syn}}),\end{eqnarray} \tag{ 28 }$$

where ${{\rm{\Omega }}}_{0}\equiv {\rm{\Omega }}(t=0)$ is the initial spin frequency, Ω_eq is the equilibrium spin frequency as $t\to \infty$, and the total synchronization time of the planet is defined to be

$$\begin{eqnarray}&&{\tau }_{\mathrm{syn}}\equiv \displaystyle \frac{{\tau }_{\mathrm{syn},m}{\tau }_{\mathrm{syn},* }}{{\tau }_{\mathrm{syn},m}+{\tau }_{\mathrm{syn},* }}.\end{eqnarray} \tag{ 29 }$$

The easiest way to estimate Ω_eq is to just set $d{\rm{\Omega }}/{dt}=0$ with Equation (10), and then solve for Ω, resulting in.

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{\mathrm{eq}}\equiv \displaystyle \frac{{\tau }_{\mathrm{syn},m}{n}_{* }+{\tau }_{\mathrm{syn},* }{n}_{m}}{{\tau }_{\mathrm{syn},m}+{\tau }_{\mathrm{syn},* }}.\end{eqnarray} \tag{ 30 }$$

Equation (28) shows that the planet's spin just exponentially decays on a timescale τ_syn to the final equilibrium spin of Ω_eq.

In general, this solution is not exactly correct because a_m can also potentially evolve on similar timescales as shown by Equation (26). This in turn changes ${\tau }_{\mathrm{syn},* }$ and ${\tau }_{\mathrm{syn},m}$, so that they are not constant as assumed to derive Equation (28). An example where this fails is for the Sun–Earth–Moon system, where the migration of the Moon cannot be ignored.

Nevertheless, the full evolutions of the star–planet–moon system will demonstrate that the concept of Ω_eq plays an important role. In particular, note that from Equation (30) it is apparent that Ω_eq cannot exactly equal either n_* or n_m. This means that there will be tidal forcing on the planet by both the star and the moon, and tidal heating will always play some role in a star–planet–moon system.

Now that the main background has been covered, I consider the full evolution of the star–planet–moon system. This is solved by numerically integrating forward in time Equations (10), (13), and (14), which are three coupled differential equations in the dependent variables Ω, a_*, and a_m, respectively.

5.1. Example Evolution

It is helpful to first just focus on one fiducial example that exemplifies the main features on the solutions, which is presented in Figure 3. For this specific case, I use M_* = 0.5 M_⊙, M_p = M_⊕, and ${M}_{m}={M}_{{\rm{M}}}$ (in all further examples, M_p and M_m use these values for simplicity). The initial orbital separations are ${a}_{* }\,=0.3\,\mathrm{au}$ and ${a}_{m}=0.5\,{a}_{{\rm{M}}}$. Finally, the initial spin of the planet must be chosen, in this case I use ${P}_{0}=2\pi /{{\rm{\Omega }}}_{0}=7\,\mathrm{hr}$ (note that the Earth is generally thought to have had an initial spin period of roughly 6 hrs).

Zoom In Zoom Out Reset image size

From this example, a number of important features are seen that will inform our more detailed parameter survey below. First, in the upper panel of Figure 3, I summarize the key periods of the system. These are the actual spin period of the planet P = 2π/Ω (solid red line), the equilibrium spin period related to the analytic solutions from Section 4.2 ${P}_{\mathrm{eq}}=2\pi /{{\rm{\Omega }}}_{\mathrm{eq}}$ (dashed turquoise line), $4\pi /{\hat{\omega }}_{* }$ (dotted blue line), and $4\pi /{\hat{\omega }}_{m}$ (long dashed purple line). These latter two periods are the periods of the tidal forcing from the star and moon, respectively, with an extra factor of two in each case to cancel the factor of two from the two tidal bulges. Also, for these periods the absolute value is plotted because they can be either negative or positive. At early times we see that $P\approx 4\pi /{\hat{\omega }}_{* }\approx 4\pi /{\hat{\omega }}_{m}$. This is because the fast initial spin of the planet dominates setting the tidal frequencies. Furthermore, both of the toques are negative, and the planet is spinning down. At a time of ≈2 × 10⁸ years, the planet has spun down sufficiently that ${\rm{\Omega }}\lt {n}_{m}$, which causes the torque from the moon to instead want to spin the planet up. This switch in the sign of the moon's torque can be seen in the plot of $4\pi /{\hat{\omega }}_{m}$, which shows a cusp at this time. Now the combined torques of the moon (spinning the planet up) and the star (spinning the planet down) push the planet's spin toward P ≈ P_eq.

Even though the planet has reached this spin equilibrium, the story is not over because this equilibrium is not stable. As shown by Equation (30), if spin equilibrium is reached, then the planet is not tidally locked with either the moon or the star, and thus the tidal forces still persist. Given the relative frequencies, in this situation the ordering is ${n}_{m}\gt {{\rm{\Omega }}}_{\mathrm{eq}}\gt {n}_{* }$, and therefore the planet will be spun down by the star but spun up by the moon. Thus, we find that an equilibrium is reached, but the equilibrium is not stable as summarized in Figure 4. As ${n}_{m}\gt {{\rm{\Omega }}}_{\mathrm{eq}}$, the moon torques up the planet. This then moves the moon closer to the planet to conserve angular momentum. With the new orbital separation, once again ${n}_{m}\gt {{\rm{\Omega }}}_{\mathrm{eq}}$, and thus the process repeats until the moon moves close enough to the planet to be disrupted.

Zoom In Zoom Out Reset image size

This is seen in the middle panel of Figure 3, which plots the evolution of the moon's orbital separation (red solid line). At early times, a_m increases as the planet spins down and donates its angular momentum to the moon's orbit until leveling off when P ≈ P_eq. This does not last though because of the unstable situation and eventually a_m comes crashing back to the planet. Also plotted in Figure 3 is the critical radius for the moon to be lost to the parent star ${a}_{\mathrm{crit},m}$ (black dashed line) given by Equation (1). In this particular case, ${a}_{m}\lt {a}_{\mathrm{crit},m}$ for the entire evolution, but if this critical radius is ever exceeded, one should expect the moon to be tidally stripped from the planet.

Finally, in the bottom panel of Figure 3, I summarize the tidal heating rates. At early times, this is dominated by the star, and at late times the moon. But in either case, it must remain nonzero because the combined tides always make sure the planet is asynchronous with the orbits of both the star and moon. Also plotted is the heat flux of $\approx 4\times {10}^{20}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$ that is currently emanating through the Earth (black dashed line, Dye 2012) due to a combination of radioactivity, latent heat, and heat left over from the Earth's formation. This demonstrates that the tidal heating is actually similar or exceeds this value for a timescale of ≈10⁸ years, and thus one might expect such a planet to have tectonic activity similar to the Earth (or perhaps even Io, which has an even greater heating rate) during this time. Although compared to the timescale necessary for life to develop, this appears too short (albeit, such a timescale is obviously very uncertain).

5.2. Is the Moon Disrupted or Lost?

In the example from the previous section and in Figure 3, it appears the moon will be forced to migrate into the planet rather than expelled by the star. But this example also showed that the moon was fairly close to being removed by the star if only the moon migrated out a little further (compare the red solid line and black dashed lines in the middle panel of Figure 3).

To better explore what controls the fate of the moon, in Figure 5 we plot the evolution of the star–planet–moon system for a variety of different initial spins for the planet. The initial spins vary from from periods of 3 hr (blue curves) to 13 hr (red curves). This demonstrates that the initial spin of the planet can have a dramatic affect on the fate of the moon. The reason is that the moon's migration is driven by the extraction of angular momentum from the planet's spin, and the more spin the planet has, then the further out the moon will migrate. This dependency on P₀ is not emphasized in the work of Sasaki et al. (2012), which used a constant value of P₀ = 6 hr for Earth-like planets (although see some of the discussion in Sasaki & Barnes 2014).

Zoom In Zoom Out Reset image size

Whether the moon is disrupted or lost can be addressed with some simple arguments. The initial angular momentum in the planet's spin is $2\pi \lambda {M}_{p}{R}_{p}^{2}/{P}_{0}$. If all of this angular momentum goes into the moon's orbit, with angular momentum ${M}_{m}{({{GM}}_{p}{a}_{m})}^{1/2}$, then the orbital separation of the moon would be

$$\begin{eqnarray}&&{a}_{m}=\displaystyle \frac{4{\pi }^{2}{\lambda }^{2}{M}_{p}{R}_{p}^{4}}{{{GM}}_{m}^{2}{P}_{0}^{2}}.\end{eqnarray} \tag{ 31 }$$

Equating this to ${a}_{\mathrm{crit},m}$ given by Equation (1) and solving for P₀ provides the critical initial period for the planet

$$\begin{eqnarray}{P}_{\mathrm{crit},0} & = & \displaystyle \frac{2\pi \lambda {M}_{p}^{1/2}{R}_{p}^{2}}{{({{Gfa}}_{* })}^{1/2}{M}_{m}}{\left(\displaystyle \frac{3{M}_{* }}{{M}_{p}}\right)}^{1/6},\\ & = & 5.7{\left(\displaystyle \frac{{M}_{p}}{{M}_{\oplus }}\right)}^{1/3}{\left(\displaystyle \frac{{M}_{m}}{{M}_{{\rm{M}}}}\right)}^{-1}{\left(\displaystyle \frac{{R}_{p}}{{R}_{\oplus }}\right)}^{2}\\ & & \times {\left(\displaystyle \frac{{a}_{* }}{0.3\mathrm{au}}\right)}^{-1/2}{\left(\displaystyle \frac{{M}_{* }}{0..5{M}_{\odot }}\right)}^{1/6}\mathrm{hr}.\end{eqnarray} \tag{ 32 }$$

This matches fairly closely the critical value for the initial period found from the evolutions in Figure 5.

In detail, this critical period does not exactly hold. This is because it depends on how much angular momentum the moon initially has as well. For example, in Figure 6 for each of the evolutions the initial spin period of the planet is set to P₀ = 9 hr (in Figure 5, this P₀ was found to be sufficiently long that the moon is expected to migrate back into the planet). But varying the moon's initial orbital separation from $0.2{a}_{{\rm{M}}}$ to $0.7{a}_{{\rm{M}}}$ demonstrates that if the moon is initially sufficiently far, then it will be easier to strip. This effect is not captured by Equation (32).

Zoom In Zoom Out Reset image size

Another physical mechanism that will play a role in determining the fate of the moon is the role of tidal resonances. Of particular interest may be the evection resonance, where the moon's perihelion precession rate becomes equal to the orbital period of the planet (e.g., Touma & Wisdom 1998). This can alter the moon's orbit dramatically, but it depends in detail on the rate of the orbital evolution and the width of the resonance. Due to these complications, I save a more detailed study of these processes for future work.

5.3. Fate of the Disrupted Moon

In the cases where the moon is forced to tidally migrate back toward the planet, there are likely two potential outcomes: (1) the moon is tidally disrupted, or (2) the moon directly impacts the planet. Here, it is argued that at least for a system similar to the Earth–Moon, tidal disruption is more likely.

The moon will migrate inward until its radius hits the Roche lobe, i.e., the equipotential surface where material will no longer be gravitationally bound to the moon. As long as ${M}_{m}/{M}_{p}\ll 1$, this can be approximated by the condition (Frank et al. 2002)

$$\begin{eqnarray}&&{R}_{m}=0.462{\left(1+\displaystyle \frac{{M}_{p}}{{M}_{m}}\right)}^{-1/3}{a}_{t}\approx 0.462{\left(\displaystyle \frac{{M}_{m}}{{M}_{p}}\right)}^{1/3}{a}_{t},\end{eqnarray} \tag{ 33 }$$

where a_t denotes the semimajor axis when tidal disruption occurs. At the moment of disruption, the L1 Lagrange point is located a distance ${d}_{t}\approx 0.7{({M}_{m}/{M}_{p})}^{1/3}{a}_{t}$ from the center of mass of the moon. Thus, for tidal disruption to occur, then ${a}_{t}\gtrsim {R}_{p}+{d}_{t}$, otherwise the planet and moon are too close together and direct impact will occur instead. Rewriting this inequality by making use of Equation (33),

$$\begin{eqnarray}&&2.16{\left(\displaystyle \frac{{M}_{p}}{{M}_{m}}\right)}^{1/3}\gtrsim \displaystyle \frac{{R}_{p}}{{R}_{m}}+1.51.\end{eqnarray} \tag{ 34 }$$

Simply using the properties of the Earth and Moon, the left side is 9.35 while the right side is 5.19, thus the result is tidal disruption. Another way to understand this is to multiple both sides of this expression by ${R}_{m}/{R}_{p}$, which results in

$$\begin{eqnarray}&&{\left(\displaystyle \frac{\langle {\rho }_{p}\rangle }{\langle {\rho }_{m}\rangle }\right)}^{1/3}\gtrsim 0.46+0.70\displaystyle \frac{{R}_{m}}{{R}_{p}},\end{eqnarray} \tag{ 35 }$$

where $\langle {\rho }_{p}\rangle$ and $\langle {\rho }_{m}\rangle$ are the average density of the planet and moon, respectively. Because the right-hand side of Equation (35) has a value between 0.46 and 1.16 (because ${R}_{m}\lt {R}_{p}$), this demonstrates that as long as the density of the planet is similar or greater than the moon, then tidal disruption is expected. Conversely, direct impact only occurs when the density of the planet is much less than the moon (similar arguments were presented in Metzger et al. 2012, in a somewhat different context).

Following disruption, the moon will likely form some sort of ring-like structure around the planet. This may indicate that some rocky planets around low-mass stars should be expected to have circumplanetary rings. To get some idea of what such rings may look like, consider that at the moment of disruption at a separation a_t, the orbital angular momentum of the moon is

$$\begin{eqnarray}{J}_{t} & = & {[G({M}_{p}+{M}_{m}){a}_{t}]}^{1/2}{M}_{m}\\ & \approx & 1.47{[G({M}_{p}+{M}_{m}){R}_{m}]}^{1/2}{M}_{p}^{1/6}{M}_{m}^{5/6}.\end{eqnarray} \tag{ 36 }$$

where for the second expression I have used Equation (33). For the values of the Earth–Moon system, this gives an angular momentum of ${J}_{t}\approx 6\times {10}^{40}\,\mathrm{erg}\,{\rm{s}}$. This same amount of angular momentum should be roughly stored in the resulting rings, which would be given by

$$\begin{eqnarray}&&{J}_{r}\approx {({{GM}}_{p}{R}_{r})}^{1/2}{M}_{r},\end{eqnarray} \tag{ 37 }$$

where R_r and M_r are the radius and mass of the rings, respectively. The actual mass that goes into the rings would depend on the details of the ring dynamics, but even if ≈10% of the moon's mass went into the rings it would imply a radius of ${R}_{r}\approx 2\times {10}^{11}\,\mathrm{cm}$ (and an even larger radius if the ring mass is smaller). Such a large radius cannot be maintained for the rings. Material interior to the so-called fluid Roche limit will remain in rings, where the fluid Roche limit is given by (Murray & Dermott 1999)

$$\begin{eqnarray}&&{R}_{\mathrm{FRL}}=2.46{R}_{p}{\left(\displaystyle \frac{\langle {\rho }_{p}\rangle }{\langle \rho \rangle }\right)}^{1/3},\end{eqnarray} \tag{ 38 }$$

where $\langle \rho \rangle$ is the average density of the material that makes up the rings. This results in typical ring radii of ≈2 × 10⁹ cm. Searching for such rings (using the methods described in, for example, Barnes & Fortney 2004; Ohta et al. 2009; Zuluaga et al. 2015) may verify that processes as discussed here took place in a specific system. Exterior to the fluid Roche limit, material can coalesce into a new moon that is some fraction of the mass of the moon before. This new moon could then migrate inward again and the process repeat. This suggests that the planet may go through phases where it alternatively has a moon or rings, which has actually been suggested to be the case for Mars and Phobos (Hesselbrock & Minton 2017). More work is needed to better understand the evolution and duty cycle of such rings.

5.4. Exploration of Additional Evolutions

The above discussions spell out some of the general features of the evolution, but there are many parameters that can be varied for a star–planet–moon system. Thus, here I highlight some other example evolutions to provide some sense of the diversity of potential results. In each example, I consider a range of values for P₀, as this factor has proven to be key in determining a moon's fate.

In the above examples, I set ${a}_{* }=0.3\,\mathrm{au}$ for the initial separation so that the planet would be within the habitable zone, but what if the planet is much further out? In Figure 7, I consider such a case with ${a}_{* }=0.6\,\mathrm{au}$. The two main results of this change are that (1) the moon avoids being tidally stripped by the star for a larger range of initial spin periods for the planet, and (2) it takes a much longer timescale for the moon to migrate back into the planet. For the first case, this is simply because the ability of the star to strip the moon depends most strongly on the separation, as shown by Equation (1). For the second case, this is because the tidal timescales depend on a high power of the separation. In addition, tidal heating also stays strong for a longer period of time, exceeding the current heat flux on Earth for up to ∼10⁹ years.

Zoom In Zoom Out Reset image size

In Figure 8, I consider the case of a smaller star with ${M}_{* }=0.2\,{M}_{\odot }$, which requires an initial separation of a_* = 0.1 au for the planet to be near the habitable zone. Although the mass of the star is smaller, the planet's orbital separation plays a much stronger role in setting the tidal interaction timescales. Thus, the evolution occurs much faster, in this case with the fate of the moon being decided within ∼10⁶ years. In comparison, in Figure 9, just by tripling the initial separation to a_* = 0.3 au, now the evolution can occur for 10¹⁰ years or longer.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In this work, I have investigated the tidal interactions of star–planet–moon systems. In comparison to the broader work of Sasaki et al. (2012), this study has a more specific focus on Earth–Moon-like systems around low-mass stars, which is motivated by recent surveys finding planets in the habitable zones of these stars. Furthermore, because these habitable zones are relatively close to the star, the tidal interactions between the planet and star are naturally comparable to the interactions between the planet and moon. I especially highlight the role of the initial spin period of the planet P₀ in determining the fate of the moon, and I use a constant τ formalism rather than a quality factor for assessing the impact of the tides.

Solving for the time evolution of these systems, my main conclusions are as follows.

• 1.  
  The combined tidal interactions cause the moon to eventually be stripped by the star or migrate back toward the planet.
• 2.  
  Which of these fates befall the moon depends sensitively on the initial spin period of the planet P₀, with a small spin making it more likely for the moon to be stripped.
• 3.  
  In cases where the moon migrates into the planet, the moon will be tidally disrupted rather than directly impact the planet because of the relatively similar densities of the rocky planet and moon. This may produce rings around rocky planets.
• 4.  
  The combined tidal interactions force the planet to always spin asynchronously with respect to the both the moon and star, often generating an amount of tidal heating similar to the current heat flux coming up through the Earth for up to ∼10⁹ years.
• 5.  
  The overall evolution of this system until the time the moon is stripped or disrupted depends on the initial separation of the planet and star, and can vary greatly from less than ${10}^{6}\,\mathrm{years}$ to greater than 10¹⁰ years.

In the future, as extrasolar moons are inevitably discovered, the formalism presented here can be used to better understand the lifetime and fate of these moons (although in some cases, additional planetary bodies may also impact the moon, something outside the scope of this work). If instead the presence of rings around rocky planets orbiting low-mass stars is found, it would provide evidence that processes as described here have occurred in specific systems. Alternatively, if no rings or moons are ever present, it could indicate that the moons are stripped because of the short initial spin when the planet is formed, providing insight into the planet formation process. Further calculations are needed to understand the details of these rings, how they evolve, and how long they should be present.

It is interesting to ponder the implications of the long timescale found for the loss of the moon through tidal disruption or stripping. One might assume that the fate of a moon would be determined relatively early during the formation process of the solar system and planet. But in fact, without fine tuning the parameters much, a moon could orbit a planet for well over 10¹⁰ years before being tidally disrupted or lost from the planet completely! If such a planet could harbor life, this would presumably be sufficiently long for an advanced civilization to develop, only to be subject to a catastrophic event. Nevertheless, as we have seen over the last century, technology advances quickly. And although our species still seems to be sorting out issues more local to home, one might hope that in a relatively short amount of time in comparison to astrophysical timescales, an advanced society may be able to overcome such a unique challenge.

I thank Johanna Teske for suggestions on a previous draft of this manuscript, and Jason Barnes for answering my questions about his work. I also thank Konstantin Batygin, Alexandre Correia, Michael Efroimsky, Jim Fuller, and Valeri Makarov for their feedback.