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

Moons orbiting exoplanets ( “ exomoons ” ) may hold clues about planet formation, migration, and habitability. In this work, we investigate the plausibility of exomoons orbiting the temperate ( T eq = 294 K ) giant ( R = 9.2 R ⊕ ) planet HIP 41378 f, which has been shown to have a low apparent bulk density of 0.09 g cm − 3 and a ﬂ at near-infrared transmission spectrum, hinting that it may possess circumplanetary rings. Given this planet ’ s long orbital period ( P ≈ 1.5 yr ) , it has been suggested that it may also host a large exomoon. Here, we analyze the orbital stability of a hypothetical exomoon with a satellite-to-planet mass ratio of 0.0123 orbiting HIP 41378 f. Combining a new software package, astroQTpy , with REBOUND and EqTide , we conduct a series of N -body and tidal migration simulations, demonstrating that satellites up to this size are largely stable against dynamical escape and collisions. We simulate the expected transit signal from this hypothetical exomoon and show that current transit observations likely cannot constrain the presence of exomoons orbiting HIP 41378 f, though future observations may be capable of detecting exomoons in other systems. Finally, we model the combined transmission spectrum of HIP 41378 f and a hypothetical moon with a low-metallicity atmosphere and show that the total effective spectrum would be contaminated at the ∼ 10 ppm level. Our work not only demonstrates the feasibility of exomoons orbiting HIP 41378 f but also shows that large exomoons may be a source of uncertainty in future high-precision measurements of exoplanet systems.

1. INTRODUCTION 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).Moonplanet 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 (> 5, 000), 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 aide in the interpretation of future high-precision observations.For example, numerical N-body simulations have set limits on exomoon eccentricities and semi-major 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;Quarles et al. 2020;Jagtap et al. 2021;Sucerquia et al. 2019), 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, 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 groundbased 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 HARPS, HARPS-N, HIRES, and PFS spectrographs, and 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, private communication).A summary of the most recent planet parameters from Santerne et al. (2019) is shown in Table 1 1 .
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 −3 ; 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 (Ohno & Tanaka 2021;Gao & Zhang 2020;Wang & Dai 2019).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 HIP 41378 f's apparently low density 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 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 −3 ) 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 (χ 2 r ≈ 1.03), high-altitude photochemical hazes (χ 2 r ≈ 0.97), or a high-metallicity atmosphere (χ 2 r ≈ 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, HIP 41378 f's large separation from its host star relative to the population of hot Jupiters (P ≲ 10 d) 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 highprecision facilities such as ground-based 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 follow: 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.

