Stability and Detectability of Exomoons Orbiting HIP 41378 f, a Temperate Jovian Planet with an Anomalously Low Apparent Density

The discovery and characterization of extrasolar moons (i.e., "exomoons") will play an important role in developing a more complete understanding of planetary systems. As remnants of planet formation, moons may provide useful information about the formation and evolution of exoplanets (e.g., Sasaki et al. 2010; Heller et al. 2014; Spalding et al. 2016; Ronnet et al. 2018). Moon-planet interactions may also be crucial for the long-term stability of exoplanet habitability; for example, moons may stabilize planetary obliquity (Laskar et al. 1993; Lissauer et al. 2012) and drive tidal heating (e.g., Piro 2018a), thereby significantly affecting planetary climate. Moreover, sufficiently heated exomoons may themselves be capable of hosting life, providing an alternative path to habitability beyond conventional habitable-zone terrestrial planets (e.g., Heller & Barnes 2013; Heller et al. 2014).

Despite the large number of confirmed exoplanets surveyed to date (>5000), there has yet to be a confirmed detection of an exomoon, largely due to the minuscule signals that putative moons imprint on measurements of their exoplanet hosts. Though tentative evidence for exomoons has recently been claimed (Teachey & Kipping 2018; Kipping et al. 2022), robustly detecting exomoons is likely beyond the reach of current technology. Nonetheless, theoretical constraints on exomoon formation and evolution can provide valuable insight into plausible satellite-hosting exoplanets, and aid in the interpretation of future high-precision observations. For example, numerical N-body simulations have set limits on exomoon eccentricities and semimajor axes that are plausibly stable over long timescales (e.g., Domingos et al. 2006; Payne et al. 2013; Rosario-Franco et al. 2020). Moreover, previous studies have investigated how tidal interactions between bodies can result in the radial migration of an exomoon, leading to its eventual escape or tidal disruption by the planet (e.g., Barnes & O'Brien 2002; Sasaki et al. 2012; Sasaki & Barnes 2014; Piro 2018a; Sucerquia et al. 2019; Quarles et al. 2020; Jagtap et al. 2021), as well as the effects of planetary migration on exomoon stability (e.g., Alvarado-Montes et al. 2017; Sucerquia et al. 2020).

In this work, we investigate the plausibility of exomoons orbiting the temperate low-mass Jovian planet HIP 41378 f. HIP 41378 is a bright (V ≈ 8.9), slightly evolved late F-type star (T_eff = 6320 K, $\mathrm{log}g=4.294$, [Fe/H] = − 0.10, M_⋆ = 1.16 M_⊙; Santerne et al. 2019) hosting five known transiting exoplanets (Vanderburg et al. 2016a; Berardo et al. 2019). The system was initially observed over the course of about 75 days during Campaign 5 (C5) of the extended Kepler mission (K2; Borucki et al. 2010; Howell et al. 2014), after which the two innermost planets were classified as sub-Neptunes (R_p,b ≈ 2.9 R_⊕, R_p,c ≈ 2.6 R_⊕) orbiting in a near 2:1 resonance at approximately 15.6 and 31.7 days (Vanderburg et al. 2016a). However, only single transits were initially observed for the three outer planets (R_p,d ≈ 4.0 R_⊕, R_p,e ≈ 5.5 R_⊕, R_p,f ≈ 10 R_⊕), and hence their orbital periods and ephemerides could not be well constrained (Vanderburg et al. 2016a).

More precise limits on the orbits of planets d and f were later derived after one additional transit was observed for each planet during C18 of the K2 mission (Becker et al. 2019; Berardo et al. 2019). A subsequent third transit observation of planet f with the ground-based Next Generation Transit Survey (NGTS; Wheatley et al. 2018) confirmed its orbital period and suggested the presence of transit timing variations (TTVs; Bryant et al. 2021). Furthermore, radial velocity monitoring of the host star over hundreds of epochs was used to constrain the masses of all of the planets in the system, providing new insights into their physical properties. Santerne et al. (2019) conducted a joint analysis using the two campaigns of K2 photometry together with 464 epochs of radial velocity measurements from the High Accuracy Radial velocity Planet Searcher (HARPS) spectrograph at the 3.6-m ESO telescope in Chile, the HARPS-North spectrograph at the 3.58-m TNG telescope on the Canary Islands, Spain, the High Resolution Echelle Spectrometer (HIRES) at the 10-m Keck telescope on Mauna Kea, Hawaii, and the Prime Focus Spectrograph at the 8.2-m Subaru telescope on Mauna Kea, Hawaii. They reported updated planet radii, masses, bulk densities, and orbital parameters for the five transiting planets, as well as a tentative detection of a sixth nontransiting planet. However, we note that an updated analysis of the data likely rules out the hypothesis of a tentative sixth planet in the HIP 41378 system (A. Santerne 2023, private communication). A summary of the most recent planet parameters from Santerne et al. (2019) is shown in Table 1 ²⁴

Table 1. HIP 41378 System Parameters from Santerne et al. (2019)

Parameter	Median and 68.3% Credible Interval
Planet Parameters	HIP 41378 b	HIP 41378 c	HIP 41378 d
Orbital period, P (days)	15.57208 ± 2 × 10−5	31.70603 ± 6 × 10−5	278.3618 ± 5 × 10−4
Eccentricity, e	0.07 ± 0.06	${0.04}_{-0.03}^{+0.04}$	0.06 ± 0.06
Semimajor axis, a (au)	0.1283 ± 1.5 × 10−3	0.2061 ± 2.4 × 10−3	0.88 ± 0.01
Inclination, i (°)	88.75 ± 0.13	${88.477}_{-0.061}^{+0.035}$	89.80 ± 0.02
Radius, Rp (R⊕)	2.595 ± 0.036	2.727 ± 0.060	3.54 ± 0.06
Mass, Mp (M⊕)	6.89 ± 0.88	4.4 ± 1.1	<4.6 b
Bulk density, ρp (g cm−3)	2.17 ± 0.28	1.19 ± 0.30	<0.56 b
Equilibrium temperature a , Teq (K)	${959}_{-5}^{+9}$	${757}_{-4}^{+7}$	${367}_{-2}^{+3}$
Insolation flux, S (S⊕)	${140}_{-3}^{+5}$	${54}_{-1}^{+2}$	${3.01}_{-0.06}^{+0.11}$
Planet parameters	HIP 41378 e	HIP 41378 f
Orbital period, P (days)	369 ± 10	542.07975 ± 1.4 × 10−4
Eccentricity, e	0.14 ± 0.09	${0.004}_{-0.003}^{+0.009}$
Semimajor axis, a (au)	${1.06}_{-0.02}^{+0.03}$	1.37 ± 0.02
Inclination, i (°)	${89.84}_{-0.03}^{+0.07}$	${89.971}_{-0.008}^{+0.01}$
Radius, Rp (R⊕)	4.92 ± 0.09	9.2 ± 0.1
Mass, Mp (M⊕)	12 ± 5	12 ± 3
Bulk density, ρp (g cm−3)	0.55 ± 0.23	0.09 ± 0.02
Equilibrium temperature*, Teq (K)	335 ± 4	${294}_{-1}^{+3}$
Insolation flux, S (S⊕)	2.1 ± 0.1	${1.24}_{-0.02}^{+0.05}$
Stellar parameters	HIP 41378
Effective temperature, Teff (K)	${6320}_{-30}^{+60}$
Surface gravity, $\mathrm{log}g$ (cgs)	4.294 ± 0.006
Metallicity, [Fe/H] (dex)	−0.10 ± 0.07
Mass, M⋆ (M⊙)	1.16 ± 0.04
Radius, R⋆ (R⊙)	1.273 ± 0.015
Age, τ⋆ (Gyr)	3.1 ± 0.6
Distance from Earth, D (pc)	103 ± 2

Notes.

^a Assumes zero albedo and full heat redistribution. ^b 95% credible upper limit.

Download table as:  ASCIITypeset image

The largest and outermost planet in the system, HIP 41378 f, orbits near the inner edge of the system's habitable zone with a period of about 1.5 years. Mysteriously, the planet has an inferred radius of 9.2 R_⊕ but it has a mass of only 12 M_⊕, making it one of the lowest bulk density planets currently known (ρ = 0.09 g cm⁻³; Santerne et al. 2019). It has been proposed that such low bulk densities in exoplanets could be explained by atmospheric outflows entraining optically thick photochemical hazes, which act to increase the planet's apparent radius in transit (Wang & Dai 2019; Gao & Zhang 2020; Ohno & Tanaka 2021). However, this mechanism may not be compatible with HIP 41378 f because of the planet's relatively high mass and cool equilibrium temperature (e.g., see Equation (34) of Ohno & Tanaka 2021). Moreover, Belkovski et al. (2022) showed that HIP 41378 f would require a very high envelope mass fraction (≳75%) to remain in hydrostatic balance, in tension with the core accretion paradigm of planet formation.

An alternative explanation for the apparently low density of HIP 41378 f is that the planet itself is significantly smaller, but is accompanied by a system of optically thick circumplanetary rings (Zuluaga et al. 2015; Akinsanmi et al. 2020; Piro & Vissapragada 2020). Such a ring system, if viewed at a sufficiently high inclination, would effectively increase the observed transit depth while allowing the planet to have a bulk density more typical of sub-Neptunes. For example, Akinsanmi et al. (2020) demonstrated that observations of HIP 41378 f are consistent with a 3.7 R_⊕ (ρ = 1.2 g cm⁻³) planet with opaque rings extending from 1.05 to 2.6 planet radii. In addition, the presence of optically thick rings potentially causes a featureless transmission spectrum (Ohno & Fortney 2022; Ohno et al. 2022). Interestingly, our recent observations of HIP 41378 f taken with the Widefield Camera 3 (WFC3) aboard the Hubble Space Telescope (HST) reveal a featureless near-infrared (1.1–1.7 μm) transmission spectrum (HST-GO Program 16267; Alam et al. 2022). Within the precision of these observations, we currently cannot distinguish between circumplanetary rings (${\chi }_{r}^{2}\approx 1.03$), high-altitude photochemical hazes (${\chi }_{r}^{2}\approx 0.97$), or a high-metallicity atmosphere (${\chi }_{r}^{2}\approx 1.84$). Future spectral observations at mid-infrared wavelengths (e.g., with JWST/MIRI) may be able to distinguish between these scenarios (Alam et al. 2022).

The possibility of a circumplanetary ring system around HIP 41378 f naturally raises the question of whether the planet could also host moons. For example, Saillenfest et al. (2023) concluded that the migration of a former moon is a viable formation pathway for the proposed ring and tilt of HIP 41378 f. More generally, moons and circumplanetary rings are ubiquitous in the outer solar system, and earlier studies have suggested that rings are originally formed from moons (e.g., Canup 2010), or vice versa (e.g., Crida & Charnoz 2012).

Because the Hill radius scales linearly with the distance of a planet from its host star, the large separation of HIP 41378 f from its host star relative to the population of hot Jupiters (P ≲ 10 days) makes it a favorable environment for exomoons. While HIP 41378 f itself is likely not habitable due to its high gas fraction, its temperate orbital separation raises the possibility that it hosts habitable exomoons. We note that the habitability of icy moons in the outer solar system is still an open question. Though the detailed level of characterization required to assess exomoon habitability, or even confirm exomoon detections at all, is beyond current instrumental sensitivity, it may be prudent to consider the extent to which the presence of exomoons could manifest as a source of uncertainty in measurements of cool gaseous exoplanets, especially as we look toward future high-precision facilities such as ground-based extremely large telescopes (ELTs) and large space-based UV/optical/IR observatories (e.g., Dalba et al. 2015; The LUVOIR Team 2019).

Here we assess the stability of exomoons orbiting HIP 41378 f using a suite of numerical N-body simulations and models of tidal evolution. Then, from the results of our stability analysis, we simulate the observable impact of a large exomoon on the white light curve and infrared transmission spectrum of HIP 41378 f. The rest of this paper is organized as follows: in Section 2 we describe our theoretical methods used to constrain the allowable range of parameters for a stable moon orbiting HIP 41378 f. In Section 3, we present the results from our analysis and discuss the future observational implications. We conclude with a summary in Section 4.

We first consider whether it is feasible that HIP 41378 f could host an exomoon over the system's lifetime. Here we outline a theoretical approach to investigating this possibility. First, we apply a suite of dynamic simulations to assess the long-term orbital stability of satellites around HIP 41378 f. We then separately consider the effects of tidal friction on the orbits of satellites and whether tidal migration can further destabilize the system. Throughout this section, we assume that HIP 41378 f has properties consistent with measurements from Santerne et al. (2019), and we do not self-consistently include circumplanetary rings in our simulations.

2.1.  N-body Simulations

2.1.1. Three-body Simulation

Because HIP 41378 is a multiplanet system, gravitational interactions between neighboring planets may affect the long-term stability of a satellite orbiting planet f. Before considering the effects of additional planets, we start with an idealized three-body system with only one satellite orbiting planet f, which is in orbit around HIP 41378. To simulate the dynamical evolution of the system, we use the REBOUND N-body code (Rein & Liu 2012) and the 15th-order "IAS15" integrator with adaptive time stepping (Rein & Spiegel 2015).

In our simulations, we fix the stellar mass and the mass and semimajor axis of planet f to the values from Santerne et al. (2019), provided in Table 1. Given a vast range of possible satellite configurations and properties, which would be computationally prohibitive to fully explore, we chose to consider a satellite with the largest satellite-to-planet mass ratio observed in the solar system (i.e., the Moon–Earth system). This choice allows us to explore a physically motivated scenario that maximizes the potential signal produced by the satellite—larger satellites may be exceedingly rare, based on the population of hundreds of moons in the solar system, and smaller moons are far less likely to be detected. Therefore, as an extreme case, we consider a satellite of HIP 41378 f with a mass of M_s = 0.15 M_⊕ to be consistent with the satellite-to-planet mass ratio for the Earth-Moon system (M_M/M_⊕ ≈ 0.0123). We prescribe the satellite an Earth-like rocky composition (ρ_s = 5.5 g cm⁻³), yielding a radius of R_s = 0.53 R_⊕ (approximately the size of Mars).

For simplicity, the moon is initiated on a circular orbit around planet f, then, following the results of previous studies estimating stability limits of hierarchical systems (e.g., Rosario-Franco et al. 2020; Jagtap et al. 2021; Quarles et al. 2021), each simulation is run for a duration of up to 10⁵ yr. A simulation is halted and marked as "unstable" if either of the following criteria are met:

• 1.  
  the satellite collides with the planet, ²⁵ or
• 2.  
  the satellite escapes the gravitational influence of the planet by crossing outside of the planet's Hill radius, defined as

  $$\begin{eqnarray}&&{R}_{{\rm{H}}}={a}_{{\rm{p}}}{\left(\displaystyle \frac{{M}_{{\rm{p}}}}{3{M}_{\star }}\right)}^{1/3},\end{eqnarray} \tag{ 1 }$$

  where a_p is the planet's semimajor axis, and M_p and M_⋆ are the masses of the planet and star, respectively.

As a test case, we first reproduced the results of Domingos et al. (2006) and Rosario-Franco et al. (2020) by assessing the stability of the hypothetical satellite as a function of the initial planet-moon semimajor axis, a_s, and the eccentricity, e_f, of the host planet's orbit. Because the eccentricity of planet f is well constrained by observations (see Table 1), the main purpose of this setup was to benchmark our N-body code and confirm that it generates expected results—we provide a more thorough description of this procedure in Appendix A. In summary, we qualitatively recover a similar trend in the stability of the exomoon as a function of (a_s, e_f) as previous studies, and the stability limit we recover for circular planetary orbits, a_crit ≈ 0.41R_H, is consistent with published values (Rosario-Franco et al. 2020).

Next, we investigate how the satellite's orbital stability depends on its initial inclination. Our setup is similar to that of our benchmark simulations and Domingos et al. (2006) and Rosario-Franco et al. (2020), but here we fix the eccentricity of planet f to its measured 95% confidence upper limit of e_f = 0.035 (Santerne et al. 2019), while varying the satellite's orbital inclination, i_s, and semimajor axis, a_s. We also tested a zero-eccentricity case and found similar results (in general, lower-eccentricity systems are expected to be more stable regardless of satellite inclination; see Appendix A).

We constrained our simulations to vary i_s from $\cos {i}_{{\rm{s}}}=0$ (aligned with the planet's reference plane) to $\cos {i}_{{\rm{s}}}=1$ (perpendicular to the planet's reference plane), and a_s from 0.1R_H to 0.8R_H. Here, R_H is the Hill radius defined in Equation (1) (R_H ≈ 0.03 au ≈76R_p). The planet and satellite are both initialized with their ascending node longitudes and pericenter arguments set to zero (Ω = ω = 0°), and the mean anomaly of the planet is also set to zero (${{ \mathcal M }}_{{\rm{f}}}=0^\circ$).

Then, to efficiently explore the moon's stability within our defined parameter space, we run a suite of N-body simulations using a quadtree data structure implemented in astroQTpy ²⁶ (Harada 2023), which is based on the transit injection and recovery quadtree algorithm described by Ment & Charbonneau (2023). We initially create a 4 × 4 stability grid by evenly subdividing the parameter space in $\cos {i}_{{\rm{s}}}$ and a_s. Within each grid cell, we run 25 independent N-body simulations ²⁷ with the initial conditions $\cos {i}_{{\rm{s}}}$ and a_s drawn from uniform distributions. The satellite's initial mean anomaly, ${{ \mathcal M }}_{{\rm{s}}}$, is also chosen randomly from a uniform distribution for each simulation (see Table 2). Each simulation is run forward for up to 10⁵ yr, and then the fractional stability of each cell, f_stab, is calculated as the fraction of simulations (out of 25) that satisfy our stability criteria defined above.

Table 2. Parameter Distributions Used for N-body Simulations

Parameter	Distribution	Distribution
	(3-body)	(4-body)
f
Mass, Mf (M⊕)	12.3	${ \mathcal T }(12.3,3.1,0.0,\infty )$
Orbital period, Pf (days)	542.0798	${ \mathcal N }(542.0798,0.0001)$
Eccentricity, ef	0.035	${ \mathcal T }(0.005,0.005,0.0,1.0)$
Mean anomaly, ${{ \mathcal M }}_{{\rm{f}}}$	${ \mathcal U }(0,2\pi )$	${ \mathcal U }(0,2\pi )$
Argument of periastron, ωf (°)	0.0	${ \mathcal T }(231,120,0.0,360)$
Inclination, if (°)	90.0	${ \mathcal T }(89.97,0.01,0.0,90.0)$
Planet e
Mass, Me (M⊕)	N/A	${ \mathcal T }(12,5,0.0,\infty )$
Orbital period, Pe (days)	N/A	[344, 394)
Eccentricity, ee	N/A	[0.0, 0.3)
Mean anomaly, ${{ \mathcal M }}_{{\rm{e}}}$	N/A	${ \mathcal U }(0,2\pi )$
Argument of periastron, ωe (°)	N/A	${ \mathcal T }(114,55,0,360)$
Inclination, ie (°)	N/A	${ \mathcal T }(89.84,0.05,0.0,90.0)$
Satellite
Mass, Ms (M⊕)	0.15	0.15
Semimajor axis, as (RH)	[0.1, 0.8)	${ \mathcal U }(0.1,0.4)$
Eccentricity, es	0.0	0.0
Mean anomaly, ${{ \mathcal M }}_{{\rm{s}}}$	${ \mathcal U }(0,2\pi )$	${ \mathcal U }(0,2\pi )$
Argument of periastron, ωs (°)	0.0	0.0
a Cosine of inclination, $\cos ({i}_{{\rm{s}}})$	[0, 1)	${ \mathcal U }(0,1)$

Note. ${ \mathcal N }(\mu ,\sigma )$ is a normal distribution with mean μ and standard deviation σ; ${ \mathcal U }(a,b)$ is a uniform distribution bounded between a and b; ${ \mathcal T }(\mu ,\sigma ,a,b)$ is a truncated normal distribution with mean μ and standard deviation σ and bounded between a and b.

^a The inclination of the satellite's orbit, i_s, is defined from the reference plane of the planet's orbit.

Download table as:  ASCIITypeset image

In this setup, each grid cell corresponds to a child node of the quadtree which occupies (1/4) × (1/4) of the total initial parameter space $\{\cos {i}_{{\rm{s}}}\}\times \{{a}_{{\rm{s}}}\}$. The grid is then iteratively refined to higher resolution by splitting each node into four equal child nodes where the following criteria are satisfied:

• 1.  
  the fractional stability of the node differs from that of a neighboring node by at least 20%, ²⁸ and
• 2.  
  each of the resulting child nodes occupies at least (1/32) × (1/32) of the total parameter space.

After a node has split, all of the completed simulations stored within that node are distributed to its four child nodes according to their $\cos {i}_{{\rm{s}}}$ and a_s values, and each child node is subsequently filled by running more N-body simulations until each node contains a maximum of 25 simulations. This process is repeated until the splitting criteria are no longer true for any nodes. Therefore, the total number of N-body simulations run for the completed stability map is between 4 × 4 × 25 = 400 and 32 × 32 × 25 = 25, 600. Figure 1 shows the final stability map for this setup, which we will discuss further in Section 3.

Zoom In Zoom Out Reset image size

2.1.2. Four-body Simulation

We next run an additional set of N-body simulations, now including the second-outermost planet in the system, HIP 41378 e (which is separated from planet f by ∼10 mutual Hill radii). It is beyond the scope of this work to simulate the global stability of the system, so we do not include planet b (P = 15.57 days), planet c (P = 31.71 days), or planet d (P = 278.36 days) in our simulations in order to reduce the complexity of the problem and lower computational costs. For the purpose of studying the stability of a satellite around planet f, this is a valid assumption given that the force of gravity on planet f due to the three inner planets at closest approach is orders of magnitude smaller than that of planet e (planets b, c, and d are separated from planet f by approximately 65, 61, and 19 mutual Hill radii, respectively). Moreover, the transit timing variations of planet f are thought to be dominated by a 2:3 period commensurability with planet e (Bryant et al. 2021).

For this case, our setup is similar to that of our suite of three-body simulations, but we now evaluate the fractional stability as a function of planet e's orbital period, P_e, and eccentricity, e_e. Furthermore, to account for observational uncertainties in the masses and orbits of planets e and f, we randomly draw other system parameters that are consistent with the results of Santerne et al. (2019). For each trial simulation we randomly draw the mass, argument of pericenter, and inclination of planets e and f, and the period and eccentricity of planet f from normal distributions with standard deviations consistent with the uncertainties given in Table 1.

For the satellite of planet f, we assume a circular orbit and randomly draw the cosine of the inclination (with respect to the planet's orbital plane) from a uniform distribution between 0 and 1, and the initial semimajor axis between 0.1 and 0.4R_H (i.e., within the stability limits for a low-eccentricity host planet derived from our benchmark simulations). We assume the same fixed stellar mass and satellite mass as before, and initialize each object with a random mean anomaly between 0° and 360°. A summary of the distributions used to select parameters for each iteration is provided in Table 2. Figure 2 shows an example of 20 random realizations of the system setup for fixed values of P_e and e_e.

Zoom In Zoom Out Reset image size

With the addition of planet e to our simulations, we update the stability criteria defined in Section 2.1.1, such that a system is also flagged as "unstable" if planet e enters the Hill sphere of planet f, or if any planet's orbit crosses the orbit of another planet. The latter criterion places an upper limit on the eccentricity of planet e for a given orbital period due to the intersection of its periastron distance with the expected orbit of planet d (note, however, that we do not explicitly include planet d in the N-body simulations). Assuming that planet d has a circular orbit, ²⁹ this limit is set by

$$\begin{eqnarray}&&{P}_{{\rm{e}},\mathrm{crit}}({e}_{{\rm{e}}})\simeq \displaystyle \frac{2\pi }{\sqrt{{{GM}}_{\star }}}{\left(\displaystyle \frac{{a}_{{\rm{d}}}}{1-{e}_{{\rm{e}}}}\right)}^{3/2},\end{eqnarray} \tag{ 2 }$$

where a_d is the semimajor axis of planet d and G is the Newtonian gravitational constant. Note that for the measured eccentricity of planet e (e_e = 0.14; Santerne et al. 2019), P_e,crit ≈ 351 days, which is only within ∼1.8σ of the measured period of planet e (P_e = 369 ± 10 days; Santerne et al. 2019).

We explore the fractional stability of the system in the parameter space around this limit over a timescale of 10⁵ yr, selecting a range of P_e and e_e consistent with their observed values within ∼2σ (see Table 2; Santerne et al. 2019). Again, we start by evenly dividing the entire (P_e, e_e) parameter space into a 4 × 4 grid, and computing f_stab in each cell by running a set of 25 independent trials with the set of parameters described above. We then proceed by iteratively subdividing the quadtree according to the criteria defined in Section 2.1.1 and computing more N-body simulations until no more cells can be split. The resulting stability map is shown in Figure 3.

Zoom In Zoom Out Reset image size

2.2. Tidal Migration

In addition to gravitational interactions with other planets, tidal interactions can lead to the secular evolution of a moon's orbit. Here we consider the possible effects of tides on the orbital stability of moons around HIP 41378 f. Here, we implement an equilibrium tide (ET) model, which assumes a simple parameterization to capture the relevant tidal processes without considering the detailed physics of the interiors of the bodies or 3D effects (see Barnes 2017, and references therein). In summary, the ET model assumes that the gravitational force of a satellite creates an elongated bulge in the planet in a state of equilibrium, whose long axis is misaligned with the line connecting the planet and moon centers of mass due to dissipative processes within the planet. This misalignment, or "lag" between the passage of the tidal bulge and the satellite, creates torques that lead to the secular evolution of the spin and orbital angular momenta of the two bodies.

As detailed in Barnes (2017), two distinct parameterizations of the tidal lag are well established in the literature: the so-called "constant time lag" (CTL) model, in which the time interval between the passage of the moon and the tidal bulge is fixed, and the "constant phase lag" (CPL) model in which the phase between the moon and the tidal bulge is constant, independent of frequency. In both frameworks, the energy dissipation and tidal bulge are coupled via two parameters: the tidal Love number of degree 2, k₂, which accounts for the internal forces of the deformed body (where k₂ = 0 describes a perfectly rigid body and k₂ = 3/2 a perfectly fluid body), and a parameter representing the offset between the tidal bulge axis and the line connecting the centers of mass. In the CTL model, the latter parameter is referred to as the "tidal time lag," τ, which is inversely related to the tidal quality factor Q in the CPT model. These parameters effectively describe the efficiency of the tidal response between the two bodies—larger values of τ (smaller Q) imply faster tidal evolution. While both CTL and CPL generally produce qualitatively similar results for the solar system, one may be more suitable than the other depending on the specific scenario. In this work, we choose to implement the CTL model, as it allows us to explore the impact of planetary rotation rate on the tidal evolution of the moon.

In both ET model frameworks, the tidal evolution of two bodies is described by a set of six differential equations and six free parameters: the semimajor axis a, the eccentricity e, the rotation rates of the primary and secondary bodies Ω_i , and their two obliquities ψ_i (where i = 1,2 correspond to the primary and secondary bodies). Following Barnes (2017) and references therein (Mignard 1979; Hut 1981; Greenberg 2009; Leconte et al. 2010; Heller et al. 2011), the evolution for the CTL model is described by the following set of differential equations:

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{de}}{{dt}} & = & \displaystyle \frac{11{ae}}{2{{GM}}_{1}{M}_{2}}\\ & & \times \,\displaystyle \sum _{i=1}^{2}{Z}_{i}^{\mathrm{ctl}}\left(\cos ({\psi }_{i})\displaystyle \frac{{f}_{4}(e)}{{\beta }^{10}(e)}\displaystyle \frac{{{\rm{\Omega }}}_{i}}{n}-\displaystyle \frac{18}{11}\displaystyle \frac{{f}_{3}(e)}{{\beta }^{13}(e)}\right),\end{array}\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{da}}{{dt}} & = & \displaystyle \frac{2{a}^{2}}{{{GM}}_{1}{M}_{2}}\\ & & \times \,\displaystyle \sum _{i=1}^{2}{Z}_{i}^{\mathrm{ctl}}\left(\cos ({\psi }_{i})\displaystyle \frac{{f}_{2}(e)}{{\beta }^{12}(e)}\displaystyle \frac{{{\rm{\Omega }}}_{i}}{n}-\displaystyle \frac{{f}_{1}(e)}{{\beta }^{15}(e)}\right),\end{array}\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{d{{\rm{\Omega }}}_{i}}{{dt}} & = & \displaystyle \frac{{Z}_{i}^{\mathrm{ctl}}}{2{M}_{i}{r}_{g,i}^{2}{R}_{i}^{2}n}\\ & & \times \,\left(2\cos ({\psi }_{i})\displaystyle \frac{{f}_{2}(e)}{{\beta }^{12}(e)}-\left[1+{\cos }^{2}({\psi }_{i})\right]\displaystyle \frac{{f}_{5}(e)}{{\beta }^{9}(e)}\displaystyle \frac{{{\rm{\Omega }}}_{i}}{n}\right),\end{array}\end{eqnarray} \tag{ 5 }$$

