# Oblique rings from migrating exomoons: A possible origin for long-period exoplanets with enlarged radii

Melaine Saillenfest, Sophia Sulis, Paul Charpentier, Alexandre Santerne

► To cite this version:

Melaine Saillenfest, Sophia Sulis, Paul Charpentier, Alexandre Santerne. Oblique rings from migrating exomoons: A possible origin for long-period exoplanets with enlarged radii. Astronomy and Astrophysics - A&A, 2023, 675, pp.A174. 10.1051/0004-6361/202346745. hal-04128031

HAL Id: hal-04128031

https://hal.science/hal-04128031

Submitted on 25 Aug 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Oblique rings from migrating exomoons: A possible origin for long-period exoplanets with enlarged radii

Melaine Saillenfest<fnref n="1" />, Sophia Sulis<fnref n="2" />, Paul Charpentier<fnref n="3" />, and Alexandre Santerne<fnref n="2" />

Received 26 April 2023 / Accepted 7 June 2023

## ABSTRACT

*Context.* The extremely low density of several long-period exoplanets in mature systems is still unexplained – with HIP 41378 f being archetypical of this category. It has been proposed that such planets could actually have normal densities but be surrounded by a ring observed approximately face on, mimicking the transit depth of a puffy planet. This configuration would imply that the equator of the planet is nearly perpendicular to its orbit plane, which is at odds with the formation process of gas giants. Yet, in the context of the Solar System planets, it has recently been shown that after gigayears of evolution, the tidal migration of a moon can naturally lead to a very tilted planet with a ring.

*Aims.* As exomoons are expected to be ubiquitous around giant exoplanets, this mechanism may be responsible for the anomalous radii of some observed exoplanets. In preparation for the future discoveries of the PLATO mission, we present a simple method for checking the plausibility of this mechanism for a given exoplanet.

*Methods.* Analytical formulas give the probability density function of the relevant precession harmonics of the planet. For each harmonic, simple criteria set the moon mass and other properties required for the mechanism to operate.

*Results.* We applied this methodology to HIP 41378 f, and we show that in order to reproduce the observed configuration, a hypothetical former moon should have had a moon-to-planet mass ratio of a few times $10^{-4}$ (i.e. roughly the mass of our Moon) and have migrated over a distance of a few planet's radii on a gigayear timescale. These orders of magnitude match the properties of moons expected to exist around gaseous exoplanets.

*Conclusions.* We conclude that the migration of a former moon is a viable formation pathway for the proposed ring and tilt of HIP 41378 f. This example strengthens the ring hypothesis and motivates its application to other promising targets.

**Key words.** planets and satellites: dynamical evolution and stability – planets and satellites: rings – planets and satellites: individual: HIP 41378 f

## 1. Introduction

The so-called super-puff exoplanets have moderate masses (typically $\le 15 M_\oplus$) but surprisingly large radii ($\gtrsim 6 R_\oplus$), giving them extremely low bulk densities ($\le 0.3 \text{ g cm}^{-3}$; see e.g. Lee & Chiang 2016). Although relatively rare, super-puffs form a growing class of exoplanets. Among the puffiest exoplanets with the longest orbital periods, we can cite the iconic HIP 41378 f, Kepler-87 c, Kepler-79 d, Kepler-177 c, and Kepler-51 b, c, and d. Super-puffs must be distinguished from inflated hot Jupiters, which show a correlation between stellar irradiation and radius inflation (see e.g. Laughlin et al. 2011; Lopez & Fortney 2016). This correlation indicates that hot Jupiters have extended atmospheres connected in some way to their close proximity to the star (see e.g. Burrows et al. 2000; Chabrier & Baraffe 2007; Batygin et al. 2011; Grunblatt et al. 2016). A similar conclusion can be reached for short-period sub-Neptunes (Pu & Valencia 2017; Millholland 2019), but not for distant super-puffs, because they have much cooler equilibrium temperatures and undergo negligible star-planet tidal dissipation.

Initiated by the preprint of Santerne et al. (2019), the low density of exoplanet HIP 41378 f, in particular, immediately raised much discussion. HIP 41378 f is mature ($2.1^{+0.4}_{-0.3}$ Gyr; Lund et al. 2019) and has a long period (542 days) and low equilibrium

temperature (300 K). Its low density ($0.09 \pm 0.02 \text{ g cm}^{-3}$) puts this planet among the puffiest exoplanets known to date. Even though other super-puffs are known, most of them are likely young and/or have shorter periods (Lee & Chiang 2016). Instead of a radius inflation, Akinsanmi et al. (2020) propose that HIP 41378 f could be a standard Neptune-sized planet surrounded by an inclined opaque ring that would mimic the transit depth of an inflated planet. As no significant distortion is visible in the transit ingress and egress of HIP 41378 f, the hypothetical ring should be optically thick and seen roughly face on. This configuration would imply that the obliquity of the planet<fnref n="1" /> is nearly 90°.

The ring hypothesis was investigated by Piro & Vissapragada (2020) for other super-puff exoplanets. Good candidates are Kepler-87 c, Kepler-79 d, and Kepler-177 c, even though their moderate temperatures – as that of HIP 41378 f – do not allow for water ice to exist around them. Therefore, unlike Saturn's ring, their rings would need to be composed of porous

<footnote n="1">¹ IMCCE, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Université, Université de Lille, 75014 Paris, France
e-mail: melaine.saillenfest@obspm.fr</footnote>

<footnote n="2">² Université Aix-Marseille, CNRS, CNES, LAM, Marseille, France</footnote>

<footnote n="3">³ Université de Toulouse, CNRS, IRAP, Toulouse, France</footnote>

<footnote n="1">¹ Not to be confused with the stellar obliquity (i.e. the angle between the spin axis of the star and the orbit pole of a given planet). Throughout this article, the term obliquity is exclusively used for the planetary obliquity (i.e. the angle between the spin axis of the planet and its orbit pole).</footnote>
rocky particles. According to the results of Piro & Vissapragada (2020), HIP 41378 f is currently the best candidate for a ring. Its long period would protect a ring against destructive irradiation levels and a strong warp due to the stellar torque; it also results in negligible star-planet tidal dissipation, which means that no particular mechanism would be required for the planet to maintain a large obliquity<fnref n="2" />. The low eccentricity of HIP 41378 f also guarantees a small level of orbital perturbations for the ring particles.

In order to determine the planets' atmospheric properties and test the ring hypothesis, near-infrared transmission spectra have been acquired for Kepler-51 b and d (Libby-Roberts et al. 2020), Kepler-79 d (Chachan et al. 2020), and HIP 41378 f (Alam et al. 2022). These spectra ended up being featureless, ruling out clear, low-metallicity atmospheres. The ring hypothesis is therefore not contradicted for these planets, but flat spectra can also be produced by high-altitude hazes or high-metallicity atmospheres. In fact, convincing atmospheric models have been put forward for Kepler-51 b and d, as well as Kepler-79 d (see also Wang & Dai 2019; Gao & Zhang 2020; Ohno & Tanaka 2021). Interestingly, these models of extended atmospheres appear to be inapplicable to HIP 41378 f as it is too massive ($M = 12 \pm 3 M_\oplus$), too cold, and too old.

The question of the possible physical composition of HIP 41378 f was explicitly tackled by Belkovski et al. (2022). The authors show that photoevaporation is not nearly enough to explain the extreme density disparity between planet f and other planets in the system. Moreover, the observed mass and radius of HIP 41378 f would require an envelope-to-core mass fraction larger than 75% together with a high entropy (e.g. produced by recent collisions). Such a massive envelope is unlikely from the perspective of planetary formation, as it would require runaway gas accretion to have started precisely during the dissipation of the gas disc, and planet HIP 41378 f may not be massive enough anyway to have triggered runaway accretion.

Hence, the ring hypothesis appears to be favoured for HIP 41378 f, and it may apply as well to a restricted number of other observed super-puffs. Tidal rings are confined below the Roche limit, very close to their host planets. As such, they are strongly coupled to the centrifugal bulge of the planets, and they directly materialise their equatorial planes. In order to produce a substantial increase in a planet's transit depth (i.e. a very noticeable super-puff), its ring must be oriented roughly in the sky plane. This means that the planet's spin axis must point roughly along the observer's direction; its obliquity is therefore $\epsilon \approx 90^\circ$ as proposed by Akinsanmi et al. (2020). Such an exotic configuration may seem questionable from a formation point of view. Because of the angular momentum acquired during gas accretion, gaseous planets are expected to form with low obliquities. The obliquities of the Solar System giant planets are therefore interpreted as strong tracers of their dynamical evolution, and much effort is put into understanding their origin (see e.g. Tremaine 1991; Ward & Hamilton 2004; Hamilton & Ward 2004; Boué et al. 2009; Boué & Laskar 2010; Morbidelli et al. 2012; Vokrouhlický & Nesvorný 2015; Rogoszinski & Hamilton 2020, 2021; Saillenfest et al. 2020, 2021a, 2022; Salmon & Canup 2022; Rufu & Canup 2022; Wisdom et al. 2022). In this context, the ring hypothesis for super-puffs would greatly benefit

from an underlying mechanism that may be responsible for their unusual configuration. The existence of such a mechanism would not certify whether a given planet does possess a ring or not, but it would show whether known dynamical processes are able to (or are even likely to) produce the proposed configuration.

In the Solar System, a substantial tidal migration of moons has recently been observed to be at play around gaseous planets (see Lainey et al. 2009, 2017, 2020; Jacobson 2022) – even though it involves mechanisms of energy dissipation that are vastly different from those responsible for the well known rapid migration of our Moon (see e.g. Farhat et al. 2022). These results have strong implications for the orbital dynamics of moons around gaseous planets, but also for the gigayear-timescale dynamics of planetary spin axes. Indeed, moons affect the spin-axis precession rate of planets in a way that is intimately related to their distance (see e.g. Boué & Laskar 2006). The migration of a moon is therefore accompanied with a variation in the planet's spin-axis precession rate. In turn, this variation can drive the planet into a so-called secular spin-orbit resonance, that is, a resonance between the planet's spin-axis precession and one harmonic of its orbital nodal precession. As a matter of fact, this kind of resonances abound in multi-planetary systems. Provided that a planet has a substantially massive migrating moon, it may therefore be guaranteed to encounter one of these resonances sooner or later during its evolution. Once captured in resonance, the still ongoing migration of the moon produces a gradual tilting of the planet's spin axis (unless, as for the Earth, resonances are so numerous that they overlap massively; see Laskar & Robutel 1993; Néron de Surgy & Laskar 1997). This phenomenon is probably responsible for the $27^\circ$ obliquity of Saturn (Saillenfest et al. 2021a,b; Wisdom et al. 2022), and it is predicted to happen to Jupiter in the future (Saillenfest et al. 2020). It may also have played a role in the tilting of Uranus (Saillenfest et al. 2022).