METHODS
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 dynamical 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.

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 semi-major 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 satellitelarger 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 −3 ), 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;Quarles et al. 2021;Jagtap et al. 2021), each simulation is run for a duration of up to 10 5 yr.A simulation is halted and marked as "unstable" if either of the following criteria are met: • the satellite collides with the planet2 , or • the satellite escapes the gravitational influence of the planet by crossing outside of the planet's Hill radius, defined as where a p is the planet's semi-major 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 semi-major 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 planet f's eccentricity 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 semi-major axis, a s .We also tested a zeroeccentricity case and found similar results (in general, lower-eccentricity systems are expected be more stable regardless of satellite inclination; see Appendix A).We constrained our simulations to vary i s from cos i s = 0 (aligned with the planet's reference plane) to cos i 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 (M f = 0 • ).
Then, to efficiently explore the moon's stability within our defined parameter space, we run a suite of Nbody simulations using a quadtree data structure im-plemented in astroQTpy3 (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 s and a s .Within each grid cell, we run 25 independent N-body simulations4 with the initial conditions cos i s and a s drawn from uniform distributions.The satellite's initial mean anomaly, M s , is also chosen randomly from a uniform distribution for each simulation (see Table 2).Each simulation is run forward for up to 10 5 yr, and then the fractional sta- Note-N (µ, σ) is a normal distribution with mean µ and standard deviation σ; U (a, b) is a uniform distribution bounded between a and b; T (µ, σ, 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, is, is defined from the reference plane of the planet's orbit.
bility of each cell, f stab , is calculated as the fraction of simulations (out of 25) which satisfy our stability criteria defined above.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 s } × {a 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: • the fractional stability of the node differs from that of a neighboring node by at least 20% 5 , and • 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 s and a s values, and each child node is subsequently filled by running more N-body 5 We replicated our simulations using a more stringent threshold of 10%, but found that this had no major impact on our results.
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.

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.57d), planet c (P = 31.71d), or planet d (P = 278.36d) 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 semi-major 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 .
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 orbit6 , this limit is set by where a d is the semi-major 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 d, which is only within ∼1.8σ of the measured period of planet e (P e = 369 ± 10 d; Santerne et al. 2019).
We explore the fractional stability of the system in the parameter space around this limit over a timescale of 10 5 yrs, 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.

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 which 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 litera- ture: 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 2 , which accounts for the internal forces of the deformed body (where k 2 = 0 describes a perfectly rigid body and k 2 = 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 val-ues 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 semi-major 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: and 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 ctl i and ξ i are given by 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: 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 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 hours7 , 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., Hellard et al. 2019;Akinsanmi 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 2 = 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 9 yr (the approximate age of the system; Santerne et al. 2019), starting the planet with its currently measured semi-major axis and orbital eccentricity.We find that the relative change in semi-major 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 Earthcomposition 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 semi-major 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 tidallysynchronous, circular orbit.We calculate a suite of models, varying the initial semi-major 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 Nbody 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 semi-major axis of the moon to the stability limit from dynamical simulations to determine whether tidal migration can significantly destabilize the system.Our results are shown in Figure 4.

RESULTS AND DISCUSSION
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 then simulating idealized observations from the PLATO mission.Finally, we explore the potential impact of large moons with atmospheres on the observed transmission spectrum of HIP 41378 f.

Stability Limits from Dynamical Simulations
Our three-body simulations expanded upon the work of Rosario-Franco et al. (2020), but in greater detail for the HIP 41378 system.While Rosario-Franco et al. (2020) studied the stability of a Neptune-size moon orbiting a Kepler-1625b analog as a function of planet eccentricity and moon semi-major axis, here we simulated HIP 41378 f with a 0.15 M ⊕ Mars-size satellite, varying the semi-major axis, a s , inclination, i s and initial mean anomaly, M s , of the satellite.Ignoring the other planets in the system and assuming a fixed eccentricity of 0.035 (95% confidence upper limit; Santerne et al. 2019) for planet f, we defined the quantity f stab as the fraction of 25 simulations within each quadtree grid cell that are stable over 10 5 yr. Figure 1 shows the results of our simulations, where f stab is quantified by color.
From the benchmark three-body simulations described in Appendix A, the outer stability limit for a satellite on a co-planar orbit around a planet with an eccentricity of 0.035 is a crit ≈ 0.39 R H .Our results shown in Figure 1 show that this limit depends only weakly on the inclination of the satellite for cos i s ≲ 0.8.At higher inclinations however, Kozai-Lidov oscillations can play a significant role in the dynamical evolution of the system, forcing the destabilization of the satellite over time.For these high-inclination systems, a trade off can occur between the satellite's inclination and eccentricity in order to conserve angular momentum, such that initially circular orbits become highly eccentric over time (i.e., bringing the satellite's pericenter distance closer to the planet, and the apocenter distance farther away).As shown in Figure 1, this process starts to become important for cos i s ≳ 0.8, and dominates for cos i s ≳ 0.95 where there are no stable orbits.
The likelihood of these different inclination scenarios occurring in nature depends strongly on the moon's formation pathway.For moons formed within a circumplanetary disk, as is thought to be the case for the Galilean moon system around Jupiter, the inclination of the moon system relative to the planet spin axis should be low because of momentum conservation from the rotating circumplanetary disk.However, other formation scenarios, such as a kinetic impacts or dynamical capture, may result in moons with more misaligned inclinations, which we've demonstrated are less stable over long timescales, especially for closer-in orbits.We note, however, that tidal interactions may lead to significant outward migration of such close-in moons, thereby increasing their overall stability.We will discuss the effects of tidal migration in greater detail in Section 3.2.
We then used the results from this three-body case and from our benchmark test in Appendix A to inform .Each model assumes a different constant time lag τ for the planet, indicated by color (where τJ = 0.00766 s is roughly Jupiter's value).The fiducial planet rotation period is assumed to be 10 hours (solid lines), but we also consider a slow rotation case (Prot = 24 hours; dashed lines) and a fast rotation case (Prot = 3 hours; dot-dashed lines).The gray hatched regions highlight regions where the moon tends to become dynamically unstable, as ≳ 0.41 RH) (see Figure 1).Note that the ratio of the final-to-initial moon semi-major axis approaches unity as the starting semi-major axis increases in each obliquity case we considered.This is consistent with the very strong dependence of τmig on the semi-major axis in Equation 9 for the zero-obliquity case (i.e., τmig ∝ a 8 ).Note also that for the retrograde (ψp = 180 • ) case, moons with initial orbits close to the planet tend to migrate inward and ultimately collide with the planet, depending on the assumed planetary rotation and tidal lag constant.
our second suite of simulations, which more realistically represented the HIP 41378 system with the inclusion of planet e.Since the observed orbit of planet f is nearly circular, for each simulation we randomly selected the initial semi-major axis of its satellite from a uniform distribution between 0.1R H and 0.4R H , the approximate stability limit for low eccentricity planet orbits.The mass of the satellite was assumed to be the same as in our initial simulations (0.15 M ⊕ ), and the initial mean anomaly and inclination were drawn from uniform distributions.The other initial parameters for the satellite and planets e and f were chosen randomly from the distributions given in Table 2.As before, we determined f stab by taking the fraction of 25 simulations within a given quadtree cell spanning {e e } and {P e } that survived over 10 5 yr.Our results are shown in Figure 3, where again f stab is indicated by color.We also show the upper limit on e e as a function of P e (given by Equation 2), which is constrained by the (assumed fixed) orbit of planet d; and the observed values of e e and P e from Santerne et al. (2019).
The quasi-stable regions (left of the blue dashed line in Figure 3) are dominated by at least two independent effects.First, a trend of decreasing quasi-stability toward larger e e is set by gravitational interactions between planet e and planet f's satellite.At higher eccentricities, the apastron distance of planet e becomes larger, reducing the separation between its orbit and the orbit of planet f.Therefore, when planet e has a high eccentricity, there is a greater likelihood that it will gravitationally perturb planet f's satellite, leading to an unstable system.By the same physical reasoning, we also see a trend of decreasing quasi-stability along the period axis-for fixed eccentricity, the system becomes less stable as the period of planet e increases.In this case, the apastron distance also increases with period, according to Kepler's third law, which again subse-quently increases the likelihood of planet e gravitationally perturbing planet f's satellite.
Second, the quasi-stable regions in Figure 3 are also influenced by the randomly distributed parameters in each simulation, which are independent of the orbit of planet e.These manifest as random noise in the stability map.For example, in cases where the mass of planet f is ≳ 2σ less than the "best-fit" measured mass, the planet's Hill radius becomes sufficiently small such that moons with larger initial orbital separations can more easily be stripped away.Moreover, moons that start out with orbits highly inclined to the reference plane of the planet (especially those with smaller semi-major axes) are more likely to be unstable due to Kozai-Lidov oscillations.
Because the observed period and eccentricity of planet e reside in a transitional quasi-stable region of the map, we cannot robustly conclude whether a moon orbiting planet f could be stable given the uncertainties in planet e's measured period and eccentricity.Nonetheless, improved measurements of planet e's orbit may yet reveal that the system actually resides in a stable region of Figure 3, for instance, if its eccentricity is determined to be smaller than currently thought.If this turns out to be the case, then we would expect that most moons within ∼ 0.4R H of planet f should remain stable as long as their orbits are not highly inclined relative to the planet's reference plane.This result highlights the importance of modeling gravitational effects from other planets in multi-planet systems in the context of exomoon stability (and by extension, habitability), as well as the need for improved radial velocity measurements to measure planet masses in multi-planet systems.

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 2 and τ .Though both k 2 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 2 , τ , r g , and Ω for the planet, as described earlier in Section 2. We subsequently tested cases where τ was two 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 2 , τ , 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 semi-major 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 semi-major 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 semi-major axis in Figure 4, as well as how the final semi-major axis of the satellite varies as a function of its initial semi-major 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 semi-major 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 τ mig ∼ a/ ȧ (where ȧ 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: 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 1 (e), and f 2 (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 semi-major axis to increase on faster timescales than smaller values of τ , and shorter planet rotation periods (larger Ω) also lead to faster evolution of the semi-major 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 semi-major 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 which 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 two 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 dynamical interactions.Additionally, if HIP 41378 f has a smaller radius than assumed here, then the migration timescale could drastically increase as well (τ mig ∝ R −5 p ).This could be the case, for ex-ample, if HIP 41378 f is indeed a smaller size 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.

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 andKepler-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 challening), 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.

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 opti-mized for computational speed and accuracy.First, for the planet-only model, we fixed the barycentric orbital period, semi-major 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 1 and q 2 (Kipping 2013), and an offset ∆T 0 in the barycentric time of transit center T 0 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(Buchner , 2019) ) implemented in the UltraNest8 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 0 of ±0.1 days9 .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 1 and q 2 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 multi-modal (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.
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  2022) (points with errorbars) and best fit planet-only Pandora model (red dashed line).The black lines show simulated light curves with satellites injected at arbitrary orbital phase for ease of comparison-the solid line shows the signal from a 0.53 R⊕ moon, while the dash-dotted line shows the signal from an Earth-size moon.Both simulated moons are assumed to have orbits aligned with the planet's reference plane and an orbital period of 30 days.A zoomed-in view near the transit center is shown in the inset panel, demonstrating the small effect produced by the presence of exomoons-the maximal flux difference between the moon-free and 0.53 R⊕ moon models (∼14 ppm; i.e., the transit depth of a 0.53 R⊕ body transiting a 1.3 R⊙ star) is highlighted in cyan.Posterior distributions for the planet-only model parameters are shown in Figure 8. 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); since 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 one hour 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 0 and set Gaussian priors on the R p /R ⋆ , q 1 , and q 2 .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 planetonly 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 PLAnetary Transits and Oscillation of stars (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 min), 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 that was randomly drawn from a Gaussian distribution with a standard deviation of 49 ppm (i.e., 20 ppm/hr scaled to 10-min 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.

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 (χ 2 r ≈ 1.05; Alam et al. 2022)10 .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 highmetallicity clear atmosphere, or an optically thick ring system (Alam et al. 2022;Akinsanmi et al. 2020).Since 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 spectrum11 .
We computed model transmission spectra and evaluated the likelihood of each one using the PLATON pack-age (Zhang et al. 2019(Zhang et al. , 2020)).First, we calculated the model for the planet based on the results of Alam et al. (2022), assuming a 300x 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.
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: Here, c s = (kT /µ) 1/2 is the isothermal sound speed, where k is the Boltzmann constant and T is the isothermal temperature; r s = GM sat /(2c 2 s ) is the sonic radius; and ρ 0 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 2x that of Jupiter and 1.5x that of Neptune.This corresponds to a metallicity of approximately 120x 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, then combined this with the planet-only model as follows: where (R p /R ⋆ ) 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 2 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 (χ 2 r ≈ 1.55) or the flat line (χ 2 r ≈ 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 (χ 2 r ≈ 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., Lustig-Yaeger et al. 2023;Coulombe 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, 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(Rackham et al. , 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.

CONCLUSION
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 5 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 planet f's gravitational sphere of influence 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 s for most (low) moon inclinations.However, Kozai-Lidov oscillations can start to become an important destabilizing influence for cos i s ≳ 0.8, and disallow stable orbits for the highest inclinations (cos i s ≳ 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 < cos i s ≲ 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 dynamical modeling of exomoon stability in multi-planet 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 of the order 15ppm.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 which 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 possess 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 modelling 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 −3 ), 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 −3 .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 19 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 a 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. 2022), 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 phasefolded 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(Noyola et al. , 2016;;Narang et al. 2022), 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 semi-major 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 make 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 stateof-the-art observational facilities such as ground-based extremely large telescopes (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 −9 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.
where a crit is implicitly in terms of the fraction of the Hill radius, c 1 is the stability limit for circular orbits, and c 2 is a slope parameter.
From their N-body simulations, the authors determined that c 1 = 0.4061 ± 0.0028 and c 2 = 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 semi-major axis; we assumed a semi-major 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.

B. ADDITIONAL FIGURES
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 planet-only transit models.q 2 = 0.1916 +0.0437 0.0327

Figure 1 .
Figure1.Three-body orbital stability map for a system consisting of the star HIP 41378, planet f, and planet f's hypothetical satellite, as a function of the satellite's semi-major axis as (expressed in terms of the Hill radius, RH) and the cosine of the satellite's inclination cos is (with respect to the planet's reference plane).The color of each grid cell corresponds to the fraction f stab of the 25 simulations contained inside the cell (shown as points) which survive for 10 5 yr.Absolutely unstable regions are represented by the lightest areas, while absolutely stable regions are shown by the darkest red regions.Note that moons with high orbital inclinations (cos is ≳ 0.8) tend to become unstable due to Kozai-Lidov oscillations, and there are no stable orbits for cos is ≳ 0.95.

Figure 2 .
Figure2.Example two-dimensional projections of the initial positions and orbits of each object in our four-body simulations, including planet e. Panel (a) shows the top-down and edge-on projections of 20 random realizations of the system, with the initial orbits (lines) and positions (points) of planet e shown in orange, and the initial orbits and positions of planet f shown in blue.The black "⋆" symbol shows the position of the host star, and the dashed black line shows the (assumed) fixed orbital semi-major axis for planet d.The small red boxes indicate the limits of the area plotted in panel (b), which shows a zoomed-in version of the system highlighting the initial position and orbit of the satellite, marked "s" (red), orbiting planet f.

Figure 3 .
Figure3.Four-body orbital stability map for a system including the star HIP 41378, planets e and f, and a satellite orbiting planet f, as a function of planet e's initial orbital period (Pe) and eccentricity (ee).As in Figure1, the color of each grid cell shows the fraction of 25 simulations which are stable over 10 5 yr.Here, the dashed blue line shows Pe,crit(ee), the theoretical upper limit on planet e's eccentricity vs. period due to the orbit of planet d (see Equation2)-the region to the right of this line is unstable by definition because planet e would cross the orbit of planet d.The measured period and eccentricity of planet e fromSanterne et al. (2019) are shown by the green point with errorbars representing 1-σ uncertainties.For reference, the period of planet e corresponding to a 2:3 mean motion resonance (MMR) with planet f is shown by the horizontal dashed line.

Figure 4 .
Figure4.Equilibrium tide evolution models for a 0.15 M⊕(0.53 R⊕) satellite orbiting HIP 41378 f computed with the CTL mode in EqTide(Barnes 2017) for different planet obliquities.Top row: Moon semi-major axis as a function of time over the lifetime of the system, starting from an initial orbital separation of 4Rp.Bottom row: Final vs. initial moon semi-major axis after a timescale of 3.1 Gyr, the current estimated age of the system(Santerne et al. 2019).The columns from left to right show the results for different assumed planet obliquities (moon inclinations): ψp = 0 • (aligned prograde orbit), ψp = 90 • (polar orbit), and ψp = 180 • (aligned retrograde orbit).Each model assumes a different constant time lag τ for the planet, indicated by color (where τJ = 0.00766 s is roughly Jupiter's value).The fiducial planet rotation period is assumed to be 10 hours (solid lines), but we also consider a slow rotation case (Prot = 24 hours; dashed lines) and a fast rotation case (Prot = 3 hours; dot-dashed lines).The gray hatched regions highlight regions where the moon tends to become dynamically unstable, as ≳ 0.41 RH) (see Figure1).Note that the ratio of the final-to-initial moon semi-major axis approaches unity as the starting semi-major axis increases in each obliquity case we considered.This is consistent with the very strong dependence of τmig on the semi-major axis in Equation 9 for the zero-obliquity case (i.e., τmig ∝ a 8 ).Note also that for the retrograde (ψp = 180 • ) case, moons with initial orbits close to the planet tend to migrate inward and ultimately collide with the planet, depending on the assumed planetary rotation and tidal lag constant.

Figure 5 .
Figure 5. HST white light curve of HIP 41378 f from Alam et al. (2022) (points with errorbars) and best fit planet-only Pandora model (red dashed line).The black lines show simulated light curves with satellites injected at arbitrary orbital phase for ease of comparison-the solid line shows the signal from a 0.53 R⊕ moon, while the dash-dotted line shows the signal from an Earth-size moon.Both simulated moons are assumed to have orbits aligned with the planet's reference plane and an orbital period of 30 days.A zoomed-in view near the transit center is shown in the inset panel, demonstrating the small effect produced by the presence of exomoons-the maximal flux difference between the moon-free and 0.53 R⊕ moon models (∼14 ppm; i.e., the transit depth of a 0.53 R⊕ body transiting a 1.3 R⊙ star) is highlighted in cyan.Posterior distributions for the planet-only model parameters are shown in Figure8.

Figure 6 .
Figure 6.Top: Transmission spectrum of HIP 41378 f (points with error bars;Alam et al. 2022), compared to 1D forward models and a flat line.For the planet-only model (black line), we assume a cloud-free 300x solar metallicity atmosphere and an isothermal (T = 294 K) temperature profile; the combined planet and moon model (blue line) is described in the main text.The fainter lines show the computed unbinned (R ∼ 1000) PLATON spectra, while the colored points show each model binned to the same resolution as the WFC3 observations.Middle: HIP 41378 f spectrum and models shown in the full PLATON wavelength range, including wavelengths that are accessible with JWST.The darker bold lines are binned by a factor of 30 to facilitate easier comparison.Bottom: Residual signal between the observed spectrum and the planet-only spectrum model, and between the planet-only spectrum model and the planet-moon model, over the same wavelength range.

Figure 7 .
Figure 7. Three-body orbital stability map for a system consisting of HIP 41378, planet f, and planet f's satellite, as a function of planetary eccentricity, e f , and satellite semi-major, axis as (expressed in terms of the Hill radius, RH).The color of each grid cell corresponds to the fraction f stab of 25 simulations within each cell (indicated by the black points) that survive over 10 5 yr.The stability limit determined by Rosario-Franco et al. (2020), shown as the dashed line, indicated good agreement with our simulations at low eccentricities.

Figure 8 .Figure 9 .
Figure 8. HST white light curve posterior distributions for the planet-only model shown in Figure 5.

Table 2 .
Parameter distributions used for N-body simulations.