# Consequences of dynamically unstable moons in extrasolar systems

Bradley M. S. Hansen,1,
1Mani L. Bhaumik Institute for Theoretical Physics, Department of Physics and
Astronomy, University of California, Los Angeles, CA, 90095, USA
<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Moons orbiting rocky exoplanets in compact orbits about other stars experience
an accelerated tidal evolution, and can either merge with their parent planet
or reach the limit of dynamical instability within a Hubble time. We review
the parameter space over which moons become unbound, including the effects of
atmospheric tides on the planetary spin. We find that such tides can change
the final outcome from merger to escape, albeit over a limited parameter
space. We also follow the further evolution of unbound moons, and demonstrate
that the overwhelmingly most likely long-term outcome is that the unbound moon
returns to collide with its original parent planet. The dust released by such
a collision is estimated to reach optical depths $\sim 10^{-3}$, exhibit
characteristic temperatures of a few hundred degrees Kelvin, and last for a
few thousand years. These properties make such events an attractive model for
the emerging class of middle-aged main sequence stars that are observed to
show transient clouds of warm dust. Furthermore, a late collision between a
planet and a returning moon on a hyperbolic orbit may sterilise an otherwise
habitable planet.

###### keywords:

planets and satellites: detection – planets and satellites: dynamical
evolution and stability – planets and satellites: terrestrial planets –
planet-star interactions – infrared: planetary systems – (stars:) planetary
systems‘

##  1 Introduction

The planets of the Solar System are accompanied by a variety of satellites –
from the moons and rings of the giant planets to the compact moons of Earth
and Mars (Peale, 1999). The fact that such structures are relatively common in
the Solar system suggests that at least some fraction of the planets formed
around other stars must be, in turn, accompanied by their own moons and rings.
These satellites also offer insights into the processes that sculpt planetary
systems – a subject of great interest at present as observations push our
knowledge well past the parameter space occupied by the Solar system planets.
Although of great theoretical interest, the low mass ratios and small sizes
imply that most exoplanetary moon systems are likely to remain unobservable
for some time. Studies of transitting giant planets (Simon et al., 2007;
Kipping, 2009; Awiphan & Kerins, 2013; Lewis, 2013; Heller et al., 2014;
Hippke, 2015; Agol et al., 2015) are now beginning to probe the upper
plausible mass range of exomoons (Teachey & Kipping, 2018) although it is not
clear whether such masses can be explained by traditional processes (Ochiai et
al., 2014; Hamers & Portegies Zwart, 2018; Hansen, 2019; Cilibrasi et al.,
2021).

The moons of rocky planets are of particular interest. The formation of
compact, low-mass, planetary systems is still a question of active discussion.
If the late stages of planetary assembly occur in situ (Hansen & Murray, 2012,
2013; Chiang & Laughlin, 2013), then moon formation is expected to be a
natural consequence of the giant impacts that occur during the final stages of
planetary system clearing. On the other hand, if the planets migrate a
significant distance inwards to their current configurations (Ida & Lin, 2010;
Izidoro et al., 2017), then the additional torques experienced by the system
may render systems unstable (Spalding et al., 2016). In this sense, the
presence of moons in compact, extrasolar planetary systems may serve as a
marker of a late stage of in situ giant impacts. The moons of terrestrial
planets may also play a role in preserving planetary habitability. The large
lever arm of a moon orbit is potentially a stabilising force against chaotic
variations in planetary obliquity caused by resonance between planetary
precession frequencies and the secular eigenmodes of planetary systems (Laskar
et al., 1993). However, at least in the case of the Earth, it may not be
absolutely necessary (Lissauer et al., 2012).

An important consideration in extending the study of satellites to other,
compact, planetary systems is that the tidal evolution of moon systems
proceeds at a faster pace for planets that are close to their parent stars
(Barnes & O’Brien, 2002; Sasaki et al., 2012), although this depends somewhat
on the parameters of the tidal model (Piro, 2018; Tokadjian & Piro, 2020). At
some level this question is also moot, because terrestrial class planets are
themselves difficult to detect, and so the detection of the moons of
extrasolar terrestrial in place is still somewhat beyond current capabilities.
However, if these systems are short lived, there is the possibility that we
may observe the destruction of such objects, and thereby understand something
about their frequency of occurrence.

In particular, moons that spiral out to the point at which the circumplanetary
orbit becomes unstable will enter an orbit around the star with approximately
the same semi-major axis as the parent planet. Although no longer bound to the
original planet, they will still be bound to the host star. As such, they will
therefore undergo repeated close passages and scatter off the planet. For a
moon freed from a massive planet on scales of several AU, this may result in
the ejection of the moon from the system. However, for the kinds of compact
systems susceptible to tidal evolution, the potential well of the star is too
deep for planetary scattering for ejection to be a significant loss channel.
As such, nearly all of the freed exomoons will eventually collide with the
original host planet and we anticipate a dramatic release of dust would
result, and the formation of a temporary infra-red excess around the parent
star. Such a scenario may offer an explanation for the class of objects known
as EDD (Balog et al., 2009; Melis, 2016; Kral et al., 2017; Moór et al., 2021;
Melis et al., 2021) hereafter, which are middle aged stars (100 Myr or older)
that show infrared excesses well in excess of that expected from the
traditional models of planetary system assembly. It is the quantitative
aspects of this scenario that we wish to investigate here.

In § 2 we review the evolution of plausibe moons in extrasolar rocky planet
systems, and examine the expected age range over which they go unstable. In
addition to the traditional stellar tides discussed by previous authors
(Barnes & O’Brien, 2002; Sasaki et al., 2012; Tokadjian & Piro, 2020), we also
consider the effects of atmospheric tides, which have been postulated (Gold &
Soter, 1969; Correia & Laskar, 2001; Leconte et al., 2015; Cunha et al., 2015;
Auclair-Desrotour et al., 2017) to provide asynchronous equilibria around
lower mass stellar hosts. In § 3 we discuss the post-ejection evolution of
unstable moons, and consider their ultimate fates. In § 4 we discuss the
potential observational signatures of such fates.

##  2 Orbital Evolution of Moons

The basic picture for a terrestrial-type planet+moon system is that the moon
forms as the result of a collision during the late stages of planet assembly.
A vapor disk produced in the collision condenses into solid particles at the
appropriate Roche radius and these bodies rapidly co-alesce into a moon
(Stevenson, 1987; Canup et al., 2021). The collision also leaves the planet
with significant spin and the long-term evolution of the moon is determined by
the tidal coupling between the planet and the moon, with three limiting
outcomes available (Counselman, 1973). If the planet spin is too low, the
planet will extract angular momentum from the orbit and the moon will spiral
in to merge with the planet. If the planet spin is high enough, angular
momentum will be transferred to the moon and it will spiral outwards,
eventually to achieve a synchronous rotation or to spiral outwards towards
escape. Within the context of a planet orbiting a star, this outwards
evolution is truncated when the separation approaches the Hill sphere and the
orbit becomes dynamically unstable.

We now know of many planetary systems where the planets are much closer to the
host star than the planets in our own Solar system. In such cases, an
additional influence must be accounted for. The host star also raises a tide
on the planet, and so dissipation of the planetary tide also acts to transfer
angular momentum from the spin of the planet to its orbit around the star.
This acts to spin the planet down, and so can cause the outwards orbital
evolution to reverse, increasing the range of parameter space leading to
inspiral. In this scenario, achieving a synchronous state is also unlikely,
because synchronism between the planetary spin, planetary orbital frequency
and moon orbital frequency would place the synchronous orbit at the
approximate location of the Hill sphere.

###  2.1 Orbital Evolution Equations

There have been several prior calculations of the orbital evolution of moons
due to tides (Barnes & O’Brien, 2002; Sasaki et al., 2012; Piro, 2018;
Tokadjian & Piro, 2020), and they all give the same qualitative behaviour, but
can differ in quantitative ways, primarily because of variations in the
assumption of tidal strengths. As in prior works, the orbital evolution of the
moon is determined by the gravitational torque acting to transfer angular
momentum between the satellite orbit and the planetary spin

$\frac{d}{dt}\left(M_{m}(GM_{p}a_{m})^{1/2}+I_{m}n_{m}\right)=-N_{m}$ (1)

where $M_{m}$, $a_{m}$, $I_{m}$ and $n_{m}$ represent the moon mass, semi-
major axis, moment of inertia and orbital frequency respectively. The torque
$N_{m}$ is driven by dissipation in the planet and is given by