When the planet's obliquity reaches $\epsilon \ge 70^\circ$, however, regular moons are known to be unstable in some range of distance (Tremaine et al. 2009). Interestingly, the migration of a single moon makes the system converge to this unstable zone, putting a dramatic end to the tilting process (see Saillenfest & Lari 2021; Saillenfest et al. 2022). At this point, the moon may be ejected or be destructed below the planet's Roche limit, eventually forming a tidal disc of debris. In the latter case, the final state of the system is a ringed planet with very high obliquity. This final state recalls the exotic configuration proposed for super-puff exoplanets. It would therefore be valuable to determine whether this mechanism could apply to them and provide a plausible dynamical background to the ring hypothesis.

In this article, we aim to present a generic methodology to assess whether the migrating-moon mechanism can realistically produce a tilted ring around a given exoplanet. Even though the number of known distant super-puffs is small today, the future PLATO mission (Rauer et al. 2014, 2016) will considerably increase our knowledge of the population of long-period exoplanets – including their masses through an intensive radial-velocity follow-up. In this context, we need efficient methods for a routine characterisation of the newly discovered planets and identification of the most interesting targets for follow up. For this reason, we design our methodology to be applicable even if the minimum amount of information about the planetary system is available (masses, periods, and sky-plane inclinations).

The article is organised as follows. In Sect. 2, we recall the basics of the tilting mechanism. In Sect. 3, we compute the probability density function of the dominant orbital precession frequencies of a planet, and we present an example of application to the super-puff exoplanet HIP 41378 f. From these

<footnote n="2">² High-obliquity equilibrium states also exist for short-period planets (Millholland & Laughlin 2019; Millholland & Spalding 2020); however, because of tidal despinning and obliquity damping, their obliquity needs to be continuously forced through dynamical interactions involving several planets (see also Su & Lai 2022a,b).</footnote>
results, we estimate in Sect. 4 the mass and migration rate that a moon around this planet would need in order to trigger the full tilting mechanism. In Sect. 5, we check that the resonance is large enough to enable an adiabatic capture and tilting, and we illustrate this mechanism with numerical simulations. We then discuss our results in Sect. 6 and conclude in Sect. 7.

## 2. Basic mechanism

As shown by Saillenfest & Lari (2021), the tilting of a planet from a low obliquity $\epsilon$ up to $\epsilon \approx 90^\circ$ can be achieved on a gigayear timescale via the tidal migration of a moon. This process occurs through the adiabatic drift of the system along the centre of a secular spin-orbit resonance. In this section, we recall the physical quantities involved and the conditions required to trigger this process.

We write $I$ the orbital inclination of the planet and $\Omega$ its longitude of ascending node. We decompose the inclination dynamics of the planet in a quasi-periodic series truncated to $N$ terms:

$$ \zeta = \sin \frac{I}{2} \exp(i\Omega) = \sum_{j=1}^{N} S_j \exp[i\phi_j(t)], \quad (1) $$

where $S_j$ is a positive real constant, and $\phi_j(t) = v_j t + \phi_j^{(0)}$ evolves linearly over time $t$ with frequency $v_j$. Resonance capture from a low obliquity is possible only for resonances with a harmonic having a negative frequency $v_j$ such that $|v_j| \ge p$, where

$$ p = \frac{3}{2} \frac{G M_{\star}}{a^3 (1 - e^2)^{3/2}} \frac{J_2}{\omega \lambda} \quad (2) $$

is the characteristic spin-axis precession rate of the planet. In this expression, $G$ is the gravitational constant, $M_{\star}$ is the mass of the star, $a$ and $e$ are the semi-major axis and eccentricity of the planet on its orbit around the star, $J_2$ is the second zonal gravity coefficient of the planet, $\omega$ is its spin rate, and $\lambda$ is its normalised polar moment of inertia. The parameters $J_2$ and $\lambda$ must be defined through the same normalising radius $R$ (which is generally chosen as the equatorial radius of the planet).

The influence of a regular moon on the long-term spin-axis dynamics of the planet can be quantified by its non-dimensional ‘mass parameter’ $\eta$ defined by

$$ \eta = \frac{1}{2} \frac{m}{M} \frac{r_{\mathrm{M}}^2}{J_2 R^2}, \quad (3) $$

where $m$ is the mass of the moon, $M$ is the mass of the planet, and $r_{\mathrm{M}}$ is the following characteristic length<fnref n="3" />:

$$ r_{\mathrm{M}}^{5} = 2 \frac{M}{M_{\star}} J_{2} R^{2} a^{3} (1 - e^{2})^{3/2}. \quad (4) $$

Under the hypothesis that the moon’s mass ratio $m/M$ is small (which does not necessarily imply that $\eta$ is small), Saillenfest & Lari (2021) show that all resonances with a nodal harmonic having a negative frequency $v_j$ verifying

$$ p \le |v_j| \le p \frac{\eta}{2} \quad (5) $$

<figure><img src="image_2.png" /><figcaption>Fig. 1. Level curves of the planet’s spin-axis precession rate (adapted from Saillenfest et al. 2022 for a moon mass-parameter $\eta = 20$). If the planet is trapped in a secular spin-orbit resonance, the system evolves along one of these curves as the moon migrates. The condition in Eq. (5) corresponds to the pink region; the left and right inequalities are the green and red curves, respectively. In the blue area, the moon is unstable.</figcaption></figure>

can allow the planet’s obliquity to grow from $\epsilon = 0^\circ$ to $\epsilon = 90^\circ$. This condition is illustrated in Fig. 1. Knowing the harmonics $v_j$ of the planet’s orbital precession, Eq. (5) allows one to compute the minimum mass required for the moon to produce the tilting. As resonances converge to an unstable region, the moon is ultimately lost at the end of the tilting process (see Fig. 1).

When Eq. (5) is verified, the adiabatic capture and tilting of the planet within a given resonance requires an adequate hierarchy of timescales. First, we introduce the timescale $\tau$ of secular oscillations of the moon around its equilibrium ‘Laplace plane’ (see Tremaine et al. 2009) as $\tau = 2\pi/\kappa$, where

$$ \kappa^2 = \frac{9 M_{\star}}{4 M} \frac{r_{\mathrm{M}}^3}{a^3 (1 - e^2)^3} \frac{G M_{\star}}{a^3}. \quad (6) $$

An adiabatic capture in resonance requires that $\tau$ is much shorter than the spin-axis precession timescale of the planet $T = 2\pi/p$; this conditions is generally well verified in practice.

Then, a given observed planet may have been adiabatically tilted via a resonance only if the timescale $T_{\mathrm{lib}}$ of libration inside the resonance is much smaller than the age of the system. For a given secular spin-orbit resonance, the value of $T_{\mathrm{lib}}$ near the resonance centre can be computed as $T_{\mathrm{lib}} = 2\pi/\mu$, where