and

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{d{\psi }_{i}}{{dt}} & = & \displaystyle \frac{{Z}_{i}^{\mathrm{ctl}}\sin ({\psi }_{i})}{2{M}_{i}{r}_{g,i}^{2}{R}_{i}^{2}n{{\rm{\Omega }}}_{i}}\\ & & \times \,\left(\left[\cos ({\psi }_{i})-\displaystyle \frac{{\xi }_{i}}{\beta }\right]\displaystyle \frac{{f}_{5}(e)}{{\beta }^{9}(e)}\displaystyle \frac{{{\rm{\Omega }}}_{i}}{n}-2\displaystyle \frac{{f}_{2}(e)}{{\beta }^{12}(e)}\right),\end{array}\end{eqnarray} \tag{ 6 }$$

where t is time, M_i are the masses of the two bodies, R_i are their radii, and n is the mean motion. The quantity r_g is the radius of gyration, which is related to the moment of inertia by $I=M{({r}_{g}R)}^{2}$. The above equations are averaged over the orbital period and therefore represent the mean variation of the orbital parameters.

The factors ${Z}_{i}^{\mathrm{ctl}}$ and ξ_i are given by

$$\begin{eqnarray}&&{Z}_{i}^{\mathrm{ctl}}\equiv 3{G}^{2}{k}_{2,i}{M}_{j}^{2}({M}_{i}+{M}_{j})\displaystyle \frac{{R}_{i}^{5}}{{a}^{9}}{\tau }_{i},\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{\xi }_{i}\equiv \displaystyle \frac{{r}_{g,i}^{2}{R}_{i}^{2}{{\rm{\Omega }}}_{i}{an}}{{{GM}}_{j}},\end{eqnarray} \tag{ 8 }$$

