Tidal Dissipation in Satellites Prevents Hill Sphere Escape

The transit method is a promising means to detect exomoons, but few candidates have been identified. For planets close to their stars, the dynamical interaction between a satellite's orbit and the star must be important in their evolution. Satellites beyond synchronous orbit spiral out due to the tide raised on their planet, and it has been assumed that they would likely escape the Hill sphere. Here we follow the evolution with a three-body code that accounts for tidal dissipation within both the planet and the satellite. We show that tidal dissipation in satellites often keeps them bound to their planet, making exomoons more observable than previously thought. The probability of escape depends on the ratio of tidal quality factors of the planet and satellite; when this ratio exceeds 0.5, escape is usually avoided. Instead, the satellite moves to an equilibrium in which the spin angular momentum of the planet is not transferred into the orbit of the satellite, but is transferred into the orbit of the planet itself. While the planet continues spinning faster than the satellite orbits, the satellite maintains a semi-major axis of approximately 0.41 Hill radii. These states are accompanied with modest satellite eccentricity near 0.1 and are found to be stable over long timescales.


INTRODUCTION
Satellites of planets in the Solar System are common, and the giant planets all host extensive systems of satellites.Their existence, composition, and orbital states have all contributed to theories of the origin of the planets themselves (Peale 1999).With the discovery of more than 5000 exoplanets, it is natural to ask whether those bodies have so-called exomoons.
The most straightforward way to detect exomoons to date involves both the planet and the exomoon transiting in front of the star (Sartoretti & Schneider 1999).Even if the moon itself is not seen transiting, the timing variations it induces on the planet may be seen (Kipping 2009).The first precise transit curves of the first transiting planet were thus searched for exomoons (Brown et al. 2001).
With a lack of any definitive detections so far, theory becomes more important, both to reexamine what we really expect for the frequency of exomoons, and perhaps to provide some guidance on where the search could be most fruitful going forward.
The most immediate reason we expect exomoons is that moons are nearly ubiquitous in the Solar System, and there are theoretical reasons to expect moons around exoplanets may be even more detectable than straight analogues to the Solar System satellites (Hamers & ★ E-mail: akisare@uchicago.eduPortegies Zwart 2018; Moraes & Vieira Neto 2020).Since satellite formation has numerous pathways, it is hard to develop confident predictions.
Many have done tidal calculations for moons in the Solar System and hypothetical moons around exoplanets.Once a moon has formed, Barnes & O'Brien (2002) calculated the tidal evolution, showing that planets which are tidally locked with their stars will cause their satellites to decay, but rapidly-spinning planets could push their satellites out of the Hill sphere.In fact, Counselman (1973) had already speculated that this type of evolution may have removed primordial moons from Mercury and Venus.Both Sasaki et al. (2012) and Piro (2018) also took into account the slowing of a planet's rotation by the tide raised by a moon.Cassidy et al. (2009) looked at the problem even closer, to see how the tidal response of the moon might change the dynamics.We pick up the theoretical story there.
The Hill sphere is the region around a planet where the gravity of a planet dominates the motion of test particles (and real moons, to a good approximation).Its size is  H =  p ( p /3 * ) 1/3 , where  p is the planet's semi-major axis and  p and  * are the mass of the planet and star, respectively.At this distance from the planet lie the Lagrange  1 and  2 equilibrium points.These are unstable equilibria though, and a natural satellite placed at them would eventually fall into either astrocentric or planetocentric orbits.This radius limits the motion of a satellite, but stable motion does not usually extend to that distance.Rather, prograde orbits have 0.5 H as a typical maximum distance of stability (Shen & Tremaine 2008), and orbits initialised further are able to escape beyond  1 and  2 .The evection resonance, in which the apsidal period of the satellite equals the planet's orbit around the star, is thought to be the mechanism that finally disrupts a satellite orbit (Grishin et al. 2017).In a rotating reference frame, the fictitious force known as the Coriolis force enhances the binding of retrograde satellites and destabilizes prograde ones (Innanen 1980).For the purposes of the n-body code, the center of  p is treated as the origin.Relative sizes and distances are not to scale.
Though contained within a zero-velocity curve, stable orbits can quickly transition between circular and eccentric (Hamilton & Burns 1991).These works did not include tides, so we suspected these orbits with fluctuating eccentricities may evolve due to dissipation.We set out to find the outcome.
In the following sections, we describe our numerical methods ( § 2), display the results of a suite of n-body simulations ( § 3), derive analytically the stability conditions ( § 4), and explain how the results affect the big picture of exomoon evolution and even detection ( § 5).

