© The Authors 2023

1 Introduction

The so-called super-puff exoplanets have moderate masses (typically ≲15 M_⊕) but surprisingly large radii (≳6 R_⊕), giving them extremely low bulk densities (≲0.3 g cm⁻³; 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 (***(eq1)*****${2.1}_{-0.3}^{+0.4}$ Gyr; Lund et al. 2019) and has a long period (542 days) and low equilibrium temperature (300 K). Its low density (0.09 ± 0.02 g cm⁻³) 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¹ 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 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². 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 ± 3 M_⊕), 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 ε ≈ 90° 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° 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 ε ≳ 70°, 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 ε up to ε ≈ 90° 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 Ω its longitude of ascending node. We decompose the inclination dynamics of the planet in a quasi-periodic series truncated to N terms: $$\zeta =\text{sin}\frac{I}{2}\exp \left(i\text{\hspace{0.17em}Ω}\right)={\displaystyle \sum _{j=1}^N{S}_j\exp \left[i{\varphi }_j\left(t\right)\right]},$$(1)

where S_j is a positive real constant, and ***(eq3)*****${\varphi }_j\left(t\right)=v_{j}t+{\varphi }_j^{\left(0\right)}$ 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| ⩾ p, where $$p=\frac{3}{2}\frac{\mathcal{𝒢}{M}_{\star }}{a^3{\left(1-e^2\right)}^{3/2}}\frac{J_2}{\omega \lambda }$$(2)

is the characteristic spin-axis precession rate of the planet. In this expression, 𝒢 is the gravitational constant, M_⋆ 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₂ is the second zonal gravity coefficient of the planet, ω is its spin rate, and λ is its normalised polar moment of inertia. The parameters J₂ and λ 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’ η defined by $$\eta =\frac{1}{2}\frac{m}{M}\frac{r_{\text{M}}^2}{J_2{R}^2},$$(3)

where m is the mass of the moon, M is the mass of the planet, and r_M is the following characteristic length³: $$r_{\text{M}}^2=2\frac{M}{M_{\star }}{J}_2{R}^2{a}^3{\left(1-e^2\right)}^{3/2}.$$(4)

Under the hypothesis that the moon’s mass ratio m/M is small (which does not necessarily imply that η is small), Saillenfest & Lari (2021) show that all resonances with a nodal harmonic having a negative frequency v_j verifying $$p⩽\left|v_j\right|⩽p\frac{\eta }{2}$$(5)

can allow the planet’s obliquity to grow from ε = 0° to ε = 90°. 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 τ of secular oscillations of the moon around its equilibrium ‘Laplace plane’ (see Tremaine et al. 2009) as τ = 2π/κ, where $${\kappa }^2=\frac{9}{4}\frac{M_{\star }}{M}\frac{r_{\text{M}}^3}{a^3{\left(1-e^2\right)}^3}\frac{\mathcal{𝒢}{M}_{\star }}{a^3}.$$(6)

An adiabatic capture in resonance requires that τ is much shorter than the spin-axis precession timescale of the planet T = 2π/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_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_lib near the resonance centre can be computed as T_lib = 2π/μ, where $${\mu }^2={\left(p\prime \right)}^2\left(\frac{{\beta }^2}{{\sin }^2{\epsilon }_0}+\beta \sin {\epsilon }_0\right).$$(7)

In this expression, ε₀ 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 γ = ρ₁/p′ and ß = ρ₂/p′, where $$\begin{array}{l}{\rho }_1=-\left(v_k-2{\displaystyle \sum _{j=1}^N{v}_j{S}_j^2}\right),\hfill \\ {\rho }_2=-S_k\left(2v_k+v_k{S}_k^2-2{\displaystyle \sum _{j=1}^N{v}_j{S}_j^2}\right),\hfill \end{array}$$(8)