$N_{m}=\frac{3k_{2}}{2Q}\frac{GM_{m}^{2}R_{p}^{2}}{a_{m}^{6}}b_{m}(n_{m},\Omega_{p})$
(2)

where $k_{2}$ and $Q$ are the tidal Love number and quality factor
respectively. The function $b_{m}$ encodes the frequency dependance of the
tidal response and will be discussed below.

The orbital evolution of the moon is coupled to the spin evolution of the
planet, and so $N_{m}$ will also act on the spin of the planet. In addition,
planets that orbit close to their host stars will also experience a
gravitational tide due to the influence of the host star, leading to the
operation of a second torque (also driven by dissipation in the planet)

$N_{p}=\frac{3k_{2}}{2Q}\frac{GM_{*}^{2}R_{p}^{2}}{a_{p}^{6}}b_{p}(n_{p},\Omega_{p})$
(3)

where $M_{*}$ is the mass of the host star and $R_{p}$, $a_{p}$, $n_{p}$ are
the radius, semi-major axis and orbital frequency for the planet. Similarly,
the function $b_{p}$ encodes the frequency response.

In addition to the gravitational tide, planets with thick atmospheres may
experience an atmospheric tide, driven by the heating of the atmosphere by
stellar irradiation, and regulated by the forces exerted by the atmosphere on
the solid body rotation of the planet. This torque has long been held
responsible for the peculiar retrograde spin of Venus in our solar system
(Ingersoll & Dobrovolskis, 1978; Dobrovolskis & Ingersoll, 1980; Correia &
Laskar, 2001) and may also apply more generally for planets in close orbits
(Leconte et al., 2015; Cunha et al., 2015; Auclair-Desrotour et al., 2017). We
therefore include an atmospheric torque on the planetary spin

$N_{a}=-\frac{3}{2}q_{a}\frac{3M_{*}R_{p}^{3}}{5\bar{\rho}a_{p}^{3}}b_{a}(n_{p},\Omega_{p})$
(4)

where $\bar{\rho}$ is the planetary mean density, $q_{a}$ is the amplitude of
pressure field at the base of the planetary atmosphere, and $b_{a}$ represents
the frequency dependance. The sum of the torques $N_{m}$, $N_{p}$ and $N_{a}$
will determine the evolution of the planetary spin, which will, in turn, drive
the orbital evolution of the moon.

One important consideration is the amplitude and frequency dependance of the
various torques, which is determined by the physical mechanisms that underlie
the dissipation in the planet. For an Earth-like planet, the dissipation is
believed to result from the dissipation in the boundary layers of marginal
seas (Jeffreys, 1921) and resulting from pelagic turbulence (Bell, 1975;
Egbert & Ray, 2000), with a relatively weak frequency response. These
dissipation mechanisms require the existence of oceans and dry planets may
have weaker dissipation, as estimated for Mars from Lainey et al. (2007). On
the other hand, the presence of substantial Hydrogen envelopes observed on
many extrasolar planets may enable the propagation of gravity waves, causing
dissipation from the conversion of baroclinic waves into barotropic waves on
topographic features (Bell, 1975; Balmforth et al., 2002). In the light of the
various potential mechanisms, we will adopt a broad parameterisation,
calibrated in terms of the tidal Q parameter ($Q_{\oplus}=30$ serving as a
benchmark). We will assume $Q$ is frequency independent, but changes sign with
$n_{p}-\Omega_{p}$ or $n_{m}-\Omega_{p}$. Many previous treatments assume that
this occurs as a step function at zero frequency. However, studies of the
frequency response of solid bodies (Efroimsky, 2012) suggest that, at low
enough frequencies, the tidal response is driven by the internal viscosity of
the body. As such, we will employ a smoother transition (an arctangent
function), parameterised by the viscous timescale of Earth-like planets. For
the atmospheric tide model, we adopt the results of Leconte et al. (2015),
which are determined from global circulation models under different levels of
irradiation and atmospheric pressures. Our default model is for an atmospheric
surface pressure of 1 bar, and a level of irradiation appropriate to the inner
habitable zone. Appendix A describes the atmospheric torque in quantitative
terms and demonstrates the kinds of spin equilibria that can result.

In normalised form, our equation for the evolution of the moon orbit is

$\frac{dx}{d\tau}=\frac{1}{x^{7}}\left(yx^{3/2}-1\right)\bar{b}_{m}$ (5)

where $x=a_{m}/R_{\oplus}$, $y=\Omega_{p}/\Omega_{k}$, $\tau=t/t_{0}$ and

$\displaystyle\Omega_{k}$ $\displaystyle=$
$\displaystyle\left(\frac{GM_{\oplus}}{R_{\oplus}^{3}}\right)$ (6)
$\displaystyle t_{0}$ $\displaystyle=$ $\displaystyle 0.068{\rm
yrs}\frac{Q}{30}\left(\frac{M_{m}/M_{p}}{0.0123}\right)^{-1}\left(\frac{M_{p}}{M_{\oplus}}\right)^{-1/2}\left(\frac{R_{p}}{R_{\oplus}}\right)^{-5}$
(7)

with $\alpha_{p}$ is the moment of inertia constant for the planet. The
function $\bar{b}_{m}$ is

$\bar{b}_{m}=\frac{2}{\pi}\arctan\left(\gamma_{0}(x^{-3/2}-y)\right)$ (8)

where $\gamma_{0}=10\Omega_{k}\tau_{M}=10^{5}$ for $\tau_{M}=10^{10}s$ (the
characteristic viscous timescale – the Maxwell time). This function ensures a
smooth transition from positive to negative torques as the forcing frequency
drops below $1/\tau_{M}$.

Our equation for the evolution of the stellar spin is more complicated, with
three different torques contributing.

$\displaystyle\frac{dy}{d\tau}$ $\displaystyle=$
$\displaystyle\frac{\gamma_{1}}{x^{15/2}}\left(yx^{3/2}-1\right)\bar{b}_{m}+\gamma_{2}\bar{b}_{p}+$
(10)
$\displaystyle\gamma_{3}\frac{\left(y-y_{p}\right)}{1+\gamma_{4}\left(y-y_{p}\right)^{2}}$

and

$\displaystyle\gamma_{1}$ $\displaystyle=$
$\displaystyle\frac{1}{\alpha_{p}}\frac{M_{m}}{M_{p}}\left(\frac{R_{p}}{R_{\oplus}}\right)^{-2}$
(11) $\displaystyle\gamma_{2}$ $\displaystyle=$
$\displaystyle\frac{1}{\alpha_{p}}\frac{M_{*}^{2}}{M_{m}M_{p}}\left(\frac{R_{\oplus}}{a_{p}}\right)^{6}\left(\frac{R_{p}}{R_{\oplus}}\right)^{-2}$
(12) $\displaystyle\gamma_{3}$ $\displaystyle=$
$\displaystyle\frac{4\pi}{5}\frac{Q\bar{q}_{0}}{\alpha
k_{2}\bar{\omega}_{0}}\frac{M_{*}}{M_{p}}\frac{M_{\oplus}}{M_{m}}\frac{R_{\oplus}}{R_{p}}\left(\frac{a_{p}}{R_{\oplus}}\right)^{-3}$
(13) $\displaystyle\gamma_{4}$ $\displaystyle=$ $\displaystyle
1/\bar{\omega}_{0}^{2}=(\Omega_{k}/\omega_{0})^{2}$ (14)
$\displaystyle\bar{q}_{0}$ $\displaystyle=$
$\displaystyle\frac{q_{0}R_{\oplus}^{4}}{GM_{\oplus}^{2}}$ (15)
$\displaystyle\bar{b}_{p}$ $\displaystyle=$
$\displaystyle\frac{2}{\pi}\arctan\left(\gamma_{0}(y_{p}-y)\right)$ (16)

The quantities $q_{0}$ and $\omega_{0}$ are taken from Leconte et al. (2015).
Moons are assumed to form at just outside the Roche radius (taken to be $x=3$
here) and the outer limit for dynamical stability is taken to be approximately
half of the Hill sphere radius $R_{H}=a_{p}(M_{p}/3M_{*})^{1/3}$. This factor
is based on numerous prior studies (Domingos et al., 2006; Donnison, 2010) and
confirmed with our own calculations in § 3. Thus, we truncate our outward
evolution because dynamical instability occurs for orbits with or
$x>x_{i}=119(a_{p}/1AU)(M_{p}/M_{\oplus})^{1/3}(M_{*}/M_{\odot})^{-1/3}$.