where k_2,i are the second-order tidal Love numbers, τ_i are the tidal time lags, and the subscripts i and j refer to the two bodies. The remaining quantities are given by:

$$\begin{eqnarray*}\begin{array}{rcl}\beta (e) & = & \sqrt{1-{e}^{2}},\\ {f}_{1}(e) & = & 1+\displaystyle \frac{31}{2}{e}^{2}+\displaystyle \frac{255}{8}{e}^{4}+\displaystyle \frac{185}{16}{e}^{6}+\displaystyle \frac{25}{64}{e}^{8},\\ {f}_{2}(e) & = & 1+\displaystyle \frac{15}{2}{e}^{2}+\displaystyle \frac{45}{8}{e}^{4}+\displaystyle \frac{5}{16}{e}^{6},\\ {f}_{3}(e) & = & 1+\displaystyle \frac{15}{4}{e}^{2}+\displaystyle \frac{15}{8}{e}^{4}+\displaystyle \frac{5}{64}{e}^{6},\\ {f}_{4}(e) & = & 1+\displaystyle \frac{3}{2}{e}^{2}+\displaystyle \frac{1}{8}{e}^{4},\\ {f}_{5}(e) & = & 1+3{e}^{2}+\displaystyle \frac{3}{8}{e}^{4}.\end{array}\end{eqnarray*}$$

