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

Caleb K. Harada, Courtney D. Dressing, Munazza K. Alam, James Kirk, Mercedes López-Morales, Kazumasa Ohno, Babatunde Akinsanmi, Susana C. C. Barros, Lars A. Buchhave, A. Collier Cameron, et al.

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

To cite this version:

Caleb K. Harada, Courtney D. Dressing, Munazza K. Alam, James Kirk, Mercedes López-Morales, et al. Stability and Detectability of Exomoons Orbiting HIP 41378 f, a Temperate Jovian Planet with an Anomalously Low Apparent Density. The Astronomical Journal, 2023, 166, 10.3847/1538-3881/ad011c. insu-04479038

HAL Id: insu-04479038

https://insu.hal.science/insu-04479038

Submitted on 10 Mar 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_1.png" /></figure>

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

Caleb K. Harada<sup>1,22</sup>①, Courtney D. Dressing<sup>1</sup>①, Munazza K. Alam<sup>2</sup>①, James Kirk<sup>3</sup>①, Mercedes López-Morales<sup>4</sup>①, Kazumasa Ohno<sup>5,6</sup>①, Babatunde Akinsanmi<sup>7</sup>①, Susana C. C. Barros<sup>8,9</sup>①, Lars A. Buchhave<sup>10</sup>①, A. Collier Cameron<sup>11</sup>①, Ian J. M. Crossfield<sup>12</sup>, Fei Dai<sup>13,14,23</sup>①, Peter Gao<sup>2</sup>①, Steven Giacalone<sup>1</sup>①, Salomé Grouffal<sup>15</sup>①, Jorge Lillo-Box<sup>16</sup>①, Andrew W. Mayo<sup>1,17</sup>①, Annelies Mortier<sup>18</sup>①, Alexandre Santerne<sup>15</sup>①, Nuno C. Santos<sup>8,9</sup>①, Sérgio G. Sousa<sup>8</sup>①, Emma V. Turtelboom<sup>1</sup>①, Andrew Vanderburg<sup>19</sup>①, and Peter J. Wheatley<sup>20,21</sup>①