###  2.2 Solutions

Our baseline calculation assumes a replica of the Earth–Moon system, orbiting
the Sun at different values of the semi-major axis. The moon evolution begins
at $x=3$, with an initial spin $y=y(0)$, which then evolves according
equations (5) and (10). The system is either evolved until the moon spirals in
to $x=3$ again (corresponding to tidal disruption), spirals out to the
dynamical instability limit, or reaches a system age of $10^{10}$ years.

####  2.2.1 Gravitational Tides only

If we omit the atmospheric tides, our system of equations describes a version
of the calculations previously calculated (Barnes & O’Brien, 2002; Sasaki et
al., 2012; Tokadjian & Piro, 2020). The simplest evolution results in the case
where the stellar tide is also weak enough to be irrelevant. In that case, the
equations have a a simple analytic solution – corresponding to basic angular
momentum conservation – given by

$y=y(0)-2\gamma_{1}\left(x^{1/2}-x(0)^{1/2}\right).$ (17)

where we assume $x(0)=3$ here. In the limit $yx^{3/2}\gg 1$, and large $y(0)$,
we can also estimate a timescale it will take to spiral out to a particular
separation $x$, namely

$\tau\sim\frac{2}{13y(0)}x^{13/2}$ (18)

and so we can estimate an outspiral time

$\displaystyle T_{esc}$ $\displaystyle=$ $\displaystyle 5.4\times 10^{11}{\rm
yrs}\left(\frac{a}{1AU}\right)^{13/2}\left(\frac{0.0123M_{\oplus}}{M_{m}/M_{p}}\right)\frac{Q^{\prime}}{30}$
(19)
$\displaystyle\left(\frac{M_{p}}{M_{\oplus}}\right)^{5/3}\left(\frac{R_{p}}{R_{\oplus}}\right)^{-5}\left(\frac{M_{*}}{M_{\odot}}\right)^{-13/6}\left(\frac{y(0)}{0.6}\right)^{-1}.$

Thus, within the context of the Earth-Moon system, the moon would survive for
$a>0.48$AU (in the absence of stellar tides).

The addition of the stellar tides bleeds angular momentum from the system and
causes the orbital evolution to reverse if $yx^{3/2}<1$. Figure 1 shows the
nature of the solutions for the gravitational tides only case, for an Earth-
like configuration located at three different semi-major axes (1 AU, 0.7 AU
and 0.4 AU). The initial spin ($y(0)=0.6$) and tidal $Q$ ($Q=30$) were chosen
to yield present-day Earth-Moon conditions (the solid point on the right). The
equations have been integrated to $10^{10}$ years, but the solid part of each
curve extends only to an age of 4.5 Gyr.

Figure 1: The black, green and blue curves show the evolution of the planetary
spin and moon semi-major axis for three systems with the Earth–Moon
parameters, but different heliocentric semi-major axes. The solid part of each
curve extends up to an age of 4.5 Gyr. The dotted part extends this to 10 Gyr.
The short dashed line indicates the locus of spin–orbit synchronism. The
horizontal long dashed lines indicate the spin at which the planet is
synchronous with its orbital frequency around the Sun (for each case). The
solid point indicates the parameters of the present Earth–Moon system.

We see that the Earth–Moon system has largely followed the simple solution in
which angular momentum conservation is the principal determinant of the
orbital properties. At a closer distance of a=0.7 AU (approximately the
distance of Venus), we see that the stellar tides are capable of spinning the
planet down to synchronism with the heliocentric orbit, and that the moon has
begun to spiral in and will eventually approach the synchronous orbit,
although it will not be tidally disrupted within 10 Gyr. This is the Type II
case described by Sasaki et al. (2012). At an even closer orbit (0.4 AU) we
see that the tidal evolution proceeds all the way to inspiral and merger.

####  2.2.2 Including Atmospheric Tides

The addition of atmospheric tides to the equations changes the spin dynamics
of the planet, because the atmospheric tides are of the opposite sign than the
corresponding gravitational tide, and allow for the possibility of
asynchronous equilibria (Correia & Laskar, 2001; Leconte et al., 2015; Cunha
et al., 2015; Auclair-Desrotour et al., 2017). Figure 2 shows the consequence
of introducing atmospheric tides into the dynamics of the planet+moon system.
Once again, we choose masses appropriate to the Terrestrial system, but all
curves are calculated at $a=0.7$AU.

The dotted curve in Figure 2 shows the locus for $dx/dt=0$, i.e. the change in
orbital evolution from inspiral to outspiral. The solid curve shows the locus
of $dy/dt=0$ for the same equivalent case shown in Figure 1 – Q=30 and no
atmospheric tides included. The red curve shows the consequence of adding our
atmospheric tide model to the system. Although there is no dramatic
topological difference, the spin reversal locus is moved to larger x, and this
can have important consequences for the system evolution.

Figure 2: The dotted curve shows the locus for which the sign of the orbital
evolution reverses, i.e. $dx/dt=0$. The solid black curve shows the parameters
for which the equivalent spin evolution reverses ($dy/dt=0$), for the case of
Q=30 and no atmospheric tides (i.e. the green curve in Figure 1). If we
include our atmospheric tide model, we get the solid red curve for $dy/dt=0$.
With the same atmospheric tide but a larger Q=100, the sign reversal locus now
shows multiple branches (the solid blue curves) corresponding to the
appearance of new, asynchronous, spin equilibria. The vertical dot-dashed line
shows the edge of the dynamically unstable region.

An even more dramatic consequence occurs if we increase Q to Q=100 (weakening
the gravitational tide and therefore enhancing the importance of the
atmospheric tide). In this case, the new asynchronous equilibria observed in
Leconte et al. (2015) – from which we derived our atmospheric tide model –
appear, and are shown by the blue curves. Of particular interest is the branch
that appears above the synchronous curve. Since the evolutionary curves start
at large y and move down as the moon spirals out, they will encounter thi sign
reversal first, before the reversal in the direction of the orbital motion. As
a consequence, the injection of angular momentum into the system by
atmospheric tides will refresh the outward orbital evolution. This is shown in
Figure 3, which shows the temporal evolution in each case for the initial
conditions used in Figure 1.

Figure 3: The black curve shows the same evolution as the green curve in
Figure 1, with Q=30 and no atmospheric tides. In the lower panel we show the
temporal evolution of the spin in each case. The green curve in this figure
shows the evolution if we keep Q=30 but include the atmospheric tide. We see
that the character of the evolution is fundamentally changed, as the moon
continues to spiral out once atmospheric tides become important. The cyan
curve shows the case for which Q=100 and atmospheric tides are included. The
solid curve extends to $10^{10}$ years, and the dotted curve extends beyond
this. We see that the moon will eventually spiral out to escape, driven by the
new, asynchronous equilibrium.

The inclusion of the atmospheric tide does indeed refresh the outwards orbital
evolution. The action of the atmospheric tide is to weaken the influence of
the tidal forces that spin the planet down, and so this means there is more
angular momentum available to drive the planet outwards. This effect becomes
stronger if the gravitational tide is weaker (Q is larger) but the evolution
is slower as well because the timescale depends on Q. This is explored further
in Appendix A.

It is important to note that, despite the increased complexity, there are
still no globally stable equilibria within the Hill sphere. Even in cases
where $dy/dt=0$ in Figure 2, there are no intersections with the $dx/dt=0$
curve, so that the orbital evolution of the moon will always pull the system
out of (or along) the spin equilibria.

###  2.3 The influence of tidal strength

The uncertainty in the nature and strength of the tidal forces at play means
that there is some uncertainty in the exact outcomes in certain parts of
parameter space. To understand the plausible range of variation, we have
calculated the evolution of an Earth–Moon analogue system (in terms of mass)
for a range of semi-major axis and spin, at both strong (Q=10) and weak
(Q=100) dissipation levels, both with and without atmospheric torques
parameterised in the models of Leconte et al. (2015). Figure 4 shows the
resulting distribution of outcomes.

Figure 4: The black points indicate systems where the moon survives for 10
Gyr. The red points indicate systems where the moon spirals out to the point
of dynamical instability, while the blue points indicate those systems where
the moon ultimately spirals in and is tidally disrupted.