$$ \mu^2 = (p')^2 \left( \frac{\beta^2}{\sin^2 \epsilon_0} + \beta \sin \epsilon_0 \right). \quad (7) $$

In this expression, $\epsilon_0$ is the planet’s obliquity at the resonance centre and $p'$ is a modified version of $p$ that takes into account the presence of the planet’s moon (see Saillenfest & Lari 2021). We define the non-dimensional variables $\gamma = \rho_1/p'$ and $\beta = \rho_2/p'$, where

$$ \begin{aligned} \rho_1 &= -\left( v_k - 2 \sum_{j=1}^{N} v_j S_j^2 \right), \\ \rho_2 &= -S_k \left( 2v_k + v_k S_k^2 - 2 \sum_{j=1}^{N} v_j S_j^2 \right), \end{aligned} \quad (8) $$

and $k$ is the index in Eq. (1) of the considered resonance. $T_{\mathrm{lib}}$ depends on the distance of the moon through $p'$ and $\epsilon_0$. However, an upper bound for $T_{\mathrm{lib}}$ is obtained at the time of resonance

<footnote n="3">³ $r_{\mathrm{M}}$ is called ‘mid-point radius’ by Saillenfest & Lari (2021). It is sometimes defined as the Laplace radius in other publications, either with or without the leading factor 2.</footnote>
capture, for which $\gamma^{2/3} + \beta^{2/3} = 1$ (see Henrard & Murigande 1987; Saillenfest et al. 2019). In this case, $p'$ is equal to $p' = (\rho_1^{2/3} + \rho_2^{2/3})^{3/2}$, and the planet's obliquity at the centre of the resonance is

$$ \cos \varepsilon_0 = \gamma - \gamma^{1/3} + \sqrt{\gamma^2 + \gamma^{2/3} - \gamma^{4/3}}. \quad (9) $$

Thanks to these expressions, we can compute $T_{\text{lib}}$ from Eq. (7) as a mere function of the planet's orbital dynamics in Eq. (1).

### 3. Orbital precession modes of the planet

To apply this mechanism to a given planet, we need to know its orbital precession spectrum, which depends on planet-planet mutual interactions. However, the masses and orbital elements of exoplanets are generally not well known. For given parameters and their uncertainties, the most simple way to explore the variety of possible long-term orbital solutions is to use the Lagrange–Laplace system (see e.g. Murray & Dermott 1999).

#### 3.1. The Lagrange–Laplace proper modes

The Lagrange–Laplace system is a secular theory at second order in eccentricity and inclination. As such, it assumes that all eccentricities and inclinations are small and it neglects the long-term influence of mean-motion resonances. Small mutual inclinations are indeed strongly favoured in multi-planetary systems in which most planets are observed to transit their star. This is the case of HIP 41378, around which the transits of five planets are observed (Vanderburg et al. 2016). Eccentricities are also expected to be small in multi-planetary systems for stability reasons. Moreover, according to the statistical distribution of multi-planetary systems (Xie et al. 2016) and to theoretical arguments about chaotic diffusion (which leads to the statistical equipartition of angular momentum deficit; see Laskar & Petit 2017), planets having small mutual inclinations tend to have small eccentricities, and vice versa. Hence, the use of the Lagrange–Laplace theory is generally justified in this regard for multi-planetary systems. Neglecting the long-term effect of mean-motion resonances may seem more questionable, as many pairs of exoplanets are observed to be close to important resonances (see e.g. Fabrycky et al. 2014). Yet, the strongest mean-motion resonances in planetary systems – and those enabling smooth captures – are of eccentricity type. As such, they mainly affect eccentricities. Here, instead, we are only interested in the inclination degree of freedom of the planets because it is by far the main driver of their long-term spin-axis dynamics. The planets' eccentricity dynamics only enter into play at order three and beyond (see Saillenfest et al. 2019), so mean-motion resonances can safely be ignored in this analysis.

As above, we describe the nodal precession and inclination dynamics of a planet *k* in the planetary system through the complex variable

$$ \zeta_k = \sin \frac{I_k}{2} \exp(i \Omega_k), \quad (10) $$

where $I_k$ is the orbital inclination of planet *k* and $\Omega_k$ is its longitude of ascending node. The Lagrange–Laplace system gives the linear equation of motion

$$ \frac{d\zeta}{dt} = iB\zeta, \quad (11) $$

in which $\zeta$ is the vector containing the $\zeta_k$ variable of all planets and $B$ is a constant matrix that only depends on the masses and semi-major axes of the planets (see e.g. Laskar & Robutel 1995). The solution of this equation for a given planet *k* has the form of a quasi-periodic series as in Eq. (1):

$$ \zeta_k(t) = \sum_{j=1}^{N_p} S_j \exp \left[ i \left( \nu_j t + \phi_j^{(0)} \right) \right], \quad (12) $$

where the number of terms *N* is equal to the number of planets $N_p$ in the system. Equation (12) is a linear combination of proper modes whose frequencies $\nu_j$ are the eigenvalues of the matrix $B$. As $B$ only depends on the masses and semi-major axes of the planets, this is also the case of the frequencies $\nu_j$. Because of the conservation of total angular momentum, one of the frequencies $\nu_j$ is identically equal to zero; the related constant term in Eq. (12) gives the orientation of the system's invariant plane.

Thanks to the fast computation of the solution of the Lagrange–Laplace system (which amounts to a mere matrix inversion), millions of trials can be performed at virtually no cost. In order to explore the distribution of possible values for the frequencies $\nu_j$, the first step is to draw the masses and semi-major axes of the $N_p$ planets from their respective statistical distributions – which represent our knowledge of their values. A similar approach was followed by Becker & Adams (2016) in their study of the compact multi-planetary systems observed by Kepler. Each sequence of masses and semi-major axes for the $N_p$ planets represents a possible realisation of the planetary system. In case the mass of a given planet has not been measured, a broad distribution of mass can be adopted (e.g. a uniform distribution in a given interval, or a law drawn for an assumed mass-radius relationship; see below). From a large number of realisations of the planetary system, a histogram for each frequency $\nu_j$ can be computed. These histograms define the possible locations of secular spin-orbit resonances given our current knowledge of the planetary system.

In practice, the largest values of the Lagrange–Laplace matrix $B$ in Eq. (11) are often located along its diagonal (meaning that the planetary system is only weakly coupled); this implies that each planet *k* has its own dominant proper mode, which appears in Eq. (12) as the term with largest amplitude. The frequency of the dominant proper mode of planet *k* is usually noted $s_k$. In the context of the Lagrange–Laplace approximation, the quasi-periodic series in Eq. (12) contains exactly $N = N_p$ terms and the frequencies $\nu_j$ are each equal to one of the $s_k$. More generally, the orbital evolution of any planet in a stable system can be written as in Eq. (12), but where *N* tends to infinity and each harmonic $\nu_j$ is a linear combination of the fundamental frequencies of the system (see Sect. 2). The first few strongest harmonics of the series are however proper modes given by the Lagrange–Laplace approximation; hence, the analysis presented here can be thought of as the dominant component of a more general theory.

While building the histogram for each proper mode $s_k$, a complication may arise. Indeed, if the masses and semi-major axes of the planets have large uncertainties, the distributions of the various frequencies may overlap. In this case, identifying each eigenvalue $\nu_j$ of matrix $B$ as the correct proper mode $s_k$ requires some caution. As the hierarchy of proper modes depends on the planetary system considered, a specific identification process is required. As an example, we subsequently present the case of the HIP 41378 system.
<table><caption>**Table 1.** Parameters of planets in the HIP 41378 system used in this article.</caption><thead><tr><th>k</th><th>Name</th><th>M<sub>k</sub> (M<sub>⊕</sub>)</th><th>P<sub>k</sub> (day)</th><th>a<sub>k</sub> (au)</th><th>I<sub>k</sub> (°)</th></tr></thead><tbody><tr><td>1</td><td>b</td><td>6.89 ± 0.88</td><td>15.57208 ± 0.00002</td><td>0.1283 ± 0.0015</td><td>88.75 ± 0.13</td></tr><tr><td>2</td><td>c</td><td>4.4 ± 1.1</td><td>31.70603 ± 0.00006</td><td>0.2061 ± 0.0024</td><td>88.48 ± 0.07</td></tr><tr><td>3</td><td>g</td><td>7.0 ± 1.5</td><td>62.06 ± 0.32</td><td>0.3227 ± 0.0036</td><td>88</td></tr><tr><td>4</td><td>d</td><td>12.7 ± 6.0</td><td>278.3618 ± 0.0005</td><td>0.88 ± 0.01</td><td>89.80 ± 0.02</td></tr><tr><td>5</td><td>e</td><td>12 ± 5</td><td>369 ± 10</td><td>1.06 ± 0.03</td><td>89.84 ± 0.07</td></tr><tr><td>6</td><td>f</td><td>12 ± 3</td><td>542.0798 ± 0.0002</td><td>1.37 ± 0.02</td><td>89.97 ± 0.01</td></tr></tbody></table>

**Notes.** Parameters are the mass $M_k$, period $P_k$, semi-major axis $a_k$, and inclination $I_k$ with respect to the sky plane. By convention, inclinations are given with values $I_k \le 90^\circ$, but transit and radial velocity data cannot discriminate between an inclination value $I_k$ or $180^\circ - I_k$. Uncertainties are assumed to be Gaussian; quoted intervals are $1\sigma$. Parameters come from the preprint by Santerne et al. (2019) except the mass of planet d (see text). The star mass is taken to be $1.16 \pm 0.04 M_\odot$ for consistency with the other parameters of Santerne et al. (2019). The inclination of the non-transiting planet g is assumed to be $88^\circ$ as in Santerne et al. (2019).

## 3.2. Application to the HIP 41378 system

HIP 41378 is bright F-type star<fnref n="4" /> which harbours at least five planets called b, c, d, e, and f (Vanderburg et al. 2016). Dynamical analysis reveals that planets b and c are slightly off the 2:1 mean-motion resonance, similarly to many *Kepler* planets (see e.g. Fabrycky et al. 2014). A tentative detection of a sixth planet – planet g – is reported in the preprint of Santerne et al. (2019) close to the 2:1 mean-motion resonance with planet c. As of today, only planets b, c, and f have been observed during successive transits (see Vanderburg et al. 2016; Becker et al. 2019; Bryant et al. 2021; Alam et al. 2022) and unambiguously detected in radial velocity (Santerne et al. 2019). Therefore, only planets b, c, and f have secured periods and masses.

Two transits of planet d have been observed by the K2 mission but they are separated by a three-year observation gap, leading to a discrete set of possible periods (Becker et al. 2019). From stability considerations, and thanks to additional observations by TESS, this discrete set is further reduced to only two likely values (278 and 371 days; see Berardo et al. 2019; Lund et al. 2019; Grouffal et al. 2022). In contrast, only one transit of planet e has been observed so far, so its period suffers from large uncertainties. The best period estimate for planet e is $260^{+160}_{-60}$ days (Lund et al. 2019). The period of 369 ± 10 days obtained by Santerne et al. (2019) is compatible with this estimate, and it results in a mass of $12 \pm 5 M_\oplus$ for planet e. The mass of planet d, however, is unknown.

As explained in Sect. 1, planet HIP 41378 f is a paradigmatic case of distant super-puff. Its period is about 542 days, and it has a radius of $9.2 \pm 0.1 R_\oplus$ and mass $12 \pm 3 M_\oplus$, giving it a bulk density of $0.09 \pm 0.02 \text{ g cm}^{-3}$ (Santerne et al. 2019). Under the ring hypothesis, current data suggests a planet with radius $R = 3.7 \pm 0.3 R_\oplus$ surrounded by a ring with radius $2.6 \pm 0.2 R$ and inclination $25 \pm 4^\circ$ from the sky plane (Akinsanmi et al. 2020). This new planetary radius yields a bulk planet density of $1.2 \pm 0.4 \text{ g cm}^{-3}$, similar to that of Uranus. The hypothetical equatorial ring provides an indirect measure of the obliquity of the planet, namely<fnref n="5" /> $\varepsilon = 92 \pm 7^\circ$.

In order to compute the orbital precession modes of HIP 41378 f, our choice of prior for the masses and semi-major

axes of the planets must reflect our partial knowledge of the HIP 41378 system. We sort the planets by increasing orbital periods such that the indexes $k = (1, 2, 3, 4, 5, 6)$ correspond to the planets (b, c, g, d, e, f). We assume all masses and semi-major axes to have Gaussian distributions centred on the best-fit values of Santerne et al. (2019) given in Table 1. Planet d needs a specific treatment: even though its period has tentatively been confirmed by Grouffal et al. (2022), it has still not been detected by the radial velocity method, so its mass is highly uncertain. We choose to remain as agnostic as possible as regards its mass, and draw it from a Gaussian fit to the mass-radius distribution of all known exoplanets having a radius between 3 and $4 R_\oplus$ and a mass measurement. From the NASA Exoplanet Archive<fnref n="6" /> on date 2022-11-23, we obtain a central mass value of $12.7 M_\oplus$ and a standard deviation of $6.0 M_\oplus$. The high tail of this distribution may not be compatible with radial velocity measurements; yet, this broad interval gives us confidence that the actual mass of planet d is contained in our analysis. The low tail of the distribution (from which we cut the portion $< 0.1 M_\oplus$) corresponds to cases in which planet d barely exists at all. The system may also contain additional massive planets that have not been discovered yet. Hence, we stress that the analysis below represents our current knowledge of the system and it may need to be revisited in the future.

For the HIP 41378 system as considered in Table 1, a look at the diagonal and off-diagonal values in the Lagrange–Laplace matrix B reveals a peculiar hierarchical configuration. The system is composed of two weakly coupled subsystems: i) the inner subsystem (planets 1–2–3) is characterised by planets 1 and 3 interacting with each other and affecting the motion of the low-mass planet 2; and ii) the outer subsystem (planets 4–5–6) is made of the two strongly coupled planets 4 and 5, interacting as a whole with planet 6.

This peculiar hierarchy can be visualised by solving the Lagrange–Laplace system a first time using reasonable values for the parameters. The exact values of the parameters do not matter for now; this first step only serves as a guide to identify the frequencies and choose an adequate naming convention. Figure 2 shows an example obtained from the nominal masses and semi-major axes of the planets. We name the proper frequencies according to their qualitative role in the dynamics: $s_1$ is the precession frequency of planets 1 and 3 about their total angular momentum vector; $s_2$ is the precession frequency of the low-mass planet 2 under the action of planets 1 and 3; $s_3$ is the slow