Unlike hot Jupiters (P ≲ 10 days), which are expected to be in tidally synchronous rotation states with their host stars, we expect the orbital separation of HIP 41378 f to be large enough such that tidal effects from its host star are negligible. We test this assumption using a CTL equilibrium tide model as implemented in the open-source software package EqTide (Barnes 2017).

We assign the host star a mass, radius, and rotation period consistent with Santerne et al. (2019), and a solar radius of gyration, r_g,⋆ ≈ 0.26 (e.g., Claret & Gimenez 1989). We assume a Love number of k_2,⋆ = 0.5 and a constant time lag of τ_⋆ = 0.01 s, following Barnes (2017). For HIP 41378 f, we use the planet mass and radius reported in Santerne et al. (2019) and assume a rotation period of 10 hr, ³⁰ similar to that of the solar system gas giants. For the radius of gyration and Love number (which are generally not measurable for exoplanets except in special cases, e.g., Akinsanmi et al. 2019; Hellard et al. 2019; Barros et al. 2022), we assume analogous values to Jupiter's, and use recent measurements of Jupiter from the Juno spacecraft, r_g,f = 0.52 (Ni 2018) and k_2,f = 0.565 (Idini & Stevenson 2021). Note that Saturn has similar properties, but a slightly smaller Love number of k₂ = 0.390 (Lainey et al. 2017), which would not significantly affect our results. For the constant time lag, we also assume an analogous value to solar system gas giants, τ_f = 0.00766 s (e.g., Tokadjian & Piro 2020). For reference, note that τ ≈ 0.766 s for Neptune-like planets and τ ≈ 638 s for rocky planets (Tokadjian & Piro 2020).

Under these assumptions, we evolve the CTL model forward for a duration of 3.1 × 10⁹ yr (the approximate age of the system; Santerne et al. 2019), starting the planet with its currently measured semimajor axis and orbital eccentricity. We find that the relative change in the semimajor axis for planet f is ≪1 ppm, and the change in rotation period is less than 0.1%. Therefore, tidal interactions between planet f and its host star are negligible and we can safely ignore the tidal influence of the star for the purposes of this study.

Next, ignoring the influence of tides from the host star and any other planets, we investigate the tidal stability of a hypothetical moon orbiting planet f. As in Section 2.1, we consider the extreme case of an Earth-composition satellite with the same mass ratio to planet f as the Moon–Earth system (M_s = 0.15 M_⊕, R_s = 0.53 R_⊕). Again, this choice is motivated by the largest moon mass ratio in the solar system and the fact that smaller moons would be much harder to measure. For the tidal properties, we arbitrarily assume that the satellite has a gyration radius similar to the measured value for the solar system's largest moon r_g,s = 0.56 (Ganymede; Showman & Malhotra 1999), and a tidal Love number and constant time lag consistent with measurement of rocky planets in the solar system, k_2,s = 0.3 and τ_s = 638 s (e.g., Tokadjian & Piro 2020, and references therein). The mass and radius of the planet are again assumed to be the values reported in Santerne et al. (2019), and the gyration radius and Love number are set to the respective Jovian values (Ni 2018; Idini & Stevenson 2021).

Because the obliquity, constant time lag, and rotation period of the HIP 41378 f are almost entirely unconstrained, especially given the low bulk density inferred from mass and radius measurements, we explore how changing these parameters affects the final semimajor axis of the moon after a duration of 3.1 Gyr. Again using the EqTide code (Barnes 2017), we evolve the system forward in time starting the satellite on a tidally synchronous, circular orbit. We calculate a suite of models, varying the initial semimajor axis of the moon from 4R_p (≈0.0525R_H) to 0.4R_H (i.e., the approximate stability limit for a circular orbit of planet f from our N-body simulations). As an initial case, we choose τ_f to be consistent with the estimates for the solar system gas giants, τ_J = 0.00766 s (Tokadjian & Piro 2020) and a rotation period similar to that of Jupiter, P_rot,f = 10 hr. Then, we investigate the effect of planetary rotation rate by running the same grid of models assuming a slow planet rotation (P_rot,f = 24 hr) and a fast planet rotation (P_rot,f = 3 hr). For each case, we also test two extreme cases of the constant time lag by setting τ_f to 0.01 and 100 times the value assumed for Jupiter τ_J. We then repeated each of these simulations assuming three different planet obliquities, ψ_p: 0°, 90°, and 180°. Lastly, after running each model for 3.1 Gyr, we compare the final semimajor axis of the moon to the stability limit from dynamic simulations to determine whether tidal migration can significantly destabilize the system. Our results are shown in Figure 4.

Zoom In Zoom Out Reset image size

In this section, we discuss the results of our numerical N-body and tidal migration simulations. We then discuss the observability of large exomoons orbiting HIP 41378 f in transit photometry by comparing theoretical finely sampled light curves to observations from HST and K2 and simulating idealized observations from the PLAnetary Transits and Oscillation of stars (PLATO) mission. Finally, we explore the potential impact of large moons with atmospheres on the observed transmission spectrum of HIP 41378 f.

3.1. Stability Limits from Dynamical Simulations

3.2. Limits from Tidal Evolution

Tidal forces between two rotating bodies can lead to secular changes in their rotation and orbits. For example, a well-studied case in our solar system is the gradual recession of Earth's Moon (and subsequent lengthening of the Earth's day) over Gyr timescales. We investigated how such tidal migration may influence the long-term orbital stability of a hypothetical moon orbiting HIP 41378 f.

In our ET simulations outlined in Section 2, the secular evolution was parameterized by two intrinsic properties of the orbiting bodies, k₂ and τ. Though both k₂ and τ are completely unknown for HIP 41378 f (as are the moment of inertia and spin rate), there are constraints on these parameters for many solar system bodies (though τ and Q have uncertainties spanning orders of magnitude). We therefore assumed that the properties of HIP 41378 f and any hypothetical moon are analogous to solar system bodies. As the fiducial case, we assumed Jupiter-like values for k₂, τ, r_g , and Ω for the planet, as described earlier in Section 2. We subsequently tested cases where τ was 2 orders of magnitude larger and smaller than the estimated Jovian value (with the larger value representing a planet with Neptune-like properties; e.g., Tokadjian & Piro 2020), and the rotation period was changed to either 3 hr ("fast") or 24 hr ("slow"). Finally, we tested the effect of planet obliquity by repeating our simulations for ψ_p = 0°, ψ_p = 90°, and ψ_p = 180°. For the satellite, we assumed tidally synchronous rotation and k₂, τ, and r_g consistent with rocky solar system bodies.