For high enough initial spin, there is sufficient angular momentum for moons
to become dynamically unstable for semi-major axes $<0.7$AU (with some
dependance on the tidal parameters). For lower values of the spin, the sapping
of planetary spin by tides causes inspiral and merger. Inclusion of the
atmospheric tide increases the range over which escape is possible, although
the effects of atmospheric tides is weaker if gravitational tides are
stronger. Furthermore, the variation in the strength of atmospheric tide due
to atmospheric pressure (as parameterised in Leconte et al. (2015)) yield
small shifts in the boundaries between outcomes, but do not introduce
substantial qualitative changes.

###  2.4 Moon/Planet Mass ratio

The angular momentum balance in a moon–planet system is intimately tied to the
mass ratio between the two bodies. Figure 5 shows the effect of varying the
moon mass while keeping the planet mass fixed. We show two representative
cases, corresponding to the upper left and lower right panels in Figure 4.
These represent the two extermes of tidal interaction. The case of Q=10 and no
atmospheric tides shows systems with efficient tidal interactions, while the
case of Q=100 and atmospheric tides included represent a less efficient tidal
evolution. In this case, all initial spins were set to $y(0)=0.6$, so we probe
how the outcomes vary as a function of semi-major axis and moon/planet mass
ratio.

Figure 5: The upper left panel shows the eventual outcome for moons of mass
$M_{m}$, orbiting a $1M_{\oplus}$ planet around a $1M_{\odot}$ star, at the
stated semi-major axis. In all cases, the initial planetary spin $y(0)=0.6$.
As before, red represents dynamical instability, blue indicates inspiral and
merger, and black represents a moon that survives for 10 Gyr. In this panel,
only gravitational tides are included and are parameterised by Q=10. In the
lower left panel, we show equivalent models except that the gravitational
tides are parameterised by $Q=100$ and atmospheric tides are included. The
four open magenta circles represent four models whose evolution is shown in
the right hand panel, labelled by the moon mass in Lunar masses. In two cases,
the moons reach dynamical instability and in two others, they result in
inspiral and merger.

We see that dynamical instability occurs only for moon/planet mass ratios
similar to those of the terrestrial system. In the case of very massive moons,
this is easy to understand, since the initial angular momentum budget in the
spin of the planet is fixed for all of these models. There is simply not
enough angular momentum to drive a massive moon out to Hill sphere distances,
and the planet eventually reaches synchronism and drags the moon inwards. Less
intuitively, we find that moons that are too small also move inwards. This
results from the fact that small mass moons are very sensitive to changes in
the angular momentum budget due to planetary spin-down. As we see from the red
curve in the right panel of Figure 5, this evolution initially begins at the
shallowest slope of all, but eventually dips down faster as the stellar tides
reduce the planetary spin. We see that the moon mass is too small to halt the
evolution at the synchronous curve and the spin drops down all the way to the
stellar synchronism line, only rising again when the moon spirals in close
enough to make its tidal influence significant.

###  2.5 Planetary Mass Range

The potential observability of the catastrophic destruction of a moon will be
enhanced if the moon is larger. Thus, we may ask which combination of
plausible parameters yields the largest moons that can escape. We have already
seen that we cannot increase the mass of the moon too much, relative to the
planet, before the angular momentum required to reach the instability limit
becomes prohibitive (§ 2.4). However, we can try to scale up the mass of the
planet too.

Although the known planets span several orders of magnitude in mass, our
interest here is primarily in rocky planets. The moons of gas giant or Neptune
class planets – those with substantial gaseous envelopes – may still
experience collisions with returning escaped moons, but the impact of a rocky
body into a gaseous atmosphere is likely to excavate little debris that will
later be observed as a dust population. Instead, we might expect a brief
optical flash from the impact, but the bulk of the impact energy is expected
to be absorbed by the planet.

Figure 6: The red points show moons that eventually evolved to the point of
dynamical instability, while the blue points evolved to the point of inspiral
and merger. All systems orbit a planet of mass $M_{p}=10M_{\oplus}$ and
$R_{p}=1.78R_{\oplus}$, that is hosted by a 1 $M_{\odot}$ star. Filled circles
reach their endpoints at ages $>0.1$ Gyr, while open circles end at earlier
ages. The density of points gets lower at larger separations, as many of the
moons at these distances survive for a Hubble time and are therefore not
plotted here.

It is therefore only for rocky bodies that we expect a significant excavation
of debris. The wide variety of exoplanet properties doesn’t lend itself to a
clear boundary between rocky and non-rocky bodies, with the observed mass-
radius distribution suggesting a continuum of Hydrogen mass fractions ranging
from negligible to dominant. A significant Hydrogen envelope contribution will
also affect the radius, which is the strongest parameter dependance in
equation (19). Inspection of the mass–radius relation for known extrasolar
planet population (e.g Weiss & Marcy, 2014; Martinez et al., 2019; Swain et
al., 2019), restricted to planets with well measured radii and masses,
suggests rocky planets can extend up to masses $\sim 10M_{\oplus}$ and radii
$\sim 1.78R_{\oplus}$.

Thus, we wish to examine the outcomes for moons orbiting massive rocky
planets. Figure 6 shows the outcomes for a series of moon systems around a
planet with $M_{p}=10M_{\oplus}$ and $R_{p}=1.78R_{\oplus}$, assuming $Q=100$
and including the atmospheric tide model. We allow the semi-major axis to vary
between 0.2–0.8 AU, the moon mass to vary between 1 and 10 Lunar masses, and
the initial spin $y(0)$ from 0.2–0.9. Figure 6 shows all cases which result in
dynamical instability (red) and in merger (blue)

These results demonstrate that we can get dynamical instability for moons up
to ten times the mass of the moon, although these require close-in planets and
the moons end their evolution quite quickly. If we also require that the
dynamical instability occur at late times ($>100$ Myr), then the maximum moon
mass is more like $\sim 0.07M_{\oplus}$.

###  2.6 Stellar Mass range

The mass of the host star also plays an important role in setting both the
Hill radius and the strength of the tide acting on the star. At fixed semi-
major axis, a more massive star exerts a stronger gravitational influence and
drives tidal and dynamical evolution on a faster timescale. However, lower
mass stars are also less luminous, and so the semi-major axis at fixed
luminosity (such as would be relevant if one were to track the habitable zone
as a function of mass – Kopparapu et al. (2013)) decreases rapidly. The end
result is that an Earth-analogue, located in the habitable zone of stars of
lower mass than the Sun, must lie closer to the star and the tidal forces lead
to an acceleration in their evolution (Piro, 2018). Figure 7 demonstrates this
for the cases of an approximately Earth-like tidal dissipation (Q=30 and no
atmospheric tide) and also for a case where the atmospheric tide is important
(Q=100). In each case, the semi-major axis was chosen to match Earth-like
levels of irradiation, given the host star – using the models of Kopparapu et
al. (2013).

Figure 7: The upper panel shows the evolution of an Earth–moon system located
around stars of various mass, with semi-major axes chosen so as to produce the
same level of irradiation as the Solar case. The colours correspond to semi-
major axes if 1.0 (black), 0.9 (red), 0.8 (green), 0.7 (cyan), 0.6 (blue), 0.5
(magenta) and 0.4 AU (yellow). The dotted line indicates a synchronous moon
orbit. For 0.7–1 $M_{\odot}$, the moon survives to the present day, while it
spirals in and merges at 0.6 AU. Interior to this the moon reaches the point
of dynamical instability. The system age at merger/instability is listed in
Myr. The lower panel is equivalent except that the planetary tidal model now
uses Q=100 and incorporates atmospheric tides.

We see that, for the more massive host stars, the tidal forces are weak enough
that the moon can survive for 10 Gyr without disruption. At the lower end of
this mass range, the planet does pass through a phase of synchronism with its
orbital period, and the moon begins to spiral inwards. For host stars $\sim
0.6M_{\odot}$, this inspiral is able to proceed to disruption and merger
within 2 Gyr. At smaller masses, the moon becomes dynamically unstable before
synchronism is reached. This is because the dynamical threshold moves inwards
quite rapidly as the host star mass decreases, because the habitable zone
distance is a strong function of mass due to the dependance on stellar
luminosity. The timescale to instability also decreases as the host star mass
decreases.

This distinction in timescale is relevant for the observational identification
of mergers or collisions between moons and planets. Planetary systems undergo
an extended period of scattering and disruption during the process of clearing
out the debris from planetary assembly, and any signature on timescales $<100$
Myr are likely to be confused with the production of dust during the standard
debris disk clearing stage. However, those systems that take several hundred
million years, or more, to complete their evolution will be easily
distinguished, as middle-aged stars are not expected to produce significant
amounts of dust.

