© The Authors 2023

1 Introduction

A common convention regarding the rotation of close-in exoplanets is that their tidal de-spinning is fully defined by the tides generated on these planets by their host stars. It is for this reason that close-in gas giants are assumed to be synchronised or pseudo-synchronised. We pose the question of whether a moon orbiting a planet can change the outcome of the planet’s tidal de-spinning and make the planet end up in a non-synchronous (perhaps, even retrograde) spin state. Specifically, we explore the possibility of a planet being synchronised not by the star but by a moon. If this scenario is plausible for some inner planets, these planets’ rotation will be much faster than the speed at which they orbit the star. These planets then avoid spiraling into their host stars, and maintain finite orbital eccentricities¹. An internally synchronised planet–moon pair continues to shrink because of the tides caused by the star, but the rate of this process is orders of magnitude slower than the orbital decay of a moonless planet. Indeed, in an internally synchronised or pseudo-synchronised planet-moon system, the tidal torque acting on the planet from the star is counterbalanced by the tidal torque acting on the planet from the moon. The time-averaged total torque on the planet in a stable equilibrium is nil, but both these components continue to dissipate energy. This energy is taken from the orbital kinetic energy of the planet and the moon, respectively. The supply of orbital energy is orders of magnitude greater than the rotational energy of the planet, which makes the decay process very slow. A decisive test of these models will be possible once more accurate photometric measurements of the secondary eclipses of hot Jupiters become available (i.e. the transits of the planets behind the disks of the stars in the upper conjunction), where a non-synchronous rotation is expected to produce a measurable phase shift with the peak brightness of the out-of-eclipse light curve. A statistically significant shift towards the evening terminator was detected for the WASP-12b planet by Owens et al. (2021). On a fast-rotating planet, the hot spot generated by the stellar irradiation is not stationary but rather leads the substellar meridian.

The dynamical evolution of exomoons in star-planet systems has been investigated in a number of publications (e.g., Sasaki et al. 2012; Barnes & O’Brien 2002; Guimarães & Valio 2018). The commonly shared conclusion is that the chances of exo-moons surviving for an extended duration of time in known planetary systems are low. Considering three basic scenarios for a moon in such systems, with only one of these scenarios allowing the moon to survive for at least 1 Gyr, Dobos et al. (2021) conclude that the survival rate is close to nil for planets with orbital periods of 10 days or less, and that it gradually increases to 70% for periods of 300 days. The median orbital period of all detected and tentative exoplanets is about 11.8 days, and only 10% of them have periods above 442 days, so we should expect to find very few exoplanets with moons.

In the current work, our goal is to demonstrate that there are additional pathways for exomoons to survive, even around fairly close-in planets. Furthermore, exomoons can help their parent planets survive for longer times in the close vicinity of stars. As an example of such previously ignored possibilities, moons that are initially on retrograde orbits with respect to the spin of a slowly rotating planet will generate a tidal torque on the planet that will be opposed to the torque caused by the star. Possible in principle, this particular scenario is unlikely to play a large role in real systems, because the capture of external moons onto retrograde orbits is relatively rare². Starting with solid principles of tidal dynamics, we aim to show the existence of mechanisms plausible for a long-term coexistence of moons with their planets (most of the considered planets being dangerously close to their stars). Along with Earths and super-Earths, our study addresses Jupiters. It is commonly assumed that close-in low-viscosity giant planets are synchronised by their host stars. Although the tidal quality of Jupiters is likely to be as high as that of our own Jupiter (Efroimsky & Makarov 2022), the dissipation of orbital energy inside a synchronised planet is almost nullified because of the circularisation process bringing down the planet’s eccentricity. The fate of the planet is then defined by the rotation rate of the host star, which is often slower than the mean motion. Slowly but relentlessly, the tidal break will remove the orbital energy, and the planet will plunge into the star. We however show here that a sufficiently massive moon can arrest or even reverse this process for a significantly long time, until the moon’s orbital energy is depleted.

The nomenclature used in this paper is listed in Table 1.

Fig. 1 Reduced Hill radii of 580 known inner exoplanets versus their radii in units of Jupiter’s radius. The diagonal line indicates the equality of the two radii.

2 The Hill sphere and the reduced Hill sphere

The Hill sphere is a region where the motion of moons is defined predominantly by the gravity of the planet, and less by that of the star. Analysis by Hamilton & Burns (1992) set its radius at

$$r_{\text{H}}=a_{\text{p}}\left(1-e_{\text{p}}\right){\left(\frac{M_{\text{p}}}{3M_{*}}\right)}^{1/3},$$(1)

where a_p and e_p are the planet’s semi-major axis and eccentricity, while M_p and M_* are the masses of the planet and the star, correspondingly.

However, numerical analysis by Astakhov et al. (2003) and Domingos et al. (2006) indicated that the orbits too closely approaching this radius become unstable in the long term, while those confined to a smaller domain (which we term ‘the reduced Hill sphere’) remain stable. The radius of the reduced Hill sphere (‘the reduced Hill radius’) for a circular (e_m = 0) orbit of the moon is

$$r_{\text{H}}^{\prime }=B{a}_{\text{p}}\left(1-b{e}_{\text{p}}\right){\left(\frac{M_{\text{p}}}{3M_{*}}\right)}^{1/3}$$(2)

and depends on the direction of the orbit:

$$B=0.49,b=1.03,\text{\hspace{1em}for\hspace{0.17em}prograde-orbiting\hspace{0.17em}moons;}$$(3a)

$$B=0.93,b=1.08,\text{\hspace{1em}for\hspace{0.17em}retrograde-orbiting\hspace{0.17em}moons.}$$(3b)

Figure 1 shows the distribution of calculated reduced Hill radii for a sample of 580 inner exoplanets with the necessary estimated parameters entering Eq. (2) provided in the NASA Exoplanet Archive³ (assuming prograde orbits). These values are plotted against the estimated radii of the planets themselves. One should expect the Hill spheres to be outside of the planets’ radii, but this seems to not always be the case. The planets orbiting HAT-P-32, HAT-P-65, HATS-27, HATS-69, HIP 65 A, WASP-121, and WASP-19 appear to be in violation of this condition. These are bloated Jupiters with orbital periods between 0.79 and 4.64 days orbiting solar-type stars. However, caution should be exercised with the archival data for specific exoplan-ets, where generous upper limits are sometimes given for orbital eccentricity instead of true determinations. This concerns the planets HAT-P-65, HATS-27, and HATS-69. Still, the existing models suggest that these giant planets, being so close to the Hill radii, should be shedding their outer layers, probably forming a gaseous disk around the stars. The origin of such critical proximity of inflated hot Jupiters to their host stars is not known.

Equations (2)–(3) also yield an inequality to be utilised below:

$$a_{\text{m}}<r_{\text{H}}^{\prime }\approx B{a}_{\text{p}}{\left(\frac{M_{\text{p}}}{3M_{*}}\right)}^{1/3},$$(4)

a_m and a_p being the semi-major axes of the moon and the planet.

3 The Roche radius

The Roche radius (or the Roche limit) is the closest distance at which a moon can get to its planet without being pulled apart by tides. This parameter sets the lower limit on the survival niche for exomoons. For a spherical host, the Roche radius is given by

$$r_{\text{R}}=A{R}_{\text{m}}{\left(\frac{M_{\text{p}}}{M_{\text{m}}}\right)}^{1/3},$$(5)

where M_m and R_m are the mass and radius of the moon, while M_p is the mass of the planet.

For a strengthless moon, the dimensionless parameter would assume the value A ≈ 2.46 borrowed from a calculation by Chandrasekhar (1987) for an incompressible fluid body. The shear strength of realistic rubble piles enables them, however, to survive at smaller radii, with A ≈ 2.2 (Leinhardt et al. 2012)⁴.

If the moon was not captured by the planet but formed alongside, the parameter A assumes the value 2.2 appropriate for rubble piles, because freshly accreted layers of a moon in formation are rubble. So we agree to use the expression

$$r_{\text{R}}=2.2R_{\text{m}}{\left(\frac{M_{\text{p}}}{M_{\text{m}}}\right)}^{1/3}.$$(6)

4 The synchronous value of the semi-major axis

Defined by the synchronicity condition ${\dot{\theta }}_{\text{p}}=n_{\text{m}}$, the synchronous value $a_{\text{m}}^{\left(\text{s}\right)}$ of the moon’s semi-major axis is

$$a_{\text{m}}^{\left(\text{s}\right)}=\left(\frac{G\left(M_{\text{p}}+M_{\text{m}}\right)}{{\dot{\theta }}_p^2}\right).$$(7)

where G = 6.6743 × 10⁻¹¹ m³ kg⁻¹ s⁻² is Newton’s gravitational constant, while θ_p and ${\dot{\theta }}_{\text{p}}$ are the sidereal rotation angle and rotation rate of the planet. Commonly, $a_{\text{m}}^{\left(\text{s}\right)}$ is referred to as the synchronous radius.

Bear in mind that this radius is defined only for moons in orbits prograde relative to the planet’s spin. For a retrograde moon, the semi-diurnal tidal bulge it creates on the planet will always be pulling the moon back in its orbital motion, and thus will be working to bring it down⁵.

5 A planet synchronised by its moon

As we shall see below, a moon whose mass is about one percent of the planet’s mass can in some situations synchronise the planet. For this to happen, two inequalities must be obeyed.