and k is the index in Eq. (1) of the considered resonance. T_lib depends on the distance of the moon through p′ and ε₀. However, an upper bound for T_lib is obtained at the time of resonance capture, for which γ^2/3 + β^2/3 = 1 (see Henrard & Murigande 1987; Saillenfest et al. 2019). In this case, p′ is equal to ***(eq11)*****$p^{\prime }={\left({\rho }_1^{2/3}+{\rho }_2^{2/3}\right)}^{3/2}$, and the planet’s obliquity at the centre of the resonance is $$\cos {\epsilon }_0=\gamma -{\gamma }^{1/3}+\sqrt{{\gamma }^2+{\gamma }^{2/3}-{\gamma }^{4/3}}.$$(9)

Thanks to these expressions, we can compute T_lib from Eq. (7) as a mere function of the planet’s orbital dynamics in Eq. (1).

Fig. 1 Level curves of the planet’s spin-axis precession rate (adapted from Saillenfest et al. 2022 for a moon mass-parameter η = 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.

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 \left(i{\text{Ω}}_k\right),$$(10)

where I_k is the orbital inclination of planet k and Ω_k is its longitude of ascending node. The Lagrange-Laplace system gives the linear equation of motion $$\frac{\text{d}\zeta }{\text{d}t}=iB\zeta ,$$(11)

in which ζ is the vector containing the ζ_k variable of all planets and Β 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\left(t\right)={\displaystyle \sum _{j=1}^{N_{\text{p}}}{S}_j\exp \left[i\left(v_{j}t+{\varphi }_j^{\left(0\right)}\right)\right]},$$(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 v_j are the eigenvalues of the matrix Β. As Β only depends on the masses and semi-major axes of the planets, this is also the case of the frequencies v_j. Because of the conservation of total angular momentum, one of the frequencies v_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 v_j, the first step is to draw the masses and semimajor 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 v_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 Β 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 v_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 v_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 v_j of matrix Β 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 HIP41378 system

HIP 41378 is bright F-type star⁴ 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 ***(eq16)*****${260}_{-60}^{+160}$ 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 ± 5 M_⊕ 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 ± 0.1 R_⊕ and mass 12 ± 3 M_⊕, giving it a bulk density of 0.09 ± 0.02 g cm⁻³ (Santerne et al. 2019). Under the ring hypothesis, current data suggests a planet with radius R = 3.7 ± 0.3 R_⊕ surrounded by a ring with radius 2.6 ± 0.2 R and inclination 25 ± 4° from the sky plane (Akinsanmi et al. 2020). This new planetary radius yields a bulk planet density of 1.2 ± 0.4 g cm⁻³, similar to that of Uranus. The hypothetical equatorial ring provides an indirect measure of the obliquity of the planet, namely⁵ ε = 92 ± 7°.

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_⊕ and a mass measurement. From the NASA Exoplanet Archive⁶ on date 2022-11-23, we obtain a central mass value of 12.7 M_⊕ and a standard deviation of 6.0 M_⊕. 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_⊕) 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₁ is the precession frequency of planets 1 and 3 about their total angular momentum vector; s₂ is the precession frequency of the low-mass planet 2 under the action of planets 1 and 3; s₃ is the slow rigid precession of the inner and outer subsystems (planets 1–2–3 and 4–5–6); s₄ is the precession frequency of planets 4 and 5 about their total angular momentum vector; s₅ is identically zero; s₆ 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⁶ 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 Β 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 ζ_k(t = 0) designed to magnify the specific term we are looking for. For instance, the frequency s₃ would appear as strongly dominant for all planets if we set ζ_k(t = 0) = 0 for k = {1, 2, 3} and ***(eq17)*****${\zeta }_k\left(t=0\right)=\sqrt{2}/2$ for k = {4, 5, 6}. Then, one may identify s₂ as the dominant term in the solution of planet 2 by setting ***(eq18)*****${\zeta }_2\left(t=0\right)=\sqrt{2}/2$ and ζ_k(t = 0) = 0 for k ≠ 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⁶ realisations of the system. Frequency s₄ has a broad distribution due to the large uncertainties in the masses of planets d and e. Frequency s₃, 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₃ results to be essentially due to the uncertainty in the mass of planet d (see Appendix A).

As illustrated in Fig. 2, frequency s₃ 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₆. 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₆ has a relatively peaked distribution. Its most probable value is −136″ yr⁻¹, with 68.3% occurrences within [−181, −97]″ yr⁻¹, 95.4% occurrences within [−241, −65]″ yr⁻¹, and 99.7% occurrences within [−405, −35]″ yr⁻¹. As shown in Appendix A, the value of s₆ is essentially set by the mass of the perturbing planet e, with a Spearman correlation coefficient ρ_S ≈ −0.8. The value of

s₆ is only weakly (|ρ_S| ≲ 0.3) correlated with the parameters of planet f itself. This low correlation allows us to investigate different values for the frequency s₆ 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).

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 ζk are set from the nominal inclinations Ik (all assumed to be Ik ⩽ 90°) and random longitudes of node Ωk in a 0.2°-wide interval; this choice is commented in Sect. 5.