For simplicity, we assumed that each simulated moon started out with a circular orbit. We ran each simulation for 3.1 Gyr, varying the initial semimajor axis of the moon from 4R_p (≈0.0525R_H) to 30.5R_p (≈0.4R_H) such that the moon's initial orbital period was always longer than the planet's rotation period—this ensured that no inward migration would occur (for prograde satellites), in which case the moon could ultimately collide with the planet (e.g., for the 24 hr rotation case, the moon would migrate inward instead of outward if its semimajor axis were ≲1.6R_p ). For each rotation rate, constant time lag, and obliquity case, we show an example of the time evolution of the satellite's semimajor axis in Figure 4, as well as how the final semimajor axis of the satellite varies as a function of its initial semimajor axis.

Our results demonstrate that tidal migration does not significantly alter the stability of a moon orbiting HIP 41378 f over the system's lifetime, except if the moon is in a retrograde orbit relative to the planet's spin (ψ_p = 180°) or if τ or Ω of the planet is very large. In the former case, the moon can be forced to migrate inward sufficiently to collide with the planet, depending on the assumed tidal parameters and the planet's rotation rate (outward migration is not allowed here, so these results are independent of the N-body stability limit). Note that for the ψ_p = 90° case, some inward migration may occur depending on the initial conditions, but not sufficiently to cause a collision—however, as we showed with our N-body simulations, Kozai–Lidov oscillations play a significant role in the dynamical stability of high-inclination orbits. In fact, we expect that Kozai–Lidov oscillations would not permit a stable polar orbit (ψ_p = 90°) scenario to begin with. We note also that, while we attempted to adequately represent the full dynamics of the system here, in nature we would expect both secular tidal migration and Kozai–Lidov oscillations to affect the stability of exomoons in a complex and coupled way. It is beyond the scope of our work here to self-consistently model both tides and orbital dynamics simultaneously.

For the prograde orbit case (ψ_p = 0°), because the secular evolution of the moon's semimajor axis determines whether the moon can remain dynamically stable, it is useful for our discussion to introduce a timescale associated with the evolution of a, which we can write as ${\tau }_{\mathrm{mig}}\sim a/\dot{a}$ (where $\dot{a}$ is the time derivative of a). Under the assumptions that ψ_i = e = 0, M_f ≫ M_s, and Ω_s/n = 1 (i.e., the moon is tidally locked), we can use the time derivative from Equation (4) to show that the migration timescale scales as follows:

$$\begin{eqnarray}&&{\tau }_{\mathrm{mig}}\propto \displaystyle \frac{{a}^{8}}{{k}_{2,{\rm{p}}}{\tau }_{{\rm{p}}}{M}_{{\rm{s}}}{R}_{{\rm{p}}}^{5}}{\left(\displaystyle \frac{{{\rm{\Omega }}}_{{\rm{p}}}}{n}-1\right)}^{-1},\end{eqnarray} \tag{ 9 }$$

where the subscript "p" indicates the planet (primary) and subscript "s" indicates the satellite (secondary). Note that in Equation (4), the assumption of setting e = 0 implies that the β(e), f₁(e), and f₂(e) terms become unity.

This scaling relation is consistent with our numerical results shown in Figure 4 for ψ_p = 0°: larger values of τ cause the semimajor axis to increase on faster timescales than smaller values of τ, and shorter planet rotation periods (larger Ω) also lead to faster evolution of the semimajor axis. Therefore, we note that tides could potentially destabilize a moon with a prograde orbit if the planet has a very large τ value (small Q) or is rotating quickly, which could increase the semimajor axis sufficiently after 3.1 Gyr to be beyond the empirically determined stability limit (see the blue dashed line in the right panel of Figure 4). In fact, if τ is sufficiently large, any moons that formed close to the stability limit could migrate outward far enough to be stripped away from the planet by dynamical interactions. This scenario would require that τ be roughly 2 orders of magnitude larger than the value estimated for the solar system giants. However, we cannot rule out this possibility because τ remains virtually unconstrained, especially for planets like HIP 41378 f, whose structure remains unknown with current measurements, and for which we have no close analogs in the solar system.

Furthermore, we can gain insights regarding different moon scenarios from the scaling relation in Equation (9). For example, this relation shows that less massive moons (which are presumably more common than the large moon we simulated here) migrate on a longer timescale, which means they would be less likely to be stripped away by dynamic interactions. Additionally, if HIP 41378 f has a smaller radius than assumed here, then the migration timescale could drastically increase as well (${\tau }_{\mathrm{mig}}\propto {R}_{{\rm{p}}}^{-5}$). This could be the case, for example, if HIP 41378 f is indeed a smaller-sized planet surrounded by optically thick circumplanetary rings (although modeling such a scenario is beyond the scope of this work). This demonstrates that the cases we have chosen to simulate here probe near the upper limit of plausible moon configurations: decreasing the mass of the moon or shrinking the size of the planet would only serve to increase the chances of the moon's survival.

3.3. Exomoon Observability

Despite significant challenges in robustly detecting moons orbiting exoplanets, a number of efforts in recent years have attempted to uncover the first exomoon signals. Notably, the "Hunt for Exomoons with Kepler" (HEK; Kipping et al. 2012) systematically analyzed the light curves of a subset of all Kepler planet candidates most amenable to hosting a moon. Using Bayesian multimodal nested sampling with a photodynamical forward model (LUNA; Kipping 2011) to generate combined planet-moon transit light curves, Kipping et al. (2013b) constrained the satellite-to-planet mass ratios for seven potential satellite-hosting exoplanets, but did not find compelling evidence for an exomoon around any. Subsequent HEK searches also resulted in null exomoon detections in dozens of Kepler systems (Kipping et al. 2015), but set upper mass limits on satellites orbiting the habitable-zone planet Kepler-22b (Kipping et al. 2013a) in addition to eight planets in M-dwarf systems (Kipping et al. 2014). Finally, HEK has also constrained a tentative upper limit on the occurrence rate of Galilean-size moons around planets orbiting between 0.1 and 1 au (Teachey et al. 2018).

Although these surveys have not yet robustly confirmed any exomoons, they have detected at least two exomoon candidates. Teachey & Kipping (2018) reported evidence from transit observations and timing variations of the first exomoon candidate Kepler-1625b I, consistent with a Neptune-size satellite orbiting a Jovian-like planet. More recently, Kipping et al. (2022) found tentative evidence supporting the existence of the exomoon candidate Kepler-1708b I, consistent with a 2.6 R_⊕ satellite. However, the veracity of either detection is the subject of ongoing discussion in the literature (e.g., Kreidberg et al. 2019; Teachey et al. 2020; Cassese & Kipping 2022).

Based on a suite of N-body and ET simulations in this work, we have shown that it is feasible for HIP 41378 f to host a relatively large exomoon. Naturally, this leads to the question of whether such a moon could be detected with current and future facilities. While direct photodynamical searches for exomoons may be promising for other systems (e.g., Kepler-1625b and Kepler-1708b; Teachey & Kipping 2018; Kipping et al. 2022), several factors severely complicate such an effort for HIP 41378 f given the existing data and the properties of the system.

For example, unlike the original Kepler mission, the extended K2 mission experienced significant systematic noise correlated with the spacecraft pointing (Howell et al. 2014). Though it is possible to correct some correlated noise algorithmically (e.g., Vanderburg et al. 2016b), residual systematics can persist at levels comparable to an exomoon signal, and stellar granulation and p-mode oscillations may persist as important noise sources. Moreover, for HST observations, systematic noise and gaps in the phase/time coverage introduced by the orbit of the telescope around Earth must be treated with caution, as these effects may also complicate detecting an exomoon in the light curve (e.g., Kreidberg et al. 2019). For multiplanet systems such as HIP 41378, one must also be wary of additional uncertainties in the inferred masses and orbits of other planets in the system, which can in turn lead to uncertain TTVs that affect the exomoon host (e.g., planet e). Without a firm grasp of the planet–planet-induced TTVs in this system, it is prohibitively difficult to disentangle potential TTVs caused by an orbiting moon, though in principle it may be possible for other systems (Kipping 2021).

While these factors likely preclude a meaningful direct photodynamical search for exomoons orbiting HIP 41378 f with current observations, it is still worth considering the effects that a large exomoon may have on transit observations. That is, although it is currently very difficult to detect exomoons (and robust confirmation is even more challenging), they are likely included in extant and future observations of their exoplanet hosts. As our observations become more precise with newer technology (e.g., JWST and ELTs) it is necessary to consider exomoons as a source of uncertainty in our inferences about exoplanets, especially as we begin to probe the atmospheres of smaller planets in search of signatures of life. In this section we consider an idealized analog of HIP 41378 f to explore how large exomoons can affect transit observations.

3.3.1. White Light Curve

To demonstrate the high signal-to-noise ratio required to detect a transiting exomoon, we modeled the fully detrended HST white light curve of HIP 41378 f from Alam et al. (2022) and injected the signal from a hypothetical moon orbiting the planet. We computed model light curves and evaluated the likelihood of each model using the Pandora package (Hippke & Heller 2022), an open-source photodynamical exomoon transit code optimized for computational speed and accuracy. First, for the planet-only model, we fixed the barycentric orbital period, semimajor axis, and inclination to the current literature values for HIP 41378 f (Santerne et al. 2019), and fit for the planet-to-star size ratio R_p /R_⋆, quadratic limb-darkening coefficients q₁ and q₂ (Kipping 2013), and an offset ΔT₀ in the barycentric time of transit center T₀ to allow for small deviations from the transit center reported in Alam et al. (2022). We assumed a circular orbit for the planet. To effectively compute a planet-only model using Pandora, we tuned the satellite mass and radius to negligible values (see Hippke & Heller 2022).