First, a synchronised moon’s apocentre, $a_{\text{m}}^{\left(\text{s}\right)}\left(1+e_{\text{m}}\right)$, must end up below the reduced Hill radius. For small e_m, we write this simply as

$$a_{\text{m}}^{\left(\text{s}\right)}<r_{\text{H}}^{\prime }.$$(8)

Combined with formulae (2), (3a), and (7), this inequality indicates that the planet’s mean motion must be about five times smaller than its rotation rate:

$$n_{\text{p}}<0.198{\dot{\theta }}_{\text{p}}{\left(1-1.03e_{\text{p}}\right)}^{3/2}.$$(9)

Second, a synchronised moon’s pericentre, $a_{\text{m}}^{\left(\text{s}\right)}\left(1-e_{\text{m}}\right)$, must be higher than the planet’s Roche radius r_R. For low e_m, it is sufficient to simply compare $a_{\text{m}}^{\left(\text{s}\right)}$with r_R. By inserting Eqs. (6) and (7) into the inequality

$$r_{\text{R}}<a_{\text{m}}^{\left(\text{s}\right)}$$(10)

and approximating the moon’s density with that of our own Moon, we find (see Appendix A) that the planet’s rotation rate in the end state, when it is synchronised by the moon, must obey

$${\dot{\theta }}_{\text{p}}<1.07{\text{hr}}^{-1}.$$(11)

This sets a lower bound on the length of the sidereal day of the planet:

$$T_{\text{p}}=\frac{2\pi }{{\dot{\theta }}_{\text{p}}}>5.87\text{\hspace{0.17em}hr.}$$(12)

The sidereal year of almost all known exoplanets satisfies the same restriction⁶:

$$P_{\text{p}}=\frac{2\pi }{n_{\text{p}}}>5.87\text{\hspace{0.17em}hr.}$$(13)

Bear in mind that we are addressing the planets synchronised by their moons, not stars. So the coincidence of the low bounds (12) and (13) in no way implies the equality of ${\dot{\theta }}_{\text{p}}$ and n_p.

A close-in planet lacking a moon or having a moon incapable of synchronising the planet either will be tidally synchronized by the star or will end up in a higher spin-orbit resonance, if it has a sufficient eccentricity (Makarov et al. 2018). According to Eq. (150) from Boué & Efroimsky (2019), the tides in a synchronised planet will be working to reduce the planet’s semimajor axis. The planet will therefore spiral in and get engulfed by the star.

As an interesting aside, lower bound (13) on the planet period is remarkably close to the shortest rotation periods of field solar-type stars (Schmitt & Mittag 2020). This suggests that a moonless planet, if it migrates too close to the star, is not capable of synchronising the star rotation, and therefore gets destroyed. Indeed, had such a planet been able to orbit faster than 5.87 hr and to synchronise the star, then stars with such short rotation periods would have been observed – which is not the case.

On the other hand, a planet with a period of ~6 h can hardly harbour a moon, because the gap between the planet’s radius and reduced Hill radius (Sect. 2) vanishes in most cases or becomes quite small. We estimate that for our working sample of 580 inner exoplanets, the reduced Hill radius computed for the critical value of orbital period is above the actual planet’s radius for just 26 systems, with only 5 of them having a niche for exomoons that is wider than 1 R_p.

6 The niche for exomoons

For a moon to survive on an orbit around a planet, two conditions must be fulfilled. They may come out stronger than (8) and (10), because those were set on the end-state – while here we intend to address a moon’s parameters prior to synchronisation.

(1) Its apocentre must be confined to the reduced Hill sphere:

$$a_{\text{m}}\left(1+e_{\text{m}}\right)<r_{\text{H}}^{\prime };$$(14)

(2) Its pericentre must stay above the Roche sphere:

$$a_{\text{m}}\left(1-e_{\text{m}}\right)>r_{\text{R}}.$$(15)

The following sequence of inequalities ensues:

$$r_{\text{H}}^{\prime }>a_{\text{m}}>r_{\text{R}}.$$(16)

Most of the known exoplanets have $r_{\text{H}}^{\prime }>r_{\text{R}}$, but not all. We find 31 planets that possibly violate this rule, and, therefore, cannot harbour any satellites. Figure 2 shows them as dots lying below the diagonal line of equal radii. These planets have ultrashort periods between 0.447 day (TOI-561 b) and 3.313 days (HATS-10 b), with more than half (18) being shorter than 1 d. Their eccentricities are often unknown, with only upper limits given in the literature. For example, e < 0.519 for HATS-69 b. In our calculations involving the reduced Hill radius, this introduces additional dispersion and bias. The outlying planets are likely to have small eccentricities below the sensitivity threshold of the transit detection method. The planet HATS-69 b, for example, orbits an inconspicuous solar twin that does not have any known stellar companions or detectable outer planetary companions (Hord et al. 2021), in agreement with the eccentricity migration model. There are no alternative mechanisms to excite the eccentricity of close planets apart from outer massive per-turbers. Only one of the planets in this category, HIP 65 A b (with e = 0), resides in a wide stellar binary with a smaller red dwarf companion, which apparently failed to perturb its orbit because of the significant separation > 230 AU.

Additional useful constraints on the orbital parameters of the planet-moon system can be obtained from the brackets for the moon separation. By Kepler’s third law for the planet-moon and star-planet pairs, we have:

$$\frac{M_{\text{p}}{P}_{\text{m}}^2}{M_{\text{*}}{P}_{\text{p}}^2}\simeq \frac{a_{\text{m}}^3}{a_{\text{p}}^3}.$$(17)

Combining this approximate equality with inequality (4) and with the expression (2) for the reduced Hill radius, we obtain

$$\frac{P_{\text{m}}}{P_{\text{p}}}<0.198{\left(1-1.03e_{\text{p}}\right)}^{3/2}.$$(18)

This result is stronger than formula (4), because it pertains to a prograde moon at any stage of its evolution, not only at the synchronous end state (if that state is achieved).

Thus, prograde moons cannot have longer orbital periods than about 0.2 P_p. This condition is obviously satisfied for the Earth-Moon system. We note that this upper bound on the orbital period ratio is independent of the three masses.

On the other hand, utilising the lower bound on the moon separation a_m > r_R, Eq. (6), and Kepler’s third law, we derive the lower bound on the orbital period of the moon in its survival niche:

$$P_{\text{m}}>\frac{5.55}{\sqrt{G{\rho }_{\text{m}}}}=\left[3.26\text{\hspace{0.17em}hr}\right]\sqrt{{\rho }_{\text{Moon}}/{\rho }_{\text{m}}},$$(19)

where p_m is the average density of the moon and ρ_Moon is the average density of the Moon. This constraint depends only on the average density.

Fig. 2 Reduced Hill radii of 580 known inner exoplanets versus their estimated Roche limits in units of Jupiter’s radius. The dashed diagonal line marks the locus of equal radii where the dynamical niche for satellites vanishes. The planets that fall below this line cannot have any stable satellites.

7 Mutual synchronisation of a planet with its moon

In this work, we consider the case of a planet mutually synchronised with its moon. This configuration can be attained within two different scenarios, the first of which can be initiated by various means.

7.1 Scenario 1

Mutual synchronism is achieved through mutual receding accompanied with a slow-down of both bodies’ rotation. For this process to begin, the moon should find itself above the syn-chronicity orbit of the planet, and should not have too large an eccentricity⁷. This initial state can be prepared by three different methods:

(1.1) the moon is formed above the synchronicity orbit of the planet;

(1.2) the moon is captured above the synchronicity orbit of the planet;

(1.3) a rapidly rotating moon is formed or captured slightly below synchronism. If the initial spin is higher than the mean motion, and the separation from the synchronous orbit is small, the tides in the moon can overpower those in the planet and can push the moon above the synchronous radius.

Whichever of the three mechanisms places the moon above the planet’s synchronous orbit, the moon begins receding away, and the planet begins slowing down its spin.

7.2 Scenario 2

Mutual synchronism is achieved for a moon in a prograde orbit below the synchronous radius, or a moon in a retrograde orbit. In both cases, the tides in the planet work to shrink the moon’s orbit. The tides in the moon, if its spin is synchronised with the orbital motion, work in the same direction when e_m > 0. As the moon descends onto the planet, the planet’s spin is either accelerated, for a prograde moon, or decelerated, for a retrograde moon. In the latter case, the planet may stop rotating in the sidereal frame -and then start rotating in the opposite direction. If the planet’s angular acceleration rate is higher than the rate of moon’s orbit decay, the planet’s spin can be equalised with the orbital rotation before the moon reaches the Roche radius.

Our overall objective is to understand how the star-induced tidal torque is acting on an internally synchronised planet-moon system.

8 Tidal torques

To simplify the terminology, here we often employ the adjectives ‘solar’ and ‘lunar’ in application to the properties of the star and the moon. For example, instead of the clumsy ‘stargenerated tidal torque’ or ‘moon-generated tidal torque’, we say ‘solar torque’ or ‘lunar torque’ (with the word ‘tidal’ omitted, because no other torques are to be included).

8.1 Setting

Our study will address a two-dimensional configuration: the lunar orbit coincides with the ecliptic, while both the planetary and lunar obliquities onto the equator of the planet are zero. This will enable us to consider only the polar components of tidal torques.