<footnote n="4">⁴ Also known as K2-93 and EPIC 211311380.</footnote>

<footnote n="5">⁵ The ring obtained by Akinsanmi et al. (2020) is inclined by $i_r = 25 \pm 4^\circ$ from the sky plane and rotated by $\theta = 95 \pm 17^\circ$ from the transit direction. The spin-orbit obliquity $\varepsilon$ of the planet is given by $\cos \varepsilon = \cos I \cos i_r + \sin I \sin i_r \cos \theta$, where $I = 89.97 \pm 0.01^\circ$ is the orbital inclination of HIP 41378 f.</footnote>

<footnote n="6">⁶ https://exoplanetarchive.ipac.caltech.edu</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Example of inclination evolution of the six planets in the HIP 41378 system. In this example, the Lagrange–Laplace equation is solved using the nominal masses and semi-major axes of all planets given in Table 1. The initial conditions $\zeta_k$ are set from the nominal inclinations $I_k$ (all assumed to be $I_k \le 90^\circ$) and random longitudes of node $\Omega_k$ in a 0.2°-wide interval; this choice is commented in Sect. 5.</figcaption></figure>

rigid precession of the inner and outer subsystems (planets 1–2–3 and 4–5–6); $s_4$ is the precession frequency of planets 4 and 5 about their total angular momentum vector; $s_5$ is identically zero; $s_6$ is the precession frequency of planet 6 and planets 4–5 about their total angular momentum vector. We stress that all precession modes actually appear in the dynamics of all planets (see Eq. (12)), but this qualitative description gives us a good idea of the relative importance of each term in the orbital evolution of each planet.

In order to compute the probability density function of each frequency $s_j$ given our current knowledge of the planetary system, we drew $10^6$ realisations of the star’s mass and planets’ masses and semi-major axes. For each of these realisations, we computed the eigenvalues of the Lagrange–Laplace matrix B and identified them to the frequencies $s_j$ according to their qualitative role described above. In practice, this identification can be made by choosing fictitious initial conditions $\zeta_k(t=0)$ designed to magnify the specific term we are looking for. For instance, the frequency $s_3$ would appear as strongly dominant for all planets if we set $\zeta_k(t=0)=0$ for $k=\{1,2,3\}$ and $\zeta_k(t=0) = \sqrt{2}/2$ for $k=\{4,5,6\}$. Then, one may identify $s_2$ as the dominant term in the solution of planet 2 by setting $\zeta_2(t=0) = \sqrt{2}/2$ and $\zeta_k(t=0)=0$ for $k \ne 2$, etc. This way, all frequencies can be correctly identified one by one. Moreover, we remind the reader that the frequencies $s_j$ only depend on the masses and semi-major axes of the planets, so they do not depend on the fictitious initial conditions chosen here, and they are not plagued with our ignorance of the actual orientations of the planets’ orbital planes.

Figure 3 shows the frequency distribution for each inclination proper mode obtained from our $10^6$ realisations of the system. Frequency $s_4$ has a broad distribution due to the large uncertainties in the masses of planets d and e. Frequency $s_3$, on the contrary, is very peaked, which means that the hierarchy of the two subsystems is a robust property of the HIP 41378 system – unless it contains additional massive planets yet to be discovered. In order to quantify the relative importance of each parameter in the value of each frequency, a correlation analysis can be performed on our large sample of realisations. Here, the small spread in frequency $s_3$ results to be essentially due to the uncertainty in the mass of planet d (see Appendix A).

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Probability density of the inclination proper modes of the HIP 41378 system. Histograms are built from $10^6$ realisations of the Lagrange–Laplace system with the mass and semi-major axis uncertainties in Table 1. The histogram for frequency $s_4$ has a long tail extending beyond the left border of the figure (with 99.7% occurrences above –1850" yr$^{-1}$ and 95.4% above –1000" yr$^{-1}$). The histogram for frequency $s_3$ peaks above the top border of the figure. Frequency $s_5$ is identically equal to zero from the conservation of angular momentum. The upper axis shows the minimum moon mass needed for HIP 41378 f to be fully tilted through a resonance with a given frequency value (see Sect. 4).</figcaption></figure>

As illustrated in Fig. 2, frequency $s_3$ is expected to have a strong contribution in the motion of all planets. Next to it, the dominant inclination proper mode of planet f has frequency $s_6$. This frequency would produce a strong (if not the strongest) secular spin-orbit resonance for this planet. Figure 3 shows that despite observational uncertainties, frequency $s_6$ has a relatively peaked distribution. Its most probable value is –136" yr$^{-1}$, with 68.3% occurrences within [–181, –97]" yr$^{-1}$, 95.4% occurrences within [–241, –65]" yr$^{-1}$, and 99.7% occurrences within [–405, –35]" yr$^{-1}$. As shown in Appendix A, the value of $s_6$ is essentially set by the mass of the perturbing planet e, with a Spearman correlation coefficient $\rho_s \approx –0.8$. The value of
$s_6$ is only weakly ($|\rho_S| \le 0.3$) correlated with the parameters of planet $f$ itself. This low correlation allows us to investigate different values for the frequency $s_6$ independently of the mass and semi-major axis of planet $f$ (that we fix, from now on and in the rest of the article, to their nominal values in Table 1).

## 4. Properties of the hypothetical former moon

Knowing the dominant harmonics in the orbital precession of a planet, Eq. (5) gives the conditions required to tilt the planet and form a ring through the tidal migration and disruption of a moon. In addition to the mass and orbital elements of the planet, Eq. (5) depends on the planet's normalising radius $R$, its oblateness coefficient $J_2$, and the product $\omega\lambda$. For a given super-puff exoplanet, we may assume that the anomalous planet's density is entirely due to the existence of a ring; therefore, the value of $R$ can be chosen so as to produce a conventional bulk density (e.g. that of Uranus or Neptune). In the specific case of HIP 41378 f, Akinsanmi et al. (2020) show that, under the ring hypothesis, its true radius would be $3.7^{+0.3}_{-0.2} R_⊕$. Hence, we adopt the value $R = 3.7 R_⊕$ below as our normalising radius.

For given values of the parameters $J_2$ and $\omega\lambda$, Eq. (5) provides a direct relation between the frequency $v_j$ of the resonance and the minimum mass $m_{min}$ of the former moon. Even though $J_2$ and $\omega\lambda$ are completely unknown for exoplanets, we know that they are related, and in first approximation $J_2 \propto \omega^2$ (planets spinning faster are more flattened; see e.g. Chandrasekhar 1969). For a given moon mass $m$, the condition $|v_j| \le p\eta/2$ in Eq. (5) corresponds to a power law $J_2 \propto \omega^{5/2}$. Because of the coincidental near match between these two exponents (2 and 5/2), our total ignorance of $J_2$ and $\omega\lambda$ does not affect much our estimate of $m_{min}$: we may just set $J_2$ and $\omega\lambda$ to realistic values (e.g. obtained from the Solar System planets) and be assured to obtain relevant results – unless the planet has a particularly exotic internal structure which violates $J_2 \propto \omega^2$. This property is verified in Appendix B in the case of planet HIP 41378 f. As the mass and radius proposed by Akinsanmi et al. (2020) for HIP 41378 f are relatively close to those of Uranus, we choose to apply Eq. (5) using the parameters $J_2$ and $\omega\lambda$ of Uranus (see e.g. Yoder 1995).

Independently of the resonance considered, Eq. (5) can be fulfilled only if the mass parameter $\eta$ of the moon is $\eta \ge 2$. Using the $J_2$ value of Uranus, this condition translates into $m/M \ge 1.2 \times 10^{-4}$. This is the minimum mass ever that the former moon of HIP 41378 f should have had. For a larger moon, the minimum mass $m_{min}$ needed to tilt the planet is proportional to the frequency $v_j$ of the considered resonance. The top horizontal axis in Fig. 3 shows the values of $m_{min}$ computed from Eq. (5) using the parameters $J_2$ and $\lambda\omega$ of Uranus (the tics start at $1.2 \times 10^{-4}$ and go from right to left).

The characteristic spin-axis precession rate of HIP 41378 f computed from Eq. (2) is $p \approx 25'' \text{yr}^{-1}$. According to the left inequality in Eq. (5), this value almost certainly rules out a resonance with frequency $s_3$, because frequency $s_3$ sharply peaks at $s_3 = -15.7'' \text{yr}^{-1}$ (see Fig. 3). The fact that $p > |s_3|$ means that the $s_3$ resonance is located in the green portion of Fig. 1; therefore no capture from a low obliquity is possible in this resonance whatever the mass of the moon. Frequency $s_6$, on the contrary, is the closest resonance reachable by HIP 41378 f. This resonance is expected to be strong for planet $f$, if not the strongest (see Sect. 3). Figure 3 shows that a capture and full tilting within the $s_6$ resonance requires a moon with minimum mass ratio ranging between about $2 \times 10^{-4}$ and $10 \times 10^{-4}$. This corresponds to an absolute mass ranging roughly between Triton's mass and

the mass of our Moon, respectively. More precisely, when the parameters $J_2$ and $\omega\lambda$ of Uranus are assumed for HIP 41378 f, the value of frequency $s_6 = -136_{-269}^{+101''} \text{yr}^{-1}$ obtained in Sect. 3 translates into a minimum moon mass $m_{min}/M = 6_{-5}^{+13} \times 10^{-4}$ ($3\sigma$ uncertainty).

This mass range seems realistic when viewed in the context of the regular moons of the Solar System giant planets. For comparison, the moon-to-planet mass ratio of Titan is $2 \times 10^{-4}$, and the summed masses of the largest moons of Jupiter and Uranus yield ratios of about $2 \times 10^{-4}$ and $1 \times 10^{-4}$, respectively. This similarity among planets motivated the work of Canup & Ward (2006), who found that the formation mechanism of moons around the Solar System giant planets may naturally lead to a common mass scaling, with final mass ratios of a few times $10^{-4}$. Yet, these results do not rule out the existence of larger moons, either because of differing external conditions during their formation, or because of different formation processes (see e.g. the discussion by Saillenfest et al. 2022).

In order to fully incline the planet starting from a low obliquity, the distance that the migrating moon needs to cover depends on the resonance considered, but Fig. 1 shows that one can expect in general a migration from $a_m \approx 0.5 R_M$ to $1 R_M$. Using the $J_2$ value of Uranus, Eq. (4) gives a characteristic length $r_M \approx 11 R_M$ for planet HIP 41378 f, which implies that the moon would need to migrate from roughly 5 to $10 R_M$. Given that $r_M$ is proportional to $J_2^{1/5}$, other realistic values of $J_2$ may change these distances by a small amount (see discussion in Appendix B).