METHODS
Our method involves direct N-body simulations, using the "direct code" of Mardling & Lin (2002), described in their Section 2. The direct code uses a Bulirsch-Stoer algorithm with four integration steps per circular orbit of the innermost body.Zollinger et al. (2017) noted that Mardling & Lin's equation (7) contains a minor error; the term  3  34 within the first set of brackets should be changed to  4  34 .However, the equations of motion in the code provided to us by Mardling do not contain this error.The hierarchical coordinate system is illustrated by Figure 1.

Variables
Table 1 defines many variables used throughout this paper.Quantities with a subscript "p" reference the planet, starred quantities reference the stellar body, and quantities without either generally refer to the satellite unless otherwise clarified.We describe the orbital elements of the inner pair as those of the satellite around the planet, and the orbital elements of the outer orbit as those of the planet around the star.However, the coordinates system used by the code is in the frame centred on the planet.When plotting results, rather than using Keplerian elements, which were found to be very poor at qualitatively describing the orbits this far out in the Hill sphere, we use  = ( max +  min )/2 and  = ( max −  min )/(2). min and  max were determined by splitting the positional data into 625 intervals with 32 data points per interval, calculating the minimum and maximum radial position in each interval and then interpolating between these values using a piece-wise cubic Hermite interpolating polynomial.Additionally, the satellite is defined as having escaped once  >  H .
We limit our explorations to prograde, coplanar, Moon-like satellites, with zero-obliquity, Earth-like planets.Previous work has used an analytical approach (Cassidy et al. 2009) to determine how a star perturbs a tidally-evolving satellite.Our numerical approach is able to see new effects, even though the model is kept simple.The planet was placed in a low-eccentricity orbit, also for simplicity.The parameters we varied were the -factors of the planet and satellite.The -factor controls the strength of tidal dissipation in a body and depends on its interior structure and composition (Goldreich & Soter 1966;Greenberg 2009).The influence that the -factor has on the stability of satellites close to the Hill sphere is therefore important in predicting the properties of any exomoon systems that may be discovered.

Normalization
The semi-major axis is normalized such that  = / H .Note that the magnitude of  H hardly changes in our simulations.We further define a normalized time  = / dyn where  dyn is the dynamical timescale: The value of  dyn is the theoretical timescale for the satellite to migrate outwards to its present orbit due purely to tidal dissipation in the planet.† Finally, the spin frequencies Ω and Ω p are normalized by the satellite's mean motion such that  = Ω/ and  p = Ω p /.For  = 1, the satellite is rotating synchronously.For  p = 1, the planet is rotating synchronously with the satellite.

Momentum Conservation
The total angular momentum of the system is given by: where  p =  * ( p + )/( * +  p + ) and  =  p /( p + ).
The four terms are, in order, the orbital angular momentum of the planet around the star  p , the orbital angular momentum of the satellite around the planet , the spin angular momentum of the planet  p , and the spin angular momentum of the satellite .The total angular momentum  should be conserved.To test how well the code conserves angular momentum, we ran a simulation with  * = 0 for simplicity, and  p =  = 100.Over a period of 5 Myr, the fractional error in the total angular momentum | |/ grew approximately linearly to ∼0.1%.The proportion of orbital angular momentum / grew from 50% to 60% at the expense of  p /, a value much greater than the numerical error.The only component with a contribution smaller than the error was the spin angular momentum of the satellite / ∼ 0.001%.However, as  is entirely coupled to  and  p , its value is not compromised by the numerical error.