A justification for the neglect of the lunar orbit inclination on the planet’s orbit comes from the assumption that a significant initial inclination would make the hierarchical triple star-planet-moon system subject to strong Lidov-Kozai interactions, which would increase interchangeably the inclinations and eccentricity of the inner pair (Innanen et al. 1997; Ford et al. 2000; Veras & Ford 2010; Batygin et al. 2011; Correia et al. 2011; Veras et al. 2018). As discussed in this paper, the dynamical niche for satellites is quite narrow for most known exoplanets. A large increase in eccentricity would push the reduced Hill radius limit down, possibly eliminating the niche altogether per Eq. (2). On the other hand, close-in moons can be disrupted if the pericentre distance becomes smaller than the Roche radius. This factor further diminishes the expected frequency of exomoons. The surviving moons are thus likely to have orbits nearly coplanar with the ecliptic.

Depending on the balance of the solar and lunar torques acting on the close planet, their latitudinal components should relatively quickly align the planet’s rotation axis with either the planet-star or the planet-moon orbit axis (see Sect. 9). In the latter case, the moons may be long-term viable near close planets. A limited nodal precession and nutation of the planet caused by the residual latitudinal component of the tidal torque created by the star adds to tidal dissipation in the planet and to the ensuing orbital evolution. Complex three-dimensional dynamical interactions in such nearly coplanar three-body systems are best investigated by numerical simulations, which are outside the scope of this study.

8.2 Secular part of the polar torque

A derivation of the secular polar tidal torque from the Darwin (1880) and Kaula (1964) theory is provided in Efroimsky (2012, Eq. (109)), see also Boué & Efroimsky (2019, Eqs. (123), (132), and (191)):

$$\begin{array}{c}\mathcal{T}=2\frac{G{M^{\prime }}^2}{a}{{\displaystyle \sum _{l=2}^{\infty }\left(\frac{R_{\text{p}}}{a}\right)}}^{2l+1}{\displaystyle \sum _{m=0}^l\frac{\left(l-m\right)!}{\left(l+m\right)!}m}\\ {\displaystyle \sum _{p=0}^l{F}_{lmp}^2\left(i\right){\displaystyle \sum _{q=-\infty }^{\infty }{G}_{lpq}^2}\left(e\right)K_{\text{l}}\left({\omega }_{\text{lmpq}}\right),}\end{array}$$(20)

where M′ is the mass of the perturber, G is the Newton gravity constant, while R_p is the radius of the planet. The expression involves the inclination functions F_lmp(i) and the eccentricity functions related to the Hansen coefficients by $G_{lpq}\left(e\right)=X_{l-p+q}^{-\left(l+1\right),l-2p}\left(e\right)$, see e.g., Murray & Dermott (1999, p. 232). The quantities

$$\begin{array}{l}{\omega }_{\text{lmpq}}\equiv \left(l-2p\right)\dot{\omega }+\left(l-2p+q\right)n+m\left(\stackrel{}{\text{Ω˙}}-{\dot{\theta }}_{\text{p}}\right)\hfill \\ \approx \left(l-2p+q\right)n-m{\dot{\theta }}_{\text{p}}\hfill \end{array}$$(21)

are the Fourier tidal modes over which both the potential of the perturber and the additional tidal potential of the body are expanded. These modes’ absolute values,

$${\chi }_{\text{lmpq}}=\left|{\omega }_{\text{lmpq}}\right|,$$(22)

are the actual physical forcing frequencies excited in the planet by the perturber⁸.

The order-l quality functions of the perturbed body are

$$K_{\text{l}}\left({\omega }_{\text{lmpq}}\right)\equiv k_{\text{l}}\left({\omega }_{\text{lmpq}}\right)\sin {ϵ}_{\text{l}}\left({\omega }_{\text{lmpq}}\right),$$(23a)

where both the Love numbers k₁ and the phase lags ϵ₁ are functions of the Fourier tidal modes ω_lmpq. They can also be written as⁹

$$K_{\text{l}}\left({\omega }_{\text{lmpq}}\right)\equiv \frac{K_{\text{l}}\left({\omega }_{\text{lmpq}}\right)}{K_{\text{l}}\left({\omega }_{\text{lmpq}}\right)}\mathrm{sign}\left({\omega }_{\text{lmpq}}\right),$$(23b)

where Q₁ are the tidal quality factors introduced via

$$Q_{\text{l}}^{-1}\left({\omega }_{\text{lmpq}}\right)=\left|\sin {ϵ}_{\text{l}}\left({\omega }_{\text{lmpq}}\right)\right|.$$(24)

Phase lags ϵ₁(ω_lmpq) being odd functions, their absolute values $Q_1^{-1}\left({\omega }_{\text{lmpq}}\right)$ are obviously even. Even are also the Love numbers k₁(ω_lmpq). So, altogether, the quality functions are odd. A function K₁(ω_lmpq) changes its sign when the system is crossing the Impq spin–orbit resonance, one defined by ω_lmpq = 0. On this crossing, the passage of the function K₁(ω_lmpq) through the zero value goes rapidly but smoothly (see Eq. (32)).

Since the functions k₁(ω_lmpq) and $Q_1^{-1}\left({\omega }_{\text{lmpq}}\right)$ are even, we always can regard both the Love numbers and quality factors as functions of the positive definite physical frequencies (22):

$$k_{\text{l}}\left({\omega }_{\text{lmpq}}\right)=k_{\text{l}}\left({\chi }_{\text{lmpq}}\right),Q_{\text{l}}\left({\omega }_{\text{lmpq}}\right)=Q_{\text{l}}\left({\chi }_{\text{lmpq}}\right).$$(25)

Combined with Eq. (23b), this gives:

$$K_{\text{l}}\left({\omega }_{\text{lmpq}}\right)=K_{\text{l}}\left({\chi }_{\text{lmpq}}\right)\mathrm{sign}\left({\omega }_{\text{lmpq}}\right).$$(26)

The quality functions appear not only in the expression for the tidal torque, but also in the formulae for the tidal heating rate and the tidal evolution rates of orbital elements. A detailed discussion of the generic shape of a quality function K₁(ω_lmpq) is provided in Appendix B. There, it is explained that such a function has the form of a kink, as in Fig. 3. For simple rheologies (Maxwell, Andrade), the kink has only one sharp maximum on the right and one sharp minimum on the left. These principal extrema come into being due to interplay of self-gravitation with rheological response¹⁰. For more elaborate rheological models containing Debye peaks (i.e., Sundberg-Cooper or Burgers), additional maxima – and symmetrically located minima – will emerge. Additional local extrema show up also when a multilayer structure of the body is taken into consideration¹¹. One way or another, the principal extrema shown in Fig. 3 exist for any realistic rheology and supersede in magnitude other possible peaks.

Dictated by hydrodynamical effects, tidal response of gas planets is much more complex. Nonetheless, qualitatively, the behaviour near a resonance must be similar: the largest peaks are defined, at large, by combination of material response and self-gravitation, and the transition through zero must be smooth.

As demonstrated in Appendix B for a homogeneous Maxwell sphere (in realistic situations, also for a homogeneous Burgers, Andrade, or Sundberg-Cooper sphere, see Footnote 17), the depicted in Fig. 3 two extrema of K₁(ω) are positioned at

$${\omega }_{pea{k}_l}=\pm \frac{{\tau }_M^{-1}}{1+{\mathcal{A}}_l},$$(27)

their values being

$$K_{\text{l}}^{\left(\text{peak}\right)}=\pm \frac{3}{2\left(l-1\right)}\frac{{\mathcal{A}}_l}{1+{\mathcal{A}}_l}.$$(28)

In these expressions, 𝒜_l are dimensionless effective rigidities rendered by Eq. (B.6). Proportional to the mean shear rigidity µ, these parameters serve as dimensionless measures of rheology-dominated versus gravity-dominated tidal responses, see Sects. 4.1 and 5.2 in Efroimsky (2015). For the Earth, Mars, and the Moon, the dimensionless parameter 𝒜₂ assumes the values of 2.2, 19, and 64.5, correspondingly. Assuming that an exomoon is smaller than Mars, we have¹²

$$\frac{{\mathcal{A}}_l}{\left(1+{\mathcal{A}}_l\right)}\simeq 1,\text{for moons},$$(29)

wherefrom

$$K_{\text{l}}^{\left(peak\right)}\simeq \pm \frac{3}{4\left(l-1\right)},\text{for moons},$$(30)

an observation to be utilised down the road.

Equation (28) also renders an overall limit on the values of a quality function:

$$\left|K_{\text{l}}\right|<\frac{3}{4\left(l-1\right)}$$(31)

Within the narrow inter-peak interval, the functions K_l(ω) are near-linear:

$$\left|\omega \right|<\left|{\omega }_{pea{k}_l}\right|\Rightarrow K_{\text{l}}\left(\omega \right)\simeq \frac{3}{2\left(l-1\right)}\frac{{𝒜}_l}{1+{𝒜}_l}\frac{\omega }{\left|{\omega }_{pea{k}_l}\right|},$$(32)

and fall off as ω⁻¹ outside it:

