11institutetext: IMCCE, Observatoire de Paris, PSL Research University, CNRS,
Sorbonne Université, Université de Lille, 75014 Paris, France
11email<EMAIL_ADDRESS>22institutetext: Université Aix
Marseille, CNRS, CNES, LAM, Marseille, France 33institutetext: Université de
Toulouse, CNRS, IRAP, Toulouse, France

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

Melaine Saillenfest 11 Sophia Sulis 22 Paul Charpentier 33 Alexandre Santerne
22

(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 $\lesssim
15$ $M_{\oplus}$) but surprisingly large radii ($\gtrsim 6$ $R_{\oplus}$),
giving them extremely low bulk densities ($\lesssim 0.3$ 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$ 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 planet111Not 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). is nearly $90^{\circ}$.

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 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 obliquity222High-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, 2022b, a).. 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 $\varepsilon\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 $\varepsilon\gtrsim 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
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 $\varepsilon$ up to $\varepsilon\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)]\,,$
(1)

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

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

is the characteristic spin-axis precession rate of the planet. In this
expression, $\mathcal{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}}\,,$ (3)

where $m$ is the mass of the moon, $M$ is the mass of the planet, and
$r_{\mathrm{M}}$ is the following characteristic length333$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$.:

$r_{\mathrm{M}}^{5}=2\frac{M}{M_{\star}}J_{2}R^{2}\leavevmode\nobreak\
a^{3}(1-e^{2})^{3/2}\,.$ (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 $\nu_{j}$
verifying

$p\leqslant|\nu_{j}|\leqslant p\frac{\eta}{2}$ (5)

can allow the planet’s obliquity to grow from $\varepsilon=0^{\circ}$ to
$\varepsilon=90^{\circ}$. This condition is illustrated in Fig. 1. Knowing the
harmonics $\nu_{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).

Figure 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.

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}{4}\frac{M_{\star}}{M}\frac{r_{\mathrm{M}}^{3}}{a^{3}(1-e^{2})^{3}}\frac{\mathcal{G}M_{\star}}{a^{3}}\,.$
(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^{\prime})^{2}\left(\frac{\beta^{2}}{\sin^{2}\varepsilon_{0}}+\beta\sin\varepsilon_{0}\right)\,.$
(7)

In this expression, $\varepsilon_{0}$ is the planet’s obliquity at the
resonance centre and $p^{\prime}$ 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^{\prime}$ and
$\beta=\rho_{2}/p^{\prime}$, where

$\displaystyle\rho_{1}$
$\displaystyle=-\left(\nu_{k}-2\sum_{j=1}^{N}\nu_{j}S_{j}^{2}\right)\,,$ (8)
$\displaystyle\rho_{2}$
$\displaystyle=-S_{k}\left(2\nu_{k}+\nu_{k}S_{k}^{2}-2\sum_{j=1}^{N}\nu_{j}S_{j}^{2}\right)\,,$

and $k$ is the index in Eq. (1) of the considered resonance.
$T_{\mathrm{lib}}$ depends on the distance of the moon through $p^{\prime}$
and $\varepsilon_{0}$. However, an upper bound for $T_{\mathrm{lib}}$ is
obtained at the time of resonance capture, for which
$\gamma^{2/3}+\beta^{2/3}=1$ (see Henrard & Murigande, 1987; Saillenfest et
al., 2019). In this case, $p^{\prime}$ is equal to
$p^{\prime}=(\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}}\,.$
(9)