Fig. 3 Probability density of the inclination proper modes of the HIP 41378 system. Histograms are built from 106 realisations of the Lagrange-Laplace system with the mass and semi-major axis uncertainties in Table 1. The histogram for frequency s4 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 s3 peaks above the top border of the figure. Frequency s5 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).

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₂, and the product ωλ. 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 ***(eq19)*****${3.7}_{-0.2}^{+0.3}$ R_⊕. Hence, we adopt the value R = 3.7 R_⊕ below as our normalising radius.

For given values of the parameters J₂ and ωλ, 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₂ and ωλ are completely unknown for exoplanets, we know that they are related, and in first approximation J₂ ∝ ω² (planets spinning faster are more flattened; see e.g. Chandrasekhar 1969). For a given moon mass m, the condition |v_j| ⩽ pη/2 in Eq. (5) corresponds to a power law J₂ ∝ ω^5/2. Because of the coincidental near match between these two exponents (2 and 5/2), our total ignorance of J₂ and ωλ does not affect much our estimate of m_min: we may just set J₂ and ωλ 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₂ ∝ ω². 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₂ and ωλ of Uranus (see e.g. Yoder 1995).

Independently of the resonance considered, Eq. (5) can be fulfilled only if the mass parameter η of the moon is η ⩾ 2. Using the J₂ value of Uranus, this condition translates into m/M ⩾ 1.2 × 10⁻⁴. 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₂ and λω of Uranus (the tics start at 1.2 × 10⁻⁴ and go from right to left).