$$\left|\omega \right|>\left|{\omega }_{pea{k}_l}\right|\Rightarrow K_{\text{l}}\left(\omega \right)\simeq \frac{3}{2\left(l-1\right)}\frac{{𝒜}_l}{1+{𝒜}_l}\frac{\left|{\omega }_{pea{k}_l}\right|}{\omega },$$(33)

an expression to be employed shortly.

Fig. 3 Typical shape of the quality function K1(ω)=k1(ω) sin ϵ1(ω), where ω is a shortened notation for the tidal Fourier mode ωlmpq. (Reprinted from Icarus, Vol. 241, Noyelles et al., Spin-orbit evolution of Mercury revisited, pages 26–44, Copyright (2014), with permission from Elsevier.)

8.3 Approximation valid under no synchronism

The secular part of the polar torque may be approximated with

$$\begin{array}{l}𝒯={𝒯}_{l=2}+O\left(ϵ{\left(R/a\right)}^7\right)\hfill \\ ={𝒯}_{\left(\text{lmp}\right)=\left(220\right)}+O\left(ϵi^2\right)+O\left(ϵ{\left(R/a\right)}^7\right),\hfill \end{array}$$(34)

where ϵ is a typical value of a phase lag (Efroimsky 2012, Eq. 114). Details on this are provided in Appendix C, both for the solar torque acting on the planet and for the lunar torque acting on the planet. When the planet is not synchronised with either orbital motion, the leading parts of the lunar and solar torques are

$${𝒯}_{\left(\text{lmp}\right)=\left(2200\right)}^{\left(m\right)}=\frac{3}{2}\frac{G{M}_{\text{m}}^2}{a_{\text{m}}}{\left(\frac{R_{\text{p}}}{a_{\text{m}}}\right)}^5{K}_2\left(2n_{\text{m}}-2{\dot{\theta }}_{\text{p}}\right)$$(35)

and

$${𝒯}_{\left(\text{lmp}\right)=\left(2200\right)}^{(\ast )}=\frac{3}{2}\frac{G{M}_{\ast }^2}{a_{\text{p}}}{\left(\frac{R_{\text{p}}}{a_{\text{p}}}\right)}^5{K}_2\left(2n_{\text{p}}-2{\dot{\theta }}_{\text{p}}\right),$$(36)

correspondingly. Hence, in neglect of O(e²), the planetary rotation is described by

$$\begin{array}{l}{\ddot{\theta }}_{\text{p}}=\frac{3}{2}G{\xi }^{-1}{M}_{\text{p}}^{-1}{R}_{\text{p}}^3[M_{\text{m}}^2{a}_{\text{m}}^{-6}{K}_2\left(2n_{\text{m}}-2{\dot{\theta }}_{\text{p}}\right)\hfill \\ +M_{\ast }^2{a}_{\text{p}}^{-6}{K}_2\left(2n_{\text{p}}-2{\dot{\theta }}_{\text{p}}\right)],\hfill \end{array}$$(37)

ξ being the prefactor of the moment of inertia $\xi M_{\text{p}}{R}_{\text{p}}^2$ of the planet. The lunar and stellar inputs in expression (37) relate as

$$\frac{{𝒯}_{\left(\text{lmp}\right)=2200}^{\left(m\right)}}{{𝒯}_{\left(\text{lmp}\right)=2200}^{(\ast )}}={\left(\frac{M_{\text{m}}}{M_{\ast }}\right)}^2{\left(\frac{a_{\text{p}}}{a_{\text{m}}}\right)}^6\frac{K_2\left(2n_{\text{m}}-2{\dot{\theta }}_{\text{p}}\right)}{K_2\left(2n_{\text{p}}-2{\dot{\theta }}_{\text{p}}\right)}$$(38)

If the moon is massive and close enough, it can synchronise the planet. As the system is approaching synchronicity $\left(\text{i.e.}{\dot{\theta }}_{\text{p}}\to n_{\text{m}}\right)$, the quality function $K_2\left(2n_{\text{m}}-2{\dot{\theta }}_{\text{p}}\right)$ in the numerator of (38) goes through a sharp peak and rapidly approaches zero (see Fig. 3). The solar torque will then try to drive the planet out of its synchronism with the moon. Whether it will succeed in doing this will depend on the ratio between the peak value of the lunar torque and a non-peak value of the solar torque.

9 Two pathways to a synchronous planet-moon system

As we explained in Sect. 7, mutual synchronism of the planet and the moon can be attained within two distinct dynamical scenarios having different initial conditions.

9.1 Scenario 1: the moon is initially on a wide prograde orbit

The first pathway to the 1:1 inner resonance is for the moons that are initially above the synchronous radius. By construction, the initial mean motion of such a moon falls short of the initial rotation rate of the planet: $n_{\text{m}}<{\dot{\theta }}_{\text{p}}$. Both the tides raised in the planet by the moon and the tides raised in the planet by the star are working in the same direction – to slow down the spin of the planet and, consequently, to increase its syn-chronicity radius. The moon is receding from the planet (and, possibly, is increasing its eccentricity e_m). The planet-moon 1:1 spin-orbit resonance can be reached if the growing synchronous radius catches up with the expanding orbit while the moon is still within the reduced Hill sphere. We now show that the vicinity of this resonance may be stable under certain conditions, which are plausible.

9.1.1 Approach to synchronism

Under the exact planet–moon synchronicity, i.e., for $n_{\text{m}}={\dot{\theta }}_{\text{p}}$, the e⁰ part of the lunar torque (C.1) vanishes, while the O(e²) part of this torque survives and becomes leading. It writes as:

$${𝒯}_{\text{moon-synchr}}^{\left(m\right)}=12e_{\text{m}}^2\frac{G{M}_{\text{m}}^2}{a_{\text{m}}}{\left(\frac{R}{a_{\text{m}}}\right)}^5{K}_2\left(n_{\text{m}}\right)+O\left(e_{\text{m}}^4ϵ\right),$$(39)

the stellar torque still being given by expression (36).

Quadratic in eccentricity, torque (39) may appear to be much smaller than the torque exerted by the star. However, it is not the small torque at the exact resonance but the maximum torque value at a slightly shifted peak of the tidal quality function that matters in the subtle mechanism of resonance capture.

In the vicinity of the 1:1 resonance, the quadrupole tidal quality of the planet, K₂(ω₂₂₀₀), is an odd function of the semi-diurnal tidal mode

$${\omega }_{2200}=2\left(n_{\text{m}}-{\dot{\theta }}_{\text{p}}\right),$$(40)

as in Fig. 3, with two symmetric extrema bracketing the point of resonance crossing. Were it not for the star, the 1:1 planet-moon spin-orbit resonance would imply $n_{\text{m}}={\dot{\theta }}_{\text{p}}$, and consequently ω₂₂₀₀ = 0, and therefore ${𝒯}_{\left(\text{lmp}\right)=\left(2200\right)}^{\left(m\right)}=0$ according to Eq. (35).

In a triple system under consideration, the tides generated by the star in the planet are working to drive the planet out of its spinorbit resonance with the moon. Aiming to synchronise the planet not with the moon but with the star, the solar tides are slowing down the spin of the planet. For this reason, the tidal mode ω₂₂₀₀ is not exactly zero but acquires a slightly positive value. The value of the quality function K₂ begins to climb, from the left, the steep positive peak in Fig. 3. Therefrom emerges a positive lunar torque compensating for the negative solar torque – and the equilibrium stays stable.

For the planet to leave said equilibrium, that is, to break out of the resonance, the tidal torque from the star has to overpower the peak tidal torque from the moon. The maximal tidal quality value achieved at the peak frequency depends on the effective rigidity 𝒜_l via Eq. (28). For l = 2, the tidal quality of the planet cannot exceed 3/4, for a homogeneous planet. For our purposes, it will be justified to use in the following estimates the maximal possible value of the torque exerted on the planet by the moon. Using this value, we write the ratio of the torques (given by Eqs. (30) and (36), correspondingly) wherewith the star and the moon are acting on the planet:

$$\frac{\left|{𝒯}^{(\ast )}\right|}{{𝒯}^{\left(m,\text{peak}\right)}}\simeq \frac{4}{3}{\left(\frac{M_{\ast }}{M_{\text{m}}}\right)}^2{\left(\frac{a_{\text{m}}}{a_{\text{p}}}\right)}^6\left|K_2\left(2n_{\text{p}}-2{\dot{\theta }}_{\text{p}}\right)\right|$$(41)

For the moon to remain locked in the 1:1 spin–orbit resonance, this ratio must be less than unity.

The equilibrium state of a synchronised planet–moon system depends on the composition of the planet. To illustrate this, we resort to Fig. 4. There, the blue solid line depicts the leading, semi-diurnal component (35) of the lunar tidal torque applied to the planet. Proportional to the semi-diurnal quality function $K_2\left({\omega }_{2200}\right)=K_2\left(2\left(n_{\text{m}}-{\dot{\theta }}_{\text{p}}\right)\right)$, this torque naturally comes out kink-shaped. Bear in mind, though, that in Fig. 4 this torque is given as a function of the difference ${\dot{\theta }}_{\text{p}}-n_{\text{m}}$, not of the semidiurnal frequency ${\omega }_{2200}=2\left(n_{\text{m}}-{\dot{\theta }}_{\text{p}}\right)$. The red solid line in Fig. 4 is the opposite torque with which the planet is acting on the moon’s orbit.