Then, to derive posterior probability distributions for the planet-only fit, we used the nested sampling Monte Carlo algorithm MLFriends (Buchner 2016, 2019) implemented in the UltraNest ³¹ package (Buchner 2021). We imposed a Gaussian prior on R_p /R_⋆ based on the results of Alam et al. (2022), and a conservative uniform prior on ΔT₀ of±0.1 days. ³² For the limb-darkening coefficients, we estimated initial values with the ExoTiC-LD package (Wakeford & Grant 2022) using the 3D stellar models from Magic et al. (2015) and stellar parameters from Santerne et al. (2019), then used these values to define Gaussian priors on q₁ and q₂ with a standard deviation of 0.1. Using a Gaussian likelihood function, we then ran the UltraNest ReactiveNestedSampler with 800 live points, and a step sampler with 2000 steps, until convergence was achieved. We chose this particular nested sampling approach because it can be easily extended to conduct efficient retrievals of the full set of planet and moon parameters, whose posterior probability distributions can be multimodal (see Hippke & Heller 2022). Figure 5 shows the detrended white light curve from Alam et al. (2022) along with the best-fit transit model from our nested sampling analysis. The full posterior probability distributions from the analysis are shown in Figure 8 in the appendix.

Zoom In Zoom Out Reset image size

To investigate how a transiting exomoon would change the white light curve, we then calculated Pandora transit models with hypothetical 0.53 R_⊕ and 1 R_⊕ satellites injected on an arbitrary 30 day orbit (a_s ≈ 11 R_p ≈ 0.14 R_H) around HIP 41378 f. We calculated these models with Pandora assuming the best-fit parameters determined from the planet-only fit, but we set the moon radius to 0.53 R_⊕ or 1 R_⊕. In both cases, we arbitrarily selected the moon's orbital phase for ease of comparison and ignored the TTV caused by the moon by leaving the satellite mass set to be negligible. For comparison, these moon-injected transit models are shown along with the best-fit planet-only model in Figure 5. The light curve models in Figure 5 clearly demonstrate the need for much higher precision (≲15 ppm) than is available in current HST observations in order to photometrically detect exomoons around HIP 41378 f.

We note that in an independent analysis of the HST/WFC3 transit of HIP 41378 f, Edwards et al. (2022) obtained an rms in the HST transit light curve of 125 ppm, compared to the 485 ppm rms of Alam et al. (2022); as this precision is still ≫15 ppm, our conclusions here are independent of the choice of HST reduction. This high level of precision could, however, be within reach of JWST for certain systems—for example, Coulombe et al. (2023) showed that the JWST NIRISS/SOSS white light curve of WASP-18b bins down to ∼5 ppm over 1 h timescales, and Lustig-Yaeger et al. (2023) demonstrated a similar result for LHS 475b using JWST NIRSpec/BOTS.

For completeness, we repeated the above procedure using observations from Campaigns 5 and 18 of the K2 mission (Vanderburg et al. 2016a; Santerne et al. 2019), which were detrended with K2SFF (Vanderburg et al. 2016b). For each K2 light curve, we removed low-frequency variability by fitting a basis spline to the out-of-transit flux and dividing it out of the data (as in Berardo et al. 2019). We then repeated the procedure used to fit a planet-only model to the HST data, but we used the results from Santerne et al. (2019) to define T₀ and set Gaussian priors on the R_p /R_⋆, q₁, and q₂. The K2 data and best-fit planet-only model are shown in Figure 9 in the appendix, along with the light curves with injected hypothetical satellites with radii of 0.53 R_⊕ and 1 R_⊕. The posterior probability distributions for the planet-only fit are shown in Figure 10 in the appendix. Similar to the HST observations, the K2 data are not precise enough to distinguish between planet-only and planet-moon transit models, and residual systematic variations at the ≳20 ppm level can mimic the small signal produced by a moon.

Finally, we simulated an idealized single transit observation of HIP 41378 f from the future PLATO mission, expected to launch in 2026 (Rauer et al. 2014). We considered two scenarios: one in which only planet f transits a bright, photometrically quiet star with properties of HIP 41378, and one in which the planet is orbited by an Earth-size satellite with the same configuration as previously described. We assumed a cadence of 600 s (10 minutes), as will be used for PLATO light curves, and a noise level of 20 ppm/hr for a m_V ∼ 9 star in PLATO's P1 sample (though this number will actually depend on the number of cameras used to observe the target; Rauer et al. 2014). Ignoring astrophysical variability, we added white noise to our simulated Pandora light curves; the noise was randomly drawn from a Gaussian distribution with a standard deviation of 49 ppm (i.e., 20 ppm/hr scaled to a 10 minute cadence assuming Poisson noise). The transit models and simulated PLATO observations are shown in the bottom panel of Figure 9 in the appendix. Even in these idealized PLATO simulations, the photometric noise floor makes it extremely difficult to distinguish between the moon and moon-free scenarios (even in the absence of systematic uncertainties). Therefore, any variations on the order of 10 ppm due to instrumental and astrophysical systematics must be very well understood in order to robustly attempt to constrain the presence of an exomoon from transit observations alone.

3.3.2. Implications for Transmission Spectroscopy

As observational capabilities continue to improve in the era of JWST and beyond, it is worth considering the extent to which exomoons may affect the observed wavelength dependence of transit depth, or transmission spectrum, of temperate and potentially habitable worlds. Current HST WFC3 observations of HIP 41378 f's atmosphere are consistent with a featureless transmission spectrum (${\chi }_{r}^{2}\approx 1.05;$ Alam et al. 2022). ³³ Given the low density of the planet, its flat spectrum could be explained by a variety of possible scenarios, including an extended hazy atmosphere, a high-metallicity clear atmosphere, or an optically thick ring system (Akinsanmi et al. 2020; Alam et al. 2022). As these scenarios are all consistent with the observations, for simplicity we chose here to model HIP 41378 f as if it had a high-metallicity clear atmosphere, and then assessed how an exomoon with a substantial atmosphere would alter its transmission spectrum. ³⁴

We computed model transmission spectra and evaluated the likelihood of each one using the PLATON package (Zhang et al. 2019, 2020). First, we calculated the model for the planet based on the results of Alam et al. (2022), assuming a 300× solar metallicity atmosphere and an isothermal (T = 294 K) T–P profile. Following Alam et al. (2022), we fit the model to the observed spectrum by binning the model predictions to the wavelength channels of the observations and performing a least-squares fit. Keeping all other parameters fixed, we fit for a constant vertical offset in R_p/R_⋆ to preserve the shape of the spectrum. The observed spectrum and best-fit 300x solar metallicity model are shown in Figure 6.

Zoom In Zoom Out Reset image size

We then computed a model for the combined transmission spectrum of a planet and a moon, assuming bodies both have substantial atmospheres and are in an optimal transit geometry (i.e., we see the maximum contribution of a potential moon's atmosphere). Following our stability simulations, we considered the scenario of a 0.53 R_⊕ (0.15 M_⊕) moon. For the hypothetical moon, we assumed the same isothermal T–P profile (T = 294 K) but imposed a metallicity that would result in a maximally puffy atmosphere in order to explore an upper limit on the exomoon's signal in the transmission spectrum.

We determined the lowest feasible atmospheric metallicity by considering the timescale associated with atmospheric boil-off due to an isothermal Parker wind (see, for example, Owen & Wu 2016; Wang & Dai 2019). We set the boil-off timescale, τ_boil, equal to the system age (∼3.1 Gyr; Santerne et al. 2019), and solved for the mean molecular weight, μ, using the following set of equations:

$$\begin{eqnarray}&&{\tau }_{\mathrm{boil}}\sim \displaystyle \frac{{M}_{\mathrm{atmo}}}{{\dot{M}}_{\mathrm{Parker}}},\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{Parker}}=4\pi {r}_{{\rm{s}}}^{2}{c}_{{\rm{s}}}{\rho }_{0}\exp \left(\displaystyle \frac{3}{2}-\displaystyle \frac{2{r}_{{\rm{s}}}}{{R}_{\mathrm{sat}}}\right).\end{eqnarray} \tag{ 11 }$$

Here, c_s = (kT/μ)^1/2 is the isothermal sound speed, where k is the Boltzmann constant and T is the isothermal temperature; ${r}_{{\rm{s}}}={{GM}}_{\mathrm{sat}}/(2{c}_{{\rm{s}}}^{2})$ is the sonic radius; and ρ₀ is the atmospheric mass density at the surface.

Assuming a surface pressure of 1 bar for the 0.53 R_⊕ moon's atmosphere (comparable to the atmospheric surface pressure on Titan), we found a minimum mean molecular weight of μ ∼ 4.1, approximately 2× that of Jupiter and 1.5× that of Neptune. This corresponds to a metallicity of approximately 120× solar. We note, however, that this estimation is optimistic, as we do not consider other mass loss processes, such as photoevaporation. Importantly, we also note that we do not attempt to self-consistently model the moon's atmospheric or interior structure, nor how the moon acquired an atmosphere with this composition in the first place. The results presented here are simply intended as a thought experiment for a maximally detectable exomoon, rather than a robust constraint on plausible exomoon atmospheres.

Under the assumption that the hypothetical moon would have had a viable formation pathway to achieve such a low-metallicity atmosphere, we computed its spectrum with PLATON and then combined this with the planet-only model as follows:

$$\begin{eqnarray}&&{({R}_{{\rm{p}}}/{R}_{\star })}_{\mathrm{eff}}=\sqrt{{\delta }_{{\rm{p}}}+{\delta }_{{\rm{s}}}},\end{eqnarray} \tag{ 12 }$$

where ${({R}_{{\rm{p}}}/{R}_{\star })}_{\mathrm{eff}}$ is the effective planet-to-star radius ratio due to the planet and moon, and δ_p and δ_s are the wavelength-dependent transit depths due to the planet and moon, respectively. Similar to before, we then performed a least-squares fit to the observed spectrum from Alam et al. (2022), preserving the shape of the total spectrum by fitting only for a constant offset in R_p/R_⋆. Both models are shown alongside the observations in Figure 6, along with the residual R_p/R_⋆ between the planet-only and combined planet-moon models.