The characteristic spin-axis precession rate of HIP 41378 f computed from Eq. (2) is p ≈ 25″ yr⁻¹. According to the left inequality in Eq. (5), this value almost certainly rules out a resonance with frequency s₃, because frequency s₃ sharply peaks at s₃ = −15.7″ yr⁻¹ (see Fig. 3). The fact that p > |s₃| means that the s₃ 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₆, 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₆ resonance requires a moon with minimum mass ratio ranging between about 2 × 10⁻⁴ and 10 × 10⁻⁴. 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₂ and ωλ of Uranus are assumed for HIP 41378 f, the value of frequency ***(eq20)*****$s_6=-{136}_{-269}^{+101}″$ yr⁻¹ obtained in Sect. 3 translates into a minimum moon mass ***(eq21)*****$m_{\min }/M=6_{-5}^{+13}\times {10}^{-4}$ (3σ 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 × 10⁻⁴, and the summed masses of the largest moons of Jupiter and Uranus yield ratios of about 2 × 10⁻⁴ and 1 × 10⁻⁴, 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⁻⁴. 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 ≈ 0.5 r_M to 1 r_M. Using the J₂ value of Uranus, Eq. (4) gives a characteristic length r_M ≈ 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_M is proportional to ***(eq22)*****$J_2^{1/5}$, other realistic values of J₂ may change these distances by a small amount (see discussion in Appendix B).

The HIP 41378 system is ***(eq23)*****${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 cm yr⁻¹ 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₂/Q ≈ 3 × 10⁻⁵ for a moon mass m/M = 2 × 10⁻⁴, and higher than k₂/Q ≈ 6 × 10⁻⁶ for a moon mass m/M = 10⁻³. For comparison, the value measured for Jupiter’s satellite Io is k₂/Q = 1.102 ± 0.203 × 10⁻⁵ (Lainey et al. 2009), and the value measured globally for Saturn’s main satellites is k₂/Q = 1.59 ± 0.74 × 10⁻⁴ (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₂ and ωλ, 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₆ 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 Ω_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 Ω_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 Ω_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 inner-most planet may possibly transit 100% of the time, even if we set all the Ω_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 HIP41378 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 Ω_k in the sky plane. As shown in Appendix C, allowing for just a few degrees dispersion in Ω_k can increase the amplitude Sj 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 Ω_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 Ω_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₆ resonance angle when the separatrix appears is T_lib ≈ 547 000 years. This value is much smaller than the age of the system (***(eq24)*****${2.1}_{-0.3}^{+0.4}$ 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₆ 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 × 10⁻⁴ (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₂/Q ≈ 10⁻⁵. 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⁻⁴, with random argument of pericentre and longitude of ascending node. As expected, Fig. 4 shows that the adiabatic capture and tilting in resonance s6 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 τ defined in Eq. (6). Here, one obtains τ ≈ 100 yr, which means that the eccentricity increase is extremely fast compared to the planet’s spin-axis precession timescale (T ≈ 52000 yr), to the oscillations of the planet inside the resonance (T_lib ≈ 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σ with the obliquity ε = 92 ± 7° 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°.

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 × 10−4. The moon migrates away at constant rate emulating a tidal parameter k2/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 sk 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.

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₆ 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).

Fig. 5 Probability of underestimating the density of a transiting exo-planet 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.53 (exact face-on configuration). In case it possesses a ring, planet HIP 41378 f would have q ≈ 13 (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⁸ (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⁶ 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₂ and ωλ of the planet, but thanks to the approximate relation J₂ ∝ ω², 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 ~ 10⁻⁴ 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₂/Q for the planet. Expected values are of the order of 10⁻⁵ to 10⁻⁴ 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 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 Ω_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 Ω_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 ≳ 2 × 10⁻⁴ to 10 × 10⁻⁴, 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₂/Q larger than 2 × 10⁻⁵ (for the smallest possible moon) or larger than about 6 × 10⁻⁶ (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 (***(eq26)*****${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 ε reaches values in the interval [70°, 110°], and its moon becomes unstable (Tremaine et al. 2009; Saillenfest & Lari 2021). Due to the short instability timescale of the exomoon (τ ≈ 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 ε ≈ 90^°, 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.

Appendix A Correlations between parameters and orbital precession frequencies

Fig. A.1 Scatter plot of the proper frequency s3 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 ρS as computed from 106 realisations.

Fig. A.2 Same as Fig. A.l, but for frequency s6.

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 Β (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₃ and s₆ 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₃ and s₆ as a function of all parameters.

From these scatter plots, it is visible that the value of s₃ mostly depends on M₄, while the value of s₆ mostly depends on M₅. The strength of these correlations can be quantified by Spearman’s correlation coefficient ρ_S. The coefficients obtained are given in Figs. A.1 and A.2 for each parameter. The correlations between s₃ and M₄ (ρ_S = −0.85) and between s₆ and M₅ (ρ_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₂ and ωλ of the planet. In this section, we consider them as free parameters and study their influence on our results.

The values of J₂ and ωλ are related and depend on the interior properties of the planet. In the simplest case of a homogeneous planet, J₂ and ω are linked through the law of Maclaurin’s ellipsoid (see e.g. Chandrasekhar 1969), while λ = 2/5. This law simplifies to J₂ ∝ ω² for a nearly spherical planet, and it can be rewritten as $$J_2\approx \frac{25}{8}\frac{R^3}{𝒢M}{\left(\omega \lambda \right)}^2.$$(B1)

Fig. B.1 Minimum moon mass needed for HIP 41378 f to undergo the full tilting mechanism through a resonance with s6. Each panel corresponds to a given possible value of s6 (labels). The minimum moon mass mmin (colour scale) is obtained from Eq. (5) as a function of the unknown parameters J2 and ωλ of HIP 41378 f. The white curves highlight three mass levels, namely mmin/M = 3 × 10−4, 6 × 10−4, and 10 × 10−4. The grey region is forbidden by the left inequality in Eq. (5). The black dotted line shows the relation between J2 and ωλ obtained for a homogeneous planet (see Eq. B.1). For reference, the black dots show the parameters (ωλ, J2) of the Solar System giant planets (Yoder 1995).

Even though planets are not homogeneous, this approximate relation gives an idea of where to look for realistic combinations of parameters in the plane (ωλ, J₂) 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 v_j in the orbital precession spectrum of a planet, the conditions p ⩽ |v_j| in Eq. (5) corresponds to a straight line J₂ ∝ ωλ, whereas the condition |v_j| ⩽ pη/2 corresponds to a power law J₂ ∝ (ωλ)^5/2. Both conditions can be visualised in Fig. B.1 for planet HIP 41378 f using three different values of the frequency v_j (most probable value of s₆ 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₂ and ωλ (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₂ and ωλ. Whatever realistic values are chosen, we obtain a minimum moon mass m_min/M ≈ 6 × 10⁻⁴ for the most probable value of s₆ (Fig. B.1b), with a dispersion at 95.4% ranging from m_min/M ≈ 3 × 10⁻⁴ (Fig. B.1a) to m_min/M ≈ 1 × 10⁻³ (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₂, but not on ωλ (see Eq. 4). The distance covered by the migrating moon depends therefore on the value chosen for J₂. Yet, this dependence is not very steep (***(eq29)*****$r_{\text{M}}\propto J_2^{1/5}$). Moreover, extreme values of J₂, 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₂ (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 Ψ_ik of two planets i and k can be written as $$\cos {\text{Ψ}}_{ik}=\cos I_i\cos I_k+\sin I_i\sin I_k\cos \left({\text{Ω}}_i-{\text{Ω}}_k\right),$$(C.1)

where I and Ω 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 Ω_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 Ψ_ik = |I_i − I_k| reached for Ω_i − Ω_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); Ψ_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 randomly the longitudes of nodes of all planets in a given interval ΔΩ₀ and building a histogram of the amplitudes S_j obtained. The result is shown in Fig. C.1 for an interval ΔΩ₀ ranging from 0° to 2°. As expected, the terms s₁ and s₂ have very small amplitudes; they contribute negligibly to the dynamics of planet HIP 41378 f. The dominant terms of the dynamics are s₃ and s₆, with the qualitative roles described in Sect. 3. The amplitudes of all terms are generally minimum for ΔΩ₀ = 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₄ term, for which a dispersion as small as ΔΩ₀ = 1° can make the amplitude increase by a factor 1000 or so.

Fig. C.1 Possible values of the amplitudes Sj for the inclination terms of planet HIP 41378 f. The frequency of each term is labelled. The amplitude of the zero-frequency term s5 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 Ik are assumed to be Ik ⩽ 90°. The unknown longitudes of node Ωk of the planets in the sky plane are drawn from a uniform random distribution in an interval ΔΩ0 common for all planets. For a given value of the range ΔΩ0 (vertical axis), a histogram of each amplitude is built from 106 realisations of the Lagrange-Laplace system (colour scale). A bin is coloured grey if no occurrence is found among our 106 realisations. The scale is the same for all graphs to highlight the disparity of amplitude between the different terms.

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. ΔΩ₀ = 0), then five planets or more transit about 30% of the time. If we allow a dispersion of ΔΩ₀ = 1°, this fraction can be reduced to less than 5%.

Fig. C.2 Probability of observing five transiting planets in the HIP 41378 system over the planets’ precession cycles. The unknown longitudes of node Ωk of the planets in the sky plane are drawn from a uniform random distribution in an interval ΔΩ0 common for all planets. For a given value of the range ΔΩ0 (vertical axis), a histogram of the transit probability is built from 106 realisations of the Lagrange–Laplace system (colour scale). A bin is coloured grey if no occurrence is found among our 106 realisations.

References

All Tables

All Figures

Fig. 1 Level curves of the planet’s spin-axis precession rate (adapted from Saillenfest et al. 2022 for a moon mass-parameter η = 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.

In the text

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 ζk are set from the nominal inclinations Ik (all assumed to be Ik ⩽ 90°) and random longitudes of node Ωk in a 0.2°-wide interval; this choice is commented in Sect. 5.

In the text

Fig. 3 Probability density of the inclination proper modes of the HIP 41378 system. Histograms are built from 106 realisations of the Lagrange-Laplace system with the mass and semi-major axis uncertainties in Table 1. The histogram for frequency s4 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 s3 peaks above the top border of the figure. Frequency s5 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).

In the text

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 × 10−4. The moon migrates away at constant rate emulating a tidal parameter k2/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 sk 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.

In the text

Fig. 5 Probability of underestimating the density of a transiting exo-planet 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.53 (exact face-on configuration). In case it possesses a ring, planet HIP 41378 f would have q ≈ 13 (Akinsanmi et al. 2020).

In the text

	Fig. A.1 Scatter plot of the proper frequency s3 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 ρS as computed from 106 realisations.
In the text

	Fig. A.2 Same as Fig. A.l, but for frequency s6.
In the text

Fig. B.1 Minimum moon mass needed for HIP 41378 f to undergo the full tilting mechanism through a resonance with s6. Each panel corresponds to a given possible value of s6 (labels). The minimum moon mass mmin (colour scale) is obtained from Eq. (5) as a function of the unknown parameters J2 and ωλ of HIP 41378 f. The white curves highlight three mass levels, namely mmin/M = 3 × 10−4, 6 × 10−4, and 10 × 10−4. The grey region is forbidden by the left inequality in Eq. (5). The black dotted line shows the relation between J2 and ωλ obtained for a homogeneous planet (see Eq. B.1). For reference, the black dots show the parameters (ωλ, J2) of the Solar System giant planets (Yoder 1995).

In the text

Fig. C.1 Possible values of the amplitudes Sj for the inclination terms of planet HIP 41378 f. The frequency of each term is labelled. The amplitude of the zero-frequency term s5 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 Ik are assumed to be Ik ⩽ 90°. The unknown longitudes of node Ωk of the planets in the sky plane are drawn from a uniform random distribution in an interval ΔΩ0 common for all planets. For a given value of the range ΔΩ0 (vertical axis), a histogram of each amplitude is built from 106 realisations of the Lagrange-Laplace system (colour scale). A bin is coloured grey if no occurrence is found among our 106 realisations. The scale is the same for all graphs to highlight the disparity of amplitude between the different terms.

In the text

Fig. C.2 Probability of observing five transiting planets in the HIP 41378 system over the planets’ precession cycles. The unknown longitudes of node Ωk of the planets in the sky plane are drawn from a uniform random distribution in an interval ΔΩ0 common for all planets. For a given value of the range ΔΩ0 (vertical axis), a histogram of the transit probability is built from 106 realisations of the Lagrange–Laplace system (colour scale). A bin is coloured grey if no occurrence is found among our 106 realisations.

In the text