Within Scenario 1, the moon starts above synchronism, so initially we have ${\dot{\theta }}_{\text{p}}-n_{\text{m}}>0$. As the moon is receding, and the planet is slowing down its rotation, the system is moving, in Fig. 4, from the right to the left. At the point of exact 1:1 resonance, where ${\dot{\theta }}_{\text{p}}-n_{\text{m}}=0$, the tidal torque generated by the moon vanishes, while the negative torque applied by the star on the planet (the dashed blue line) is finite. For a rigid planet with a finite triaxial elongation, this negative solar torque gets compensated by a positive torque from the moon, caused by a constant tilt of the planet’s longest axis with respect to the moon direction, which resembles the constant tilt of the Moon with respect to the mean Earth direction (Rambaux & Williams 2011). This way, a planet having a permanent dynamic triaxiality (a terrestrial or icy world) can be captured into an exact spin–orbit resonance with the moon, due to the presence of a restoring triaxial torque.

Inner giant planets are more likely to be fluid, apart from a possible compact rigid core. If their effective viscosity is lower than a critical value (Makarov 2015), the planet gets locked in a pseudo-sychronous rotation, a long-term stable equilibrium where the relative frequency ${\dot{\theta }}_{\text{P}}-n_m$ changes with time very slowly, as we shall see now.

Since no restoring triaxial torque is acting on a fluid planet at the point of exact 1:1 synchronism, evolution of this planet does not stop on arrival to that resonance. Indeed, for ${\dot{\theta }}_{\text{p}}-n_{\text{m}}>0$ the net torque is negative due to the presence of the negative torque from the star. As a result of this, the angular acceleration ${\ddot{\theta }}_{\text{p}}$ is negative there. In Fig. 4, the system transcends the resonance and keeps moving leftwards – and approaches the dotted vertical line, which denotes a situation where the two torques acting on the planet are compensating each other. There, the angular acceleration of the planet becomes zero: ${\ddot{\theta }}_{\text{p}}=0$. This state is still not the end point, and gets transcended too, because a negative torque from the planet (shown by the red curve) is acting on the lunar orbit. Impelled by this torque, the moon is now slowly descending, and the value of n_m is increasing. Thus, the negative relative frequency ${\dot{\theta }}_{\text{p}}=n_{\text{m}}$ is decreasing, and the system keeps moving leftwards in the figure.

The system eventually reaches a quasi-equilibrium point residing somewhere to the left of the dotted vertical line, but to the right of the peak tidal torque. The equilibrium point is where the relative angular acceleration equalises, i.e. ${\ddot{\theta }}_{\text{p}}={\dot{n}}_m$ This equilibrium, however, will be slowly evolving. Indeed, in the interval between said peak and the dotted line, the net torque on the planet is positive, resulting in a spin-up. The net torque on the moon’s orbit is negative, so the moon is very slowly descending back onto the planet. Therefore, in Scenario 1, the moon of a fluid planet is initially rapidly receding, and then, on arrival to equilibrium, begins to slowly descend.

In this complex pseudo-equilibrium, there is one source of energy (the orbital energy of the moon) and three sinks (tidal dissipation in the planet, planet’s rotation spin-up, and the expansion of the planet’s orbit around the star). From the energy exchange aspect, the necessary condition for this equilibrium to be stable is ${\dot{E}}_{\text{kin}}/\left|{\dot{E}}_{\text{orb}}\right|<1$, that is to say, the rate of the released orbital energy should be greater than the rate of the planet’s rotational energy. This condition coincides with the requirement that the right-hand part of Eq. (50) should be less than 1. The deviation from the exact synchronicity is expected to slowly diminish due to the decay of the moon’s orbit and the corresponding increase of the amplitude of the tidal torque.

The relative contribution of the lunar tide (as compared to the solar one) diminishes and becomes minimal near the reduced Hill radius. Using inequality (4) and the Kepler law, we obtain the condition on the torque ratio required for the moon to remain locked in the 1:1 spin–orbit resonance up to the Hill radius as

$$\frac{\left|{𝒯}^{(\ast )}\right|}{{𝒯}^{\left(m,\text{peak}\right)}}\le 0.002{\left(\frac{M_{\text{p}}}{M_{\text{m}}}\right)}^2\left|K_2\left(2n_{\text{p}}-2{\dot{\theta }}_{\text{p}}\right)\right|<1$$(42)

For realistic moons whose mass is about 1% of the planet’s mass, this condition is fulfilled if the tidal quality of the planet at the solar-tide frequency satisfies

$$\left|K_2\left(2n_{\text{p}}-2{\dot{\theta }}_{\text{p}}\right)\right|<\frac{1}{20}$$(43)

Inequality (43) is easily obeyed even by Earth-like terrestrial planets, because the rate of the planet’s rotation in our scenario is

$${\dot{\theta }}_{\text{p}}\simeq n_{\text{m}}\gg n_{\text{p}},$$(44)

and therefore the semi-diurnal tidal frequency

$$\chi ={\chi }_{2200}=\left|{\omega }_{220}\right|=\left|2n_{\text{p}}-2{\dot{\theta }}_{\text{p}}\right|\simeq 2n_{\text{m}}$$(45)

is much greater than the degree-two peak frequency

$${\chi }_{{\text{peak}}_2}=\left|{\omega }_{pea{k}_2}\right|=\frac{{\tau }_M^{-1}}{1+{𝒜}_2},$$(46)

which is close to zero. For $\chi \gg {\chi }_{\text{peak\hspace{0.17em}}2}$, the quadrupole quality function is given by expression (33), with l = 2:

$$K_2\left(\omega \right)=\frac{3}{2}\frac{{𝒯}_2}{{𝒯}_2}\frac{\left|{\omega }_{pea{k}_2}\right|}{\omega }=\frac{3}{2}\frac{{𝒯}_2}{{\left(1+{𝒯}_2\right)}^2}\frac{1}{\omega {\tau }_M}$$(47)

The third fraction in this expression, 1/(ωτ_m), assumes very small values for planets with the mean Maxwell times longer than their orbital periods around the host stars (like the Earth or Mars – planets whose Maxwell times are in the hundreds of years). This warrants the smallness of K₂(ω) of these planets.

The mean Maxwell times of hot Jupiters and hot super-Earths may be comparable to or even shorter than these planets’ orbital periods, in which case the factor 1/(ωτ_m) is not warranted to be small. For such planets, however, the effective rigidities 𝒜₂ are very large¹³. Hence, the ratio 𝒜₂/ (1 + 𝒜₂)² is vanishingly small, and the condition for resonance is still fulfilled, even at the maximum planet-moon separation equal to the reduced Hill radius.

Fig. 4 Schematic depiction of the three torques involved in the quasi-stable pseudo-synchronous equilibrium of the planet-moon system. The tidal torque wherewith the moon is acting on the planet (solid blue line) is coupled with an opposite torque wherewith the planet is acting on the orbit of the moon (red curve). The negative tidal torque on the planet from the star (dashed blue line) is weakly dependent on the relative frequency. The vertical dotted line indicates the location where the tidal torques acting on the planet are equalised.

9.1.2 Competition in speed

In Scenario 1, the 1:1 resonance, $n_{\text{m}}\simeq {\dot{\theta }}_{\text{p}}$, must be attained before the moon reaches its reduced Hill radius. Since in this scenario the moon starts above synchronism $\left({\dot{\theta }}_{\text{p}}-n_{\text{m}}>0\right)$, a configuration with ${\dot{\theta }}_{\text{p}}-n_{\text{m}}=0$ can be achieved only if $\frac{d}{dt}\left({\dot{\theta }}_{\text{p}}-n_{\text{m}}\right)<0$, that is to say, if

$${\ddot{\theta }}_{\text{p}}<{\dot{n}}_{\text{m}}$$(48)

Both these rates being negative, we can rewrite the above as

$$\left|{\ddot{\theta }}_{\text{p}}\right|>\left|{\dot{n}}_{\text{m}}\right|$$(49)

Thus, the absolute rate of the planet’s spin-down (which is equivalent to the lengthening of the day on Earth) must be higher than the absolute rate of moon’s orbital slowdown.

As demonstrated in Appendix D, in a two-body planet-moon system (i.e., in neglect of the torques from star) the ratio of the orbital expansion and the planet’s spin-down rate, for e_m = 0 and M_p M_m/(M_p + M_m) ≈ M_m, is:

$$\frac{{\dot{n}}_{\text{m}}}{{\ddot{\theta }}_{\text{p}}}\approx 3\xi \frac{M_{\text{p}}}{M_{\text{m}}}{\left(\frac{R_{\text{p}}}{a_{\text{m}}}\right)}^2$$(50)

which is independent of the tidal quality or a specific tidal model. Within the considered scenario, both the orbital angular acceleration in the numerator and the spin acceleration in the denominator of this expression are initially negative.

From Eq. (49), the necessary condition for the moon’s survival is that the value of ratio (50) must become less than 1, beginning from some value a_m. Assuming for an estimate that the factor 3ξ equals 1, we get a limit on the moon’s semi-major axis value, starting from which inequality (49) must become valid:

$$a_{\text{m}}>a_{\text{crit}}=R_{\text{p}}\sqrt{\frac{M_{\text{p}}}{M_{\text{m}}}}$$(51)