NUMERICAL RESULTS
In Figure 2(a) we show an individual system evolution with  p =  = 100, which we refer to as Simulation 1.The satellite approaches the Hill sphere until  ≃ 0.05, after which the satellite is saved from escaping from the planet and falls back inwards.After  ≃ 0.3, the satellite appears to approach an equilibrium, with  consistently fluctuating near 0.42, moderate eccentricity and a stable spin frequency.
The spin of the planet, however, continues decreasing for the entire duration of the simulation.
The angular momentum of the planetary spin is thus decreasing, yet the orbit of the satellite has stopped evolving; the angular momentum must be transferring into the planetary orbit, i.e. the orbit of the planet around the star is growing.After equilibrium was achieved at  ≃ 0.3,  p / grew linearly by approximately 0.2 ppm from this point until the end of the simulation, while  p / decreased linearly by nearly the exact same amount, confirming that the planet's spin angular momentum is indeed being transferred into its own orbit.showing that both  p and  matter for whether the satellite stays bound to the planet.Each pixel represents a combination of  and  p .All the white pixels had the satellite staying bound to the planet for the full duration of the simulation, while black pixels correspond to simulations where the satellite escaped.The hatched region indicates -factors that were not tested.The dotted line demarcates the four separate sets of simulation data that were stitched together.The upper left quadrant has the least variation in outcome, so a lower resolution of 0.02 px/ was used there, while the other three quadrants have a resolution of (43/450) px/ ≃ 0.10 px/.
We now experiment with two variables: the tidal quality factor of the planet  p , and that of the satellite .
In Figure 2(b), we let  p ≃ 403.48 and  = 100, which we call Simulation 2. In Simulation 1, which has a lower  p / than Simulation 2, the satellite very quickly approaches the Hill sphere, nearly escaping before having a crisis at  ≃ 0.05, saving itself, falling back towards the planet, circularizing and reaching an equilibrium.In Simulation 2, the satellite never dramatically falls inwards before rebounding, and starts its evolution in a state similar to that of Simulation 1 at  ≃ 0.1.
Running 5820 simulations, each with a unique combination of  p and , we find that escape occurs with complicated boundary (Figure 2(c)).We define  as the vector (  p ) which has a slope  p / and a norm  = ||.From this set of simulation data, we organised each data point into 24 bins of slope ranging from 0.1 to 0.86.For each of these bins, the number of simulations that resulted in an escape of the satellite was divided by the total number of simulations in that bin, giving an escape probability  esc in each bin.Points with  > 1000 were excluded, so as to help balance the number of data points in each bin.This is plotted as the "all " histogram in Figure 3(a).Two more bins were then created, a "high " bin where 743 <  < 1000 and a "low " bin where  ≤ 743 .743 was chosen as the boundary as this value equalises the number of data points in the high  and low  bins.The two additional histograms corresponding to these two bins are also plotted in Figure 3(a).
There is very little discrepancy between  esc in the high- region and the low- region.The sharp drop in probability in the low- region around  p / ≃ 0.1 is a consequence of lower slopes being geometrically impossible.The consistency in the escape probability between the high and low regions suggests that escape is independent of , but highly dependent on  p /. Greater relative dissipation in the planet (low  p /) yields higher probabilities of escape, while greater relative dissipation in the satellite (high  p /) yields lower probabilities.
In Figure 3(b), we show that there exists a modest relationship between  esc and  p for a given value of . = 1000 was chosen as it maximizes the number of data points, but a similar relationship exists for all values of .Although there is much noise, escape generally occurs slightly quicker when  p is lower.This is interesting, as  is already a time quantity normalized by the dynamical timescale  dyn which from (1) is linear in  p .Therefore, there is likely some weaker, higher-order dependence between  dyn and  p which we have not accounted for.An opposite relation exists between  and  esc , shown in Figure 3(c).Stronger dissipation in the satellite (low ) increases escape time, particularly at the boundary between escaping and remaining bound.
The statistical results at the end of all simulations for which the satellite did not escape is shown in Figure 4.The value of  in this plot is the normalized time at the end of the simulation.Higher  represents systems that are more dynamically evolved.This strongly suggests that, in Figure 4(a), there is a movement of orbits from the light-colored region to the dark-colored region centred around  ≃ 0.1.From Figure 4(b), we find that highly evolved orbits can stabilise into this excited eccentricity if  p / > 1. Essentially, if the satellite does not escape, systems with strong dissipation in the satellite relative to the planet will stabilise the satellite into a slightly eccentric orbit.
The final state of the system in each of these simulations is with the planet still out of equilibrium; the planet is spinning down due to tides raised by the satellite on it.The Mardling code does not raise tides due to the star, but they would also act to spin down the planet.Future work will look at the even longer-term evolution of the satellite, once the planet has spun down to the orbital frequency of the satellite.We expect that from then on, the satellite cannot escape into an astrocentric orbit, because the total angular momentum of the planet-satellite system is too low to cause the satellite to reach the Hill sphere.Instead, there would be spin-orbit locking of both the planet and the satellite.Angular momentum would continue to drain at an even slower rate, now due mainly to the tidal bulge the star raises on the planet.This would couple to the satellite's orbit, now the dominant reservoir of angular momentum, which would very slowly decay.