These late-stage events may also be relevant to the habitability of planets
themselves. After all, each system in Figure 7 was constructed to have Earth-
like gravity and Earth-like temperatures, although changes due to tidal
despinning and different host star spectral energy distributions are certainly
relevant. Nevertheless, on these timescales $>$100 Myr, a planet may have
started its biological evolution and taken the first steps towards the
construction of an Earth-analogue biosphere. However, tidally induced mergers,
or hyperbolic moon collisions (see next section) may provide a late-type
catastrophic event that could quench any start in terms of biological
evolution. For planets in the habitable zones of these stars, the presence of
a moon may eventually wipe life out, instead of preserving it.

##  3 Final Evolution

In § 2 we have reviewed the tidal evolution of the moons orbiting extrasolar
rocky planets. As other authors have noted, the tidal evolution of such
systems proceeds faster if the planet is located closer to the parent star and
many moon systems are not expected to survive to the present day. Some are
removed by tidally driven inspiral and disruption, while others spiral out to
the point of dynamical instability, at which point they are expected to become
unbound with respect to the planet. We have also included, for the first time,
the effects of atmospheric tides on the moon evolution. We find that this
weakens the tidal evolution in some systems, and can increase the fraction of
systems that become dynamically unstable.

Most prior calculations stop at this point, but, if we wish to understand the
potential observability of these systems, we must consider what happens to the
moons after they become dynamically unstable. Although they are no longer
bound to the original host planet, they are still bound to the original host
star, and we wish to follow their dynamical evolution beyond the point of
instability.

###  3.1 Dynamical Instability

A moon that reaches the dynamical instability limit will transition to a
heliocentric orbit, with a similar semi-major axis as the original host
planet. This guarantees close passages with the host planet in the future and
the cumulative effect of the gravitational scattering that result will drive
orbital evolution. For a more massive planet, with a larger semi-major axis,
such a process can eventually eject bodies from the system, or deposit them in
the Oort cloud. For a terrestrial class planet, deep in a stellar potential
well, the more likely outcome is a planetary collision, with a high
probability that it occurs with the original host planet.

To examine the details of this dynamical disruption process, we have
integrated the orbits of potentially unstable moons using the direct
integration code outlined in Hansen (2019). We begin the moon with circular,
coplanar, orbits about the planet with semi-major axes in the range 0.5–0.6
Hill radii, for a planet mass $1M_{\oplus}$, in a 1 AU circular orbit about a
Solar mass star. The moon mass is assumed to be negligible, and the orbit is
integrated for 1600 years. At this point, the moons are unbound and the
resulting heliocentric orbital elements are used to start a direct integration
of the heliocentric motion of the test particle moons using the N-body
integrator Mercury (Chambers, 1999). The orbits are then integrated until the
test particles are lost, mostly by collision with the planet. A similar
calculation of this phase of evolution was performed by Sucerquia et al.
(2019), which did find some surviving moons, but the orbital integration
lasted only for 0.5 Myr, which is far too short to assess the long-term
stability.

###  3.2 Moon–Planet Collisions

Indeed, all of our ejected moon candidates eventually collide with the
original host planet. Most collide within 10 Myr, suggesting that the post-
instability evolution is a minor contributor to the total age of the system
when the collision occurs. In order to understand the nature of these
moon/planet collisions, Figure 8 shows the relative velocity of moon and
planet during encounters, as a function of closest distance in each passage
($b$). We see that most encounters occur at relative velocities below that of
the local escape velocity. This is particularly true for $b<1R_{\oplus}$, i.e.
the terminal close passages that result in a collision. The amount of mass
excavated in a collision will be determined by the energy of the impact, and
the fact that this is not significantly above escape velocity suggests that
the mass released is likely to be less than the original moon mass.

Figure 8: Each point shows the relative velocity of the moon and planet during
a post-escape close enounter. The ordinate $b$ is the distance of close
approach in each encounter. If $b<1R_{\oplus}$, then it represents a true
physical collision. The dotted curve represents the escape velocity from a
$1M_{\oplus}$ mass at distance $b$. Open circles represent a planet with semi-
major axis of 1 AU, while the filled circles show a planet at 0.5 AU. The
sample shown here represents 20 realisations of an escaping moon for each
case.

This may, in part, be due to the coincidence that the escape velocity from the
surface of Earth is only a factor $\sim 3$ smaller than the heliocentric
orbital velocity (and so comparable to the relative velocity of two bodies
with similar semi-major axes). We therefore also show the case where the
planet is closer to the star (0.5 AU) in Figure 8. The relative velocity is
still dominated by the gravitational focussing due to the planet (as evidence
by the fact that the shape of the distribution follows the dotted curve so
closely). This is ultimately due to the fact that the moons diffuse out of the
instability region, rather than getting a significant ‘kick’, and so their
relative velocities with respect to the planet remain quite small, regardless
of the semi-major axis. In principle, it might be possible to increase the
relative velocities if we reduce the semi-major axis further, but then the
lifetime of the moon phase becomes comparable to the lifetime of the original
planetary assembly phase and so because hard to distinguish observationally.

###  3.3 Multiple Planet Systems

The parameters of Figure 8 reflect those of an isolated Earth–mass planet on a
circular orbit. Many of the real exoplanet systems contain multiple planets
and finite eccentricities. These may, in principle, induce larger relative
velocities by introducing larger eccentricities. To examine the effects of
this, let us consider an analogue of the Kepler-62 system. Although Kepler has
discovered many multi-planet systems, most are too compact for this
comparison, as they occur at semi-major axes where the moon evolutionary
timescale is too short (see Figure 6). However, Kepler-62e and Kepler-62f have
semi-major axes of 0.427 AU and 0.718 AU respectively. We take Kepler-62e to
have mass 7 $M_{\oplus}$ and radius $1.7R_{\oplus}$. In Figure 9 we show the
close passage velocities of test particle moons released from Kepler-62e. The
open circles represent the case of the single isolated planet, while the
filled circles indicate the case where Kepler-62e orbits together with two
additional planets, Kepler-62 d (0.120 AU, 10$M_{\oplus}$ and e=0.1) and
Kepler-62 f (0.718 AU, 7$M_{\oplus}$ and e=0.1). This represents the kind of
secular evolution experienced by a compact planetary system.

Figure 9: Each point shows the relative velocity of the moon and planet during
a post-escape close enounter. The ordinate $b$ is the distance of close
approach in each encounter. If $b<1.7R_{\oplus}$, then it represents a true
physical collision. The dotted curve represents the escape velocity from a
$7M_{\oplus}$ mass at distance $b$. Open circles represent a planet with the
properties of Kepler-62e orbiting the host star alone. The filled circles show
the evolution of the same escaping moons in a system that now contains three
planets – Kepler-62d, Kepler-62e and Kepler-62f, with the parameters given in
the text.

Figure 9 shows that there is some more spread in velocities in the multi-
planet case, but not sufficient to make a qualitative difference in the
behaviour of collisions. The lifetimes of the ejected moons are also shorter
in the multi-planet system, as the precession induced by the additional
contributions to the gravitational potential shorten the intervals between
close encounters. Our fundamental conclusion is that almost all moons that are
lost through dynamical instability will end up colliding with their original
host planet, at velocities of similar order of magnitude to that of the escape
velocity from the planet.

###  3.4 Merger

In cases where the stellar tides come to dominate the evolution, the planet
reverses its migration and starts to spiral inwards, spinning the planet up.
The eventual outcome of this process is a merger or tidal disruption. The most
likely outcome for this process is tidal disruption, because moons generated
by impacts tend to be made of lower density crustal material from the host
planet. This disruption may result in the formation of a set of rings (Piro,
2018) or a set of smaller moons (Hesselbrock & Minton, 2017). Such an episode
is likely to be far less dramatic in terms of dust production than the high
speed impact associated with a collision.

##  4 Observability

The essential consequence of the previous sections is that moons formed around
rocky planets on scales $\sim 0.4$–0.8 AU are likely to become dynamically
unstable within a Hubble time. When such moons are lost from their parent
systems, they are overwhelmingly likely to collide again with their parent
planets, on a much shorter timescale than their tidal evolution. Thus, the
most likely observational consequence is the event associated with the
collision.

Collisions between rocky bodies are expected to produce debris that will
collide and grind down to dust, which absorbs and reprocesses stellar light,
leading to an infrared excess. This is a generic feature of young planet
forming systems, as collisions between rocky bodies are integral part of the
planet formation process. However, such excesses become increasingly rare with
advancing stellar age, as planetary systems mostly complete their assembly
within 100 Myr. Thus, excesses on timescales of Gyr or greater are an
indication of some additional processes causing collisions.