If condition (49) never gets obeyed, the moon will leave the reduced Hill sphere before synchronising the planet, and will eventually be lost. On the other hand, the fulfilment of condition (49) since the birth of the moon is not necessary for synchronisation to happen. Right after accretion, the moon must be located above synchronism (as required by Scenario 1), but it may still be so close to the planet, that the right-hand side of Eq. (50) is larger than unity. In this situation, while both ${\dot{n}}_{\text{m}}$ and ${\ddot{\theta }}_{\text{p}}$ are negative, the moon’s orbital recession is initially going faster than the planet’s despinning. In the course of recession, a_m is growing, and, according to Eq. (50), the evolution of the difference ${\dot{\theta }}_{\text{p}}-n_{\text{m}}$ may reverse. If, beginning from some value of a_m, inequality (49) gets satisfied, the moons acquires a chance to synchronise its host planet.

A moon of 1% relative mass should be initially separated by more than ten planet’s radii. The Moon, for example, is well outside this critical radius, and the lengthening of the day is faster than the expansion rate of the lunar orbit. So the Moon will eventually synchronise the Earth¹⁴. Equation (50) does not take into account the contribu_.t_.ion of the solar tidal brake, which in this scenario will make $\left|{\ddot{\theta }}_{\text{p}}\right|$ larger and push the lower bound of a_m slightly lower.

9.1.3 Scenario 1 is realisable only for small planets

Replacing, in Eq. (51), a_m with the reduced Hill radius $r_{\text{H}}^{\prime }$, we estimate the minimum relative mass of the moon, M_m/M_p, required to achieve the critical radius within the reduced Hill sphere:

$$\frac{M_{\text{m}}}{M_{\text{p}}}>{\left(\frac{R_{\text{p}}}{r{\prime }_{\text{H}}}\right)}^2$$(52)

Figure 5 shows the results for 580 inner exoplanets from our working sample, which have all the required values in the database. The horizontal lines are shown for two specific M_m/M_p values, 0.01 and 0.05, as benchmarks of regular and massive moons. Only those planets that are located above these lines can possibly have moons of the corresponding mass. We note that the majority of hot Jupiters lie below the limit even for the most imaginably massive moons, because for such Jupiters inequality (52) is violated. Hence, they cannot be synchronised by moons, and such systems should be long-term unstable. The situation is different with smaller planets. Generally, we count 58 planets that can be synchronised by a 0.01 M_p moon (10% of the general sample) and 150 (26%) by a 0.05 M_p moon.

Fig. 5 Reduced Hill radii of 580 known inner exoplanets versus their radii in units of Jupiter’s radius. Dashed horizontal lines mark the lower boundaries of plot areas where these exoplanets can in principle have moons with the smallest masses as labelled in the graph. For each specific |exoplanet shown with a dot in this plot, the mass of a survivable moon should be above the dashed line (above the upper line, for putative moons of Mm = 0.01 Mp, or above the lower line, for putative moons of Mm = 0.05 Mp).

9.2 Scenario 2: the moon is initially below the synchronous radius

We now consider an initial configuration, in which the moon is sandwiched between the Roche radius and the synchronous radius. The moon moves fast with respect to the planet’s surface, and the tide raised by it lags the centre line. The tide brakes the orbital momentum, and the moon’s orbit decays. It appears that the moon is doomed to reach the Roche limit and disintegrate. However, the same tidal interaction makes the planet’s rotation accelerate. The outcome of these two parallel processes depends on the ratio of the characteristic times of orbital decay and planet’s spin-up.

In this scenario, the system from the start satisfies ${\dot{\theta }}_{\text{p}}<n_{\text{m}}$, where, by distinction from the previous scenario, both rates are now positive. Equation (50) is still valid, and this ratio must be always less than 1, in order for the moon to have a chance to synchronise the planet and survive. There is no possibility of survival for the moon if this ratio is initially greater than 1, because as the orbit shrinks, the ratio of the accelerations can only increase – so the moon will simply spiral onto the Roche limit, with the planet still rotating slower than the moon’s orbital motion. Thus, a significant initial separation between the planet and the moon (e.g., greater than ~10 planet’s radii for M_m/M_p = 0.01; cf. Eq. (51)) is required. The initial separation is limited from above by the reduced Hill radius. Relatively few inner planets have sufficiently large reduced Hill radii, cf. Figs. 2 and 5. Thus, most of the known exoplanets, if they have moons located below the synchronous radius, will destroy these moons before their rotation becomes synchronised.

The major difference between Scenario 1 and Scenario 2 is that in the former the tidal torques from the star and from the moon, initially, are both negative, whereas in the latter, they are likely to be counter-directed. The torque from the rapidly moving moon is initially positive in Scenario 2, and the torque from the star is negative, unless the planet initially rotates slower than its mean orbital motion. Equation (38) is valid, given that the torque from the star and the corresponding quality function are negative. The planet-moon system approaches the point of equilibrium from the left in Fig. 4, and that point is located to the left of the peak tidal torque. The general necessary condition is that the torque from the moon overpowers the torque from the star. The exact relation between the two torques in Scenario 2 cannot be estimated without knowing the planet’s rate of rotation and its rheological properties. If we ignore the third fraction on the right-hand side of Eq. (38) for a very rough estimate, the minimum mass of the moon can be estimated by using condition (19), assuming that the moon’s average density equals that of the Moon, and requiring that the ratio of the torques should be greater than 1 in absolute value:

$$M_{\text{m}}>M_{\text{p}}{\left[\frac{3.26\text{\hspace{0.17em}hr}}{P_{\text{p}}}\right]}^2$$(53)

This turns out to be a relatively loose requirement even for the close-in exoplanets. Out of 1366 planets with available estimates of P_p and M_p in the NASA Archive, 735 (54%) have lower bounds on M_m/M_Moon within 1, and 834 (61%) within 5. There is a caveat in this consideration, though. The outcome of the competition in speed also depends on the initial rate of planet’s rotation ${\dot{\theta }}_{\text{p}}$ and its tidal properties. If the initial rotation is very slow, ${\dot{\theta }}_{\text{p}}<n_{\text{p}}$, or even retrograde, the planet also has to cross the 1:1 spin–orbit resonance with the star before it can be synchronised by the moon. To avoid capture into this resonance, the quality function in the numerator in Eq. (38), which is approximately $K_2\left(2n_{\text{m}}-2n_{\text{p}}\right)$, should not be much smaller than the peak quality of the tide caused by the star. A wide range of fairly complex scenarios emerge for planets of terrestrial composition. The most important and poorly known parameters are the effective viscosity of the mantle, how close the mantle is to solidus (Makarov et al. 2018), and the possible contribution of the Andrade-type dissipation mechanisms.

If the braking torque from the star overpowers the torque from a small or distant moon, the planet gets synchronised by the star. The moon will then inevitably spiral in to the Roche radius and disintegrate. The characteristic times for this process strongly depends on a number of planet’s parameters including its rheology. Somewhat counterintuitively, Earth-like rigid planets may be less efficient in dissipating tidal energy because of the large difference between n_p and n_m, i.e., the high tidal frequency. The time of orbital decay for such systems (and the life time of their satellites) is also dependent on the residual eccentricity of the moon.

10 Discussion

We have demonstrated, on very general principles, that exo-moons can orbit most of known exoplanets only within a relatively thin shell of a few to several Jupiter radii (Figs. 1 and 2).

The dynamics of triple star-planet-moon systems is determined by the interplay of direct gravitational interaction, on the one hand, and the tidally caused exchange of spin and orbital angular momenta, on the other. For various initial configurations, this interplay leads to a multitude of possible scenarios of dynamical evolution (Dobos et al. 2021; Tokadjian & Piro 2020, 2022). Most of them, however, predict a short-lived moon destined either to spiral in to the planet’s Roche radius or to retreat from the planet and be ejected by the gravitational pull of the star. In this study, we focus on the remaining possibilities for hypothetical exomoons to survive over longer timescales. These possibilities are critically dependent on the ability of the moon to synchronise the planet, that is to say, to overpower the tidal action from the star on the planet. This scenario is not far-fetched, for it is actually epitomised by our own Moon – which is in the process of synchronising the Earth. The settings for most of the known exoplanets are much less accommodating, however. These planets are closer to their hosts and are larger than the Solar System planets. On the one hand, in order to exert a strong torque on the planet, a moon should be sufficiently close to it and sufficiently massive. On the other hand, a very close moon of a rapidly rotating planet would undergo fast orbital expansion (cf. Eq. (50)) with a likely outcome of being driven outside the Hill sphere before the planet becomes synchronised. This places an additional lower bound on the relative mass of the moon, M_m/M_p, as expressed by Eq. (52) and depicted in Fig. 5. For each specific exoplanet shown with a dot in this plot, the mass of a surviv-able moon should be well above the dashed lines (above the upper line, for putative moons of M_m = 0.01 M_p, or above the lower line, for putative moons of M_m = 0.05 M_p). Only small fractions of the available sample satisfy this condition: 10% at M_m/M_p > 0.01, and 26% at M_m/M_p > 0.05. The majority of discovered exoplanets are massive, and the median limiting mass 0.05 M_p equals 917 M_Moon, while 0.01 M_p equals to 183 M_Moon. It is doubtful that such massive satellites are common, or that any even exist. Limiting our consideration to satellites less massive than 7 M_Moon, we estimate that only 10% of all known planets satisfy the M_m/M_p > 0.01 requirement. We therefore predict that exoplanet systems stabilised by their moons should be quite rare in the current compendium, perhaps a few per hundred.