Qualitatively, the addition of a moon with a substantial atmosphere results in a superposition of the moon and planet transmission spectra. Given our assumptions, the main differences between the planet-only spectrum and combined planet-moon spectrum are slightly enhanced H₂O absorption features (<50 ppm excess in R_p/R_⋆) and a relative deficit in mid-IR and optical absorption (<100 ppm difference in R_p/R_⋆) for the combined planet and moon models. We quantified the quality of each model fit, along with a flat line fit, using the reduced least-squares statistic. Like Alam et al. (2022), we cannot distinguish between the planet-only 300x solar metallicity model (${\chi }_{r}^{2}\approx 1.55$) or the flat line (${\chi }_{r}^{2}\approx 1.03$), given the data. Here we have also shown that the addition of the 0.53 R_⊕ exomoon with a μ ∼ 4.1 atmosphere is indistinguishable compared to the planet-only model and the flat line (${\chi }_{r}^{2}\approx 1.61$). It is worth noting that the independent analysis of Edwards et al. (2022) resulted in a different transmission spectrum from the HST/WFC3 data—their spectrum shows an upward slope from blue to red which is difficult to reproduce with equilibrium chemistry models and would lead to a qualitatively poorer fit to the forward models in Figure 6. While we do not investigate the nature of the slope seen in Edwards et al. (2022), we conclude that using their spectrum would not change our inferences about the inability to detect an exomoon in the HST data.

Though the low-metallicity atmosphere we considered here is likely an extreme scenario and not necessarily likely for a moon orbiting HIP 41378 f, the excess absorption features in the model transmission spectrum may suggest that other systems containing large exomoons with thick atmospheres could experience nonnegligible contamination if observed at high enough precision. This would be especially important for planets with smaller radii than we considered here. Such observations would need to be capable of resolving ∼10 ppm absorption features, a task that would be challenging even for JWST (e.g., Coulombe et al. 2023; Lustig-Yaeger et al. 2023).

For example, moons with substantial atmospheres could in principle contaminate otherwise high-fidelity spectra of terrestrial planets, introducing uncertainties in possible constraints on biosignature molecules. Moreover, the effects of a moon's atmosphere would vary in time, as the orbital phase of the moon at the time of observations would change with each transit. The temporal variability of the combined spectrum could therefore be misinterpreted as intrinsic temporal variability in the atmosphere of the planet. Finally, the spectra of planets with large moons could be harder to interpret even if the moon itself lacks a substantial atmosphere. Akin to the transit light source effect (Rackham et al. 2018, 2019), a moon could create an apparent time-varying inhomogeneity on the disk of the host star, leading to false spectral features. While the scenarios discussed thus far are only for singular large moons, we note that the presence of multiple smaller moons (like we see in the outer solar system) could also introduce additional timescales and complexity to the problem. Therefore, it may be important to consider moons as a source of uncertainty in transmission spectroscopy as advances in technology allow for more precise measurements of exoplanet atmospheres.

In this work, we investigated the stability of a hypothetical large moon (0.15 M_⊕, 0.53 R_⊕) orbiting the long-period (P ≈ 1.5 yr) 12 M_⊕ planet HIP 41378 f. Using a suite of numerical N-body simulations, we have shown that moons up to this size can survive in the system over long timescales over a wide range of system setups, including with the added complexity of multiple planets in the HIP 41378 system and tidal interactions between the moon and the host planet. Although our analysis has demonstrated the plausibility of a moon orbiting HIP 41378 f, we have shown that the detection of a relatively large Earth-size moon in the system is unlikely given current observations, but may be feasible in the near future with high-precision JWST observations. Moreover, an exomoon with a sufficiently puffy atmosphere, while highly idealized, could imprint ∼10 ppm features on the transmission spectrum of the planet (though it is unclear whether this signal could be recognized as such). Our main conclusions are summarized as follows:

• 1.  
  HIP 41378 f could feasibly host exomoons at least as massive as 0.15 M_⊕ that are stable over a timescale of 10⁵ yr, in agreement with stability limits determined by previous works (e.g., Rosario-Franco et al. 2020). Under ideal assumptions of zero-eccentricity and co-planar orbits, a 0.15 M_⊕ moon does not escape the gravitational sphere of influence of planet f or collide with the planet as long as it orbits within ∼40% of the planet's Hill radius. Our simulations demonstrate that this limit does not depend strongly on $\cos {i}_{{\rm{s}}}$ for most (low) moon inclinations. However, Kozai–Lidov oscillations can start to become an important destabilizing influence for $\cos {i}_{{\rm{s}}}\gtrsim 0.8$, and disallow stable orbits for the highest inclinations ($\cos {i}_{{\rm{s}}}\gtrsim 0.95$). The satellite mass chosen here represents an extreme scenario, analogous to the highest moon-to-planet mass ratio observed in the solar system.
• 2.  
  Current constraints on the mass and orbit of the inner planet HIP 41378 e do not necessarily preclude the possibility of a moon orbiting planet f. However, improved measurements of both planets' masses and orbits are needed to robustly constrain the stability of a putative exomoon in the system. If, for example, the true eccentricity of planet e is 1σ smaller than the eccentricity reported by Santerne et al. (2019), then moons orbiting within ∼40% of the Hill radius of planet f should be stable, as long as $0\lt \cos {i}_{{\rm{s}}}\lesssim 0.9$. This not only highlights the need for improved masses and orbital constraints for the planets in the HIP 41378 system but also the importance of considering multiple planets in the dynamical modeling of exomoon stability in multiplanet systems in general.
• 3.  
  For most circular prograde orbits, tidal interactions between planet f and a moon are not sufficient to cause the moon to migrate outward beyond the 0.4R_H dynamical stability limit, assuming that planet f has tidal response properties comparable to the solar system giant planets. However, under extreme circumstances where planet f is given a very large time lag (τ ≳ 0.7 s) and fast rotation rate (P ∼ 3 hr), tidal migration can lead to the moon's escape from the gravitational influence of the planet. While in principle more massive moons would also migrate faster, we do not consider such cases here because moons with higher mass ratios have not been confirmed in the solar system or in extrasolar systems. For high planet obliquities (or moon inclinations), tidal interactions alone are not sufficient to destabilize the moon, but, as mentioned in the first point, Kozai–Lidov oscillations can become significant enough to destabilize such systems. We also note that moons with retrograde orbits are more likely to collide with the planet due to inward tidal migration (independent of the dynamical stability limit), but this depends strongly on the assumed tidal parameters of the planet, which remain unconstrained.
• 4.  
  Existing observations of the HIP 41378 system from HST and K2 cannot reliably constrain the presence of exomoons, as the expected transit signal produced by a large moon (0.53 R_⊕) is only on the order of 15 ppm. Furthermore, uncertainties in the orbits and masses of the planets in the HIP 41378 system, as well as correspondingly uncertain TTVs, likely preclude the detection of a moon with planned missions such as PLATO. However, it is possible that future observations with JWST capable of precision better than ∼20 ppm on 10 minute timescales will be able to detect exomoons in other systems with more robust planet measurements. For example, the giant temperate planet PH-2b (Kepler-86b) is currently slated to be observed by JWST NIRSpec/PRISM in 2024 (GO Cycle 2 Proposal ID: 3235; PI: Fortney), and could potentially provide serendipitous observations of a transiting exomoon.
• 5.  
  Exomoons with thick atmospheres may contaminate the transmission spectra of exoplanets at the ∼10 ppm level in a time-dependent manner (well below the precision of HST observations). However, this is an optimistic estimate which assumes that a moon can acquire and retain a low mean molecular weight atmosphere over its lifetime. Nonetheless, exomoons with significant atmospheres may be sources of uncertainty in other exoplanet systems as new technology enables smaller and smaller signals to be probed.

We note that while our analysis provides conservative limits on the stability of exomoons orbiting HIP 41378 f using well-tested methods such as N-body simulations and equilibrium tide modeling, there are a few caveats that may be important to consider in the real HIP 41378 system. For example, multiple planets in the system do not have well-constrained orbital periods or masses due to the difficulty of measuring long-period planets (P ≳ 1 yr). While we attempted to incorporate some of this uncertainty in our analysis (e.g., randomly drawing masses and orbital parameters for planet e in our N-body simulations), we could not account for all the unconstrained degrees of freedom in the actual HIP 41378 system, which includes at least 5 planets in total. Moreover, we did not attempt to model the effects of multiple moons in the system, which also would have increased the complexity of our simulations.

Given the current data, another caveat is that we do not know with certainty what the true nature of HIP 41378 f is. Throughout our analysis, we assumed that the planet has a mass and radius consistent with observations (Santerne et al. 2019) and tidal properties similar to a Jovian planet. However, if in reality the planet is more similar to Neptune and it possesses a circumplanetary ring system (e.g., Akinsanmi et al. 2020; Piro & Vissapragada 2020), this could change the tidal response of the planet and therefore affect the orbital stability of moons.

Moreover, self-consistent modeling of any moon-ring interactions is beyond the scope of this work but would be an interesting subject of future simulations. Akinsanmi et al. (2020) showed that observations of HIP 41378 f are consistent with a planet of radius R_p = 3.7 R_⊕ (ρ_p = 1.4 g cm⁻³), with opaque rings extending from 1.05 to 2.6 R_⊕. If the rings disperse out to the Roche radius, it follows that the density of the ring material would be ρ_r = 1.08 g cm⁻³. In this case, assuming that the ring particles have a typical radius of about 100 cm, the total mass of the rings would be approximately 10¹⁹ kg (see Equation (4) of Piro 2018b), roughly the observed mass of Saturn's rings (Iess et al. 2019). Hence, in reality, gravitational interactions between the rings and the moon may affect their long-term stability (e.g., Nakajima et al. 2020). On a similar note, we did not attempt to model the combined observable effects of moons and rings (e.g., Ohno & Fortney 2022).