Indeed, there are reports of such signatures amongst a subset of older stars
(Rhee et al., 2008; Melis, 2016; Melis et al., 2021), where excesses of near
infrared emission are reported around stars with estimated ages of 100 Myr or
more. In order to estimate the size of the infrared excess expected from our
model, we note that the impact velocities $V_{imp}$ calculated in Figure 8 and
9 are of order $\sim 70\%$–$90\%$ of the escape velocities from the surface of
the planet. Thus, equating incoming kinetic energy of an impactor of mass
$M_{imp}$ with the amount of energy required to unbind an amount of material
$\Delta M$ from the planet,

$\frac{1}{2}M_{imp}V_{imp}^{2}\sim\alpha_{bind}\frac{GM_{p}\Delta M}{R_{p}}$
(20)

we find that

$\frac{\Delta
M}{M_{imp}}\sim\frac{1}{2\alpha_{bind}}\frac{V_{imp}^{2}R_{p}}{GM_{p}}$ (21)

which is of order unity, since $V_{imp}\sim(GM_{p}/R_{p})^{1/2}$ and
$\alpha_{bind}\sim 1/2$. Thus, an impactor like our moon should generate a
mass $\sim 0.01M_{\oplus}$, with possibly even a factor of 10 larger in the
case of more massive planet–moon systems.

If this dust is broken into particles of size $s$, it could, in principle,
cover an area

$A_{IR}\sim 1.7\times
10^{28}cm^{2}\left(\frac{M_{dust}}{0.01M_{\oplus}}\right)\left(\frac{\rho}{3g/cm^{3}}\right)^{-1}\left(\frac{s}{10\mu
m}\right)^{-1}$ (22)

which could render a star completely opaque out to semi-major axes $\sim 2$AU.
However, the generation of small dust requires that the debris released from
the collision grind down in a collisional cascade, so the amount of dust
present in the system at any instant is determined by the rate at which
material is ground down.

If the moon debris is initially broken up into pieces with characteristic
radius $R\sim 1$km, then we have approximately

$N_{rocks}\sim 2.9\times
10^{9}\frac{M_{dust}}{0.01M_{\oplus}}\left(\frac{R}{1km}\right)^{-3}$ (23)

each with a mass $\sim 2\times 10^{16}g(R/1km)^{3}$.

We can estimate the collision rate of this debris by assuming they are spread
around a torus of semi-major axis comparable to the planetary orbit and
thickness comparable to the planetary Hill sphere radius. Relative velocities
between dust grains are estimated to be $\sim 10\%$ the local circular
velocity. This yields a collision rate

$\Gamma\sim 3.9\times 10^{-4}{\rm
yr^{-1}}\frac{M_{dust}}{0.01M_{\oplus}}\frac{1km}{R}\left(\frac{a}{1AU}\right)^{-3.5}\left(\frac{M_{p}}{M_{\oplus}}\right)^{-2/3}\left(\frac{M_{*}}{1M_{\odot}}\right)^{7/6}$
(24)

Note that the collision rate increases as $R$ goes down, so we can assume that
the collisional cascade proceeds rapidly to dust size scales. This gives us a
rate at which the collisional cascade generates particles of size $s=10\mu m$,
namely

$\dot{N}\sim 1.8\times
10^{30}yr^{-1}\left(\frac{M_{dust}}{0.01M_{\oplus}}\right)^{2}\frac{1km}{R}\left(\frac{1AU}{a}\right)^{3.5}\left(\frac{M_{\oplus}}{M_{p}}\right)^{2/3}\left(\frac{10\mu
m}{s}\right)^{3}$ (25)

The amount of dust mass present in the system will be regulated by the
processes that remove dust from the system. For dust of size $s$, in orbit at
semi-major axis $a$ around a star of luminosity $L$, the timescale to spiral
inwards due to the Poynting-Robertson effect is

$T_{PR}\sim 8.4\times 10^{4}{\rm
yrs}\left(\frac{a}{1AU}\right)^{2}\frac{s}{10\mu m}\frac{L_{\odot}}{L}$ (26)

The dust can also undergo its own collisional evolution and grind down to the
point where it is removed by radiation pressure. Following the arguments of
Zuckerman & Song (Zuckerman & Song, 2012), we derive an even shorter lifetime

$T_{coll}\sim\frac{1.6\times 10^{-3}}{\tau}{\rm
yrs}\left(\frac{M_{p}}{M_{\oplus}}\right)^{1/3}\left(\frac{M_{*}}{M_{\odot}}\right)^{-5/6}\left(\frac{a}{1AU}\right)^{3/2}$
(27)

where $\tau$ is the total optical depth of the dust cloud around the star.
This estimate follows that of Zuckerman & Song (2012) except that the width of
the torus of dust is here assumed to be the Hill radius of the planet, since
the dust has a localised origin instead of being generated by a primordial
belt of planetesimals.

If we multiply the dust generation rate $\dot{N}$ by the lifetime $T_{coll}$,
we find that the overall optical depth of the dust population generated by
this collision is

$\displaystyle\tau$ $\displaystyle\sim$ $\displaystyle 2\times
10^{-3}\frac{M_{dust}}{0.01M_{\oplus}}\left(\frac{R}{1km}\right)^{-1/2}\left(\frac{a}{1AU}\right)^{-2}$
(28)
$\displaystyle\left(\frac{M_{p}}{1M_{\oplus}}\right)^{-1/6}\left(\frac{M_{*}}{M_{\odot}}\right)^{1/6}\left(\frac{s}{10\mu
m}\right)^{-1/2}.$

The expected lifetime of the infrared excess can be obtained by estimating the
time required to consume the reservoir of N objects created by the collision,
assuming the collision rate $\Gamma$. This leads to a lifetime

$T_{IR}\sim 2.4\times 10^{3}{\rm
yrs}\left(\frac{M_{dust}}{0.01M_{\oplus}}\right)^{-1}\frac{R}{1km}\left(\frac{a}{1AU}\right)^{3.5}\left(\frac{M_{p}}{M_{\oplus}}\right)^{2/3}\left(\frac{M_{*}}{1M_{\odot}}\right)^{-7/6}$
(29)

The location of the dust will also determine the temperature of reradiated
emission, by virtue of the equibrium between absorbed and emitted flux

$T_{dust}\sim
280K\frac{T_{eff}}{5780K}\left(\frac{R_{*}}{R_{\odot}}\right)^{1/2}\left(\frac{a}{1AU}\right)^{-1/2}$
(30)

For dust on scales $\sim$0.4–0.8 AU around Sun-like stars, temperatures span
300–450 K, falling into the category of ‘Warm dust’.

The age of the system when the moon becomes unstable is a strong function of
planetary semi-major axis, and can range from $<10$Myr to $>1$Gyr (see Figure
6 for example). As such, those systems for which the tidal evolution is
sufficiently rapid will likely be indistinguishable from a more traditional
scenario in which the dust is released by collisions during the actual
assembly of the planetary system. However, the age distribution of systems
with unbound moon collisions is expected to have a longer tail to Gyr-level
ages. To illustrate this, we can estimate the frequency of rocky planetary
systems as a function of semi-major axis as $df/da\sim$constant, based on the
empirical observation that the distribution with period is $df/dP\sim
P^{-1/3}$(Youdin, 2011). Given the strong dependance of $T_{esc}$ on $a$
(equation 19), we anticipate that the distribution of systems with escape time
$T$ should scale as $f(T)\propto T^{-11/13}$, which falls off much slower than
the usual exponential that is fitted to the lifetimes of young star infrared
excesses.

The features of warm temperatures, large optical depths and potentially ages
up to several Gyr make this a natural scenario for the origin of dust observed
around older Solar-type stars. The frequency of occurrence of such excesses
drops with system age, as expected for a phenomenon associated with the final
accumulation of rocky planetary systems. As such, we expect most systems with
extreme debris disks to be younger than 100 Myr. Nevertheless, there are a
handful of objects which are demonstrably much older, such as BD+20 307 (Rhee
et al., 2008; Weinberger et al., 2011), $\eta$ Corvi (Lisse et al., 2012) or
TYC 8830 410 1 (Melis et al., 2021). To estimate the frequency of such
objects, we note that the survey of Moór et al. (2021) found an upper limit of
6 candidates out of 78650 sources identified as F5–G9 sources with ages $>100$
Myr. This indicates an occurrence rate $\sim 8\times 10^{-5}$, or a
characteristic lifetime $\sim 8\times 10^{3}$ years or less. This matches well
with the lifetime estimated in equation (29).