Another set of constraints on a planet-stabilising moon can be derived from the condition of relative angular acceleration (given by Eq. (50)) being less than 1. Irrespective of the rheolog-ical properties of the planet, a moon has a chance to synchronise its host planet only if ${\dot{n}}_{\text{m}}<{\ddot{\theta }}_{\text{p}}$. Otherwise, because of the narrowness of the exomoon niche, it is likely to bump into either the reduced Hill limit or the Roche limit, depending on the initial ${\ddot{\theta }}_{\text{p}}$. The initial a_m is arbitrary in this equation, but it should be between r_R and $r_{\text{H}}^{\prime }$ Substituting a_m with these limits, and assuming 3 ξ = 1, we can estimate the corresponding limits M_ and M₊ for the mass of the moon. The results for 580 exoplanets with sufficient data in the database are presented in Fig. 6. We note that both these limits are one-sided, in that the mass should be greater than the limit, for the moon to survive. The M_p > M₋ condition is applicable in the Scenario 1 pathway, wherein the moon’s orbit is initially expanding. The M_p > M₊ condition is applicable in the Scenario 2 pathway, wherein the moon’s orbit is initially shrinking. The graph shows that the bulk of known inner exoplanets, seen as a dense cloud in the upper-right corner, require extremely massive satellites in the range of Earth-like planets for the pathways to work. We estimate the number of exo-planets with survivable moons within 5 M_Moon to be 28 (out of 580) for pathway 1, and only 2 for pathway 2. Thus, exomoons that can survive in the vicinity of known close-in planets over stellar lifetime durations should be rare.

This leads to question of how theoretical cogitations can be tested using observations. Exomoons, even unusually massive ones, cannot be imaged in the glare of the close host stars. Exoplanets harbouring synchronising moons should have finite eccentricities because of the residual tidal torque couple to the star. They may also be impervious to the tidal orbital decay, which is expected from moonless close-in planets. Depending on the balance of the tidal torques in the star-planet pair, a moon-synchronous planet may exhibit orbital expansion instead of orbital decay. These are indirect indicators, however, which allow for an alternative interpretation. The best promise is offered by precision photometry of transiting exoplanets. The light reflected by the planet modulates the general light curve due to the planetary phases. The peak flux is expected to coincide with the upper conjunction, which may be observable as a shallow secondary transit. Hot Jupiters are also expected to irradiate substantial amounts of infrared light from their surfaces. Due to thermal inertia, rapidly rotating planets would have hot spots with phase shifts leading the planet-star directions, which would result in an observable shift between the overall light curve and the eclipse in the upper conjunction. This, however, can be ambiguously interpreted as a finite eccentricity of the star-planet pair. Despite these difficulties, a search for exomoons has a tremendous impact on our understanding of cosmogony of other planetary worlds, and it should continue along both theoretical and practical routes. The two tidal-orbital pathways to long-term survival of exo-moons in the harsh conditions of known systems suggest that the fate of close-in exoplanets may be intertwined with the fate of their satellites.

Fig. 6 Lower masses of exomoons that can in principle synchronise their host planets for 580 known exoplanet systems, in units of the mass of the Moon. The M− limits are derived from the estimated reduced Hill radii and refer to the Scenario 1 pathway, when the exomoon is initially above the synchronous radius. The M+ limits are derived from the estimated Roche radii and refer to the Scenario 2 pathway, when the exomoon is initially below the synchronous radius.

Acknowledgements

The authors are grateful to James L. Hilton for useful comments and recommendations. We also thank our anonymous reviewer for a substantial and deep review, which resulted in a more consistent and focused paper.

Appendix A Comparing the Roche radius and the synchronous value of the semi-major axis

We wish to know for what critical spin rate ${\dot{\theta }}_p^{\left(crit\right)}$ the semi-major axis of the synchronous orbit of the moon about the planet, $a_{\text{m}}^{\left(s\right)}$, coincides with the planet’s Roche radius r_R. Formulae (6) and (7) entail:

$$\begin{array}{l}{r}_R<a_m^{\left(s\right)}\Rightarrow \hfill \\ {\dot{\theta }}_p<\sqrt{\frac{4}{3}\pi \frac{{\rho }_{m}G}{A^3}}\sqrt{\frac{M_p+M_m}{M_p}}\sqrt{\frac{4}{3}\pi \frac{{\rho }_{m}G}{A^3}.}\hfill \end{array}$$(A.1)

For A = 2.2, and for the density ρ_m of the moon given in kg m⁻³, this renders:

$${\dot{\theta }}_p<5.12\times {10}^{-6}\sqrt{{\rho }_m}{\text{s}}^{-1}=1.84\times {10}^{-2}\sqrt{{\rho }_m}{\text{hr}}^{-1},$$(A.2)

For the density of our Moon, ρ_m = 3340 kg m⁻³, we have

$${\dot{\theta }}_p<1.07{\text{hr}}^{-1}.$$(A.3)

Appendix B Quality functions. A typical shape

To get an idea of a typical shape of a quality function, consider a tidally perturbed homogeneous body whose rheology is defined by a complex compliance

$$\overline{J}\left(\chi \right)=ℛ\left[\overline{J}\left(\chi \right)\right]+iℐm\left[\overline{J}\left(\chi \right)\right],$$(B.1)

with

$$\chi =\chi lmpq$$(B.2)

being a concise notation for a tidal frequency ${\chi }_{lmpq}=\left|{\omega }_{lmpq}\right|$ exerted by tides.

For example, the Maxwell rheology looks like

$${}^{\left(Maxwell\right)}\overline{J}\left(\chi \right)=J+i\left(-\frac{1}{\eta \chi }\right)=J\left(1-\frac{1}{{\tau }_M\chi }\right),$$(B.3)

where η is the shear viscosity, J is the unrelaxed sheer compliance (the inverse of the unrelaxed shear rigidity µ), and τ_M= η J = η/µ is the Maxwell time. In the more realistic Andrade model, one has to include terms responsible for transient processes. The elastic, instantaneous, part of the complex compliance will, though, still be J. Neither the Maxwell nor the Andrade model are capable of taking the relaxation of the elastic part into consideration. (This relaxation is, however, an actual effect taking place mainly due to grain-boundary sliding.) So these models fail to distinguish between the unrelaxed and relaxed values of J. To fix this flaw, a more detailed model by Sundberg and Cooper should be employed. A simple review on the topic is available in Section 3 of Walterova et al. (2023) and in Section 2 of Bagheri et al. (2022).

As demonstrated, for example, in Efroimsky (2015), Eq. (40), the quality function of a homogeneous body is¹⁵

$$\begin{array}{c}{K}_l\left(\chi \right)\equiv k_l\left(\chi \right){ϵ}_l\left(\chi \right)\\ =-\frac{3}{2\left(l-1\right)}\frac{{ℬ}_lℐm\left[\overline{J}\left(\chi \right)\right]}{{\left(ℛe\left[\overline{J}\left(\chi \right)J^{-1}\right]+{ℬ}_l\right)}^2+{\left(ℐm\left[\overline{J}\left(\chi \right)\right]\right)}^2}\end{array}$$(B.4a)

$$=-\frac{3}{2\left(l-1\right)}\frac{{\mathcal{𝒜}}_lℐm\left[\overline{J}\left(\chi \right)J^{-1}\right]}{{\left(ℛe\left[\overline{J}\left(\chi \right)J^{-1}\right]+{\mathcal{𝒜}}_l\right)}^2+{\left(ℐm\left[\overline{J}\left(\chi \right)J^{-1}\right]\right)}^2},$$(B.4b)

with the dimensional constants ℬ_l; and dimensionless constants 𝒜_l defined as

$${ℬ}_l\equiv \frac{\left(2l^2+4l+3\right)}{\lg \rho R}=\frac{3\left(2l^2+4l+3\right)}{4l\pi G{\rho }^2{R}^2},$$(B.5)

$${\mathcal{𝒜}}_l\equiv {ℬ}_l{J}^{-1}={ℬ}_l\mu ,$$(B.6)

G being Newton’s gravitational constant, and ɡ, ρ, R being the surface gravity, density, and radius of the body. Be mindful that expression (20) for the secular part of the polar tidal torque contains not K_l(χ_lmpq) but the odd function Kl(ω_lmpq)=Kl(χ_lmpq) Sign(ω_lmpq).

For the Maxwell or Andrade model, the function K_l(ω_lmpq) has the shape of a kink, shown in Figure 3.¹⁶ Specifically, for the Maxwell rheology (B.3), the quality functions become

$$K_l\left(\omega \right)=\frac{3}{2\left(l-1\right)}\frac{{\tau }_M\omega {\mathcal{𝒜}}_l}{1+{\left({\tau }_M\omega \right)}^2\left(1+\right){\mathcal{𝒜}}_1^2},$$(B.7)

ω being a short notation for ω_lmpq.

From the above expression, we observe that for a Maxwell body – and, in realistic situations, for Andrade, Burgers, and Sundberg-Cooper bodies also¹⁷ –the extrema of the kink K_l(ω) are residing at

$${\omega }_{{\text{peak}}_l}=\pm \frac{{\tau }_M^{-1}}{1+{𝒜}_l},$$(B.8)