Finally, throughout our analysis, we assumed that moons are a natural outcome of planet formation, and we only considered the stability of moons once they reached a tidally locked state with planet f. That is, we did not attempt to treat exomoon formation and evolution in a self-consistent manner, and we remained agnostic to the specific formation pathway of our simulated moons.

To assess the observability of an exomoon, we chose to model transit observations because these are likely the most promising for a long-period transiting planet like HIP 41378 f. We note that aside from photodynamical transit searches for exomoons such as HEK (Kipping et al. 2012) and the future Transiting Exosatellites, Moons, and Planets in Orion (TEMPO) survey (Limbach et al. 2023), several other methods of detecting exomoons have been proposed. For example, the orbital sampling effect (OSE) technique (Heller 2014), which leverages many transit observations to achieve robust statistics, can be used to infer the size, orbital separation, and mass of transiting exomoons from a phase-folded light curve of about a dozen transits (Hippke 2015; Heller et al. 2016). However, this is not currently feasible for HIP 41378 f because the planet's long orbital period and transit duration make it difficult to acquire the necessary large number of observations.

Other novel approaches, including direct imaging and Doppler spectroscopy (Peters & Turner 2013; Agol et al. 2015; Vanderburg et al. 2018), radio emission from satellite-hosting exoplanets (Noyola et al. 2014, 2016; Narang et al. 2023), and combined planet-moon thermal phase curves (Forgan 2017), may also be fruitful probes of exomoons with future technologies. But again, HIP 41378 f is not an optimal target for such studies. The semimajor axis of planet f is too small for current direct imaging technology, with a sky-projected angular separation of ∼13 mas, and uncertainties in the other planet masses and orbits preclude a radial velocity exomoon search. Moreover, uncertainty regarding the nature of planet f makes radio emission studies challenging, and the cool equilibrium temperature of planet f (T_eq = 294 K) is not sufficient to produce an appreciable thermal phase curve.

Nonetheless, as our ability to measure exoplanets continues to improve in the future, exomoons may indeed become a more important source of uncertainty (in addition to transit light source variability, stellar spectrum noise, etc.). Further into the future, we may even be able to robustly confirm detections of exomoons and constrain the properties of their atmospheres with state-of-the-art observational facilities such as ground-based ELTs and next-generation space missions like the Habitable Worlds Observatory (HWO). For example, the Ancillary Science Case A-19 for the Large UV/Optical/Infrared Surveyor (LUVOIR) suggests spectroastrometry (i.e., measuring the astrometric shift that occurs between wavelengths with flux dominated by the exoplanet and wavelengths dominated by the exomoon) as a possible exomoon detection strategy with a 12 m class space telescope (The LUVOIR Team 2019). With contrast >10⁻⁹ and relative astrometric precision between wavelength bands of ∼0.1 mas, such an observatory could in principle detect an Earth-size exomoon orbiting a 1 au Jupiter-size planet at 10 pc (The LUVOIR Team 2019). Not only would this be a revolutionary discovery in itself, but it would also open the possibility of using exomoons as a way of studying planet formation and evolution, and even provide a potential new avenue for searching for signs of life beyond the solar system.

The authors thank Kristo Ment for helpful discussions regarding quadtree algorithms. This paper makes use of observations from the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 526555. These observations are associated with HST-GO program 16267 (PI: Dressing), and the analysis was supported by grant HST-GO-16267. This work was supported by the Programme National de Planétologie (PNP) of CNRS-INSU co-funded by CNES. This work was supported by FCT—Fundação para a Ciência e a Tecnologia through national funds and by FEDER through COMPETE2020—Programa Operacional Competitividade e Internacionalização by these grants: UIDB/04434/2020; UIDP/04434/2020, 2022.06962.PTDC. This research has also been partly funded by the Spanish State Research Agency (AEI) Project No. PID2019-107061GB-C61. A.C.C. acknowledges support from STFC consolidated grant numbers ST/R000824/1 and ST/V000861/1, and UKSA grant No. ST/R003203/1. C.K.H. acknowledges support from the National Science Foundation Graduate Research Fellowship Program under grant No. DGE 2146752. J.K. acknowledges financial support from Imperial College London through an Imperial College Research Fellowship grant. J.L.-B. acknowledges support from the Ramón y Cajal program (RYC2021-031640-I) supported by the MCIN/AEI/10.13039/501100011033 and the European Union "NextGenerationEU"/PRTR as well as partial financial support received from "la Caixa" Foundation (ID 100010434) and the European Unions Horizon 2020 research and innovation program No. 847648, with fellowship code LCF/BQ/PI20/11760023. N.C.S. acknowledges funding from the European Union (ERC, FIERCE, 101052347). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. S.G.S. acknowledges the support from FCT through Investigador FCT contract No. CEECIND/00826/2018 and POPH/FSE (EC). P.J.W. acknowledges support from the UK Science and Technology Facilities Council (STFC) under consolidated grant ST/T000406/1.

Facilities: HST (WFC3) - , K2 - .

Software: Astropy (Astropy Collaboration et al. 2013, 2018), AstroQTpy (Harada 2023), EqTide (Barnes 2017), ExoTiC-LD (Wakeford & Grant 2022), IPython (Perez & Granger 2007), Matplotlib (Hunter 2007), NumPy (van der Walt et al. 2011; Harris et al. 2020), Pandora (Hippke & Heller 2022), PLATON (Zhang et al. 2019, 2020), Rebound (Rein & Liu 2012; Rein & Spiegel 2015), UltraNest (Buchner 2021).

While in the main text we describe our N-body simulations for exploring how the semimajor axis and inclination of the satellite and the presence of planet e affect the system's stability, here we describe the set of initial simulations that were used to test our N-body code.

We benchmarked our code by following the general procedure outlined by Domingos et al. (2006) and Rosario-Franco et al. (2020). We set up our grid of simulations using a quadtree data structure implemented in astrQTpy (Harada 2023), which we described in more detail in Section 2.1.1. We constrained the initial eccentricity, e_f, of planet f between 0.0 and 0.5, and the initial semimajor axis, a_s, of the satellite between 0.25 to 0.55R_H, where R_H is the Hill radius defined in Equation (1) (R_H ≈ 0.03 au ≈76R_p). These limits were chosen in order to explore the parameter space encompassing the stability limit determined by Rosario-Franco et al. (2020). Note however that the measured eccentricity of planet f is actually constrained to ${e}_{{\rm{f}}}={0.004}_{-0.003}^{+0.009}$ (Santerne et al. 2019), as stated in the main text. As in Section 2.1.1, the planet and satellite were initialized with co-planar orbits, with the longitudes of the ascending node and arguments of pericenter set to zero (Ω = ω = 0°), and the mean anomaly of the planet set to zero (${{ \mathcal M }}_{{\rm{f}}}=0^\circ$).

We proceeded by dividing the parameter space evenly into a 4 × 4 grid and running 25 independent N-body simulations in each grid cell with (e_f, a_s) and the satellite mean anomaly, ${{ \mathcal M }}_{{\rm{s}}}$, drawn from uniform distributions. We calculated the fraction of stable systems, f_stab, in each cell as the fraction of simulations (out of 25) which satisfied the stability criteria defined in the main text over a timescale of 10⁵ yr. We then continued splitting each grid cell in the quadtree into four equal child cells and computing more N-body simulations, following the procedure described in Section 2.1.1. Figure 7 shows the final stability grid generated by our simulations, where f_stab is quantified by color. Regions where 0 < f_stab < 1 are considered quasi-stable, while regions where f_stab = 1 (0) are absolutely stable (unstable) over 10⁵ yr.

Zoom In Zoom Out Reset image size

Our results were qualitatively similar to Rosario-Franco et al. (2020) and quantitatively consistent in terms of the stability boundary for low-eccentricity orbits (i.e., the largest value of a_s for a given e_f where f_stab = 1). Rosario-Franco et al. (2020) modeled the stability boundary as

$$\begin{eqnarray}&&{a}_{\mathrm{crit}}={c}_{1}(1-{c}_{2}{e}_{f}),\end{eqnarray} \tag{ A1 }$$

where a_crit is implicitly in terms of the fraction of the Hill radius, c₁ is the stability limit for circular orbits, and c₂ is a slope parameter.

From their N-body simulations, the authors determined that c₁ = 0.4061 ± 0.0028 and c₂ = 1.1257 ± 0.0273. Their stability limit is shown in Figure 7 plotted over our stability. Our results are in good agreement for low eccentricities, with our stability limit for circular orbits being a_crit ≈ 0.41R_H. However, we note that our simulations tend to predict a slightly higher stability limit as eccentricity increases (for e_f ≳ 0.15). This discrepancy at higher eccentricities is likely due to our different choice of planetary semimajor axis; we assumed a semimajor axis of 1.37 au for HIP 41378 f, whereas Rosario-Franco et al. (2020) assumed an orbital separation of 1 au in their simulations. Because our simulated planet and its moon are farther away from the host star, the system tends to be more stable against gravitational perturbations.

Here we present additional figures that supplement the main text. Figure 8 shows the posterior distributions derived from our nested sampling analysis of the HST white light curve shown in Figure 5—maximum likelihood values from the posteriors were used to compute the planet-only model therein. Figure 9 shows additional observations of HIP 41378 f from K2 Campaigns 5 and 18, along with simulated data of a single transit event from the upcoming PLATO mission. The posterior distributions derived from nested sampling analyses of the K2 light curves are shown in Figure 10, which were used to calculate the maximum likelihood of planet-only transit models.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size