Thanks to these expressions, we can compute $T_{\mathrm{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})\,,$ (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{\mathrm{d}\boldsymbol{\zeta}}{\mathrm{d}t}=iB\,\boldsymbol{\zeta}\,,$
(11)

in which $\boldsymbol{\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_{\mathrm{p}}}S_{j}\exp\left[i\left(\nu_{j}\,t+\phi_{j}^{(0)}\right)\right]\,,$
(12)

where the number of terms $N$ is equal to the number of planets
$N_{\mathrm{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_{\mathrm{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_{\mathrm{p}}$ planets
represent 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_{\mathrm{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.

###  3.2 Application to the HIP 41378 system

HIP 41378 is bright F-type star444Also known as K2-93 and EPIC 211311380.
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\pm 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$ 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$ g cm-3, similar to that of
Uranus. The hypothetical equatorial ring provides an indirect measure of the
obliquity of the planet, namely555The ring obtained by Akinsanmi et al. (2020)
is inclined by $i_{\text{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_{\text{r}}+\sin I\sin i_{\text{r}}\cos\theta$, where $I=89.97\pm
0.01^{\circ}$ is the orbital inclination of HIP 41378 f. $\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
Archive666https://exoplanetarchive.ipac.caltech.edu 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.

Table 1: Parameters of planets in the HIP 41378 system used in this article.

$k$ | name |  $M_{k}$ ($M_{\oplus}$) |  $P_{k}$ (day) |  $a_{k}$ (au) |  $I_{k}$ (${}^{\text{o}}$)
---|---|---|---|---|---
$1$ | b | $6.89\pm 0.88$ | $15.57208\pm 0.00002$ | $0.1283\pm 0.0015$ | $88.75\pm 0.13$
$2$ | c | $4.4\pm 1.1$ | $31.70603\pm 0.00006$ | $0.2061\pm 0.0024$ | $88.48\pm 0.07$
$3$ | g | $7.0\pm 1.5$ | $62.06\pm 0.32$ | $0.3227\pm 0.0036$ | $88$
$4$ | d | $12.7\pm 6.0$ | $278.3618\pm 0.0005$ | $0.88\pm 0.01$ | $89.80\pm 0.02$
$5$ | e | $12\pm 5$ | $369\pm 10$ | $1.06\pm 0.03$ | $89.84\pm 0.07$
$6$ | f | $12\pm 3$ | $542.0798\pm 0.0002$ | $1.37\pm 0.02$ | $89.97\pm 0.01$
777Parameters 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}\leqslant 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).

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.

Figure 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}\leqslant 90^{\circ}$) and random longitudes of node
$\Omega_{k}$ in a $0.2^{\circ}$-wide interval; this choice is commented in
Sect. 5.

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 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\neq 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 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$ ′′ year-1 and
$95.4\%$ above $-1000$ ′′ year-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).

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$ ′′ year-1, with $68.3\%$ occurrences within $[-181,-97]$ ′′
year-1, $95.4\%$ occurrences within $[-241,-65]$ ′′ year-1, and $99.7\%$
occurrences within $[-405,-35]$ ′′ year-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_{\mathrm{S}}\approx-0.8$. The value of
$s_{6}$ is only weakly ($|\rho_{\mathrm{S}}|\lesssim 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.2}^{+0.3}$ $R_{\oplus}$. Hence, we adopt the value $R=3.7$
$R_{\oplus}$ 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 $\nu_{j}$ of the resonance
and the minimum mass $m_{\mathrm{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 $|\nu_{j}|\leqslant 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_{\mathrm{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\geqslant 2$. Using the $J_{2}$
value of Uranus, this condition translates into $m/M\geqslant 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_{\mathrm{min}}$ needed
to tilt the planet is proportional to the frequency $\nu_{j}$ of the
considered resonance. The top horizontal axis in Fig. 3 shows the values of
$m_{\mathrm{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^{\prime\prime}$ year-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^{\prime\prime}$ year-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^{+101}_{-269}$ ′′ year-1 obtained in Sect.
3 translates into a minimum moon mass $m_{\mathrm{min}}/M=6^{+13}_{-5}\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_{\mathrm{m}}\approx 0.5\,r_{\mathrm{M}}$ to $1\,r_{\mathrm{M}}$. Using the
$J_{2}$ value of Uranus, Eq. (4) gives a characteristic length
$r_{\mathrm{M}}\approx 11$ $R$ for planet HIP 41378 f, which implies that the
moon would need to migrate from roughly $5$ to $10$ $R$. Given that
$r_{\mathrm{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.4}_{-0.3}$ 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$ cm year-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^{\circ}$. 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^{\circ}$, 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^{\circ}$,
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.

Table 2: Solution for the long-term inclination dynamics of planet HIP 41378 f
given by the Lagrange-Laplace system.

$j$ | identification |  $\nu_{j}$ (${}^{\prime\prime}\,\text{yr}^{-1}$) | $S_{j}\times 10^{7}$ |  $\phi_{j}^{(0)}$ (${}^{\text{o}}$)
---|---|---|---|---
$1$ | $s_{5}$ | $0.000$ | $7046055$ | $0.0$
$2$ | $s_{3}$ | $-15.600$ | $18545$ | $0.0$
$3$ | $s_{6}$ | $-144.623$ | $5682$ | $0.0$
$4$ | $s_{1}$ | $-170.310$ | $972$ | $180.0$
$5$ | $s_{4}$ | $-477.109$ | $26$ | $180.0$
$6$ | $s_{2}$ | $-477.679$ | $5$ | $180.0$
888All 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^{\circ}$
(chosen to be $I_{k}\leqslant 90^{\circ}$), and all planets have the same
longitude of ascending node in the sky plane (chosen to be $\Omega_{k}=0$).

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_{\mathrm{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 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_{\mathrm{m}}=5$ $R$ to a distance $a_{\mathrm{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.

Figure 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.

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_{\mathrm{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 pericentre 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^{\circ}$. This value is roughly
compatible at $2\sigma$ with the obliquity $\varepsilon=92\pm 7^{\circ}$
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^{\circ}$ 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 about999The
closed-form expression for the border of the unstable region is
$\cos^{2}\varepsilon=(51+25\sqrt{3})/726$; see Saillenfest & Lari (2021).
$70^{\circ}$ and $110^{\circ}$.

##  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 complicate 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: 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^{\circ}$, which
is a necessary condition for observing a transiting face-on ring, but is not
sufficient: the precession phase $\psi$ 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 $\pm 40^{\circ}$ 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.

Figure 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 $\psi$; see text) is uniformly
distributed between $0$ and $2\pi$. 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^{3}$ (exact face-on configuration). In case it possesses a ring, planet
HIP 41378 f would have $q\approx 13$ (Akinsanmi et al., 2020).

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 their 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 phase101010As 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^{\circ}$ even if it was given infinite time. (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^{\circ}$ 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 $\nu_{j}$ of the main orbital precession harmonics of the planet
can be obtained through the Lagrange-Lagrange 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 $\nu_{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 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 that 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 $\nu_{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\gtrsim 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 magnitudes
smaller than the age of the system ($2.1^{+0.4}_{-0.3}$ 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 $\varepsilon$ 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 $\varepsilon\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

  * Akinsanmi et al. (2020) Akinsanmi, B., Santos, N. C., Faria, J. P., et al. 2020, A&A, 635, L8
  * Alam et al. (2022) Alam, M. K., Kirk, J., Dressing, C. D., et al. 2022, ApJ, 927, L5
  * Batygin (2018) Batygin, K. 2018, AJ, 155, 178
  * Batygin et al. (2011) Batygin, K., Stevenson, D. J., & Bodenheimer, P. H. 2011, ApJ, 738, 1
  * Becker & Adams (2016) Becker, J. C. & Adams, F. C. 2016, MNRAS, 455, 2980
  * Becker & Adams (2017) Becker, J. C. & Adams, F. C. 2017, MNRAS, 468, 549
  * Becker et al. (2019) Becker, J. C., Vanderburg, A., Rodriguez, J. E., et al. 2019, AJ, 157, 19
  * Belkovski et al. (2022) Belkovski, M., Becker, J., Howe, A., Malsky, I., & Batygin, K. 2022, AJ, 163, 277
  * Berardo et al. (2019) Berardo, D., Crossfield, I. J. M., Werner, M., et al. 2019, AJ, 157, 185
  * Boué & Laskar (2006) Boué, G. & Laskar, J. 2006, Icarus, 185, 312
  * Boué & Laskar (2010) Boué, G. & Laskar, J. 2010, ApJ, 712, L44
  * Boué et al. (2009) Boué, G., Laskar, J., & Kuchynka, P. 2009, ApJ, 702, L19
  * Bryant et al. (2021) Bryant, E. M., Bayliss, D., Santerne, A., et al. 2021, MNRAS, 504, L45
  * Burrows et al. (2000) Burrows, A., Guillot, T., Hubbard, W. B., et al. 2000, ApJ, 534, L97
  * Canup (2010) Canup, R. M. 2010, Nature, 468, 943
  * Canup & Ward (2006) Canup, R. M. & Ward, W. R. 2006, Nature, 441, 834
  * Chabrier & Baraffe (2007) Chabrier, G. & Baraffe, I. 2007, ApJ, 661, L81
  * Chachan et al. (2020) Chachan, Y., Jontof-Hutter, D., Knutson, H. A., et al. 2020, AJ, 160, 201
  * Chandrasekhar (1969) Chandrasekhar, S. 1969, Ellipsoidal figures of equilibrium (Yale University Press)
  * Correia et al. (2011) Correia, A. C. M., Laskar, J., Farago, F., & Boué, G. 2011, Celestial Mechanics and Dynamical Astronomy, 111, 105
  * Dittmann (2021) Dittmann, A. J. 2021, MNRAS, 508, 1842
  * Dong et al. (2021) Dong, J., Jiang, Y.-F., & Armitage, P. J. 2021, ApJ, 921, 54
  * Efroimsky & Lainey (2007) Efroimsky, M. & Lainey, V. 2007, Journal of Geophysical Research (Planets), 112
  * Fabrycky et al. (2014) Fabrycky, D. C., Lissauer, J. J., Ragozzine, D., et al. 2014, ApJ, 790, 146
  * Farhat et al. (2022) Farhat, M., Auclair-Desrotour, P., Boué, G., & Laskar, J. 2022, A&A, 665, L1
  * Gao & Zhang (2020) Gao, P. & Zhang, X. 2020, ApJ, 890, 93
  * Grouffal et al. (2022) Grouffal, S., Santerne, A., Bourrier, V., et al. 2022, A&A, 668, A172
  * Grunblatt et al. (2016) Grunblatt, S. K., Huber, D., Gaidos, E. J., et al. 2016, AJ, 152, 185
  * Hamilton & Ward (2004) Hamilton, D. P. & Ward, W. R. 2004, AJ, 128, 2510
  * Henrard & Murigande (1987) Henrard, J. & Murigande, C. 1987, Celestial Mechanics, 40, 345
  * Hyodo et al. (2017) Hyodo, R., Charnoz, S., Ohtsuki, K., & Genda, H. 2017, Icarus, 282, 195
  * Jacobson (2022) Jacobson, R. A. 2022, AJ, 164, 199
  * Kipping (2014) Kipping, D. M. 2014, in Proceedings of Frank N. Bash Symposium 2013: New Horizons in Astronomy — PoS(BASH 2013), Vol. 206, 012
  * Lainey et al. (2009) Lainey, V., Arlot, J.-E., Karatekin, Ö., & van Hoolst, T. 2009, Nature, 459, 957
  * Lainey et al. (2020) Lainey, V., Gomez Casajus, L., Fuller, J., et al. 2020, Nature Astronomy, 4, 1053
  * Lainey et al. (2017) Lainey, V., Jacobson, R. A., Tajeddine, R., et al. 2017, Icarus, 281, 286
  * Laskar & Petit (2017) Laskar, J. & Petit, A. C. 2017, A&A, 605, A72
  * Laskar & Robutel (1993) Laskar, J. & Robutel, P. 1993, Nature, 361, 608
  * Laskar & Robutel (1995) Laskar, J. & Robutel, P. 1995, Celestial Mechanics and Dynamical Astronomy, 62, 193
  * Laughlin et al. (2011) Laughlin, G., Crismani, M., & Adams, F. C. 2011, ApJ, 729, L7
  * Lee & Chiang (2016) Lee, E. J. & Chiang, E. 2016, ApJ, 817, 90
  * Libby-Roberts et al. (2020) Libby-Roberts, J. E., Berta-Thompson, Z. K., Désert, J.-M., et al. 2020, AJ, 159, 57
  * Lopez & Fortney (2016) Lopez, E. D. & Fortney, J. J. 2016, ApJ, 818, 4
  * Lund et al. (2019) Lund, M. N., Knudstrup, E., Silva Aguirre, V., et al. 2019, AJ, 158, 248
  * Millholland (2019) Millholland, S. 2019, ApJ, 886, 72
  * Millholland & Laughlin (2019) Millholland, S. & Laughlin, G. 2019, Nature Astronomy, 3, 424
  * Millholland & Spalding (2020) Millholland, S. C. & Spalding, C. 2020, ApJ, 905, 71
  * Morbidelli et al. (2012) Morbidelli, A., Tsiganis, K., Batygin, K., Crida, A., & Gomes, R. 2012, Icarus, 219, 737
  * Murray & Dermott (1999) Murray, C. D. & Dermott, S. F. 1999, Solar system dynamics (Cambridge University Press)
  * Néron de Surgy & Laskar (1997) Néron de Surgy, O. & Laskar, J. 1997, A&A, 318, 975
  * Ohno & Tanaka (2021) Ohno, K. & Tanaka, Y. A. 2021, ApJ, 920, 124
  * Perryman (2018) Perryman, M. 2018, The Exoplanet Handbook, 2nd edition (Cambridge University Press)
  * Piro & Vissapragada (2020) Piro, A. L. & Vissapragada, S. 2020, AJ, 159, 131
  * Pu & Valencia (2017) Pu, B. & Valencia, D. 2017, ApJ, 846, 47
  * Rauer et al. (2016) Rauer, H., Aerts, C., Cabrera, J., & PLATO Team. 2016, Astronomische Nachrichten, 337, 961
  * Rauer et al. (2014) Rauer, H., Catala, C., Aerts, C., et al. 2014, Experimental Astronomy, 38, 249
  * Rogoszinski & Hamilton (2020) Rogoszinski, Z. & Hamilton, D. P. 2020, ApJ, 888, 60
  * Rogoszinski & Hamilton (2021) Rogoszinski, Z. & Hamilton, D. P. 2021, The Planetary Science Journal, 2, 78
  * Rufu & Canup (2022) Rufu, R. & Canup, R. M. 2022, ApJ, 928, 123
  * Saillenfest & Lari (2021) Saillenfest, M. & Lari, G. 2021, A&A, 654, A83
  * Saillenfest et al. (2021a) Saillenfest, M., Lari, G., & Boué, G. 2021a, Nature Astronomy, 5, 345
  * Saillenfest et al. (2021b) Saillenfest, M., Lari, G., Boué, G., & Courtot, A. 2021b, A&A, 647, A92
  * Saillenfest et al. (2020) Saillenfest, M., Lari, G., & Courtot, A. 2020, A&A, 640, A11
  * Saillenfest et al. (2019) Saillenfest, M., Laskar, J., & Boué, G. 2019, A&A, 623, A4
  * Saillenfest et al. (2022) Saillenfest, M., Rogoszinski, Z., Lari, G., et al. 2022, A&A, 668, A108
  * Salmon & Canup (2022) Salmon, J. & Canup, R. M. 2022, ApJ, 924, 6
  * Santerne et al. (2019) Santerne, A., Malavolta, L., Kosiarek, M. R., et al. 2019, in revision in Nature Astronomy, arXiv:1911.07355
  * Su & Lai (2020) Su, Y. & Lai, D. 2020, ApJ, 903, 7
  * Su & Lai (2022a) Su, Y. & Lai, D. 2022a, MNRAS, 513, 3302
  * Su & Lai (2022b) Su, Y. & Lai, D. 2022b, MNRAS, 509, 3301
  * Tremaine (1991) Tremaine, S. 1991, Icarus, 89, 85
  * Tremaine et al. (2009) Tremaine, S., Touma, J., & Namouni, F. 2009, AJ, 137, 3706
  * Vanderburg et al. (2016) Vanderburg, A., Becker, J. C., Kristiansen, M. H., et al. 2016, ApJ, 827, L10
  * Vokrouhlický & Nesvorný (2015) Vokrouhlický, D. & Nesvorný, D. 2015, ApJ, 806, 143
  * Wang & Dai (2019) Wang, L. & Dai, F. 2019, ApJ, 873, L1
  * Ward & Hamilton (2004) Ward, W. R. & Hamilton, D. P. 2004, AJ, 128, 2501
  * Williams & Boggs (2016) Williams, J. G. & Boggs, D. H. 2016, Celestial Mechanics and Dynamical Astronomy, 126, 89
  * Wisdom et al. (2022) Wisdom, J., Dbouk, R., Militzer, B., et al. 2022, Science, 377, 1285
  * Xie et al. (2016) Xie, J.-W., Dong, S., Zhu, Z., et al. 2016, Proceedings of the National Academy of Science, 113, 11431
  * Yoder (1995) Yoder, C. F. 1995, in Global Earth Physics: A Handbook of Physical Constants, ed. T. J. Ahrens (American Geophysical Union), 1
  * Zuluaga et al. (2015) Zuluaga, J. I., Kipping, D. M., Sucerquia, M., & Alvarado, J. A. 2015, ApJ, 803, L14

##  Appendix A Correlations between parameters and orbital precession
frequencies

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
$\nu_{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 6 and 7 show the scatter
of the values of $s_{3}$ and $s_{6}$ as a function of all parameters.

Figure 6: 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_{\mathrm{S}}$ as computed from
$10^{6}$ realisations. Figure 7: Same as Fig. 6, but for frequency $s_{6}$.

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_{\mathrm{S}}$. The coefficients obtained are given in Figs.
6 and 7 for each parameter. The correlations between $s_{3}$ and $M_{4}$
($\rho_{\mathrm{S}}=-0.85$) and between $s_{6}$ and $M_{5}$
($\rho_{\mathrm{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
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}}{\mathcal{G}M}(\omega\lambda)^{2}\,.$
(13)

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 8
shows that the Solar System giant planets do fall roughly along this curve.

Figure 8: 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_{\mathrm{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_{\mathrm{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. 13).
For reference, the black dots show the parameters $(\omega\lambda,J_{2})$ of
the Solar System giant planets (Yoder 1995).

For a given frequency $\nu_{j}$ in the orbital precession spectrum of a
planet, the conditions $p\leqslant|\nu_{j}|$ in Eq. (5) corresponds to a
straight line $J_{2}\propto\omega\lambda$, whereas the condition
$|\nu_{j}|\leqslant p\eta/2$ corresponds to a power law
$J_{2}\propto(\omega\lambda)^{5/2}$. Both conditions can be visualised in Fig.
8 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. (13) 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_{\mathrm{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_{\mathrm{min}}/M\approx 6\times 10^{-4}$ for
the most probable value of $s_{6}$ (Fig. 8b), with a dispersion at $95.4\%$
ranging from $m_{\mathrm{min}}/M\approx 3\times 10^{-4}$ (Fig. 8a) to
$m_{\mathrm{min}}/M\approx 1\times 10^{-3}$ (Fig. 8c). These values are
similar to those obtained in Sect. (4) using the parameters of Uranus.

The value of characteristic length $r_{\mathrm{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_{\mathrm{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_{\mathrm{M}}$ only varies by a small
amount when considering realistic values of $J_{2}$ (see the right vertical
axis in Fig. 8).

##  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})\,,$ (14)

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. (14), 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^{\circ}$ (see Table 1 of the main text);
$\Psi_{ik}$ is minimised if $I_{i}$ and $I_{k}$ lie on the same side of
$90^{\circ}$.

Here, we quantify the influence of the unknown longitudes of nodes on the
orbital dynamics of planet HIP 41378 f by drawing 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. 9
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 9: 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}\leqslant
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.

Figure 10 shows the probability of observing five transiting planets or more
in the same experiment as in Fig. 9. 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 10: 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.