The HIP 41378 system is $2.1_{-0.3}^{+0.4}$ Gyr-old (Lund et al. 2019). As the whole tilting mechanism must have been completed before today, the required migration range for the moon can be translated into a minimum migration rate. In the case of HIP 41378 f, we obtain a migration rate of about $6 \text{ cm yr}^{-1}$ in average. This velocity is comparable to the Moon's migration rate from the Earth (Williams & Boggs 2016), and about two times less than the migration rates of Ganymede from Jupiter (Lainey et al. 2009) or Titan from Saturn (Lainey et al. 2020). In order to power this migration through tidal dissipation within the planet, classical formulas with constant parameters (see e.g. Efroimsky & Lainey 2007) imply that the planet's dissipation coefficient needs to be higher than $k_2/Q \approx 3 \times 10^{-5}$ for a moon mass $m/M = 2 \times 10^{-4}$, and higher than $k_2/Q \approx 6 \times 10^{-6}$ for a moon mass $m/M = 10^{-3}$. For comparison, the value measured for Jupiter's satellite Io is $k_2/Q = 1.102 \pm 0.203 \times 10^{-5}$ (Lainey et al. 2009), and the value measured globally for Saturn's main satellites is $k_2/Q = 1.59 \pm 0.74 \times 10^{-4}$ (Lainey et al. 2017) with a large spread for individual moons extending to much higher values (see Lainey et al. 2020; Jacobson 2022).

## 5. Adiabatic resonance capture

The analysis above shows that when assuming realistic values for the unknown parameters $J_2$ and $\omega\lambda$, the constraints obtained for the planet HIP 41378 f and its hypothetical former moon match well the properties expected for giant planets and moons (i.e. distance, mass, migration rate, and tidal dissipation), at least when viewed in the context of the Solar System. Yet, in order for a planet to be captured and adiabatically tilted within a given resonance, this resonance must be large enough. The width of secular spin-orbit resonances scales as the square root of the amplitude of the term in the orbital series (see Eq. (1)). The $s_6$ term is expected to be among the dominant terms for planet HIP 41378 f, but its amplitude may still be small, depending on the mutual inclinations between the planets' orbital planes. In
order to compute the mutual inclinations of the planets, we need their orbital inclinations $I_k$ and longitudes of ascending nodes $\Omega_k$.

As shown in Table 1, the orbital inclinations $I_k$ of transiting planets with respect to the sky plane are tightly constrained from observations, apart from the mirror degeneracy with respect to 90°. As for the longitudes of nodes $\Omega_k$ in the sky plane, they are not constrained from transit photometry, but we know that their values are likely to be close to each other. Indeed, for a given set of orbital inclinations $I_k$, mutual inclinations between the planets' orbital planes are minimum if their longitudes of node $\Omega_k$ are equal. As a general rule, low mutual inclinations minimise the planets' orbital excitation, and a low orbital excitation is expected in multi-planetary systems for stability reasons.

In systems observed by the transit method, low mutual inclinations are expected also because they maximise the probability of observing several transiting planets. Gravitational interactions produce a precession of the planets' orbital planes, possibly making some of them evolve in and out of transit configuration (see e.g. Becker & Adams 2016). Using the Lagrange–Laplace theory, it is straightforward to compute the fraction of time that a planet spends in and out of transit configuration (see e.g. Fig. 2). In the HIP 41378 system as described in Table 1, only the innermost planet may possibly transit 100% of the time, even if we set all the $\Omega_k$ values of the planets to be equal. Due to orbital precession, the probability to observe five transiting planets (as today) is 30% at best, and the probability to observe six is lower than 5%. As such, the HIP 41378 system would not be classified as ‘continually mutually transiting’ (Becker & Adams 2016, 2017).

The level of orbital excitation of a planetary system can be quantified as a function of the dispersion of their longitudes of ascending node $\Omega_k$ in the sky plane. As shown in Appendix C, allowing for just a few degrees dispersion in $\Omega_k$ can increase the amplitude $S_j$ of several modes in Eq. (12) by orders of magnitude, drastically reducing transit probabilities. In the HIP 41378 system, the level of dispersion of the planets' longitudes of node $\Omega_k$ is therefore likely to be very small, perhaps less than 1°, but their actual values are unknown.

Here, we are interested in the possibility for a planet to be captured in secular spin-orbit resonance from a low initial obliquity. In this context, the larger the resonance, the easier the capture (see e.g. Saillenfest et al. 2020); hence, we actually just need a lower bound for the resonance widths, that is, a lower bound for the amplitudes $S_j$ in Eq. (12). If we show that the resonance capture operates flawlessly for this lower bound, then we can be assured that it will operate as well or even better for the true amplitudes $S_j$. To this aim, we consider that: i) the orbital inclinations of all planets with respect to the sky plane lie on the same side of 90°, and ii) all planets have exactly the same longitude of ascending node $\Omega_k$ in the sky plane. When applied to the HIP 41378 system, this idealised system gives the solution shown in Table 2 for planet f.

In order to produce a resonance capture, the migration of the moon must be slow compared to the oscillations of the resonance angle, so that the parameter change is close to the adiabatic regime (see e.g. Su & Lai 2020). For a given resonance, the oscillation frequency near the resonance centre can be computed through Eq. (7); the frequency scales as the square root of the amplitude $S_j$. When applying Eq. (7) to HIP 41378 f by considering the orbital series in Table 2, one finds that the libration period of the $s_6$ resonance angle when the separatrix appears is $T_{\text{lib}} \approx 547\,000$ years. This value is much smaller than the age of the system ($2.1^{+0.4}_{-0.3}$ Gyr; see Lund et al. 2019). Therefore, even when considering the minimum possible width of the

<table><caption>**Table 2.** Solution for the long-term inclination dynamics of planet HIP 41378 f given by the Lagrange–Laplace system.</caption><thead><tr><th>j</th><th>Identification</th><th>ν<sub>j</sub> (" yr<sup>-1</sup>)</th><th>S<sub>j</sub> × 10<sup>7</sup></th><th>φ<sub>j</sub><sup>(0)</sup> (°)</th></tr></thead><tbody><tr><td>1</td><td>s<sub>5</sub></td><td>0.000</td><td>7 046 055</td><td>0.0</td></tr><tr><td>2</td><td>s<sub>3</sub></td><td>-15.600</td><td>18 545</td><td>0.0</td></tr><tr><td>3</td><td>s<sub>6</sub></td><td>-144.623</td><td>5682</td><td>0.0</td></tr><tr><td>4</td><td>s<sub>1</sub></td><td>-170.310</td><td>972</td><td>180.0</td></tr><tr><td>5</td><td>s<sub>4</sub></td><td>-477.109</td><td>26</td><td>180.0</td></tr><tr><td>6</td><td>s<sub>2</sub></td><td>-477.679</td><td>5</td><td>180.0</td></tr></tbody></table>

**Notes.** All planets have their nominal masses, semi-major axes, and inclinations given in Table 1. The amplitudes $S_j$ are minimised by assuming that all planets have orbital inclinations lying on the same side of 90° (chosen to be $I_k \le 90°$), and all planets have the same longitude of ascending node in the sky plane (chosen to be $\Omega_k = 0$).

resonance, the available time span is more than enough for the planet to oscillate many times within the $s_6$ resonance, allowing an adiabatic drift to occur within this resonance.

This point can be verified by performing a numerical integration of the coupled equations of motion of the planet's spin axis and the orbit of its moon. We used the same setting as Saillenfest et al. (2022): we integrated the secular equations of Correia et al. (2011) expanded at quadrupole order, and forced the orbital evolution of the planet with the quasi-periodic series in Table 2. A typical example of evolution is displayed in Fig. 4. In this example, the mass of the moon is $m/M = 7 \times 10^{-4}$ (i.e. about the mass of Jupiter's moon Europa), and we made the moon migrate outwards at a constant rate, chosen to emulate a tidal parameter $k_2/Q \approx 10^{-5}$. For such a tidal parameter, the moon is expected to migrate from a distance $a_m = 5 R$ to a distance $a_m = 10 R$ in about 1.2 Gyr. The planet was initialised with an obliquity of 0.05 rad and a random precession phase. The eccentricity of the moon and its inclination with respect to its local Laplace plane were both initialised to $10^{-4}$, with random argument of pericentre and longitude of ascending node. As expected, Fig. 4 shows that the adiabatic capture and tilting in resonance $s_6$ is guaranteed on a gigayear timescale. Due to the large separation between timescales, the obliquity oscillations of the planet inside the resonance are not even noticeable in the figure, but they build up in the curve width.

When the system reaches the unstable region, the eccentricity of the moon increases rapidly, which produces chaotic jumps in the planet's obliquity. Indeed, near the border of the unstable region, the timescale for the moon's eccentricity to be multiplied by 100 is a few times the characteristic timescale $\tau$ defined in Eq. (6). Here, one obtains $\tau \approx 100$ yr, which means that the eccentricity increase is extremely fast compared to the planet's spin-axis precession timescale ($T \approx 52\,000$ yr), to the oscillations of the planet inside the resonance ($T_{\text{lib}} \approx 547\,000$ yr), or to the tidal eccentricity damping of the moon (whose timescale is a few millions of years; see e.g. Murray & Dermott 1999).

The simulation in Fig. 4 is stopped when the moon's peri-centre goes below the Roche limit of the planet. At this point, the moon is expected to be disrupted into pieces which would rapidly reorganise into an equatorial disc confined inside the Roche limit (see e.g. Canup 2010; Hyodo et al. 2017). As the moon is lost, the planet is suddenly released from any kind of spin-orbit coupling, and its obliquity remains permanently frozen. In the example shown in Fig. 4, the final obliquity of the planet is about 77°. This value is roughly compatible at $2\sigma$ with the obliquity $\epsilon = 92 \pm 7°$.
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Example of tidal evolution of the planet HIP 41378 f and a hypothetical former moon. The mass of the moon is chosen to be $m/M = 7 \times 10^{-4}$. The moon migrates away at constant rate emulating a tidal parameter $k_2/Q = 10^{-5}$. The trajectory of the system is shown in black; it goes from the leftmost to the rightmost point in about 1.3 Gyr. The available resonances are shown in pink, with their separatrices in red; they are labelled with the frequencies $s_k$ of the corresponding modes (see Sect. 3). In this example, the resonances have the minimum possible widths according to the planets' orbital elements in Table 1. In the hatched blue region, the moon is unstable (same as Fig. 1). The top axis shows the moon distance in unit of the planetary radius.</figcaption></figure>