These late-time excesses have spurred suggestions of late time asteroid belt
collisions, comet disruptions or analogues of the Late Heavy Bombardment. Here
we suggest that the tidal evolution, dynamical instability and subsequent
collision of the unbound moon offers an alternative scenario that matches the
temperatures and optical depths observed for the EDD sample, while also
offering a natural explanation of the timescale – a feature lacking from the
other scenarios.

##  5 Discussion

Satellite systems can teach us much about the origins of planets, but are
usually well below the detection thresholds of present day astronomical
instrumentation. Direct observations are now beginning to probe the most
massive satellites orbiting extrasolar giant planets, but direct observation
of the moons of terrestrial class exoplanets will probably remain impossible
for some time. Nevertheless, indirect observation of moons or the consequences
of their presence, may provide some information on this topic. In this spirit
we have examined the evolution of moons of rocky exoplanets, subject to both
gravitational and atmospheric tides, and also followed this evolution to the
ultimate fate of the moon.

  1. 1.

As has been shown before, over much of the extant parameter space, the tidal
evolution of moons brings them to the point of orbital instability, with the
moon becoming unbound from the original host planet and entering a
heliocentric orbit about the host star. We have shown that the inclusion of
atmospheric tides aids this process by opposing the gravitational torque from
the star and thereby preserving more of the original angular momentum to drive
the outward orbital evolution of the moon.

  2. 2.

We have shown that the overwhelmingly most likely outcome of this orbital
instability is that the moon will return, in short order, to collide with the
original host planet on a near-parabolic orbit. Such impacts are expected to
generate a cloud of dust that can lead to an infrared excess. For certain
parts of the parameter space, these excesses may occur for ages of several Gyr
– much later than expected for the collisions associated with rocky planet
assembly – and so should be identifiable on the basis of their unusual time of
occurrence.

  3. 3.

Indeed, there is an emerging class of objects, entitled Extreme Debris Disks
(EDD) which exhibit all the signatures we associate with these late time,
moon–planet collisions. They occur around stars of intermediate to old age,
they show dust temperatures which place them within an AU of the host star (so
that tides will operate fast enough) and show optical depths $\tau\sim
10^{-4}$–$10^{-3}$, consistent with the amount of dust expected from the
collision of a moon-sized object. Other scenarios have been proposed to
generate the collisions required to explain EDD (Wyatt, 2008; Melis, 2016;
Melis et al., 2021), but all require some level of coincidence to trigger the
phenomenon – the onset of a late time dynamical instability (which would
require a special set of initial conditions) or the injection of a comet from
an external debris belt (which would require another perturber to operate –
(e.g. Wyatt et al., 2007)). The moon-based scenario has the distinct advantage
of inevitability, in that a moon, once formed, will evolve inexorably
according to the angular momentum budget of the system and the strength of the
tides operating on it.

  4. 4.

If EDD are the result of collisions of tidally unbound moons with planets,
then they attest to the occurrence of moons in extrasolar systems containing
rocky planets. Although such moons are a natural occurrence in formation
scenarios that include a phase of late stage in situ assembly, they are not
guaranteed in scenarios where the planets migrate inwards from more distant
locations. This migration is expected to occur during the epoch when the
protoplanetary gas disk remains and the torques exerted by the gas disk may
render such orbits unstable. The existence of moons in extrasolar planet
systems would argue that there was a late stage of planetary collisions in
these systems, either due to a genuine episode of in situ assembly (Hansen &
Murray, 2012, 2013; Chiang & Laughlin, 2013) or as the result of the dynamical
instability of compact resonant chains (Izidoro et al., 2017).

  5. 5.

The presence of the moon in the Terrestrial system is believed to aid the
habitability of the Earth by stabilising it against obliquity variations,
although some authors suggest the influence is more neutral. However, in those
systems which form a moon that later goes dynamically unstable, the formation
of moons may point in the completely opposite direction – this may ultimately
prove to be harmful for life. The return of an unbound moon on timescales
$\sim$Gyr means that the planet would experience a catastrophic impact. Such a
collision would likely sterilise a planet that had begun to form life in the
time it took for the moon to leave and return. Whether life could then start
again after this collision is uncertain, especially if such late time
collisions remove much of the water from a potentially habitable planet.

Data availability: The data underlying this article will be shared on
reasonable request to the corresponding author.

The author appreciates the constructive comments of the referee. This research
has made use of the NASA Exoplanet Archive, which is operated by the
California Institute of Technology, under contract with the National
Aeronautics and Space Administration under the Exoplanet Exploration Program.
This research has made use of NASA’s Astrophysics Data System, which is
operated by the Smithsonian Astrophysical Observatory under NASA Cooperative
Agreement 80NSSC21M0056.