the corresponding peak values being

$$K_l^{\left(\text{peak}\right)}=\pm \frac{3}{4\left(l-1\right)}\frac{{𝒜}_l}{{𝒜}_l}.$$(B.9)

In the inter-peak interval, a quality function K_l(ω) is almost linear: ¹⁸

$$\begin{array}{c}\left|\omega \right|<\left|{\omega }_{peak}\right|\Rightarrow \\ K_l\left(\omega \right)\simeq \frac{2}{2\left(l-1\right)}\frac{{𝒜}_l}{1+{𝒜}_l}\frac{\omega }{\left|{\omega }_{peak}\right|}.\end{array}$$(B.10)

Outside this interval, it falls off as the inverse ω:

$$\begin{array}{c}\left|\omega \right|>\left|{\omega }_{pea{k}_l}\right|\Rightarrow \\ K_l\left(\omega \right)\simeq \frac{3}{2\left(l-1\right)}\frac{{𝒜}_l}{1+{𝒜}_l}\frac{\left|{\omega }_{peak}\right|}{\omega }.\end{array}$$(B.11)

From Eq. (B.9), we observe that the peaks’ amplitude is independent of the viscosity value η. On the other hand, the spread between the extrema is inversely proportional to η, as follows from expression (B.8). For solid terrestrial planets, the mean viscosity values are high, wherefore the extrema are located close to zero (i.e., to the resonance ω = ω_nmpq=0), so the inter-peak interval of linear frequency-dependence of K_l is very narrow and the peaks are sharp.

For a Maxwell or Andrade body, both the kink function K_l(ω) and the similar to it function sin ${ϵ}_l\left(\omega \right)=Q_l^{-1}\left(\omega \right)$Sign( ω) have only one peak for a positive ω. The situation will, however, change if we plug into formula (B.4), and into its counterpart for sin ${ϵ}_n\left(\omega \right)=Q_n^{-1}\left(\omega \right)$ Sign( ω), a complex compliance appropriate to a more accurate description, such as the Sundberg-Cooper rheology. In that case, an additional, smaller local maximum will appear on the right slope (an a similar local minimum will emerge on the left slope). These additional peaks may have big geophysical consequences for the Moon (Walterova et al. 2023).

Appendix C Secular components of the solar and lunar tidal torques acting on the planet

The quadrupole part of the torque looks like (Efroimsky 2012):

$$\begin{array}{l}{\mathcal{T}}_{\left(lmp\right)=\left(220\right)}=\hfill \\ \frac{3}{2}\frac{GM{\prime }^2}{a}{\left(\frac{R_p}{a}\right)}^5[\frac{1}{2304}{e}^6{K}_2\left(-n-2{\dot{\theta }}_p\right)\hfill \\ +\left(\frac{1}{4}{e}^2-\frac{1}{16}{e}^4+\frac{13}{768}{e}^6\right)K_2\left(n-2{\dot{\theta }}_p\right)\hfill \\ +\left(1-5e^2+\frac{63}{8}{e}^4-\frac{155}{36}\right)K_2\left(2n-2{\dot{\theta }}_p\right)\hfill \\ +\left(\frac{49}{4}{e}^2-\frac{861}{16}{e}^4+\frac{21975}{256}{e}^6\right)K_2\left(3n-2{\dot{\theta }}_p\right)\hfill \\ +\left(\frac{289}{4}{e}^4-\frac{1955}{6}{e}^6\right)K_2\left(2n-{\dot{\theta }}_p\right)\hfill \\ +\frac{714025}{2304}{e}^4{K}_2\left(5n-3{\dot{\theta }}_p\right)\hfill \\ +O\left(e^8ϵ\right)+O\left(i^2ϵ\right),\hfill \end{array}$$(C.1)

where M′ is the mass of the perturber, ϵ is a typical value of the phase lag in the perturbed planet, and the perturber’s orbital inclination i on the planet’s equator is small enough for the order-i² terms to be dropped.

This formula will render the solar torque ${\mathcal{T}}^{(*)}$, provided we identify the mass of the perturber with that of the star: M′ = M_*, and employ the elements a = a_p, e = e_p, and use the mean motion n = n_p of the planet’s orbit.

The formula will give us the lunar torque ${\mathcal{T}}^{\left(m\right)}$, if we set M′ = M_m, a = a_m, e = e_m, and n = n_m.

Appendix D Derivation of Eq. (50)

For the planet-moon two-body system, the angular momentumis given by

$$\begin{array}{c}H=\frac{M_p{M}_m}{M_p+M_m}\sqrt{G\left(M_p+M_m\right)a_m\left(1-e_m^2\right)}\\ +C_p{\dot{\theta }}_p+C_m{\dot{\theta }}_m,\end{array}$$(D.1)

C_p and C_m being the maximal moments of inertia of the planet and the moon, correspondingly, and ${\dot{\theta }}_{\text{p}}$ and ${\dot{\theta }}_{\text{m}}$ being their rotation rates. If we neglect the eccentricity and assume that the spin angular momentum of the moon is negligible as compared to that of the planet, the above expression will get simplified to

$$H=\frac{M_p{M}_m}{M_p+M_m}\sqrt{G\left(M_p+M_m\right)a_m+C_p{\dot{\theta }}_p,}$$(D.2)

differentiation whereof will give us:

$$0=\frac{M_p{M}_m}{M_p+M_m}\sqrt{G\left(M_p+M_m\right)}\frac{d}{dt}{a}_m^{1/2}+C_p{\ddot{\theta }}_p.$$(D.3)

Casting the definition $n_{\text{m}}=\sqrt{G\left(M_p+M_m\right)a_m^{-3}}$ into the form of $a_m^{1/2}={\left(G\left(M_p+M_m\right)\right)}^{1/6}{n}_m^{-1/3}$, we derive:

$$\begin{array}{c}\frac{d}{dt}{a}_m^{1/2}=-\frac{1}{3}{\left(G\left(M_p+M_m\right)\right)}^{1/6}{n}_m^{-4/3}{\dot{n}}_m\\ =-\frac{1}{3}{\left(G\left(M_p+M_m\right)\right)}^{-1/2}{a}_m^2{\dot{n}}_m,\end{array}$$(D.4)

insertion whereof into Eq. (D.3) entails

$$0=-\frac{1}{3}\frac{M_p{M}_m}{M_p+M_m}{a}_m^2{\dot{n}}_m+C_p{\ddot{\theta }}_p.$$(D.5)

Through a dimensionless coefficient ξ, the maximal moment of inertia is conventionally expressed as $C_p=\xi M_p{R}_p^2$. Also, since en route from Eq. (D.1) to Eq. (D.2) we set $\left|C_m{\dot{\theta }}_{\text{m}}\right|\ll \left|C_p{\dot{\theta }}_p\right|$, then in about the same approximation we may put M_p M_m/(M_p+ M_m) ≈ M_m into Eq. (D.5). This gives us

$$0=-\frac{1}{3}{M}_m{a}_m^2{\dot{n}}_m+\xi M_p{R}_p^2{\ddot{\theta }}_p.$$(D.6)

References

All Tables

All Figures

	Fig. 1 Reduced Hill radii of 580 known inner exoplanets versus their radii in units of Jupiter’s radius. The diagonal line indicates the equality of the two radii.
In the text

	Fig. 2 Reduced Hill radii of 580 known inner exoplanets versus their estimated Roche limits in units of Jupiter’s radius. The dashed diagonal line marks the locus of equal radii where the dynamical niche for satellites vanishes. The planets that fall below this line cannot have any stable satellites.
In the text

	Fig. 3 Typical shape of the quality function K1(ω)=k1(ω) sin ϵ1(ω), where ω is a shortened notation for the tidal Fourier mode ωlmpq. (Reprinted from Icarus, Vol. 241, Noyelles et al., Spin-orbit evolution of Mercury revisited, pages 26–44, Copyright (2014), with permission from Elsevier.)
In the text

Fig. 4 Schematic depiction of the three torques involved in the quasi-stable pseudo-synchronous equilibrium of the planet-moon system. The tidal torque wherewith the moon is acting on the planet (solid blue line) is coupled with an opposite torque wherewith the planet is acting on the orbit of the moon (red curve). The negative tidal torque on the planet from the star (dashed blue line) is weakly dependent on the relative frequency. The vertical dotted line indicates the location where the tidal torques acting on the planet are equalised.

In the text

Fig. 5 Reduced Hill radii of 580 known inner exoplanets versus their radii in units of Jupiter’s radius. Dashed horizontal lines mark the lower boundaries of plot areas where these exoplanets can in principle have moons with the smallest masses as labelled in the graph. For each specific |exoplanet shown with a dot in this plot, the mass of a survivable moon should be above the dashed line (above the upper line, for putative moons of Mm = 0.01 Mp, or above the lower line, for putative moons of Mm = 0.05 Mp).

In the text

Fig. 6 Lower masses of exomoons that can in principle synchronise their host planets for 580 known exoplanet systems, in units of the mass of the Moon. The M− limits are derived from the estimated reduced Hill radii and refer to the Scenario 1 pathway, when the exomoon is initially above the synchronous radius. The M+ limits are derived from the estimated Roche radii and refer to the Scenario 2 pathway, when the exomoon is initially below the synchronous radius.

In the text