## Abstract

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 \text{ g cm}^{-3}$ and a flat near-infrared transmission spectrum, hinting that it may possess circumplanetary rings. Given this planet’s long orbital period ($P \approx 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.

*Unified Astronomy Thesaurus concepts: Exoplanet astronomy (486); Exoplanet dynamics (490); Exoplanet systems (484); Exoplanet tides (497); Natural satellites (Extrasolar) (483); Transits (1711); Transmission spectroscopy (2133)*

## 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). Moon-planet interactions may also be crucial for the long-term stability of exoplanet habitability; for example,

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

<footnote n="1"><sup>1</sup> Department of Astronomy, 501 Campbell Hall #3411, University of California, Berkeley, Berkeley, CA 94720, USA; charada@berkeley.edu</footnote>

<footnote n="2"><sup>2</sup> Earth & Planets Laboratory, Carnegie Institution for Science, 5241 Broad Branch Rd NW, Washington, DC 20015, USA</footnote>

<footnote n="3"><sup>3</sup> Department of Physics, Imperial College London, Prince Consort Road, London, SW7 2AZ, UK</footnote>

<footnote n="4"><sup>4</sup> Center for Astrophysics | Harvard & Smithsonian, 60 Garden St, Cambridge, MA 02138, USA</footnote>

<footnote n="5"><sup>5</sup> Division of Science, National Astronomical Observatory of Japan, Tokyo, Japan</footnote>

<footnote n="6"><sup>6</sup> Department of Astronomy and Astrophysics, University of California, Santa Cruz, 1156 High St, Santa Cruz, CA 95064, USA</footnote>

<footnote n="7"><sup>7</sup> Observatoire Astronomique de l'Université de Genève, Chemin Pegasi 51, CH-1290 Versoix, Switzerland</footnote>

<footnote n="8"><sup>8</sup> Instituto de Astrofísica e Ciências do Espaço, Universidade do Porto, CAUP, Rua das Estrelas, 4150-762 Porto, Portugal</footnote>

<footnote n="9"><sup>9</sup> Departamento de Física e Astronomia, Faculdade de Ciências, Universidade do Porto, Rua do Campo Alegre, 4169-007 Porto, Portugal</footnote>

<footnote n="10"><sup>10</sup> DTU Space, National Space Institute, Technical University of Denmark, Elektrovej 328, DK-2800 Kgs. Lyngby, Denmark</footnote>

<footnote n="11"><sup>11</sup> Centre for Exoplanet Science, SUPA School of Physics and Astronomy, University of St Andrews, North Haugh, St Andrews KY16 9SS, UK</footnote>

<footnote n="12"><sup>12</sup> Department of Physics and Astronomy, University of Kansas, Lawrence, KS, USA</footnote>

<footnote n="13"><sup>13</sup> Division of Geological and Planetary Sciences, 1200 E California Blvd, Pasadena, CA, 91125, USA</footnote>

<footnote n="14"><sup>14</sup> Department of Astronomy, California Institute of Technology, Pasadena, CA 91125, USA</footnote>

<footnote n="15"><sup>15</sup> Aix Marseille University, CNRS, CNES, LAM, Marseille, France</footnote>

<footnote n="16"><sup>16</sup> Centro de Astrobiología (CAB, CSIC-INTA), Depto. de Astrofísica, ESAC campus, 28692, Villanueva de la Cañada (Madrid), Spain</footnote>

<footnote n="17"><sup>17</sup> Centre for Star and Planet Formation, Natural History Museum of Denmark & Niels Bohr Institute, University of Copenhagen, Øster Voldgade 5-7, DK-1350 Copenhagen K., Denmark</footnote>

<footnote n="18"><sup>18</sup> School of Physics & Astronomy, University of Birmingham, Edgbaston, Birmingham B15 2TT, UK</footnote>

<footnote n="19"><sup>19</sup> Department of Physics and Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology, Cambridge, MA 02139, USA</footnote>

<footnote n="20"><sup>20</sup> Centre for Exoplanets and Habitability, University of Warwick, Gibbet Hill Rd, CV4 7AL Coventry, UK</footnote>

<footnote n="21"><sup>21</sup> Department of Physics, University of Warwick, Gibbet Hill Rd, CV4 7AL Coventry, UK</footnote>

<footnote>Received 2023 March 17; revised 2023 September 5; accepted 2023 October 5; published 2023 October 24</footnote>

<footnote n="22"><sup>22</sup> NSF Graduate Research Fellow.</footnote>

<footnote n="23"><sup>23</sup> NASA Sagan Fellow.</footnote>
<figure><img src="image_3.png" /></figure>

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

Caleb K. Harada<sup>1,22</sup>①, Courtney D. Dressing<sup>1</sup>①, Munazza K. Alam<sup>2</sup>①, James Kirk<sup>3</sup>①, Mercedes López-Morales<sup>4</sup>①, Kazumasa Ohno<sup>5,6</sup>①, Babatunde Akinsanmi<sup>7</sup>①, Susana C. C. Barros<sup>8,9</sup>①, Lars A. Buchhave<sup>10</sup>①, A. Collier Cameron<sup>11</sup>①, Ian J. M. Crossfield<sup>12</sup>, Fei Dai<sup>13,14,23</sup>①, Peter Gao<sup>2</sup>①, Steven Giacalone<sup>1</sup>①, Salomé Grouffal<sup>15</sup>①, Jorge Lillo-Box<sup>16</sup>①, Andrew W. Mayo<sup>1,17</sup>①, Annelies Mortier<sup>18</sup>①, Alexandre Santerne<sup>15</sup>①, Nuno C. Santos<sup>8,9</sup>①, Sérgio G. Sousa<sup>8</sup>①, Emma V. Turtelboom<sup>1</sup>①, Andrew Vanderburg<sup>19</sup>①, and Peter J. Wheatley<sup>20,21</sup>①

Received 2023 March 17; revised 2023 September 5; accepted 2023 October 5; published 2023 October 24

## Abstract

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 \text{ g cm}^{-3}$ and a flat near-infrared transmission spectrum, hinting that it may possess circumplanetary rings. Given this planet's long orbital period ($P \approx 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.

*Unified Astronomy Thesaurus concepts: Exoplanet astronomy (486); Exoplanet dynamics (490); Exoplanet systems (484); Exoplanet tides (497); Natural satellites (Extrasolar) (483); Transits (1711); Transmission spectroscopy (2133)*

## 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). Moon-planet interactions may also be crucial for the long-term stability of exoplanet habitability; for example,

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

<footnote n="1"><sup>1</sup> Department of Astronomy, 501 Campbell Hall #3411, University of California, Berkeley, Berkeley, CA 94720, USA; charada@berkeley.edu</footnote>

<footnote n="2"><sup>2</sup> Earth & Planets Laboratory, Carnegie Institution for Science, 5241 Broad Branch Rd NW, Washington, DC 20015, USA</footnote>

<footnote n="3"><sup>3</sup> Department of Physics, Imperial College London, Prince Consort Road, London, SW7 2AZ, UK</footnote>

<footnote n="4"><sup>4</sup> Center for Astrophysics | Harvard & Smithsonian, 60 Garden St, Cambridge, MA 02138, USA</footnote>

<footnote n="5"><sup>5</sup> Division of Science, National Astronomical Observatory of Japan, Tokyo, Japan</footnote>

<footnote n="6"><sup>6</sup> Department of Astronomy and Astrophysics, University of California, Santa Cruz, 1156 High St, Santa Cruz, CA 95064, USA</footnote>

<footnote n="7"><sup>7</sup> Observatoire Astronomique de l'Université de Genève, Chemin Pegasi 51, CH-1290 Versoix, Switzerland</footnote>

<footnote n="8"><sup>8</sup> Instituto de Astrofísica e Ciências do Espaço, Universidade do Porto, CAUP, Rua das Estrelas, 4150-762 Porto, Portugal</footnote>

<footnote n="9"><sup>9</sup> Departamento de Física e Astronomia, Faculdade de Ciências, Universidade do Porto, Rua do Campo Alegre, 4169-007 Porto, Portugal</footnote>

<footnote n="10"><sup>10</sup> DTU Space, National Space Institute, Technical University of Denmark, Elektrovej 328, DK-2800 Kgs. Lyngby, Denmark</footnote>

<footnote n="11"><sup>11</sup> Centre for Exoplanet Science, SUPA School of Physics and Astronomy, University of St Andrews, North Haugh, St Andrews KY16 9SS, UK</footnote>

<footnote n="12"><sup>12</sup> Department of Physics and Astronomy, University of Kansas, Lawrence, KS, USA</footnote>

<footnote n="13"><sup>13</sup> Division of Geological and Planetary Sciences, 1200 E California Blvd, Pasadena, CA, 91125, USA</footnote>

<footnote n="14"><sup>14</sup> Department of Astronomy, California Institute of Technology, Pasadena, CA 91125, USA</footnote>

<footnote n="15"><sup>15</sup> Aix Marseille University, CNRS, CNES, LAM, Marseille, France</footnote>

<footnote n="16"><sup>16</sup> Centro de Astrobiología (CAB, CSIC-INTA), Depto. de Astrofísica, ESAC campus, 28692, Villanueva de la Cañada (Madrid), Spain</footnote>

<footnote n="17"><sup>17</sup> Centre for Star and Planet Formation, Natural History Museum of Denmark & Niels Bohr Institute, University of Copenhagen, Øster Voldgade 5-7, DK-1350 Copenhagen K., Denmark</footnote>

<footnote n="18"><sup>18</sup> School of Physics & Astronomy, University of Birmingham, Edgbaston, Birmingham B15 2TT, UK</footnote>

<footnote n="19"><sup>19</sup> Department of Physics and Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology, Cambridge, MA 02139, USA</footnote>

<footnote n="20"><sup>20</sup> Centre for Exoplanets and Habitability, University of Warwick, Gibbet Hill Rd, CV4 7AL Coventry, UK</footnote>

<footnote n="21"><sup>21</sup> Department of Physics, University of Warwick, Gibbet Hill Rd, CV4 7AL Coventry, UK</footnote>

<footnote n="22"><sup>22</sup> NSF Graduate Research Fellow.</footnote>

<footnote n="23"><sup>23</sup> NASA Sagan Fellow.</footnote>
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 \approx 8.9$), slightly evolved late F-type star ($T_{eff} = 6320$ K, log g = 4.294, [Fe/H] = -0.10, $M_{*}$ = 1.16 $M_{\odot}$; 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} \approx 2.9 R_{\oplus}$, $R_{p,c} \approx 2.6 R_{\oplus}$) 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} \approx 4.0 R_{\oplus}$, $R_{p,e} \approx 5.5 R_{\oplus}$, $R_{p,f} \approx 10 R_{\oplus}$), 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

<table><caption>**Table 1**

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

**Notes.**

<sup>a</sup> Assumes zero albedo and full heat redistribution.

<sup>b</sup> 95% credible upper limit.

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<sup>24</sup>

The largest and outermost planet in the system, HIP 41378 f, orbits near the inner edge of the system's habitable zone with a