ANALYTIC RESULTS
Consider the the gravitational "main problem" of lunar theory (Brouwer & Clemence 1961), i.e. just three point masses and no rotational or tidal distortion, and no tidal dissipation.If the dominant mass and the barycentre of the two smaller masses follows a circular orbit, then periodic solutions of these two two smaller masses around each other exist (Henon 1969).The two families of bound, periodic, prograde orbits are called  and  ′ .Their  and  are shown in Figure 5. Family  has low eccentricity, and its nearly-elliptical orbit is a centred ellipse, with long axis forming a right angle with the direction to the star.That is, it makes two radial oscillations for every azimuthal orbit.There is a maximum distance,  ≃ 0.41, beyond which no similarly shaped orbits are stable.Instead, family  ′ curves continue to higher  at larger eccentricity, ‡ and the figure is no longer centred on the planet but has its long axis pointing either towards or away from the star.That family also ends at  ≃ 0.70.
Our computations actually achieve steady state along family  ′ (Figure 4), with an average value of  much higher than any of the family  orbits, but also not extending to as large  as the far end of family  ′ .Prior work on tidal dissipation in planets and satellites have found stable equilibria, limit cycles, and chaotic phases all to be possible states (Meyer & Wisdom 2008).We tentatively identify family  ′ as a possible stable equilibrium actually achievable with only small modification due to rotational and tidal bulges and dissipation.The presence of some unsteady behavior is from a combination of limit cycles (where no escape occurs) and chaotic phases (allowing escape).
We focus attention on how the system achieves equilibrium on family  ′ , and then determine an empirical functional form for the probability of escape.Cassidy et al. (2009) analyzed the problem of a satellite on a variational orbit (family ) and found that tidal dissipation in the satellite could considerably heat, even melt, the satellite.Meanwhile, this energy was coming from the orbit of the ★ Synchronised 1:1 with the mean motion of the satellite.
† See § 3 for the particular -factors used in Simulation 1 and 2. A total of 5820 unique combinations of  and  p were simulated.(Henon 1969).The shaded region corresponds to the boundaries of the plots in Figure 4. star around the planet-moon system, not from the orbit of the satellite around the planet.Since the former energy reservoir is enormous, the satellite could be heated on timescales of several billion years with only negligible orbital evolution of the system.However, since the energy does not come at the expense of the orbit of the satellite, and the tides on the planet continue pumping orbital energy into the satellite, this picture does not explain how the satellite can reach a steady state in our computation.