## References

  * Agol et al. (2015) Agol E., Jansen T., Lacy B., Robinson T. D., Meadows V., 2015, ApJ, 812, 5
  * Auclair-Desrotour et al. (2017) Auclair-Desrotour P., Laskar J., Mathis S., Correia A. C. M., 2017, A&A, 603, A108
  * Awiphan & Kerins (2013) Awiphan S., Kerins E., 2013, MNRAS, 432, 2549
  * Balmforth et al. (2002) Balmforth N. J., Ierley G. R., Young W. R., 2002, Journal of Physical Oceanography, 32, 2900
  * Balog et al. (2009) Balog Z., Kiss L. L., Vinkó J., Rieke G. H., Muzerolle J., Gáspár A., Young E. T., Gorlova N., 2009, ApJ, 698, 1989
  * Barnes & O’Brien (2002) Barnes J. W., O’Brien D. P., 2002, ApJ, 575, 1087
  * Bell (1975) Bell T. H. J., 1975, J. Geophys. Res., 80, 320
  * Canup et al. (2021) Canup R. M., et al., 2021, arXiv e-prints, p. arXiv:2103.02045
  * Chambers (1999) Chambers J. E., 1999, MNRAS, 304, 793
  * Chiang & Laughlin (2013) Chiang E., Laughlin G., 2013, MNRAS, 431, 3444
  * Cilibrasi et al. (2021) Cilibrasi M., Szulágyi J., Grimm S. L., Mayer L., 2021, MNRAS, 504, 5455
  * Correia & Laskar (2001) Correia A. C. M., Laskar J., 2001, Nature, 411, 767
  * Counselman (1973) Counselman Charles C. I., 1973, ApJ, 180, 307
  * Cunha et al. (2015) Cunha D., Correia A. C. M., Laskar J., 2015, International Journal of Astrobiology, 14, 233
  * Dobrovolskis & Ingersoll (1980) Dobrovolskis A. R., Ingersoll A. P., 1980, Icarus, 41, 1
  * Domingos et al. (2006) Domingos R. C., Winter O. C., Yokoyama T., 2006, MNRAS, 373, 1227
  * Donnison (2010) Donnison J. R., 2010, MNRAS, 406, 1918
  * Efroimsky (2012) Efroimsky M., 2012, ApJ, 746, 150
  * Egbert & Ray (2000) Egbert G. D., Ray R. D., 2000, Nature, 405, 775
  * Gold & Soter (1969) Gold T., Soter S., 1969, Icarus, 11, 356
  * Hamers & Portegies Zwart (2018) Hamers A. S., Portegies Zwart S. F., 2018, ApJ, 869, L27
  * Hansen (2019) Hansen B. M. S., 2019, Science Advances, 5, eaaw8665
  * Hansen & Murray (2012) Hansen B. M. S., Murray N., 2012, ApJ, 751, 158
  * Hansen & Murray (2013) Hansen B. M. S., Murray N., 2013, ApJ, 775, 53
  * Heller et al. (2014) Heller R., et al., 2014, Astrobiology, 14, 798
  * Hesselbrock & Minton (2017) Hesselbrock A. J., Minton D. A., 2017, Nature Geoscience, 10, 266
  * Hippke (2015) Hippke M., 2015, ApJ, 806, 51
  * Ida & Lin (2010) Ida S., Lin D. N. C., 2010, ApJ, 719, 810
  * Ingersoll & Dobrovolskis (1978) Ingersoll A. P., Dobrovolskis A. R., 1978, Nature, 275, 37
  * Izidoro et al. (2017) Izidoro A., Ogihara M., Raymond S. N., Morbidelli A., Pierens A., Bitsch B., Cossou C., Hersant F., 2017, MNRAS, 470, 1750
  * Jeffreys (1921) Jeffreys H., 1921, Philosophical Transactions of the Royal Society of London Series A, 221, 239
  * Kipping (2009) Kipping D. M., 2009, MNRAS, 392, 181
  * Kopparapu et al. (2013) Kopparapu R. K., et al., 2013, ApJ, 765, 131
  * Kral et al. (2017) Kral Q., et al., 2017, The Astronomical Review, 13, 69
  * Lainey et al. (2007) Lainey V., Dehant V., Pätzold M., 2007, A&A, 465, 1075
  * Laskar et al. (1993) Laskar J., Joutel F., Robutel P., 1993, Nature, 361, 615
  * Leconte et al. (2015) Leconte J., Wu H., Menou K., Murray N., 2015, Science, 347, 632
  * Lewis (2013) Lewis K. M., 2013, MNRAS, 430, 1473
  * Lissauer et al. (2012) Lissauer J. J., Barnes J. W., Chambers J. E., 2012, Icarus, 217, 77
  * Lisse et al. (2012) Lisse C. M., et al., 2012, ApJ, 747, 93
  * Martinez et al. (2019) Martinez C. F., Cunha K., Ghezzi L., Smith V. V., 2019, ApJ, 875, 29
  * Melis (2016) Melis C., 2016, in Kastner J. H., Stelzer B., Metchev S. A., eds, Vol. 314, Young Stars & Planets Near the Sun. pp 241–246 (arXiv:1509.05726), doi:10.1017/S1743921315006638
  * Melis et al. (2021) Melis C., Olofsson J., Song I., Sarkis P., Weinberger A. J., Kennedy G., Krumpe M., 2021, ApJ, 923, 90
  * Moór et al. (2021) Moór A., et al., 2021, ApJ, 910, 27
  * Ochiai et al. (2014) Ochiai H., Nagasawa M., Ida S., 2014, ApJ, 790, 92
  * Peale (1999) Peale S. J., 1999, ARA&A, 37, 533
  * Piro (2018) Piro A. L., 2018, AJ, 156, 54
  * Rhee et al. (2008) Rhee J. H., Song I., Zuckerman B., 2008, ApJ, 675, 777
  * Sasaki et al. (2012) Sasaki T., Barnes J. W., O’Brien D. P., 2012, ApJ, 754, 51
  * Simon et al. (2007) Simon A., Szatmáry K., Szabó G. M., 2007, A&A, 470, 727
  * Spalding et al. (2016) Spalding C., Batygin K., Adams F. C., 2016, ApJ, 817, 18
  * Stevenson (1987) Stevenson D. J., 1987, Annual Review of Earth and Planetary Sciences, 15, 271
  * Sucerquia et al. (2019) Sucerquia M., Alvarado-Montes J. A., Zuluaga J. I., Cuello N., Giuppone C., 2019, MNRAS, 489, 2313
  * Swain et al. (2019) Swain M. R., Estrela R., Sotin C., Roudier G. M., Zellem R. T., 2019, ApJ, 881, 117
  * Teachey & Kipping (2018) Teachey A., Kipping D. M., 2018, Science Advances, 4, eaav1784
  * Tokadjian & Piro (2020) Tokadjian A., Piro A. L., 2020, AJ, 160, 194
  * Weinberger et al. (2011) Weinberger A. J., Becklin E. E., Song I., Zuckerman B., 2011, ApJ, 726, 72
  * Weiss & Marcy (2014) Weiss L. M., Marcy G. W., 2014, ApJ, 783, L6
  * Wyatt (2008) Wyatt M. C., 2008, ARA&A, 46, 339
  * Wyatt et al. (2007) Wyatt M. C., Smith R., Greaves J. S., Beichman C. A., Bryden G., Lisse C. M., 2007, ApJ, 658, 569
  * Youdin (2011) Youdin A. N., 2011, ApJ, 742, 38
  * Zuckerman & Song (2012) Zuckerman B., Song I., 2012, ApJ, 758, 77

##  Appendix A Atmospheric Tide Model

The thermal inertia of a planetary atmosphere means that the mass
redistribution of the atmosphere, in response to the heating from the central
star, lags the true diurnal cycle. This means that the torque, applied by the
atmosphere to the planetary interior through boundary layer frictional forces,
also lags and contributes an asynchronous contribution to the total torque
(Gold & Soter, 1969; Correia & Laskar, 2001; Leconte et al., 2015; Cunha et
al., 2015; Auclair-Desrotour et al., 2017) . Although such effects are minor
for Earth-like planets, they become progressively more important for planets
with high levels of irradiation (and are thought to be the determining factor
in the retrograde spin of Venus).

To model atmospheric tides over a broad range of irradiation, we adopt the
parameterised model from Leconte et al. (2015). They performed a series of
global circulation models to quantify the strength of the atmospheric torque
as a function of surface pressure $p_{s}$ (a proxy for atmospheric mass),
composition, and incident stellar flux $F_{in}$. To these simulations they fit
an analytic model which we adopt as function $b_{a}(n_{p},\Omega_{p})$ in
equation (4), namely

$b_{a}(n_{p},\Omega_{p})=\frac{(\Omega_{p}-n_{p})/\omega_{0}}{1+\left[(\Omega_{p}-n_{p})/\omega_{0}\right]^{2}}$
(31)

which is parameterised by a characteristic frequency response $\omega_{0}$
that is a function of $p_{s}$, $F_{in}$ and composition. The strength of the
atmospheric response is also quantified by an amplitude function $q_{a}$ (see
equation (4)), which is also a function of $p_{s}$, $F_{in}$ and composition.

Leconte et al. (2015) provide a table of $q_{a}$ and $\omega_{0}$ for a series
of different simulations. Our focus is on planets similar to Earth but closer
to the Sun, so we adopt the values for $p_{s}$=1 bar, Earth-like composition,
and $F_{in}=1366W/m^{2}$. This yields $q_{0}=1180$ Pa and $2\pi/\omega_{0}=32$
days. The atmospheric torque also operates in the opposite direction to the
direct gravitational solid-body torque, opening the possibility of an
asynchronous equilibrium, as discussed by Leconte et al. (2015). Figure 10
demonstrates how this arises.

Figure 10:  The long-dashed line shows the solid body torque for a
$1M_{\oplus}$ planet, orbiting in a circular orbit at 0.7 AU from a 1
$M_{\odot}$ star. The red curve is for $Q=30$ and the blue curve is for
$Q=100$. The short-dashed line is the atmospheric torque operating on the same
planet, based on the default model described in the text. The solid line
represents the net torque operating on the planet, with the blue and red
curves corresponding to the sum of the corresponding solid body torque with
the same atmospheric torque. Where these curves cross zero (the dotted line),
we have a potential spin equilibrium. The filled solid point indicates the
first equilibrium the system would encounter as the planet spins down due to
the outspiral of a moon.

The short dashed curve represents the atmospheric torque model above. The long
dashed line indicates the equivalent solid body torque, using
$\bar{b}_{g}=arctan(26.8(y-1))$ – the numerical value of equation (16) in this
case. The red curve indicates the case of $Q=30$ and blue indicates the case
of $Q=100$. The solid curves indicate the corresponding net torques in the two
cases. The parameter choices for Figure 10 are chosen to represent the same
cases as in Figure 2.

The solid black curve in Figure 2 corresponds to the red dashed curve in
Figure 10 – the case of no atmospheric torque. The addition of the atmospheric
torque to this (resulting in the red solid line) does not change the
qualitative nature of the equilibria, but weakens the spin down torque and
allows the moon to spiral out further before the planetary spin reverses (as
seen in Figure 2). If the solid body torque is weakened (the blue case), then
new equilibria appear.

These asynchronous equilibria can have a significant effect on the moon
evolution. Shortly after the formation of the moon, the planet is rotating
rapidly and will lie off to the right of Figure 10. As the moon spirals out,
the planet spins down and moves to the left in the diagram. In the absence of
atmospheric tides, this will continue until the planetary spin is synchronous
with the orbital frequency, which lies at $\Omega_{p}/n_{p}=1$ in this
diagram. However, we see that the system shown in Figure 10 will encounter an
asynchronous equilibrium first (shown as a filled circle). As the moon
continues to spiral outwards (acting to move the system to the left of the
solid point), the atmospheric tides act to spin the planet up and move it back
to the stable equilibrium. Thus, the atmospheric tides serve to pump angular
momentum into the system, rejuvenating the outward motion of the moon, as is
seen in Figure 3.