<footnote n="24"><sup>24</sup> An additional transit of HIP 41378 f was observed in November 2022 by a coordinated ground-based effort.</footnote>
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 \approx 8.9$), slightly evolved late F-type star ($T_{eff} = 6320$ K, log g = 4.294, [Fe/H] = -0.10, $M_{*} = 1.16 M_{\odot}$; 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} \approx 2.9 R_{\oplus}$, $R_{p,c} \approx 2.6 R_{\oplus}$) 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} \approx 4.0 R_{\oplus}$, $R_{p,e} \approx 5.5 R_{\oplus}$, $R_{p,f} \approx 10 R_{\oplus}$), 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

<table><caption>**Table 1**

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

**Notes.**

<sup>a</sup> Assumes zero albedo and full heat redistribution.

<sup>b</sup> 95% credible upper limit.

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<sup>24</sup>

The largest and outermost planet in the system, HIP 41378 f, orbits near the inner edge of the system's habitable zone with a

<footnote n="24"><sup>24</sup> An additional transit of HIP 41378 f was observed in November 2022 by a coordinated ground-based effort.</footnote>
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 ($\rho = 0.09 \text{ 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 (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_⊕$ ($\rho = 1.2 \text{ 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 ($\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 \lesssim 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.

## 2. 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 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_⊕ \approx 0.0123$). We prescribe the satellite an Earth-like rocky composition ($\rho_s = 5.5 \text{ 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; Jagtap et al. 2021; Quarles et al. 2021), each
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 ($\rho = 0.09 \text{ 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 (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_⊕$ ($\rho = 1.2 \text{ 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 ($\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 \lesssim 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.

## 2. 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 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_⊕ \approx 0.0123$). We prescribe the satellite an Earth-like rocky composition ($\rho_s = 5.5 \text{ 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; Jagtap et al. 2021; Quarles 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:

1. the satellite collides with the planet,<fnref n="25" /> or

2. the satellite escapes the gravitational influence of the planet by crossing outside of the planet's Hill radius, defined as

$$R_H = a_p \left( \frac{M_p}{3M_*} \right)^{1/3}, \quad (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} \approx 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_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 \approx 0.03$ au $\approx 76R_p$). The planet and satellite are both initialized with their ascending node longitudes and pericenter arguments set to zero ($\Omega = \omega = 0^\circ$), and the mean anomaly of the planet is also set to zero ($M_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<fnref n="26" /> (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 simulations<fnref n="27" /> with the

<table><caption>**Table 2**
Parameter Distributions Used for N-body Simulations</caption><thead><tr><th>Parameter</th><th>Distribution (3-body)</th><th>Distribution (4-body)</th></tr></thead><tbody><tr><td>f</td><td></td><td></td></tr><tr><td>Mass, M<sub>f</sub> (M<sub>⊕</sub>)</td><td>12.3</td><td>T(12.3, 3.1, 0.0, ∞)</td></tr><tr><td>Orbital period, P<sub>f</sub> (days)</td><td>542.0798</td><td>N(542.0798, 0.0001)</td></tr><tr><td>Eccentricity, e<sub>f</sub></td><td>0.035</td><td>T(0.005, 0.005, 0.0, 1.0)</td></tr><tr><td>Mean anomaly, M<sub>f</sub></td><td>U(0, 2π)</td><td>U(0, 2π)</td></tr><tr><td>Argument of periastron, ω<sub>f</sub> (°)</td><td>0.0</td><td>T(231, 120, 0.0, 360)</td></tr><tr><td>Inclination, i<sub>f</sub> (°)</td><td>90.0</td><td>T(89.97, 0.01, 0.0, 90.0)</td></tr><tr><td>Planet e</td><td></td><td></td></tr><tr><td>Mass, M<sub>e</sub> (M<sub>⊕</sub>)</td><td>N/A</td><td>T(12, 5, 0.0, ∞)</td></tr><tr><td>Orbital period, P<sub>e</sub> (days)</td><td>N/A</td><td>[344, 394)</td></tr><tr><td>Eccentricity, e<sub>e</sub></td><td>N/A</td><td>[0.0, 0.3)</td></tr><tr><td>Mean anomaly, M<sub>e</sub></td><td>N/A</td><td>U(0, 2π)</td></tr><tr><td>Argument of periastron, ω<sub>e</sub> (°)</td><td>N/A</td><td>T(114, 55, 0, 360)</td></tr><tr><td>Inclination, i<sub>e</sub> (°)</td><td>N/A</td><td>T(89.84, 0.05, 0.0, 90.0)</td></tr><tr><td>Satellite</td><td></td><td></td></tr><tr><td>Mass, M<sub>s</sub> (M<sub>⊕</sub>)</td><td>0.15</td><td>0.15</td></tr><tr><td>Semimajor axis, a<sub>s</sub> (R<sub>H</sub>)</td><td>[0.1, 0.8)</td><td>U(0.1, 0.4)</td></tr><tr><td>Eccentricity, e<sub>s</sub></td><td>0.0</td><td>0.0</td></tr><tr><td>Mean anomaly, M<sub>s</sub></td><td>U(0, 2π)</td><td>U(0, 2π)</td></tr><tr><td>Argument of periastron, ω<sub>s</sub> (°)</td><td>0.0</td><td>0.0</td></tr><tr><td><sup>a</sup>Cosine of inclination, cos(i<sub>s</sub>)</td><td>[0, 1)</td><td>U(0, 1)</td></tr></tbody></table>

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

<sup>a</sup> The inclination of the satellite's orbit, $i_s$, is defined from the reference plane of the planet's orbit.

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 stability of each cell, $f_{stab}$, is calculated as the fraction of simulations (out of 25) that satisfy our stability criteria defined above.

In this setup, each grid cell corresponds to a child node of the quadtree which occupies $(1/4) \times (1/4)$ of the total initial parameter space $\{\cos i_s\} \times \{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:

1. the fractional stability of the node differs from that of a neighboring node by at least 20%,<sup>28</sup> and

2. each of the resulting child nodes occupies at least $(1/32) \times (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 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 \times 4 \times 25 = 400$

<footnote n="25">²⁵ We also considered imposing an inner stability bound defined by the planet's Roche limit radius, $R_{Roche} \approx 2.4R_s(M_p/M_*)^{1/3}$. However, due to the planet's low observed mass, the conventional definition of the Roche limit yields a nonphysical distance smaller than the planet's measured radius. Without further information about the planet's interior structure or whether it possesses a circumplanetary ring, we chose to set the inner separation bound equal to the planet's observed radius instead.</footnote>

<footnote n="26">²⁶ https://github.com/CalebHarada/astroQTpy</footnote>

<footnote n="27">²⁷ We experimented using higher numbers of simulations per grid cell but found that this has no significant impact on our results.</footnote>

<footnote n="28">²⁸ We replicated our simulations using a more stringent threshold of 10% but found that this had no major impact on our results.</footnote>
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:

1. the satellite collides with the planet,<fnref n="25" /> or

2. the satellite escapes the gravitational influence of the planet by crossing outside of the planet's Hill radius, defined as

$$R_H = a_p \left( \frac{M_p}{3M_*} \right)^{1/3}, \quad (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} \approx 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_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 \approx 0.03$ au $\approx 76R_p$). The planet and satellite are both initialized with their ascending node longitudes and pericenter arguments set to zero ($\Omega = \omega = 0^\circ$), and the mean anomaly of the planet is also set to zero ($M_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<fnref n="26" /> (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 simulations<fnref n="27" /> with the

<table><caption>**Table 2**
Parameter Distributions Used for N-body Simulations</caption><thead><tr><th>Parameter</th><th>Distribution (3-body)</th><th>Distribution (4-body)</th></tr></thead><tbody><tr><td>f</td><td></td><td></td></tr><tr><td>Mass, M<sub>f</sub> (M<sub>⊕</sub>)</td><td>12.3</td><td>T(12.3, 3.1, 0.0, ∞)</td></tr><tr><td>Orbital period, P<sub>f</sub> (days)</td><td>542.0798</td><td>N(542.0798, 0.0001)</td></tr><tr><td>Eccentricity, e<sub>f</sub></td><td>0.035</td><td>T(0.005, 0.005, 0.0, 1.0)</td></tr><tr><td>Mean anomaly, M<sub>f</sub></td><td>U(0, 2π)</td><td>U(0, 2π)</td></tr><tr><td>Argument of periastron, ω<sub>f</sub> (°)</td><td>0.0</td><td>T(231, 120, 0.0, 360)</td></tr><tr><td>Inclination, i<sub>f</sub> (°)</td><td>90.0</td><td>T(89.97, 0.01, 0.0, 90.0)</td></tr><tr><td>Planet e</td><td></td><td></td></tr><tr><td>Mass, M<sub>e</sub> (M<sub>⊕</sub>)</td><td>N/A</td><td>T(12, 5, 0.0, ∞)</td></tr><tr><td>Orbital period, P<sub>e</sub> (days)</td><td>N/A</td><td>[344, 394)</td></tr><tr><td>Eccentricity, e<sub>e</sub></td><td>N/A</td><td>[0.0, 0.3)</td></tr><tr><td>Mean anomaly, M<sub>e</sub></td><td>N/A</td><td>U(0, 2π)</td></tr><tr><td>Argument of periastron, ω<sub>e</sub> (°)</td><td>N/A</td><td>T(114, 55, 0, 360)</td></tr><tr><td>Inclination, i<sub>e</sub> (°)</td><td>N/A</td><td>T(89.84, 0.05, 0.0, 90.0)</td></tr><tr><td>Satellite</td><td></td><td></td></tr><tr><td>Mass, M<sub>s</sub> (M<sub>⊕</sub>)</td><td>0.15</td><td>0.15</td></tr><tr><td>Semimajor axis, a<sub>s</sub> (R<sub>H</sub>)</td><td>[0.1, 0.8)</td><td>U(0.1, 0.4)</td></tr><tr><td>Eccentricity, e<sub>s</sub></td><td>0.0</td><td>0.0</td></tr><tr><td>Mean anomaly, M<sub>s</sub></td><td>U(0, 2π)</td><td>U(0, 2π)</td></tr><tr><td>Argument of periastron, ω<sub>s</sub> (°)</td><td>0.0</td><td>0.0</td></tr><tr><td><sup>a</sup>Cosine of inclination, cos(i<sub>s</sub>)</td><td>[0, 1)</td><td>U(0, 1)</td></tr></tbody></table>

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

<sup>a</sup> The inclination of the satellite's orbit, $i_s$, is defined from the reference plane of the planet's orbit.

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 stability of each cell, $f_{stab}$, is calculated as the fraction of simulations (out of 25) that satisfy our stability criteria defined above.

In this setup, each grid cell corresponds to a child node of the quadtree which occupies $(1/4) \times (1/4)$ of the total initial parameter space $\{\cos i_s\} \times \{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:

1. the fractional stability of the node differs from that of a neighboring node by at least 20%,<fnref n="28" /> and

2. each of the resulting child nodes occupies at least $(1/32) \times (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 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 \times 4 \times 25 = 400$

<footnote n="25">²⁵ We also considered imposing an inner stability bound defined by the planet's Roche limit radius, $R_{Roche} \approx 2.4R_s(M_p/M_*)^{1/3}$. However, due to the planet's low observed mass, the conventional definition of the Roche limit yields a nonphysical distance smaller than the planet's measured radius. Without further information about the planet's interior structure or whether it possesses a circumplanetary ring, we chose to set the inner separation bound equal to the planet's observed radius instead.</footnote>

<footnote n="26">²⁶ https://github.com/CalebHarada/astroQTpy</footnote>

<footnote n="27">²⁷ We experimented using higher numbers of simulations per grid cell but found that this has no significant impact on our results.</footnote>

<footnote n="28">²⁸ We replicated our simulations using a more stringent threshold of 10% but found that this had no major impact on our results.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** 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 semimajor axis $a_s$ (expressed in terms of the Hill radius, $R_H$) and the cosine of the satellite's inclination $\cos i_s$ (with respect to the planet's reference plane). The color of each grid cell corresponds to the fraction $f_{\text{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 i_s \gtrsim 0.8$) tend to become unstable due to Kozai–Lidov oscillations, and there are no stable orbits for $\cos i_s \gtrsim 0.95$.</figcaption></figure>

and $32 \times 32 \times 25 = 25,600$. Figure 1 shows the final stability map for this setup, which we will discuss further in Section 3.

### 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$.

With the addition of planet *e* to our simulations, we update the stability criteria defined in Section 2.1.1, such that a system
<figure><img src="image_5.png" /><figcaption>**Figure 1.** 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 semimajor axis $a_s$ (expressed in terms of the Hill radius, $R_H$) and the cosine of the satellite's inclination $\cos i_s$ (with respect to the planet's reference plane). The color of each grid cell corresponds to the fraction $f_{\text{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 i_s \gtrsim 0.8$) tend to become unstable due to Kozai–Lidov oscillations, and there are no stable orbits for $\cos i_s \gtrsim 0.95$.</figcaption></figure>

and $32 \times 32 \times 25 = 25,600$. Figure 1 shows the final stability map for this setup, which we will discuss further in Section 3.

### 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$.

With the addition of planet e to our simulations, we update the stability criteria defined in Section 2.1.1, such that a system
<figure><img src="image_4.png" /><figcaption>**Figure 2.** 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 semimajor 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.</figcaption></figure>

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,<fnref n="29" /> this limit is set by

$$P_{e,crit}(e_e) \approx \frac{2\pi}{\sqrt{GM_*}} \left( \frac{a_d}{1 - e_e} \right)^{3/2}, \quad (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} \approx 351$ days, which is only within ~1.8σ of the measured period of planet e ($P_e = 369 \pm 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^5$ 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_{\text{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.

## 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_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,” $\tau$, 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 $\tau$ (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 $\Omega_i$, and their two obliquities $\psi_i$ (where $i = 1,2$ correspond to the primary and

<footnote n="29">²⁹ This is consistent with Santerne et al. (2019), who report a measured eccentricity for planet d of $e_d = 0.06 \pm 0.06$.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 2.** 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 semimajor 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.</figcaption></figure>

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,<fnref n="29" /> this limit is set by

$$P_{e,crit}(e_e) \approx \frac{2\pi}{\sqrt{GM_*}} \left( \frac{a_d}{1 - e_e} \right)^{3/2}, \quad (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} \approx 351$ days, which is only within ~1.8σ of the measured period of planet e ($P_e = 369 \pm 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^5$ 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_{\text{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.

## 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_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,” $\tau$, 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 $\tau$ (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 $\Omega_i$, and their two obliquities $\psi_i$ (where $i = 1,2$ correspond to the primary and

<footnote n="29">²⁹ This is consistent with Santerne et al. (2019), who report a measured eccentricity for planet d of $e_d = 0.06 \pm 0.06$.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** 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 ($P_e$) and eccentricity ($e_e$). As in Figure 1, 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 $P_{e,crit}(e_e)$, the theoretical upper limit on planet e's eccentricity vs. period due to the orbit of planet d (see Equation (2))—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 from Santerne et al. (2019) are shown by the green point with error bars representing $1\sigma$ 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.</figcaption></figure>

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:

$$ \frac{de}{dt} = \frac{11ae}{2GM_1M_2} \times \sum_{i=1}^{2} Z_i^{\text{ctl}} \left( \cos(\psi_i) \frac{f_4(e)}{\beta^{10}(e)} \frac{\Omega_i}{n} - \frac{18}{11} \frac{f_3(e)}{\beta^{13}(e)} \right), \quad (3) $$

$$ \frac{da}{dt} = \frac{2a^2}{GM_1M_2} \times \sum_{i=1}^{2} Z_i^{\text{ctl}} \left( \cos(\psi_i) \frac{f_2(e)}{\beta^{12}(e)} \frac{\Omega_i}{n} - \frac{f_1(e)}{\beta^{15}(e)} \right), \quad (4) $$

$$ \frac{d\Omega_i}{dt} = \frac{Z_i^{\text{ctl}}}{2M_i r_{g,i}^2 R_i^2 n} \times \left( 2\cos(\psi_i) \frac{f_2(e)}{\beta^{12}(e)} - [1 + \cos^2(\psi_i)] \frac{f_5(e)}{\beta^9(e)} \frac{\Omega_i}{n} \right), \quad (5) $$

and

$$ \frac{d\psi_i}{dt} = \frac{Z_i^{\text{ctl}} \sin(\psi_i)}{2M_i r_{g,i}^2 R_i^2 n \Omega_i} \times \left( \left[ \cos(\psi_i) - \frac{\xi_i}{\beta} \right] \frac{f_5(e)}{\beta^9(e)} \frac{\Omega_i}{n} - 2 \frac{f_2(e)}{\beta^{12}(e)} \right), \quad (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^{\text{ctl}}$ and $\xi_i$ are given by

$$ Z_i^{\text{ctl}} \equiv 3G^2 k_{2,i} M_j^2 (M_i + M_j) \frac{R_i^5}{a^9 \tau_i}, \quad (7) $$

$$ \xi_i \equiv \frac{r_{g,i}^2 R_i^2 \Omega_i an}{GM_j}, \quad (8) $$

where $k_{2,i}$ are the second-order tidal Love numbers, $\tau_i$ are the tidal time lags, and the subscripts $i$ and $j$ refer to the two bodies. The remaining quantities are given by:

$$ \begin{aligned} \beta(e) &= \sqrt{1 - e^2}, \\ f_1(e) &= 1 + \frac{31}{2}e^2 + \frac{255}{8}e^4 + \frac{185}{16}e^6 + \frac{25}{64}e^8, \\ f_2(e) &= 1 + \frac{15}{2}e^2 + \frac{45}{8}e^4 + \frac{5}{16}e^6, \\ f_3(e) &= 1 + \frac{15}{4}e^2 + \frac{15}{8}e^4 + \frac{5}{64}e^6, \\ f_4(e) &= 1 + \frac{3}{2}e^2 + \frac{1}{8}e^4, \\ f_5(e) &= 1 + 3e^2 + \frac{3}{8}e^4. \end{aligned} $$
<figure><img src="image_7.png" /><figcaption>**Figure 3.** 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 ($P_e$) and eccentricity ($e_e$). As in Figure 1, 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 $P_{e,crit}(e_e)$, the theoretical upper limit on planet e's eccentricity vs. period due to the orbit of planet d (see Equation (2))—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 from Santerne et al. (2019) are shown by the green point with error bars representing $1\sigma$ 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.</figcaption></figure>

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:

$$ \frac{de}{dt} = \frac{11ae}{2GM_1M_2} \times \sum_{i=1}^{2} Z_i^{\text{ctl}} \left( \cos(\psi_i) \frac{f_4(e)}{\beta^{10}(e)} \frac{\Omega_i}{n} - \frac{18}{11} \frac{f_3(e)}{\beta^{13}(e)} \right), \quad (3) $$

$$ \frac{da}{dt} = \frac{2a^2}{GM_1M_2} \times \sum_{i=1}^{2} Z_i^{\text{ctl}} \left( \cos(\psi_i) \frac{f_2(e)}{\beta^{12}(e)} \frac{\Omega_i}{n} - \frac{f_1(e)}{\beta^{15}(e)} \right), \quad (4) $$

$$ \frac{d\Omega_i}{dt} = \frac{Z_i^{\text{ctl}}}{2M_i r_{g,i}^2 R_i^2 n} \times \left( 2\cos(\psi_i) \frac{f_2(e)}{\beta^{12}(e)} - [1 + \cos^2(\psi_i)] \frac{f_5(e)}{\beta^9(e)} \frac{\Omega_i}{n} \right), \quad (5) $$

and

$$ \frac{d\psi_i}{dt} = \frac{Z_i^{\text{ctl}} \sin(\psi_i)}{2M_i r_{g,i}^2 R_i^2 n \Omega_i} \times \left( \left[ \cos(\psi_i) - \frac{\xi_i}{\beta} \right] \frac{f_5(e)}{\beta^9(e)} \frac{\Omega_i}{n} - 2 \frac{f_2(e)}{\beta^{12}(e)} \right), \quad (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^{\text{ctl}}$ and $\xi_i$ are given by

$$ Z_i^{\text{ctl}} = 3G^2 k_{2,i} M_j^2 (M_i + M_j) \frac{R_i^5}{a^9 \tau_i}, \quad (7) $$

$$ \xi_i = \frac{r_{g,i}^2 R_i^2 \Omega_i an}{GM_j}, \quad (8) $$

where $k_{2,i}$ are the second-order tidal Love numbers, $\tau_i$ are the tidal time lags, and the subscripts $i$ and $j$ refer to the two bodies. The remaining quantities are given by:

$$ \begin{aligned} \beta(e) &= \sqrt{1 - e^2}, \\ f_1(e) &= 1 + \frac{31}{2}e^2 + \frac{255}{8}e^4 + \frac{185}{16}e^6 + \frac{25}{64}e^8, \\ f_2(e) &= 1 + \frac{15}{2}e^2 + \frac{45}{8}e^4 + \frac{5}{16}e^6, \\ f_3(e) &= 1 + \frac{15}{4}e^2 + \frac{15}{8}e^4 + \frac{5}{64}e^6, \\ f_4(e) &= 1 + \frac{3}{2}e^2 + \frac{1}{8}e^4, \\ f_5(e) &= 1 + 3e^2 + \frac{3}{8}e^4. \end{aligned} $$
Unlike hot Jupiters ($P \lesssim 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,*} \approx 0.26$ (e.g., Claret & Gimenez 1989). We assume a Love number of $k_{2,*} = 0.5$ and a constant time lag of $\tau_{*} = 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,<fnref n="30" /> 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_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, $\tau_f = 0.00766$ s (e.g., Tokadjian & Piro 2020). For reference, note that $\tau \approx 0.766$ s for Neptune-like planets and $\tau \approx 638$ s for rocky planets (Tokadjian & Piro 2020).

Under these assumptions, we evolve the CTL model forward for a duration of $3.1 \times 10^9$ 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 $\ll 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 $\tau_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 \approx 0.0525 R_H$) to $0.4 R_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 $\tau_f$ to be consistent with the estimates for the solar system gas giants, $\tau_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 $\tau_f$ to 0.01 and 100 times the value assumed for Jupiter $\tau_j$. We then repeated each of these simulations assuming three different planet obliquities, $\psi_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.

### 3. 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 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

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 semimajor axis, here we simulated HIP 41378 f with a $0.15 M_⊕$ Mars-size satellite, varying the semimajor 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_{\text{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_{\text{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_{\text{crit}} \approx 0.39 R_H$. Our results shown in Figure 1 indicate that this limit depends only weakly on the inclination of the satellite for $\cos i_s \lesssim 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

<footnote n="30">³⁰ Note that faster rotation will generally lead to tidal migration on shorter timescales, which we explore later in this section.</footnote>
Unlike hot Jupiters ($P \lesssim 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,*} \approx 0.26$ (e.g., Claret & Gimenez 1989). We assume a Love number of $k_{2,*} = 0.5$ and a constant time lag of $\tau_* = 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,<fnref n="30" /> 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_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, $\tau_f = 0.00766$ s (e.g., Tokadjian & Piro 2020). For reference, note that $\tau \approx 0.766$ s for Neptune-like planets and $\tau \approx 638$ s for rocky planets (Tokadjian & Piro 2020).

Under these assumptions, we evolve the CTL model forward for a duration of $3.1 \times 10^9$ 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 $\ll 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_\oplus$, $R_s = 0.53 R_\oplus$). 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 $\tau_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 \approx 0.0525 R_H$ to $0.4 R_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 $\tau_f$ to be consistent with the estimates for the solar system gas giants, $\tau_J = 0.00766$ s (Tokadjian & Piro 2020) and a rotation period similar to that of Jupiter, $P_{\text{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_{\text{rot,f}} = 24$ hr) and a fast planet rotation ($P_{\text{rot,f}} = 3$ hr). For each case, we also test two extreme cases of the constant time lag by setting $\tau_f$ to 0.01 and 100 times the value assumed for Jupiter $\tau_J$. We then repeated each of these simulations assuming three different planet obliquities, $\psi_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.

### 3. 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 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

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 semimajor axis, here we simulated HIP 41378 f with a $0.15 M_\oplus$ Mars-size satellite, varying the semimajor 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_{\text{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_{\text{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_{\text{crit}} \approx 0.39 R_H$. Our results shown in Figure 1 indicate that this limit depends only weakly on the inclination of the satellite for $\cos i_s \lesssim 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

<footnote n="30">³⁰ Note that faster rotation will generally lead to tidal migration on shorter timescales, which we explore later in this section.</footnote>
<figure><img src="image_6.png" /></figure>

**Figure 4.** 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 semimajor axis as a function of time over the lifetime of the system, starting from an initial orbital separation of $4R_p$. Bottom row: Final vs. initial moon semimajor 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): $\psi_p = 0°$ (aligned prograde orbit), $\psi_p = 90°$ (polar orbit), and $\psi_p = 180°$ (aligned retrograde orbit). Each model assumes a different constant time lag $\tau$ for the planet, indicated by color (where $\tau_J = 0.00766$ s is roughly Jupiter's value). The fiducial planet rotation period is assumed to be 10 hr (solid lines), but we also consider a slow rotation case ($P_{\text{rot}} = 24$ hr; dashed lines) and a fast rotation case ($P_{\text{rot}} = 3$ hr; dotted–dashed lines). The gray hatched regions highlight regions where the moon tends to become dynamically unstable, $a_s \gtrsim 0.41 R_H$ (see Figure 1). Note that the ratio of the final-to-initial moon semimajor axis approaches unity as the starting semimajor axis increases in each obliquity case we considered. This is consistent with the very strong dependence of $\tau_{\text{mig}}$ on the semimajor axis in Equation (9) for the zero-obliquity case (i.e., $\tau_{\text{mig}} \propto a^8$). Note also that for the retrograde ($\psi_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.

process starts to become important for $\cos i_s \gtrsim 0.8$ and dominates for $\cos i_s \gtrsim 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's spin axis should be low because of momentum conservation from the rotating circumplanetary disk. However, other formation scenarios, such as a kinetic impact or dynamical capture, may result in moons with more misaligned inclinations, which we have 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 our second suite of simulations, which more realistically represented the HIP 41378 system with the inclusion of planet e. As the observed orbit of planet f is nearly circular, for each simulation we randomly selected the initial semimajor 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_{\text{stab}}$ by taking the fraction of 25 simulations within a given quadtree cell spanning $\{e_c\}$ and

$\{P_e\}$ that survived over $10^5$ yr. Our results are shown in Figure 3, where again $f_{\text{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 subsequently 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 $\gtrsim 2\sigma$ 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 semimajor
<figure><img src="image_8.png" /></figure>

**Figure 4.** 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 semimajor axis as a function of time over the lifetime of the system, starting from an initial orbital separation of $4R_p$. Bottom row: Final vs. initial moon semimajor 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): $\psi_p = 0°$ (aligned prograde orbit), $\psi_p = 90°$ (polar orbit), and $\psi_p = 180°$ (aligned retrograde orbit). Each model assumes a different constant time lag $\tau$ for the planet, indicated by color (where $\tau_J = 0.00766$ s is roughly Jupiter's value). The fiducial planet rotation period is assumed to be 10 hr (solid lines), but we also consider a slow rotation case ($P_{\text{rot}} = 24$ hr; dashed lines) and a fast rotation case ($P_{\text{rot}} = 3$ hr; dotted–dashed lines). The gray hatched regions highlight regions where the moon tends to become dynamically unstable, $a_s \gtrsim 0.41 R_H$ (see Figure 1). Note that the ratio of the final-to-initial moon semimajor axis approaches unity as the starting semimajor axis increases in each obliquity case we considered. This is consistent with the very strong dependence of $\tau_{\text{mig}}$ on the semimajor axis in Equation (9) for the zero-obliquity case (i.e., $\tau_{\text{mig}} \propto a^8$). Note also that for the retrograde ($\psi_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.

process starts to become important for $\cos i_s \gtrsim 0.8$ and dominates for $\cos i_s \gtrsim 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's spin axis should be low because of momentum conservation from the rotating circumplanetary disk. However, other formation scenarios, such as a kinetic impact or dynamical capture, may result in moons with more misaligned inclinations, which we have 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 our second suite of simulations, which more realistically represented the HIP 41378 system with the inclusion of planet e. As the observed orbit of planet f is nearly circular, for each simulation we randomly selected the initial semimajor 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_{\text{stab}}$ by taking the fraction of 25 simulations within a given quadtree cell spanning $\{e_c\}$ and

$\{P_e\}$ that survived over $10^5$ yr. Our results are shown in Figure 3, where again $f_{\text{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 subsequently 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 $\gtrsim 2\sigma$ 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 semimajor
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<sub>H</sub> 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 multiplanet 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 multiplanet systems.

## 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_2$ and $\tau$. Though both $k_2$ and $\tau$ 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 $\tau$ 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$, $\tau$, $r_g$, and $\Omega$ for the planet, as described earlier in Section 2. We subsequently tested cases where $\tau$ 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 $\psi_p = 0^\circ$, $\psi_p = 90^\circ$, and $\psi_p = 180^\circ$. For the satellite, we assumed tidally synchronous rotation and $k_2$, $\tau$, 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$ ($\approx 0.0525R_H$) to $30.5R_p$ ($\approx 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 $\le 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 ($\psi_p = 180^\circ$) or if $\tau$ or $\Omega$ 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 $\psi_p = 90^\circ$ 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 ($\psi_p = 90^\circ$) 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 ($\psi_p = 0^\circ$), 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_{\text{mig}} \sim a/\dot{a}$ (where $\dot{a}$ is the time derivative of $a$). Under the assumptions that $\psi_i = e = 0$, $M_f \gg M_s$, and $\Omega_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:

$$ \tau_{\text{mig}} \propto \frac{a^8}{k_{2,p}\tau_p M_s R_p^5} \left( \frac{\Omega_p}{n} - 1 \right)^{-1}, \quad (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 $\beta$ ($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 $\psi_p = 0^\circ$: larger values of $\tau$ cause the semimajor axis to increase on faster timescales than smaller values of $\tau$, and shorter planet rotation periods (larger $\Omega$) 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 $\tau$ 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 $\tau$ 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 $\tau$ be roughly 2 orders of magnitude larger than the value estimated for the solar system giants. However, we cannot rule out this possibility because $\tau$ 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
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<sub>H</sub> 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 multiplanet 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 multiplanet systems.

## 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_2$ and $\tau$. Though both $k_2$ and $\tau$ 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 $\tau$ 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$, $\tau$, $r_g$, and $\Omega$ for the planet, as described earlier in Section 2. We subsequently tested cases where $\tau$ 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 $\psi_p = 0^\circ$, $\psi_p = 90^\circ$, and $\psi_p = 180^\circ$. For the satellite, we assumed tidally synchronous rotation and $k_2$, $\tau$, 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$ ($\approx 0.0525R_H$) to $30.5R_p$ ($\approx 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 $\le 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 ($\psi_p = 180^\circ$) or if $\tau$ or $\Omega$ 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 $\psi_p = 90^\circ$ 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 ($\psi_p = 90^\circ$) 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 ($\psi_p = 0^\circ$), 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_{\text{mig}} \sim a/\dot{a}$ (where $\dot{a}$ is the time derivative of $a$). Under the assumptions that $\psi_i = e = 0$, $M_f \gg M_s$, and $\Omega_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:

$$ \tau_{\text{mig}} \propto \frac{a^8}{k_{2,p}\tau_p M_s R_p^5} \left( \frac{\Omega_p}{n} - 1 \right)^{-1}, \quad (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 $\beta$ ($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 $\psi_p = 0^\circ$: larger values of $\tau$ cause the semimajor axis to increase on faster timescales than smaller values of $\tau$, and shorter planet rotation periods (larger $\Omega$) 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 $\tau$ 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 $\tau$ 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 $\tau$ be roughly 2 orders of magnitude larger than the value estimated for the solar system giants. However, we cannot rule out this possibility because $\tau$ 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_{mig} \propto R_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_1$ and $q_2$ (Kipping 2013), and an offset $\Delta 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, 2019) implemented in the UltraNest<fnref n="31" /> 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 $\Delta T_0$ of $\pm 0.1$ days.<fnref n="32" /> 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

<footnote n="31">³¹ https://johannesbuchner.github.io/UltraNest/</footnote>

<footnote n="32">³² Alam et al. (2022) reported a much smaller uncertainty in $T_0$ of about 0.002 days.</footnote>
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_{mig} \propto R_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_1$ and $q_2$ (Kipping 2013), and an offset $\Delta 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, 2019) implemented in the UltraNest<fnref n="31" /> 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 $\Delta T_0$ of $\pm 0.1$ days.<fnref n="32" /> 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

<footnote n="31">³¹ https://johannesbuchner.github.io/UltraNest/</footnote>

<footnote n="32">³² Alam et al. (2022) reported a much smaller uncertainty in $T_0$ of about 0.002 days.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** HST white light curve of HIP 41378 f from Alam et al. (2022; points with error bars) and best-fit planet-only *Pandora* model (red dashed line). The black lines show simulated light curves with satellites injected at arbitrary orbital phases for ease of comparison—the solid line shows the signal from a 0.53 $R_⊕$ moon, while the dashed-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.</figcaption></figure>

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 Reactive Nested-Sampler* 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.

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 \approx 11 R_p \approx 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 ($\lesssim 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 $\gg 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_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 planet-only and planet-moon transit models, and residual systematic variations at the $\gtrsim 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 \sim 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
<figure><img src="image_9.png" /><figcaption>**Figure 5.** HST white light curve of HIP 41378 f from Alam et al. (2022; points with error bars) and best-fit planet-only *Pandora* model (red dashed line). The black lines show simulated light curves with satellites injected at arbitrary orbital phases for ease of comparison—the solid line shows the signal from a 0.53 $R_⊕$ moon, while the dashed-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.</figcaption></figure>

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 Reactive Nested-Sampler* 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.

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 \approx 11 R_p \approx 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 ($\lesssim 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 $\gg 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_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 planet-only and planet-moon transit models, and residual systematic variations at the $\gtrsim 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 \sim 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^2_r \approx 1.05$; Alam et al. 2022).<fnref n="33" /> 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.<fnref n="34" />

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.

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 $R_⊕$) 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, $\tau_{boil}$, equal to the system age (~3.1 Gyr; Santerne et al. 2019), and solved for the mean molecular weight, $\mu$, using the following set of equations:

$$ \tau_{\text{boil}} \sim \frac{M_{\text{atmo}}}{\dot{M}_{\text{Parker}}} \qquad (10) $$

$$ \dot{M}_{\text{Parker}} = 4\pi r_s^2 c_s \rho_0 \exp\left(\frac{3}{2} - \frac{2r_s}{R_{\text{sat}}}\right). \qquad (11) $$

Here, $c_s = (kT/\mu)^{1/2}$ is the isothermal sound speed, where k is the Boltzmann constant and T is the isothermal temperature; $r_s = GM_{\text{sat}}/(2c_s^2)$ is the sonic radius; and $\rho_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 $\mu \sim 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:

$$ (R_p/R_*)_{\text{eff}} = \sqrt{\delta_p + \delta_s}, \qquad (12) $$

where $(R_p/R_*)_{\text{eff}}$ is the effective planet-to-star radius ratio due to the planet and moon, and $\delta_p$ and $\delta_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 ($\chi^2_r \approx 1.55$) or the flat line ($\chi^2_r \approx 1.03$), given the data. Here we have also shown that the addition of the 0.53 $R_⊕$ exomoon with a $\mu \sim 4.1$ atmosphere is

<footnote n="33">³³ We note that Edwards et al. (2022) report a different transmission spectrum from their independent analysis of the HST/WFC3 data. We discuss this point later on in this section</footnote>

<footnote n="34">³⁴ Note, however, that the high-metallicity clear scenario is probably the least likely atmospheric composition for this cool giant planet (Alam et al. 2022), and it is possible that the ground truth is some combination of the scenarios we have described (e.g., the planet could have rings and a moon, rings and hazes, rings and hazes and a moon, etc.). While in principle each of these scenarios could affect the observed transmission spectrum of the system in a unique way, it is beyond the scope of this paper to model each possible ring/haze/moon permutation.</footnote>
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^2_r \approx 1.05$; Alam et al. 2022).<fnref n="33" /> 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.<fnref n="34" />

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.

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 $R_⊕$) 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, $\tau_{boil}$, equal to the system age (~3.1 Gyr; Santerne et al. 2019), and solved for the mean molecular weight, $\mu$, using the following set of equations:

$$ \tau_{\text{boil}} \sim \frac{M_{\text{atmo}}}{\dot{M}_{\text{Parker}}} \qquad (10) $$

$$ \dot{M}_{\text{Parker}} = 4\pi r_s^2 c_s \rho_0 \exp\left(\frac{3}{2} - \frac{2r_s}{R_{\text{sat}}}\right). \qquad (11) $$

Here, $c_s = (kT/\mu)^{1/2}$ is the isothermal sound speed, where k is the Boltzmann constant and T is the isothermal temperature; $r_s = GM_{\text{sat}}/(2c_s^2)$ is the sonic radius; and $\rho_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 $\mu \sim 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:

$$ (R_p/R_*)_{\text{eff}} = \sqrt{\delta_p + \delta_s}, \qquad (12) $$

where $(R_p/R_*)_{\text{eff}}$ is the effective planet-to-star radius ratio due to the planet and moon, and $\delta_p$ and $\delta_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 ($\chi^2_r \approx 1.55$) or the flat line ($\chi^2_r \approx 1.03$), given the data. Here we have also shown that the addition of the 0.53 $R_⊕$ exomoon with a $\mu \sim 4.1$ atmosphere is

<footnote n="33">³³ We note that Edwards et al. (2022) report a different transmission spectrum from their independent analysis of the HST/WFC3 data. We discuss this point later on in this section</footnote>

<footnote n="34">³⁴ Note, however, that the high-metallicity clear scenario is probably the least likely atmospheric composition for this cool giant planet (Alam et al. 2022), and it is possible that the ground truth is some combination of the scenarios we have described (e.g., the planet could have rings and a moon, rings and hazes, rings and hazes and a moon, etc.). While in principle each of these scenarios could affect the observed transmission spectrum of the system in a unique way, it is beyond the scope of this paper to model each possible ring/haze/moon permutation.</footnote>
<figure><img src="image_8.png" /><figcaption>**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 300× 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.</figcaption></figure>

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
<figure><img src="image_10.png" /><figcaption>**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 300× 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.</figcaption></figure>

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.

## 4. Conclusion

In this work, we investigated the stability of a hypothetical large moon ($0.15 M_⊕$, $0.53 R_⊕$) orbiting the long-period ($P \approx 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 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_s$ for most (low) moon inclinations. However, Kozai-Lidov oscillations can start to become an important destabilizing influence for $\cos i_s \gtrsim 0.8$, and disallow stable orbits for the highest inclinations ($\cos i_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 < \cos i_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.4 R_⊕$ 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 ($\tau \gtrsim 0.7$ s) and fast rotation rate ($P \sim 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 \gtrsim 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
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.

## 4. Conclusion

In this work, we investigated the stability of a hypothetical large moon ($0.15 M_⊕$, $0.53 R_⊕$) orbiting the long-period ($P \approx 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 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_s$ for most (low) moon inclinations. However, Kozai-Lidov oscillations can start to become an important destabilizing influence for $\cos i_s \gtrsim 0.8$, and disallow stable orbits for the highest inclinations ($\cos i_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 < \cos i_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.4 R_⊕$ 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 ($\tau \gtrsim 0.7$ s) and fast rotation rate ($P \sim 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 \gtrsim 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. Through- out 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_⊕$ ($\rho_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 $\rho_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 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^{-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.

## Acknowledgments

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 “Next-GenerationEU”/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,