proposed by Akinsanmi et al. (2020). However, we stress that the final obliquity of the planet is the result of a chaotic phase; its value strongly depends on initial conditions, on the mass of the moon, and on the widths of nearby secular spin-orbit resonances (Saillenfest et al. 2022). More massive moons and larger resonances increase the obliquity excitation of the planet during the chaotic phase. Due to chaos, obliquity values larger than 90° can be reached, but the detailed exploration of possible outcomes would require a precise knowledge of the orbital dynamics of the planet. Without this knowledge, we can only conclude that the obliquity of the planet ends up within the hatched blue region in Fig. 4, that is, between about 70° and 110°.

# 6. Discussions

## 6.1. Refining the tilting mechanism

Under the ring hypothesis, we have presented a proof of concept for producing the unusual configuration proposed for super-puff exoplanets through the tidal migration of a former moon. We have considered the effect of a single massive moon on the planet's spin axis dynamics. This does not mean that the planet only had one moon – we expect it to possibly have many – but that this big moon gathered most of the mass of the satellite system, similarly to Titan around Saturn. Now that this big moon is lost, the remaining moons (either pre-existing or formed in the debris ring) are expected to be very small and undetectable with current facilities.

The presence of several pre-existing big moons, as the Galilean satellites around Jupiter, would complicate the picture outlined here. Through their mutual gravitational perturbations, several massive moons could either inhibit or facilitate the tilting process (see Saillenfest et al. 2022). The exploration of this more complicated scenario is out of the scope of this article. More generally, additional work can refine the scenario proposed here for a

given target exoplanet, including the efficiency of ring formation, the distribution of possible final obliquities, and the combined effect of several massive moons. However, this level of detail would require an in-depth knowledge of the orbital dynamics of the planetary system.

In the case of HIP 41378 f, confirmed periods and masses are still missing for planets d, planet e, and the candidate planet g. The analysis presented here reflects our current understanding of the system, and some results may change in case of substantial modifications in the system's hierarchy. Our correlation analysis shows that planets d and g are only weakly coupled with the frequency $s_6$ of the resonance involved. A mass measurement for these planets would therefore not alter much the picture outlined above. However, substantial changes could be produced if future observations reveal a substantially different mass or period for planet e, or if the system contains an additional outer planet; the calculations presented here should therefore be updated. In this respect, the simplicity of the analytical formulas involved is a great advantage.

## 6.2. The true nature of super-puffs

Future characterisation of super-puff exoplanets is fundamental to assess the actual nature of their anomalously large radii. Unfortunately, due to the nearly face-on configuration of the proposed ring, an unambiguous detection of the ring by transit photometry or by the Rossiter-McLaughlin effect would be challenging with current instruments (Akinsanmi et al. 2020). Spectroscopic observations are much more promising. Even though the spectra of several super-puffs have been revealed to be featureless in near infrared (Libby-Roberts et al. 2020; Chachan et al. 2020; Alam et al. 2022), rings are expected to be transparent in far infrared, which would strongly reduce the transit depth of the planet. As noted by Alam et al. (2022), mid-infrared observations by the JWST would be enough to break the degeneracy between high-altitude hazes, a high-metallicity atmosphere, or the ring hypothesis. The nominal JWST mission offers only two opportunities to observe a transit of HIP 41378 f:

<footnote n="7">⁷ The closed-form expression for the border of the unstable region is $\cos^2 \varepsilon = (51 + 25\sqrt{3})/726$; see Saillenfest & Lari (2021).</footnote>
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Probability of underestimating the density of a transiting exoplanet by a factor *q* due to the presence of an opaque ring with outer radius 2.5 *R* and no inner gap. The enhanced transit depth due to the ring is supposed to be fully misinterpreted as an enlarged planetary radius. The red curve is obtained by computing the ring inclination required to divide the measured planet density by a factor *q*, and by assuming that the inclination of the ring (or equivalently, its precession angle ψ; see text) is uniformly distributed between 0 and 2π. The three planet pictures show the approximate geometries corresponding to the factor *q* in abscissa. The probability goes from 1 at *q* = 1 (exact edge-on configuration) to 0 at *q* = 2.5³ (exact face-on configuration). In case it possesses a ring, planet HIP 41378 f would have *q* ≈ 13 (Akinsanmi et al. 2020).</figcaption></figure>

October 2025 and March 2027. Considering their high scientific value, these opportunities should not be missed. In addition, the high cadence and high photometric resolution of the future PLATO mission may allow small distortions in the transit light curve to be detected (due to the non-zero inclination of the ring with respect to the sky plane and/or to a possible thin inner gap in the ring; see Akinsanmi et al. 2020).

### 6.3. The rarity of enlarged planets

Due to the generic nature of the mechanism presented here, one may wonder why in this case we do not observe many distant exoplanets with anomalously large radii. This rarity can be explained by several factors. First, the transit and radial-velocity methods are strongly biased towards the detection of short-period exoplanets (Perryman 2018). In this regard, the detection of HIP 41378 f with a period of 542 days is already an exception (the transit probability is 0.5%). In turn, the long-period planets observed in direct imaging are strongly biased towards young systems, which cannot have gone through the gigayear adiabatic tilting process described here. As of today, this leaves us with only a handful of exoplanet detections for which this mechanism may have played a role.

The second rarity factor is geometric: a strong radius enhancement able to cast suspicion requires a roughly face-on ring. The mechanism proposed here produces a final planetary obliquity more or less equal to 90°, which is a necessary condition for observing a transiting face-on ring, but is not sufficient: the precession phase ψ of the planet must also have an adequate value. For a ring with typical radius 2.5 *R*, the increase in transit depth leading to underestimating the planet density by a factor *q* > 10 requires a precession phase within ±40° of the exact face-on configuration (see e.g. Zuluaga et al. 2015). As shown in Fig. 5, this occurs about 45% of the time. This fraction is lowered if we consider the ring to have an inner optically thin gap similar to Saturn's ring.

Finally, even though the mechanism described in this article is generic, not all giant planets are expected to reach the final

instability phase in only a few gigayears. Depending on the initial configuration of its moons and the geometry of the available resonances, the planet's obliquity may only have time to increase by a few tens of degrees during its lifetime. In the Solar System, which is aged 4.5 Gyr, only Uranus may have completed the final stage today (Saillenfest et al. 2022). In contrast, Jupiter is only starting the tilting phase<fnref n="8" /> (Saillenfest et al. 2020), while Saturn is seemingly halfway in (Saillenfest et al. 2021a,b) – and it may have recently been ejected from resonance (see Wisdom et al. 2022).

Hence, even though many exoplanets are probably affected by this mechanism, the conjunction of observational biases, ring geometry, and the long timescales at play drastically reduces the probability of detecting targets as exquisite as HIP 41378 f. In this regard, the future PLATO mission is particularly promising, as its observing strategy is tailored to long-period planets, and it will be accompanied by an intensive radial-velocity follow-up to get accurate planet masses and detect possible non-transiting companions. Hopefully, the PLATO discoveries will enable us to estimate the fraction of the exoplanet population that may have gone through the mechanism described in this article.

## 7. Conclusion

The apparent enlarged radius of some long-period exoplanets may be due to the presence of a ring observed roughly face on (Piro & Vissapragada 2020; Akinsanmi et al. 2020). Despite their unconventional configuration, such hypothetical rings and the nearly 90° obliquity of their host planets can be the natural end state of former migrating moons. This mechanism involves the capture of the planet in secular spin-orbit resonance as the moon migrates away on a gigayear timescale. The planet is then gradually tilted until the moon is destabilised and may be destructed into a debris disc.

For a given exoplanet, the plausibility of this formation mechanism can be assessed through simple analytical calculations. First, we need to determine the list of secular spin-orbit resonances that may tilt the planet. The frequencies $v_j$ of the main orbital precession harmonics of the planet can be obtained through the Lagrange-Laplace theory; in this theory, orbital frequencies are the eigenvalues of a matrix which depends only on the masses and spacings of the planets contained in the system. The probability density function of each frequency can be built from numerous realisations of the system (e.g. $10^6$ or more) which are sampled according to our uncertainties on the parameters. Simple correlation analysis can then quantify the influence of each planet in the frequency values.

Then, for each frequency $v_j$, the simple formula in Eq. (5) gives the minimum mass of a moon that the planet must have in order to trigger an adequate secular spin-orbit resonance. This formula depends on the unknown parameters $J_2$ and $\omega\lambda$ of the planet, but thanks to the approximate relation $J_2 \propto \omega^2$, this lack of knowledge only weakly affects the final result. The moon-to-planet mass ratio obtained is the first plausibility check of this dynamical mechanism. Moons with mass ratio $m/M \sim 10^{-4}$ or smaller are expected to be ubiquitous around gaseous planets (see e.g. Canup & Ward 2006). Substantially larger moons cannot be categorically ruled out, but they would require non-generic formation pathways such as captures or giant

<footnote n="8">⁸ As Jupiter possesses four massive moons interacting with each other, its tilting process is somewhat different from what is presented here. Jupiter may never be able to reach an obliquity close to 90° even if it was given infinite time.</footnote>
impacts, and are therefore much less likely (see e.g. Kipping 2014).

A second consistency check is provided by the age of the planetary system considered. The Laplace radius of the planet (see Eq. (4)) sets the distance over which the moon needs to migrate to fully tilt the planet. The migration range obtained must have been covered by the moon in a smaller timespan than the age of the system. As the migration of moons is powered by tidal dissipation inside the planet, the required distance and migration timescale can be translated into a tidal parameter $k_2/Q$ for the planet. Expected values are of the order of $10^{-5}$ to $10^{-4}$ from a Solar System perspective (Lainey et al. 2009, 2017).

The last plausibility check is the consistency of timescales between the age of the planetary system and the hypothesis of an adiabatic capture into resonance. An adiabatic capture requires the libration period inside the resonance to be much shorter than the age of the system, such that many oscillations of the resonance angle may possibly have occurred during the tilting of the planet. The characteristic libration frequency is given in Eq. (7); it depends on the frequencies $v_j$ obtained above, but also on the amplitudes $S_j$ of the corresponding harmonics in the planet's orbital precession spectrum. Using the Lagrange–Laplace theory, the computation of these amplitudes requires to know the inclinations $I_k$ and longitudes of nodes $\Omega_k$ of the planets (e.g. measured in the sky plane). As the libration frequency scales as the square root of the amplitude of the resonant term, only a lower bound for $S_j$ is actually needed. This lower bound can be obtained even in case the longitudes of nodes $\Omega_k$ of the planets are unknown, allowing one to compute a maximum value for the libration period inside the resonance considered.