Equilibrium State
Important to the dynamics of satellites is the evection resonance, which occurs when the pericentre precession rate of the satellite  is commensurable with the mean motion of the planet  p .This precession is produced by stellar perturbations and by the quadrupole moment of the planetary gravitational potential.For our simulations, stellar perturbations dominate over the planetary quadrupole.However, for satellites at a substantial fraction of the Hill radius, simple methods of orbit averaging are insufficient at analytically approximating the pericentre precession rate.Henon (1969) demonstrated that the evection resonance occurs at  ≃ 0.41.Our own numerical simulations confirm this common result, with the resonance argument  =  −  p librating near zero, where  p is the orbital longitude of the planet.The dynamical equations for  and  due to the star include only factors of sin  (Yoder 1982), so  ≃ 0 implies that the stellar perturbations of these orbital elements are negligible.We can therefore direct our attention to the relevance of tides in the equilibrium state.
Our numerical simulations demonstrate that the satellite remains stable with an excited eccentricity oscillating about some equilibrium value.The stability of this eccentricity is explained by the mutual and counteracting effects of the satellite and planetary tides.We first define the following: From Yoder & Peale (1981), the equations of motion for  and  due to the mutual tides of both bodies are then The equation for the tidal change in  holds because tidal dissipation conserves the orbital angular momentum of the satellite  =  2 (1 −  2 ) 1/2 for a synchronously rotating satellite in the evection resonance, and all satellites in equilibrium at the end of our simulations have reached the evection resonance and are synchronously rotating.The evection resonance is required for the conservation of  since under this condition, the spin angular momentum of the planet is completely transferred directly into its own orbit and not into the orbit of the satellite.Letting / = / = 0 and solving for eccentricity yields To validate the stability of this equilibrium, we identify the following Lyapunov function: This choice of a Lyapunov function does not include  as an argument, since  * is independent of  and is an equilibrium value for both / and /.The derivative of the Lyapunov function is Since / ≤ 0, it follows that  * is a stable equilibrium.For the most dynamically evolved systems, which are those systems with  > 2 ≃ 0.4, the value of  * ranges from 0.08 to 0.27.This corresponds closely with the actual eccentricities attained by these systems, which range from 0.05 to 0.36.Our analysis demonstrates that evolved systems do reach an equilibrium eccentricity which acts to halt expansion of the satellite's orbit.The mutual tidal effects nullifying each other, the satellite finds itself among the  ′ orbits, which are known to be stable.Cumulatively, the equilibrium state requires an eccentricity  =  * , synchronous rotation, and the evection resonance, all of which is in accordance with our numerical results.

Probability of Escape
To determine an expression for the probability of escape, we take the magnitude of the eccentricity excitation due to capture into the evection resonance, and compare this to the maximum stable eccentricity for  ′ orbits.Yoder & Peale (1981) found that capture into a resonance excites the eccentricity by an amount Δ = 6 1/2  c , where  c is a critical eccentricity required for transitioning between circulating and librating modes in .For us,  c =  * .Meanwhile, the maximum stable eccentricity for  ′ orbits is  ≃ 0.52.Assuming  esc (Δ = 0.52) = 1, we propose an ad-hoc function which follows our numerical  esc results: In Figure 6, this function is shown as comparable to the numerical probability of escape.However, future work should test various combinations of Love numbers, -factors, masses, and radii to confirm that equation (8) holds or identify additional dependence that it misses.

CONCLUSIONS
We have computed the tidal evolution of satellite systems accounting for both tidal dissipation in the planet and in the satellite.If tidal dissipation in an exomoon is comparable or greater than that in their host planet, it can hang out at a substantial fraction of the Hill radius as its host planet's spin rate is decreasing.In Figure 2 we see that the satellite's orbit comes to an equilibrium in , , and , in which the planet's spin angular momentum is transferred into its own orbit.

Application to transiting planets
A typical planet discovered by Kepler has  ∼ 1.5 ⊕ and orbits at  ∼ 0.1 au.Depending on whether it is a sub-Neptune or a super-Earth, its  can be as high as 10 4 .Satellites would tend to be smaller, and perhaps rocky, so their  would tend to be lower.Therefore, this mechanism would likely keep the satellite from escaping, and it would keep them at a large separation from the planet ( ≃ 0.41 H ) for hundreds of millions of years.Currently, the Sun raises about half as big a tide on the Earth in comparison to the Moon, and the degree of dissipation is a function of how fast this bulge moves through the oceans, which is comparable for each because Earth's spin dominates.Therefore, if the Moon were about 2 1/3 ≃ 1.26 times further away (i.e.0.315 H ), then the tidal dissipation from the Sun and the Moon would be comparable.At the extreme end of the  family for instance, the Sun's dissipation would be stronger.Our moon would be challenging to detect with even the best photometry (Kipping et al. 2009), so it is more likely to imagine a bigger, more massive satellite that would in fact dominate the tidal despinning of the planet, relative to the star.

Limitation and future work
Due to time limitations placed on the usage of computational resources, our simulations did not follow the evolution beyond the quasi-equilibrium in which the moon is synchronous and no longer moving out, but the planet is still spinning relatively rapidly.If the spin rate of the planet keeps decreasing at the same rate as in Simulation 1, then it will become synchronous with the satellite at  ∼ 2. The dominant angular momentum transfer would be due to the direct tide the star raises on the planet, which we have not modelled.That continual leakage of angular momentum from the planet-satellite pair would cause the satellite to spiral towards the planet (Hansen 2023).Eventually it could directly impact the planet, if the planet is much less dense, or cause the satellite to be tidally shredded at its Roche limit.After such a shredding episode, the satellite could reconstitute itself as a less massive satellite (Hesselbrock & Minton 2017).At some point the satellite mass would no longer dominate the torque on the planet's spin, and the star would be able to spin-synchronise the planet.
Future work could synthesise a population of stars, planets, and moons, at a variety of starting conditions consistent with planet formation, and determine what star-planet separations and planet types would most likely have moons in the tidally-evolved state we have found:  ≃ 0.1 at  ≃ 0.41 H . Prior work, for example by Zollinger et al. (2017), Tokadjian &Piro (2020), andHansen (2023), has assumed satellites would spiral outward and escape once they reached a semi-major axis ≃0.5 H , which is what tidal dissipation in satellites can prevent.For some combinations of parameters, we anticipate the lifetimes of moons which would have otherwise escaped would be prolonged, increasing their chance of being detected.Tidal dissipation may cause such significant heating in these exomoons, that they would shine brightly in the infrared (Peters-Limbach & Turner 2013) and produce observable, volcanic signatures of sodium and potassium (Oza et al. 2019), making such exomoons even more detectable.We are encouraged that the exomoon-hunting era is not over with the results of Kepler; we anticipate moons surviving around planets with periods shorter than the Earth.Hence, TESS (Ricker et al. 2015), JWST (Gardner et al. 2006), PLATO (Rauer et al. 2022), and others have a role to play in the hunt for exomoons.

Figure 1 .
Figure1.Model components.Tidally distorted and shifted planet and satellite, with a point-mass star referred to the centre of mass of the inner system.For the purposes of the n-body code, the center of  p is treated as the origin.Relative sizes and distances are not to scale.

Figure 2 .
Figure 2. Summary of numerical results: (a) Simulation 1, with  p =  100; (b) Simulation 2, with  p ≃ 403.48 and  = 100; (c) set of simulationsshowing that both  p and  matter for whether the satellite stays bound to the planet.Each pixel represents a combination of  and  p .All the white pixels had the satellite staying bound to the planet for the full duration of the simulation, while black pixels correspond to simulations where the satellite escaped.The hatched region indicates -factors that were not tested.The dotted line demarcates the four separate sets of simulation data that were stitched together.The upper left quadrant has the least variation in outcome, so a lower resolution of 0.02 px/ was used there, while the other three quadrants have a resolution of (43/450) px/ ≃ 0.10 px/.

Figure 3 .Figure 4 .
Figure 3. Statistics of escape: (a) histogram of the probability of escape, versus the -ratio  p /, with "high " being the region where 743 <  < 1000 and "low " being the region where  ≤ 743; (b) relationship between escape time  esc and  p for  = 1000; (c) effect of  on escape time for  p = 100.

Figure 5 .
Figure5.Periodic orbits within Hill's problem, modeling the point-mass three-body problem in which one of the three masses (the star) is much larger than the other two(Henon 1969).The shaded region corresponds to the boundaries of the plots in Figure4.

Figure 6 .
Figure 6.Dotted histogram shows the probability of escape found numerically, while the solid line is the semi-analytic functional form from equation (8)

Table 1 .
Variable Definitions and their Initial Values