We applied this methodology to the planet HIP 41378 f, and obtained that all consistency checks are fulfilled. In order to tilt the planet through an adequate resonance, the hypothetical exomoon must have had a moon-to-planet mass ratio ranging from $m/M \ge 2 \times 10^{-4}$ to $10 \times 10^{-4}$, that is, a mass comparable to that of Neptune's moon Triton, Jupiter's moon Europa, or to that of our own Moon. Even though such small exomoons are very hard to detect due to the weakness of their observational signals (Kipping 2014), we expect them to be ubiquitous around giant exoplanets. Provided that the exomoon was initially formed at a distance of about 3 to 10 planetary radii (similarly to Jupiter's moons Io or Europa), its outward migration leads to a guaranteed capture of HIP 41378 f in a secular spin-orbit resonance. The migration timescale required for the moon is found to be in line with what is observed in the Solar System, with a corresponding tidal dissipation factor $k_2/Q$ larger than $2 \times 10^{-5}$ (for the smallest possible moon) or larger than about $6 \times 10^{-6}$ (for a bigger moon). Finally, the libration timescale inside the resonance is found to be orders of magnitude smaller than the age of the system ($2.1_{-0.3}^{+0.4}$ Gyr; see Lund et al. 2019), allowing for the whole tilting mechanism to have possibly occurred.

All these requirements are confirmed by an example of fully coupled numerical integration of the planet's spin axis and the moon's orbit. The planet's spin axis is gradually tilted until its obliquity $\epsilon$ reaches values in the interval $[70^\circ, 110^\circ]$, and its moon becomes unstable (Tremaine et al. 2009; Saillenfest & Lari 2021). Due to the short instability timescale of the exomoon ($\tau \approx 100$ yr in the case of HIP 41378 f), its eccentricity increase is likely to cause catastrophic events, such as collision chains between small inner moons or a tidal disruption of the moon itself when its pericentre goes below the planet's Roche limit (see e.g. Canup 2010; Hyodo et al. 2017; Wisdom et al. 2022). Hence, we argue the dynamical mechanism described here, which may be responsible for the tilting of planet HIP 41378 f to

an obliquity $\epsilon \approx 90^\circ$, can also naturally provide the material for its hypothetical ring.

We stress, however, that even though this dynamical mechanism is physically realistic for HIP 41378 f, this does not imply that it necessarily happened. Planet HIP 41378 f may have had too small and/or too distant moons for the mechanism to operate, and the anomalous transit depth and flat spectrum of this planet may still be due to a particularly tenuous atmosphere covered with high-altitude hazes (Chachan et al. 2020; Alam et al. 2022; Belkovski et al. 2022). Yet, our analysis does provide further significance to the high-obliquity ring hypothesis, by showing that such an unusual configuration is not only feasible in a physical point of view, but even expected for some fraction of exoplanets resembling HIP 41378 f – that is, for old and distant exoplanets in multi-planetary systems. As detailed above, checking the plausibility of this mechanism only requires a limited knowledge of the planetary system considered, and this methodology can be applied to other super-puff exoplanets, and in particular to the potential future discoveries of PLATO.

*Acknowledgements.* The authors thank the anonymous referee for her/his inspiring comments. This work was supported by the Programme National de Planétologie (PNP) of CNRS/INSU, co-funded by CNES.

## References

* Aksinsanmi, B., Santos, N. C., Faria, J. P., et al. 2020, A&A, 635, A8

* Alam, M. K., Kirk, J., Dressing, C. D., et al. 2022, ApJ, 927, L5

* Batygin, K. 2018, AJ, 155, 178

* Batygin, K., Stevenson, D. J., & Bodenheimer, P. H. 2011, ApJ, 738, 1

* Becker, J. C., & Adams, F. C. 2016, MNRAS, 455, 2980

* Becker, J. C., & Adams, F. C. 2017, MNRAS, 468, 549

* Becker, J. C., Vanderburg, A., Rodriguez, J. E., et al. 2019, AJ, 157, 19

* Belkovski, M., Becker, J., Howe, A., Malsky, I., & Batygin, K. 2022, AJ, 163, 277

* Berardo, D., Crossfield, I. J. M., Werner, M., et al. 2019, AJ, 157, 185

* Boué, G., & Laskar, J. 2006, Icarus, 185, 312

* Boué, G., & Laskar, J. 2010, ApJ, 712, L44

* Boué, G., Laskar, J., & Kuchynka, P. 2009, ApJ, 702, L19

* Bryant, E. M., Bayliss, D., Santerne, A., et al. 2021, MNRAS, 504, L45

* Burrows, A., Guillot, T., Hubbard, W. B., et al. 2000, ApJ, 534, L97

* Canup, R. M. 2010, Nature, 468, 943

* Canup, R. M., & Ward, W. R. 2006, Nature, 441, 834

* Chabrier, G., & Baraffe, I. 2007, ApJ, 661, L81

* Chachan, Y., Jontof-Hutter, D., Knutson, H. A., et al. 2020, AJ, 160, 201

* Chandrasekhar, S. 1969, Ellipsoidal Figures of Equilibrium (Yale University Press)

* Correia, A. C. M., Laskar, J., Farago, F., & Boué, G. 2011, Celest. Mech. Dyn. Astron., 111, 105

* Dittmann, A. J. 2021, MNRAS, 508, 1842

* Dong, J., Jiang, Y.-F., & Armitage, P. J. 2021, ApJ, 921, 54

* Efroimsky, M., & Lainey, V. 2007, J. Geophys. Res. Planets, 112, 12003

* Fabrycky, D. C., Lissauer, J. J., Ragozzine, D., et al. 2014, ApJ, 790, 146

* Farhat, M., Auclair-Desrotour, P., Boué, G., & Laskar, J. 2022, A&A, 665, A1

* Gao, P., & Zhang, X. 2020, ApJ, 890, 93

* Grunblatt, S. K., Huber, D., Gaidos, E. J., et al. 2016, AJ, 152, 185

* Grouffal, S., Santerne, A., Bourrier, V., et al. 2022, A&A, 668, A172

* Hamilton, D. P., & Ward, W. R. 2004, AJ, 128, 2510

* Henrard, J., & Murigande, C. 1987, Celest. Mech., 40, 345

* Hyodo, R., Charnoz, S., Ohtsuki, K., & Genda, H. 2017, Icarus, 282, 195

* Jacobson, R. A. 2022, AJ, 164, 199

* Kipping, D. M. 2014, in Proceedings of Frank N. Bash Symposium 2013: New Horizons in Astronomy – PoS(BASH 2013), 206, 012

* Lainey, V., Arlot, J.-E., Karatekin, Ö., & van Hoolst, T. 2009, Nature, 459, 957

* Lainey, V., Jacobson, R. A., Tajeddine, R., et al. 2017, Icarus, 281, 286

* Lainey, V., Gomez Casajus, L., Fuller, J., et al. 2020, Nat. Astron., 4, 1053

* Laskar, J., & Petit, A. C. 2017, A&A, 605, A72

* Laskar, J., & Robutel, P. 1993, Nature, 361, 608

* Laskar, J., & Robutel, P. 1995, Celest. Mech. Dyn. Astron., 62, 193

* Laughlin, G., Crismani, M., & Adams, F. C. 2011, ApJ, 729, L7

* Lee, E. J., & Chiang, E. 2016, ApJ, 817, 90
Libby-Roberts, J. E., Berta-Thompson, Z. K., Désert, J.-M., et al. 2020, AJ, 159, 57

Lopez, E. D., & Fortney, J. J. 2016, ApJ, 818, 4

Lund, M. N., Knudstrup, E., Silva Aguirre, V., et al. 2019, AJ, 158, 248

Millholland, S. 2019, ApJ, 886, 72

Millholland, S., & Laughlin, G. 2019, Nat. Astron., 3, 424

Millholland, S. C., & Spalding, C. 2020, ApJ, 905, 71

Morbidelli, A., Tsiganis, K., Batygin, K., Crida, A., & Gomes, R. 2012, Icarus, 219, 737

Murray, C. D., & Dermott, S. F. 1999, Solar System Dynamics (Cambridge University Press)

Néron de Surgery, O., & Laskar, J. 1997, A&A, 318, 975

Ohno, K., & Tanaka, Y. A. 2021, ApJ, 920, 124

Perryman, M. 2018, The Exoplanet Handbook, 2nd edn. (Cambridge University Press)

Piro, A. L., & Vissapragada, S. 2020, AJ, 159, 131

Pu, B., & Valencia, D. 2017, ApJ, 846, 47

Rauer, H., Catala, C., Aerts, C., et al. 2014, Exp. Astron., 38, 249

Rauer, H., Aerts, C., Cabrera, J., & PLATO Team. 2016, Astron. Nachr., 337, 961

Rogoszinski, Z., & Hamilton, D. P. 2020, ApJ, 888, 60

Rogoszinski, Z., & Hamilton, D. P. 2021, Planet. Sci. J., 2, 78

Rufu, R., & Canup, R. M. 2022, ApJ, 928, 123

Saillenfest, M., & Lari, G. 2021, A&A, 654, A83

Saillenfest, M., Laskar, J., & Boué, G. 2019, A&A, 623, A4

Saillenfest, M., Lari, G., & Courtot, A. 2020, A&A, 640, A11

Saillenfest, M., Lari, G., & Boué, G. 2021a, Nat. Astron., 5, 345

Saillenfest, M., Lari, G., Boué, G., & Courtot, A. 2021b, A&A, 647, A92

Saillenfest, M., Rogoszinski, Z., Lari, G., et al. 2022, A&A, 668, A108

Salmon, J., & Canup, R. M. 2022, ApJ, 924, 6

Santerne, A., Malavolta, L., Kosiarek, M. R., et al. 2019, ArXiv e-prints [arXiv:1911.07355]

Su, Y., & Lai, D. 2020, ApJ, 903, 7

Su, Y., & Lai, D. 2022a, MNRAS, 509, 3301

Su, Y., & Lai, D. 2022b, MNRAS, 513, 3302

Tremaine, S. 1991, Icarus, 89, 85

Tremaine, S., Touma, J., & Namouni, F. 2009, AJ, 137, 3706

Vanderburg, A., Becker, J. C., Kristiansen, M. H., et al. 2016, ApJ, 827, L10

Vokrouhlický, D., & Nesvorný, D. 2015, ApJ, 806, 143

Wang, L., & Dai, F. 2019, ApJ, 873, L1

Ward, W. R., & Hamilton, D. P. 2004, AJ, 128, 2501

Williams, J. G., & Boggs, D. H. 2016, Celest. Mech. Dyn. Astron., 126, 89

Wisdom, J., Dbouk, R., Militzer, B., et al. 2022, Science, 377, 1285

Xie, J.-W., Dong, S., Zhu, Z., et al. 2016, PNAS, 113, 11431

Yoder, C. F. 1995, in Global Earth Physics: A Handbook of Physical Constants, ed. T. J. Ahrens (American Geophysical Union), 1

Zuluaga, J. I., Kipping, D. M., Sucerquia, M., & Alvarado, J. A. 2015, ApJ, 803, L14
**Appendix A: Correlations between parameters and orbital precession frequencies**

<figure><img src="image_7.png" /></figure>

**Fig. A.1.** Scatter plot of the proper frequency $s_3$ and the 13 parameters involved (the stellar mass and the masses and semi-major axes of planets 1 to 6). The black dots show 150 000 realisations of the Lagrange–Laplace system with the mass and semi-major axis uncertainties from Table 1. The red label gives Spearman's correlation coefficient $\rho_S$ as computed from $10^6$ realisations.

<figure><img src="image_8.png" /></figure>

**Fig. A.2.** Same as Fig. A.1, but for frequency $s_6$.

In a long-term stable planetary system, the orbital motion of the planets can be approximated by quasi-periodic series as in Eq. (1), where the frequencies $v_j$ are integer combinations of the fundamental frequencies of the system. In the Lagrange–Laplace approximation, the fundamental frequencies $s_k$ of the series governing the inclination dynamics of the planets are the eigenvalues of the matrix B (see e.g. Murray & Dermott 1999) which only depends on the masses and semi-major axes of the planets.

The values of the frequencies $s_k$ and their role in the dynamics are intrinsic properties of the matrix B. As the fundamental frequencies reflect the gravitational couplings between the planets, some parameters contribute much more than others in the value of a given frequency; however, this contribution is not linear and it is not obvious a priori which parameters contribute the most. Here, we are mostly interested in the frequencies $s_3$ and $s_6$ as defined in Sect. 3, because they are expected to dominate the inclination dynamics of planet HIP 41378 f. Figures A.1 and A.2 show the scatter of the values of $s_3$ and $s_6$ as a function of all parameters.

From these scatter plots, it is visible that the value of $s_3$ mostly depends on $M_4$, while the value of $s_6$ mostly depends on $M_5$. The strength of these correlations can be quantified by Spearman's correlation coefficient $\rho_S$. The coefficients obtained are given in Figs. A.1 and A.2 for each parameter. The correlations between $s_3$ and $M_4$ ($\rho_S = -0.85$) and between $s_6$ and $M_5$ ($\rho_S = -0.76$) are the strongest by more than a factor of two.

**Appendix B: Influence of the unknown physical parameters of the planet**

In Sect. 4, we estimate the properties of a hypothetical former moon needed to tilt a planet from a low obliquity and create a ring of debris. The formulas, however, depend on the unknown parameters $J_2$ and $\omega\lambda$ of the planet. In this section, we consider them as free parameters and study their influence on our results.

The values of $J_2$ and $\omega\lambda$ are related and depend on the interior properties of the planet. In the simplest case of a homogeneous planet, $J_2$ and $\omega$ are linked through the law of Maclaurin's
<figure><img src="image_9.png" /><figcaption>**Fig. B.1.** Minimum moon mass needed for HIP 41378 f to undergo the full tilting mechanism through a resonance with $s_6$. Each panel corresponds to a given possible value of $s_6$ (labels). The minimum moon mass $m_{min}$ (colour scale) is obtained from Eq. (5) as a function of the unknown parameters $J_2$ and $\omega\lambda$ of HIP 41378 f. The white curves highlight three mass levels, namely $m_{min}/M = 3 \times 10^{-4}$, $6 \times 10^{-4}$, and $10 \times 10^{-4}$. The grey region is forbidden by the left inequality in Eq. (5). The black dotted line shows the relation between $J_2$ and $\omega\lambda$ obtained for a homogeneous planet (see Eq. B.1). For reference, the black dots show the parameters $(\omega\lambda, J_2)$ of the Solar System giant planets (Yoder 1995).</figcaption></figure>

ellipsoid (see e.g. Chandrasekhar 1969), while $\lambda = 2/5$. This law simplifies to $J_2 \propto \omega^2$ for a nearly spherical planet, and it can be rewritten as

$$J_2 \approx \frac{25}{8} \frac{R^3}{GM} (\omega\lambda)^2. \qquad (B.1)$$

Even though planets are not homogeneous, this approximate relation gives an idea of where to look for realistic combinations of parameters in the plane $(\omega\lambda, J_2)$ without any assumption on its composition. Figure B.1 shows that the Solar System giant planets do fall roughly along this curve.

For a given frequency $\nu_j$ in the orbital precession spectrum of a planet, the conditions $p \le |\nu_j|$ in Eq. (5) corresponds to a straight line $J_2 \propto \omega\lambda$, whereas the condition $|\nu_j| \le p\eta/2$ corresponds to a power law $J_2 \propto (\omega\lambda)^{5/2}$. Both conditions can be visualised in Fig. B.1 for planet HIP 41378 f using three different values of the frequency $\nu_j$ (most probable value of $s_6$ and 95.4% bounds; see Sect. 3). Because the exponent 2 in Eq. (B.1) is close to 5/2, the realistic combinations of $J_2$ and $\omega\lambda$ (which are located in a rough neighbourhood of the dotted curve) follow more or less the level curves of the minimum moon mass $m_{min}$. As a consequence, our estimate of the minimum moon mass is not affected much by our total ignorance of the parameters $J_2$ and $\omega\lambda$. Whatever realistic values are chosen, we obtain a minimum moon mass $m_{min}/M \approx 6 \times 10^{-4}$ for the most probable value of $s_6$ (Fig. B.1b), with a dispersion at 95.4% ranging from $m_{min}/M \approx 3 \times 10^{-4}$ (Fig. B.1a) to $m_{min}/M \approx 1 \times 10^{-3}$ (Fig. B.1c). These values are similar to those obtained in Sect. (4) using the parameters of Uranus.

The value of characteristic length $r_M$, however, depends on $J_2$, but not on $\omega\lambda$ (see Eq. 4). The distance covered by the migrating moon depends therefore on the value chosen for $J_2$. Yet, this dependence is not very steep ($r_M \propto J_2^{1/5}$). Moreover, extreme values of $J_2$, either very small or very large, can be ruled out because giant planets are expected to spin at a fraction of their breakup velocity (see e.g. Batygin 2018; Dong et al. 2021; Dittmann 2021). As a consequence, the value obtained for $r_M$ only varies by a small amount when considering realistic values of $J_2$ (see the right vertical axis in Fig. B.1).

## Appendix C: Inclination amplitudes and transit probabilities

In the Lagrange–Laplace approximation, the long-term inclination dynamics of planets are described by quasi-periodic series as in Eq. (12), in which the frequencies solely depend on the masses and spacing of the planets. The orientations of the planets’ orbital planes enter into play only in the amplitudes of the terms of the series. Apart from the zero-frequency term (which merely gives the orientation of the planetary invariant plane), the amplitudes $S_j$ in Eq. (12) depend on the mutual inclinations between the planets’ orbital planes. The mutual inclination $\Psi_{ik}$ of two planets i and k can be written as

$$\cos \Psi_{ik} = \cos I_i \cos I_k + \sin I_i \sin I_k \cos(\Omega_i - \Omega_k), \qquad (C.1)$$

where $I$ and $\Omega$ are the orbital inclination and longitude of ascending node of the planets measured with respect to a given reference plane (e.g. the sky plane). As a simple rule of thumb, the larger the mutual inclinations between each pair of planets, the larger the amplitudes $S_j$. For unknown longitudes of node $\Omega_k$, the orbital inclinations $I_k$ of a set of planets provide minimum and maximum bounds to the mutual inclinations between each pair of planets. From Eq. (C.1), the minimum mutual inclination of two planets is $\Psi_{ik} = |I_i - I_k|$ reached for $\Omega_i - \Omega_k = 0$. In practice, inclination values measured from transit data have a mirror degeneracy with respect to 90° (see Table 1 of the main text); $\Psi_{ik}$ is minimised if $I_i$ and $I_k$ lie on the same side of 90°.

Here, we quantify the influence of the unknown longitudes of nodes on the orbital dynamics of planet HIP 41378 f by drawing
<figure><img src="image_10.png" /><figcaption>**Fig. C.1.** Possible values of the amplitudes $S_j$ for the inclination terms of planet HIP 41378 f. The frequency of each term is labelled. The amplitude of the zero-frequency term $s_5$ is not shown. The masses, semi-major axes, and inclinations of all planets are set to their nominal observed values (see Table 1); all inclinations values $I_k$ are assumed to be $I_k \le 90^\circ$. The unknown longitudes of node $\Omega_k$ of the planets in the sky plane are drawn from a uniform random distribution in an interval $\Delta\Omega_0$ common for all planets. For a given value of the range $\Delta\Omega_0$ (vertical axis), a histogram of each amplitude is built from $10^6$ realisations of the Lagrange–Laplace system (colour scale). A bin is coloured grey if no occurrence is found among our $10^6$ realisations. The scale is the same for all graphs to highlight the disparity of amplitude between the different terms.</figcaption></figure>

randomly the longitudes of nodes of all planets in a given interval $\Delta\Omega_0$ and building a histogram of the amplitudes $S_j$ obtained. The result is shown in Fig. C.1 for an interval $\Delta\Omega_0$ ranging from $0^\circ$ to $2^\circ$. As expected, the terms $s_1$ and $s_2$ have very small amplitudes; they contribute negligibly to the dynamics of planet HIP 41378 f. The dominant terms of the dynamics are $s_3$ and $s_6$, with the qualitative roles described in Sect. 3. The amplitudes of all terms are generally minimum for $\Delta\Omega_0 = 0$, and they cover a wider and wider range of possibilities when we allow the planets' longitudes of nodes to be substantially distinct. This effect is particularly visible for the $s_4$ term, for which a dispersion as small as $\Delta\Omega_0 = 1^\circ$ can make the amplitude increase by a factor 1000 or so.

Figure C.2 shows the probability of observing five transiting planets or more in the same experiment as in Fig. C.1. For a given realisation of the planetary system, the transit probability of a set of planets is the fraction of time their orbits simultaneously pass in front of the star. If we assume the same longitude of node for all planets in the HIP 41378 system (i.e. $\Delta\Omega_0 = 0$), then five planets or more transit about 30% of the time. If we allow a dispersion of $\Delta\Omega_0 = 1^\circ$, this fraction can be reduced to less than 5%.

<figure><img src="image_11.png" /><figcaption>**Fig. C.2.** Probability of observing five transiting planets in the HIP 41378 system over the planets' precession cycles. The unknown longitudes of node $\Omega_k$ of the planets in the sky plane are drawn from a uniform random distribution in an interval $\Delta\Omega_0$ common for all planets. For a given value of the range $\Delta\Omega_0$ (vertical axis), a histogram of the transit probability is built from $10^6$ realisations of the Lagrange–Laplace system (colour scale). A bin is coloured grey if no occurrence is found among our $10^6$ realisations.</figcaption></figure>