Exomoon Phase Curves: Toroidal Exosphere Simulations of Exo‐Ios Orbiting 8 Exoplanets in Alkali Spectroscopy

Toroidal atmospheres and exospheres characterized at exoplanets may be fueled by volcanically active exomoons, often referred to as exo‐Ios. We study the neutral outgassing and volatile evolution of a close‐orbiting, evaporating satellite at eight candidate exoplanet‐exomoon systems WASP‐49,‐96,‐69,‐17 b, XO‐2N b, HAT‐P‐1 b, HD‐189733 b, and HD‐209458 b by developing a 3‐D test‐particle Monte Carlo simulation, Simulating the Evolution of Ring Particles Emergent from Natural Satellites. The module is coupled to dishoom, approximating the minimum mass‐flux needed to reproduce observations of alkali line profiles identified in dozens of transmission spectra. We focus on sputtered neutral sodium, limited by photoionization and radiative effects. By considering Earth‐, Io‐, and Enceladus‐like masses, we systematically simulate the imprint of a non‐hydrostatic medium (characteristic of volcanic exospheres) in density and velocity space using a novel Delaunay tesselation field estimator algorithm. Our results demonstrate how exomoons can considerably modulate gas density observations probed near exoplanet transit, depending on the orbital phase of the putative satellite at the time of observation. The density evolution, therefore, manifests on orbital timescales as “exomoon phase curves” from shadow to occultation. We find two regimes of density evolution, characteristic of a: (a) localized cloud and (b) an azimuthally symmetric exoring/torus, degenerate with an exoplanet atmosphere, ranging from ∼109.5±0.5 cm−2 to ∼1015±0.25 cm−2 at our leading candidate WASP‐69b I. In certain orbital architectures, the smallest evaporating satellite mass surprisingly generates the brightest sodium signal, fueling optimism for discovering photometrically indiscernible rocky exomoons. We suggest long baseline monitoring of alkali and SO2 systems in spectroscopy to search for the temporal and spatial variability predicted here.


Introduction
In the last few years research in exoplanetary systems has significantly increased.Methods in acquiring data have become more sophisticated, which has lead to an increasing frequency of exoplanet discoveries (Christiansen, 2022).
Little is known about these exoplanets, apart from their fundamental mass, radius, and semi-major axis measurements (M p , R p , a p ), and due to the advent of high-resolution spectroscopy, signals of neutral and ionized gas, present near transit.However it is unclear at present whether this gas is uniquely from the exoplanet itself, or a small satellite venting gas as suggested two decades ago at HD-209458 b (Johnson & Huggins, 2006).Although • Our open access module Simulating the Evolution of Ring Particles Emergent from Natural Satellites allows one to simulate the evaporation of an exomoon and its ability to form a volcanic exoring/exosphere such exomoons, natural satellites beyond our Solar System, are still undetectable in photometry searches with state-of-the-art instruments, their existence has been hypothesized to explain various inconsistencies in transit timing variations (Teachey & Kipping, 2018) and exoplanetary spectra (Oza et al., 2019).In addition to HD-209458 b, exogenic sources of sodium and potassium have now been considered at multiple exoplanet systems.At evaporative transmission spectra altitudes, alkali abundances often exceed the exoplanet's source rates, and hydrostatic exoplanet atmosphere models remain limited in their ability to explain significant line broadening and Doppler behavior seen in Na & K spectra.Besides exomoons, exogenic sources may include short-period comets, regular cometary impacts, and micrometeoritic influx.However, each of these explanations comes with its own set of challenges.For instance, unbounded comets would yield sodium signatures regardless of whether the exoplanet is in transit or not, which contradicts observations.A fundamental issue that many exogenic sources fail to address is the existence of periodic Doppler shifts exceeding ∼10 km/s.These shifts exhibit characteristics that deviate from both hydrodynamic and hydrostatic behaviors (Gebek & Oza, 2020; H. J. Hoeijmakers et al., 2020), rather resembling velocity profiles typically associated with the atmospheric sputtering of sodium (Burger et al., 2001).
The Solar System's gas giants are well known for their magnetospheres and ring systems (Burns et al., 1999;Crida et al., 2019;Rawal, 1981), in particular Saturn.Even a satellite radius of R ☽ ∼ 250 km that is, Enceladus, is able to sustain Saturn's E-ring and a hydroxide torus, fueled by hydrothermal activity venting directly to space at Ṁ ∼ 200 kg/s (Johnson et al., 2006) (by convention our star-planet-satellite system is denoted: R ⋆ -R p -R ☽ ).In 2007, the ∼20 M J companion 1SWASP J1407 b was observed to exhibit multiple and sporadic transits in continuum light, believed to be due to a large and extended (a ≳ AU) exoring system fueled by exomoons (Kenworthy & Mamajek, 2015;Mamajek et al., 2012).Upon simulating moonlet formation and subsequent evolution with several eccentricities, it has recently been determined that the exoplanet may rather be a free-floating object nevertheless with an exoring-exomoon system (Sutton et al., 2022).The search for exomoons therefore has an inextricable avenue linked to exorings.
One way to search for the imprint of exomoons and exorings in transit spectroscopy is to study tenuous gas similar to the gas giant systems in our Solar System.This gas follows an optically thin, exospheric regime, where the gas is collisionless and the spectra are dominated by the energetics and kinetics of the moon-planet systems (Gebek & Oza, 2020;Oza et al., 2019).In this study we present a Monte Carlo code to simulate the neutral clouds and tori emerging from sputtering at the surface of evaporative exomoons.We track the evolution of ejected particle distributions in time, which manifest as a cloud or torus in transit.We first present exomoon candidates in Section 2.1 and outline our assumptions on plasma (Section 2.3) and radiative processes (Section 2.4) which are essential for our Monte Carlo model Simulating the Evolution of Ring Particles Emergent from Natural Satellites (SERPENS) in Section 3. We then apply our model results to eight candidate exoplanet systems in Section 2.1, and discuss these results in Section 5. Finally, we present an avenue for future modeling work and implications for transmission spectroscopy observations of exoplanets in Section 6.

Evaporating Exomoon Candidates
More than a dozen close-orbiting gas giant systems known to host alkali metals have now been shown based on analytical estimates using the dishoom model, to be consistent with their own close-orbiting extrasolar Ios, hereafter exo-Ios, by Oza et al. (2019).We seek to further simulate select candidate exomoon systems numerically and explore their ability to host exomoons by examining particle distributions emerging from our current understanding of the systems' environment.This study focuses on eight hot Jupiter & Saturn exoplanets, characterized by rapid orbital periods of P < 10 days: HD-189733 b, HD-209458 b, HAT-P-1 b, XO-2N b and multiple systems discovered by the Wide Angle Search for Planets (WASP) consortium, namely WASP-49 b, WASP-17 b, WASP-69 b, and WASP-96 b.Although long-period exoplanets provide a more suitable environment for stable satellite orbits over Gyrs than hot Jupiters, our study of close-in exoplanets allows for a high-cadence transit search, that is extremely time-efficient.This is due to the fact that if satellites did exist, they would orbit close to the exoplanet, and have orbital timescales that are often a fraction of a day.Although Barnes and O'Brien (2002) suggested close-in gas giant satellites may have been lost, dynamical studies soon after demonstrated satellite semi-major axes ≲0.49Hill radii remain stable (Domingos et al., 2006).Indeed, tidal dissipation that nominally destabilizes the satellite orbit, if induced into the satellite itself, stabilizes the orbit for Gyrs as estimated analytically by Cassidy et al. (2009).This is now confirmed numerically by Kisare and Fabrycky (2023), with a slightly smaller satellite orbital radius ≲0.41 Hill radii.It was shown that stellar tides have the capability to exert a stabilizing effect on these orbits, while driving enormous amounts of heat into the satellite, surpassing the tidal heat flux ( Ė) observed at Io by orders of magnitude (Oza et al., 2019).The observational consequence then, is that the mass loss ( Ṁ) of these putative yet explosive exo-Ios would be the principal factor determining exomoon survival, while also providing a source to evaporative transmission spectroscopy (Gebek & Oza, 2020).
The first detection of an extrasolar atmosphere was reported in 2001 for neutral sodium at HD-209458 b (Charbonneau et al., 2002).The authors found additional spectroscopic dimming of starlight in the sodium doublet (λ D2 = 5889.95Å, λ D1 = 5895.92Å) near the exoplanet transit.Upon this remarkable observation of neutral sodium, Johnson and Huggins (2006) considered the source material provided by an orbiting exomoon fueling a gas torus and/or debris ring which would be imprinted on the exoplanet transit spectra.Although Na I has now been detected at over 28 exoplanet systems as of this writing, the latest ESPRESSO/VLT observations of HD-209458b report a null detection of Na due to the stellar sodium lines rotating into the exoplanet's rest frame (Casasayas-Barris, N. et al., 2020).Besides Na I in the optical, the next brightest spectral line is hydrogen via its ultraviolet (UV) Lyman-α band, which has the ability to provide additional plasma constraints.Based upon energetic neutral atom (ENA) production and Lyman-α observations of the planet's transit spectra, Ekenbäck et al. (2010) suggested a magnetic standoff-distance of ∼(4.2-10.5)R HD209458b (see also Holmström et al. (2008)).Models of the magnetic dipole field of tidally locked exoplanets have difficulties to explain such an extended magnetic obstacle (see e.g.M. Khodachenko et al., 2007).An increase of plasma density resulting from an evaporating exomoon and a counterpart plasma torus, not unlike Jupiter-Io, may provide an avenue.
Early ingress absorption measurements of ionized carbon at roughly ∼10 29 ions/s observed by the Hubble Space Telescope (HST) at another hot Jupiter HD-189733 b led to the inference of a plasma torus fueled by an exomoon orbiting at ∼16 R p (Ben-Jaffel & Ballester, 2014).The fact that this satellite orbit lies beyond the exoplanet's Hillsphere, that is, beyond the range of stable gravitational capture from the exoplanet, renders this exomoon candidate HD-189733 b I unlikely.However, here we examine the idea of a second candidate, HD-189733 b II, orbiting far closer, near the Hill-sphere, responsible for continuously radially outflowing neutrals and plasma (Figure 1).Na I D lines have been spectrally resolved in WASP-49 b's upper atmosphere (M.Lendl et al., 2012) using the HARPS high-resolution spectrograph at the ESO 3.6 m telescope.Analysis has shown that the presence of hot neutral sodium at high altitudes is necessary to explain the signatures while the exoplanet's atmosphere was determined to be cloud-free (A.Wyttenbach et al., 2017).A sputtered exo-Io was shown to fit the data with a tidally heated profile identical to Io's (Burger et al., 2001;Oza et al., 2019).
From the detection of alkali metals in the transmission spectra of WASP-69 b, it was inferred that the system might be hosting an evaporating exomoon.Using the upgraded Giant Metrewave Radio Telescope (uGMRT), the system was probed for radio emissions at frequencies 150 and 218 MHz (Narang et al., 2023).Although only upper limits of 3.3 mJy (ν = 150 MHz) and 0.9 mJy (ν = 218 MHz) were attained with the uGMRT, the presence of an exomoon is still conceivable.Less mass-loss, smaller exoplanet magnetic fields or strong temporal variations in emission could explain the radio-quiet behavior while still allowing the existence of an exomoon.Current exoplanet magnetosphere models led Narang et al. (2023) to suggest the magnetic field strength to be roughly B p ∼ 15 G (Narang et al., 2023;Yadav & Thorngren, 2017).The candidate exoplanet-exomoon system WASP-69 b is set to be observed in eclipse by JWST cycle-1 GTO, although a complementary infrared observation in transit would be invaluable for further exomoon analysis.
WASP-17 b lies on the boundary between hot and ultra-hot Jupiters orbiting a 6650 K, F-star.It is the third planet at which an atmospheric sodium signature was found with the help of the GIRAFFE spectrograph on the VLT (Wood et al., 2011).Wood et al. (2011) found a transit depth of 0.55 ± 0.13% at 1.5 Å and a similar transit depth of 0.58 ± 0.13% at 1.5 Å with 4.5σ confidence was later measured using the MIKE spectrograph on the 6.5-m Magellan Telescope by Zhou and Bayliss (2012).Intriguingly, recent measurements with data from the HST did not find any evidence of such a sodium signature (Alderson et al., 2022).Most recently, observations with JWST's Mid-Infrared Instrument Low Resolution Spectrometer revealed the presence of small ∼0.01 μm particles composed of SiO 2 (s), at high altitudes, in a candidate exoplanet-exomoon system WASP-17 b (Grant et al., 2023).This finding is highly suggestive of a rocky exogenic source, and may give further credence to the candidate WASP-17 b I whose companion alkali signatures are studied in this work.We seek to test the existence of an exo-Io as the source of silicates and the apparently variable neutral sodium signature, since a satellite that gets periodically occulted by its host planet provides a natural explanation for a time-variable sodium signature (N.B. this system has the shortest photoionization lifetime in our sample: τ W17bI,γ ∼ 3.4 min (Huebner & Mukherjee, 2015)).
If sodium is found in exoplanetary spectra, ≫99% of all Na atoms are confined to the Na D-lines in absorption.Therefore, interpretations on its dynamical behavior (e.g., Doppler broadening) are deduced from its line shape mostly in high-resolution spectroscopy probing the most tenuous layers of gas at line center providing a strict minimum column density.One alkali exoplanet system in particular, WASP-96 b, is claimed to possess a thicker layer of sodium, at higher pressure due to broadening observed in low-resolution by VLT/FORS2 spectrograph (Nikolov et al., 2018).Although this profile was attributed to an endogenic cloud-free atmosphere (which reveals the deeper Na feature in low-res), exogenic sodium venting from an exomoon is plausible, however at a source rate larger than an Io-     Note.The radii correspond to Earth, Io and Enceladus, respectively.Masses are calculated with a density of ρ = 3.5 g/cm 3 .mass exomoon, motivating the test case for Earth-mass satellites, as well.
The cloudy and clear atmosphere interpretations relies on T-P profiles assuming hydrostatic equilibrium which present severe degeneracies in the exact gas pressure probed (Gebek & Oza, 2020).Two other exoplanets that host sodium signatures in low-resolution spectroscopy are HAT-P-1 b and XO-2N b.We want to include these exomoon host candidates for an increase in probing and comparison between our results.
where m i is the species mass, τ i its lifetime and R * the host star's radius.
We model three different masses of each exomoon.These correspond to an Earth-sized, Io-sized and Enceladussized body with density ρ rock ≈ 3.5 g/cm 3 .
We choose three masses and radii as a parameter study grounded on the satellite size, but also as a fundamental probe for exomoon evolution.Exomoons with large mass-loss rates will erode over time depending on their initial progenitor mass.The ultimate fate of exomoons is typically one of two scenarios: either they gradually disintegrate into exorings, as discussed by for example, Barnes and O'Brien (2002) and Charnoz et al. (2018), or they undergo significant changes due to outward migration or collision with their host planet, as suggested by Alvarado-Montes et al. (2017) and Sucerquia et al. (2019).Radial migration beyond the exoplanet Hill sphere does not necessarily seal the fate of a satellite.As shown by Hansen (2022), a rogue satellite is expected to fall into the gravitational well of its host star, after which it ultimately returns again as a re-captured exomoon.When we observe a close-in exomoon, its physical state likely falls somewhere between the cases we've outlined.More extreme cases may be witnessed such as the ongoing silicate disintegration of a volcanic exomoon, or plasma torus signatures of a largely destroyed system of metal-rich moonlets.

Tidal Parameters
Incorporating the effects of tidal dissipation, recent three-body simulations (Kisare & Fabrycky, 2023) corroborate and extend the insights from Domingos et al. (2006) and Cassidy et al. (2009).These studies collectively highlight the crucial role of tidal dissipation in keeping exomoons bound to their host planets, contrary to the   previous assumption that exomoons are prone to escape their primary's Hill sphere (see e.g.Dobos et al., 2021;Sucerquia et al., 2019;Teachey, 2021;Tokadjian & Piro, 2020).Gravitational tides are, therefore, one of the most alluring aspects of exomoon science for dynamicists.That is, since the (dynamical) stability of satellites is heavily contingent on their tidal dissipation factor Q (Cassidy et al., 2009;Oza et al., 2019), an exomoon discovery will provide a constraint on the three-body problem for a close-in planet star satellite system.Q is rather difficult to constrain, but generally expected to exceed 10 9 (Cassidy et al., 2009;Goldreich & Nicholson, 1977;Ogilvie, 2014;Ogilvie & Lin, 2004;Wu, 2005), a significantly higher value than the initial assumption of a Jovian Q ∼ 10 5 representative of our solar system's gas giants (Lainey & Tobie, 2005) used by Barnes and O'Brien (2002).Since close-in exomoons are expected to deposit significantly less tidal energy into the upper atmospheres of their host planets, their tidal Qs are consequently far larger (Oza et al., 2019) rendering exomoons to be more dynamically stable than previously thought.Furthermore, the results of Kisare and Fabrycky (2023) show that if tidal dissipation in the satellite is comparable or greater than that in their host planet, it can sustain itself at a substantial fraction of the Hill radius as its host planet's spin rate is decreasing.The spin rate of the host exoplanet is of great significance when simulating sputtering from the exomoon's surface or volcanic atmosphere as it directly influences the distribution of particle velocities.In this context, we employ fundamental approximations gleaned from Io.
Of course, in simulating orbiting gas venting from an exomoon, the velocity distribution of ejected particles is crucial (cf.Section 3.1).Speeds are usually in the range of a few km s 1 with high-energy tails ranging from ∼30-100 km s 1 at Io for atmospheric sputtering and charge exchange respectively (Wilson et al., 2002).Jupiter's rotation period exceeds the Keplerian velocity of Io, such that the co-rotating plasma (the magnetic field is synchronized to the rotation of Jupiter) hits Io with high-energy ions from the back, at the plasma ram/trailing hemisphere.The hot Jupiter systems we analyze are not expected to exhibit this behavior as they are expected to be tidally locked to the star and therefore have a much longer rotation period.In this way we are able to estimate key parameters in the ejected velocity distribution by scaling to the orbital period.
The timescale for tidal locking of an exoplanet at distance a p from the star can be estimated via the initial and final rotation speeds considering the rate of change of the planetary rotation speed ω: (2) Note.Ben-Jaffel and Ballester (2014) claim the existence of HD-189733 b I at ∼16 R p which is why we assign the present putative exomoon the Roman number II. a ☽ denotes the semi-major axis of the satellite in terms of planet radii.a sync is the synchronous orbit radius at which the exomoon has a period equal to the average rotational period of the exoplanet, assuming that the planet is tidally locked.All values lie beyond the Hill sphere, which implies that the orbits are not synchronous.Ṁ☽ is an estimate for the mass-loss rate using Equation 1 in kg/s.τ γ,Na represents the lifetime of neutral sodium atoms against photoionization, which we consider to be the most prominent loss term.Both v b and v M are parameters used in the velocity distribution (Equation 8) and determine most likely velocity and maximum velocity, respectively.β is the ratio between radiation pressure acceleration on particles to gravitational acceleration.These values are interpolated from results by M. L. Khodachenko et al. (2011).
Neglecting ω f and setting ω i to Jupiter's value ω J = 1.77 ⋅ 10 4 s 1 one can rewrite Equation 2 to be (see Grießmeier et al. (2007)), where α is constant and related to the internal mass distribution, Q ′ p = 3/ (2k 2,p ) Q the modified tidal dissipation Q-value for the exoplanet defined by MacDonald (1965) and Goldreich and Soter (1966), incorporating the planet's tidal Love number k 2,p (∼0.565 for Jupiter, ∼0.32 for Saturn).Following Gu et al. (2003) we take α to be the value corresponding to a large gaseous planet with equation-of-state polytropic index κ = 1, that is α = 0.26.Doing these calculations for the systems examined in this work, substantiates the claim of tidal locking as in all cases the synchronization timescale amounts to a few hundred million years (see Table 3).

Magnetosphere and Plasma Environment
In this work we focus on neutrals as the bright Na I and K I doublet lines observed thus far in the optical are not ionized.In fact, the presence of neutrals in such close proximity to the host star is unexpected, as they are typically not anticipated to survive under these conditions.Our code is intended to be coupled with a magnetohydrodynamics module (e.g., Fatemi et al. (2022)), however here we provide some basic guiding constraints.The gas giants considered are close-in orbiting exoplanets, with a semi-major axis of ∼0.05 AU for instance, the starplanet architecture may be categorized as outskirts of the stellar corona.In order to sustain a particle cloud or torus around the exoplanet, a sufficiently strong exoplanetary magnetic field may be necessary to shield the exomoon from extreme stellar plasma flow.The (minimal) extent of the exoplanet's magnetosphere is given by the standoff distance R mag .This distance can be calculated from the balance between stellar wind ram-and thermal pressure against the planetary magnetic field pressure at the substellar point.In the case of a simple dipole magnetic field B ∝ M/ r 3 the standoff distance may then be calculated as: with planetary magnetic dipole moment M, magnetosphere form factor f 0 , stellar wind mass density ρ, incident at velocity v eff (including an aberration effect due to the planet's Kepler velocity) and stellar wind thermal pressure 2nk B T (Chapman & Ferraro, 1931).Models of the magnetic dipole moment (Busse, 1976;Mizutani et al., 1992;Sano, 1993;Stevenson, 1983) depend on the planet's angular rotation speed.Tidally locked exoplanets therefore generally exhibit smaller dipoles and thus smaller standoff distances up to the point of inability to provide efficient protection for planetary upper atmospheres against strong non-thermal mass loss essential at short orbital distances from stars.M. L. Khodachenko et al. (2011) showed that the inclusion of a magnetodisk and the magnetotail helps to solve this problem.Plasma, escaping along field lines penetrating beyond the (exoplanetary) Alfvénic surface, deforms the planetary magnetic dipole field by radially stretching the field lines.The corresponding Alfvén-sphere's radius is the boundary above which the exoplanetary dipole field is no longer strong enough to drive the plasma into rigid corotation.Here the energy density of the rotational motion of plasma ρ A ω 2 p R 2 A / 2 equals the energy density of the dipole magnetic field μ 0 M 2 / 32π 2 R 6 A ) .The resulting magnetodisk needs to have a continuous load of plasma to be supported beyond the Alfvénic surface thought to be supplied by thermal mass loss of the hot Jupiter/Saturn, but in return increases the standoff distance up to 40%-70 %.In the case of the Solar System such a magnetodisk is formed at Jupiter by Io's volcanic activity.A hot Jupiter showing low mass-loss of its atmosphere while still supporting a magnetodisk might therefore be a strong sign for the existence of an exo-Io taking in the role of the plasma source for the magnetodisk.
In order to analyze if the exomoon resides inside the Alfvén sphere of the exoplanet, we need to estimate the thermal mass loss of the planet following the considerations of M. L. Khodachenko et al. (2011).Assuming a similarity between the magnetodisk geometry of the considered systems and Jupiter, one can estimate the Alfvénic surface equatorial radius R A as where k is the angular rotation scaling power that depends on the magnetic dipole moment model used and R AJ = 19.8,characteristic to Jupiter.For the Stevenson (1983) magnetic moment model we have k = 1/2, such that the ratio Equation 5 only depends on the thermal mass-loss rate.The planetary mass-loss rate Ṁp depends on the stellar XUV flux at the orbital distance of the exoplanet, which we compute according to the host star's spectral type from Lammer, H. et al. (2009).The precise XUV flux at a system is a subject of ongoing theoretical work on photoevaporation on a population of exoplanets (Affolter, L. et al., 2023;Mordasini, 2020).According to this calculation the exomoons orbit inside the Alfvén spheres even for heating efficiencies up to 30%, such that the plasma is driven into a rigid corotation with the planet at a ☽ < R mag similar to the case of Io in the Solar System (see Table 3).
Due to the slow rotation speed of the tidally locked gas giants, the Kepler velocity of the exomoon is greater than the corotating plasma velocity in all considered cases.These considerations lead us to the assumption that photoionization is the limiting factor for the neutral sodium lifetimes and that stellar winds should not play a substantial role for neutrals in the vicinity of the exomoon.Photoionization lifetimes of sodium are calculated with rates from Huebner and Mukherjee (2015) and Huang et al. (2017).

Radiative Processes
We consider radiation pressure to substitute an additional force characterized by the radiation parameter β describing the ratio of radiation acceleration to gravitational acceleration in our simulations, that is, F rad = βF grav .In the case of dust, radiation pressure coefficients β can be calculated from grain properties, luminosity & mass of the star, and Mie scattering coefficient (Burns et al., 1979;Hayakawa & Hansen, 2023).However, when dealing with atoms and molecules, a more advanced analysis is required to calculate the radiative effects, which involves examining the absorption line profile of the species involved.Thus, both the flux of photons and the relative speed of a particle to the star determine the magnitude of the pressure.The wavelengths of incident photons are blue-or red-shifted from the most efficient absorption wavelength in the rest frame of the atom/molecule.In the case of atomic hydrogen with its relatively narrow Lyman-α line, the Doppler effect takes an important role.For the exoplanet HD-209458 b the corresponding (Lyman-α) β-value is estimated to be of order ∼4, indicating that radiation pressure is dominating over stellar gravity (β > 1) (V.Bourrier & Lecavelier des Etangs, 2013).
Neglecting relativistic Doppler shifts one can estimate β to be constant.In the case of sodium, one needs to consider the stellar flux at the line centers of the D-lines, but finding a concrete estimate for a β value is difficult.However, owing to the short neutral lifetimes, we deem an extensive study of localized radiation parameters not to be necessary (see e.g.Koskinen et al., 2010).
Radiative forces also influence the extent of the Hill sphere.A closed analytical expression for the modified Hill radius was found by Beth et al. (2016) and rests upon the position of the L 2 Lagrange point.As a consequence of revisited Hill radius calculations, they found that the L 2 hydrogen point for HD-209458b lies within the expected exobase location of the planet.Hence, no thermal energy is needed for exospheric particles to escape at low exosphere temperatures.This mass loss could contribute to plasma enrichment of a magnetodisk.Similar studies of the Neptune-like planet GJ 436b (V.Bourrier et al., 2015;Ehrenreich et al., 2015) highlight a comet-like shape explainable with the effects of radiation (Beth et al., 2016).
We include a conical shadow behind the planet that effectively removes radiation pressure on every particle inside of it.As particles enter the shadow their assumed photoionization limited lifetime is significantly higher.If only photo-ionization from stellar radiation is present, the lifetime would approach infinity.However, other loss mechanisms including charge exchange will be present.We therefore dynamically increase the lifetimes of particles to a characteristic toroidal timescale of 3 hr while they reside inside the planet's shadow (the shadowcondition calculation can be found in Appendix A4).

Modeling Methods
We have created a Python module called SERPENS (Simulating the Evolution of Ring Particles Emergent from Natural Satellites) (Meyer zu Westram, 2023) that runs the weighted particle Monte Carlo simulations.
SERPENS is publicly available on GitHub.A link to the domain is provided at the end of this study.
Parameters for the simulation are saved inside a Python file and include the integration details, how often particles are injected, what kind of species we want to simulate, as well as their properties starting from mass up to lifetimes and chemical networks.When SERPENS starts a new simulation, the first step is the initialization of the main celestial bodies in a REBOUND N-body simulation instance.REBOUND is a N-body integrator supporting a multitude of symplectic and non-symplectic integrators that was first developed by Rein and Liu (2012).The coordinate system used for simulation is depicted in Figure 3.
Neutral particles that get added to the simulation are afterward treated as massless test particles that do not interact with each other.

Creating Particles
Constant mass loss by the source object is represented by adding new particles to the simulation before each integration advance.The ejection of particles is determined by the velocity distribution of the ejecta, which is characteristic to the physical process at hand.In addition, the location of the outburst site, defined by its longitude and latitude, is determined randomly from a set of distributions.Depending on the process, each particle species may have different parameters for sampling.The current version of the code supports thermal evaporation and sputtering.
Thermal particle velocities follow a Maxwell distribution given by the χ 2 -distribution with three degrees of freedom and scale parameter v sc : with However, in this work we focus on mass loss by the mechanism of sputtering.This phenomenon describes the ejection of particles from a material or surface when it gets bombarded by high-energy particles.Multiple models have been used in sputtering, one of which was described by Smyth and Combi (1988).It combines a Boltzmann distribution core (thermal) with a high-energy tail: where D(α, v M /v b ) acts as a dimensionless normalization constant: Figure 3. Geometry and coordinate systems of the underlying simulations.The x-axis is defined by the co-alignment of exoplanet and star at t = 0.For t > 0 the exoplanet moves along its orbit characterized by the true anomaly ϕ p .The exomoon true anomaly is denoted as ϕ s ≡ ϕ ☽ .The y-axis is inside the orbital plane and perpendicular to x.The third axis z is parallel to the specific angular momentum vector of the orbital plane.The exomoon is described by the same reference frame, only shifted by the position of the exoplanet.In order to track the evolution of the exomoon, a co-rotating reference frame is of advantage.
Both Equations 8 and 9 depend on three parameters.The parameter α regulates the width of the speed distribution.
It takes a value of 3 for a classical sputtering distribution, that is, a complete collisional cascade process, and a value of 7/3 for a Thomas-Fermi modified-sputtering flux distribution, that is, the limit of a single elastic collisional ejection process (also see Smyth and Combi (1997)).v b is non-linearly related to the most probable of the flux speed distribution and depends on the exobase temperature T X (Smyth & Combi, 1988, 1997).v M describes the maximum speed that a stationary atom can acquire in an elastic collision with a plasma torus ion.Using these parameters, the largest possible value for the speed is given by v Figure 4 shows a comparison between the distributions Equations 6 and 8.
The distribution Equation 8 is not trivial and doesn't make the often used method of inverse transform the preferred choice of sampling.To generate random samples, a rejection approach is used.Adopting the Nelder-Mead method, we find the maximum max{f sp } of the sputtering distribution numerically using the Scipy package for Python.We now generate a sample value v from the target distribution V sp by instead sampling from an uniform distribution U on the support We accept the sample v if the subsequently also uniform sampled probability density value Once the speed has been sampled according to Equation 8, the velocity vector is constructed from uniform distributions in spherical coordinate angles, ranging from 0 to 2π and 0 to π/2, respectively, in order to provide a hemispheric coverage.It then gets rotated according to which latitude and longitude on the source object the particle is emitted from and anchored to that position.For our simulations we choose to have an isotropic ejection, that is, both longitude and latitude of the outburst site are sampled from a uniform distribution, as well.
After the particle is created and added to the simulation, SERPENS adds an individual signature of each particle to a hash dictionary that later allows matching of particles to their properties.

Loss Function
Reactions may happen with certain rates that can change the composition of a particle cloud, as well as the velocities of individual particles.SERPENS doesn't directly simulate reactions between test particles, but reactions defined via their reaction rate and the chemical products.At each advance, the weight of a superparticle is calculated by a decay exponential to describe depopulation: with lifetime τ i = 1/ν i for the ith reaction given the reaction rate ν i .For some constituents/systems one may have not any information about the local chemical network, in which case a single τ value has to be used as an estimate for the lifetime.
Reaction rates are currently an user input and not calculated in SERPENS directly.In addition to chemical losses all particles that leave the system get removed from the simulation after time propagation.The distance at which this should happen is a parameter of the simulation.
Collisions are detected by REBOUND after the integration (Rein & Liu, 2012).Particles that collide with a source of gravity are merged with that object and therefore also eliminated from the simulation.
In theory, one could add mass to each particle (or dust) and simulate the transfer of mass and the destruction/ growth of celestial bodies with this merging approach.However, this is not part of this work, but could be interesting for further studies about large volcanic eruptions on exomoons.

Advancing the Simulation
After particles have been added and their weights adjusted according to the time they have spent in the simulation, SERPENS passes the instance to REBOUND for integration in time.
The integrator chosen during initialization is an implementation of the Wisdom-Holman symplectic integrator (WHFast).In a nutshell, WHFast splits the problem into a Keplerian (drift) and an interaction part (kick).After doing a half-timestep Keplerian drift, WHFast applies a full timestep interaction velocity modification and completes the timestep with the other half of the Keplerian drift.A more detailed description of the WHFast algorithm can be found in the works of Rein and Tamayo (2015) and Wisdom and Holman (1991).We employ REBOUNDX to apply radiation forces to particles in the simulation, incorporating both radiation pressure and Poynting-Robertson drag (Tamayo et al., 2019).To speed up the process, SERPENS splits the simulation instance into smaller pieces that can individually and simultaneously be integrated in time using multiprocessing.As soon as the calculation has finished, these instances get recombined.Note that the splitting approach is only applicable if the particles do not directly interact with each other.
After the integration, we update the hash dictionary, make a copy of it and append it to a sub-dictionary that holds the information of each advance.Note that an advance refers to an integration in time until new particles are added to the simulation and not a single integration time-step.In our simulations an advance is typically set to 1/120 of an orbit.
Finally, we save the current state of the simulation in an archive file that holds the complete information of the system at that point in time.This allows us to later analyze the system's evolution by comparing its states at multiple points in time.

Parameter Settings
In this work, we only consider mass loss to space by atmospheric sputtering.The mass supply to the upper atmosphere is supplied by tidal and thermal heating of the surface as estimated by the dishoom model.For our smallest satellite mass (e.g., Enceladus) the source rate is equivalent to the escape rate Ṁ0 ≈ Ṁesc .For each system we consider two different particle lifetime scenarios, τ cloud and τ torus , which correspond to a photoionization-limited case and a toroidal case of 3 hr lifetime, respectively.The parameters of the speed distribution (Equation 8) of sputtered particles is essential for the morphology of the particle cloud.We assume a value of α = 7/3 and calculate rough estimates for v b and v M : 1. v M primarily determines the maximum speed and depends on the relative velocity of the surrounding plasma and the masses of colliding particles.In order to estimate an upper bound, we assume that there are no heavier plasma constituents than sodium ions, s.t.v M ≈ v rel for the sodium neutral velocity distribution.In the case of a tidally locked planet, the plasma corotation velocity is simply the orbital velocity of the exoplanet and the relative velocity v rel therefore the difference between the orbital velocities of exoplanet and exomoon.2. v b can be calculated from the exobase temperature of the exomoon.At Io this temperature ranges from 220 to 2800 K (Wong & Smyth, 2000).However, as this temperature is extremely difficult to determine even for objects in the Solar System, we resort to using the estimated surface temperature of the host exoplanet.Then, we calculate v b numerically by solving the non-linear equation Smyth and Combi (1988)) under the approximation that T X ≈ T surf .As mentioned before, in the case of Jupiter and Io the plasma corotates with a velocity much higher than Io's Keplerian velocity, whereas the tidally locked gas giants examined in this work show an opposite behavior.A hypothetical exomoon exceeds the speed of a corotating plasma.Hence, some caveats in using the velocity distribution Equation 8established for sputtering at Io may need to be examined.Cloud geometries could be reversed/mirrored to our presented simulation results as our simulations only take the magnitude of the relative plasma velocity into account and not its sign.3. Drawing from the analysis of XUV atmospheric evaporation and stellar wind interactions at HD-189733 b and HD-209458 b by V. Bourrier and Lecavelier des Etangs (2013), we utilize linear interpolation to estimate minimum values of β Na , which we use for our simulations.We interpolate between the values of β HD189733b = 2 and β HD209458b = 4 (V.Bourrier & Lecavelier des Etangs, 2013), based on the corresponding irradiance at the exoplanets.4. The semi-major axes of the exomoons are calculated as a ☽ ≲ 0.49 a Hill according to Domingos et al. (2006) following the stability analysis of Cassidy et al. (2009), Kisare and Fabrycky (2023), and Oza et al. (2019).5. Sodium lifetimes are assumed to be limited by photoionization and are adopted from Oza et al. (2019).While particles are inside the planet's shadow, their lifetimes are increased to 3 hr, reminiscent of sodium lifetimes near Jupiter's moon Io.

Densities
The main interest of the simulation is the tracking of particle and column densities over time.Extracting physical data from the distribution of point-like features, that is, the superparticles, is not trivial.The reason lies behind the fact that we, in general, do not have any information about the background particle environment.This makes conventional binning and kernel density estimations difficult, since there is no clear indication for the optimal bin size or convolution bandwidth.To overcome this obstacle, SERPENS borrows a method often used in cosmological problems (e.g., van de Weygaert et al. ( 2013)).This method is known as a Delaunay tessellation field estimation (DTFE) and was developed by Schaap and van de Weygaert (2000) based on ideas from the velocity interpolation scheme introduced by Bernardeau and van de Weygaert (1996).The Delaunay tessellation forms an ideal adaptive multi-dimensional interpolation grid from the particles populating the simulation.See Figure 5 for a schematic comparison between the DTFE and a grid-based model.The Appendix A3 provides a mathematical introduction to the concepts behind the DTFE.
Each (super-)particle's mass takes the role of a scalable weight in the density calculation Equation A3.If we track all particles, even if they represent extremely small density regions, we can later refer to the distribution at different mass-loss rates by re-evaluating the particle weights and removing all particles below a certain threshold (e.g., 1 cm 3 ) afterward.Equivalently, if a process reduces a particle's lifetime but does not change its trajectory noticeably, the morphology will remain the same as only the weight and therefore scaling will be affected.

Transit Depth
For a solid body or optically thick material, the transit depth is independent of the stellar wavelength, and is opaque to continuum radiation.Thus, from geometric considerations alone, one can derive the transit depth δ tra as first approximation to be In addition to transit, a change in the lightcurve can be seen when the exoplanet moves behind the star (occultation), where exoplanet radiation and stellar radiation reflected from the exoplanet's surface will no longer add to the flux observed.We assume that the planet is of uniform brightness, and that the observing wavelength is long enough for thermal emission to dominate.Subsequently, approximating the planet and star as blackbodies, the occultation depth is evaluated as where T p and T ⋆ denote the temperatures of planet and star (Winn, 2010).In terms of observation, the planet's temperature is wavelength-dependent, T p = T p (λ), such that one may define a brightness temperature T b as the equivalent blackbody temperature that would lead to the observed value of δ occ (λ).
When probing the exosphere of a planet, however, material is usually found to be optically thin, allowing for the identification of the absorbing gas by their imprint on the residual stellar spectrum.In that case the transit depth is dependent on the wavelength as the amount of absorption depends on the quantum transition.In our case, we assume that the exomoon's vented sodium is indeed in such an optically thin regime.The wavelength-dependent optical depth τ may then be calculated as where 〈N〉 is the average line-of-sight column density and σ abs,λ the wavelength dependent absorption crosssection.The most prominent absorption line of sodium is the D 2 -line of the Fraunhofer D-doublet.Using the Python module prometheus (Gebek & Oza, 2020), which employs equations from Draine (2011), we calculate the absorption cross-section at a bin-width of 0.2 Å typical to high-resolution observations.Results for temperatures of 1000 and 2000 K show only minor differences, such that a temperature dependency can be neglected.The resulting cross-section used for transit-depth calculations is given by σ abs = 9.84 ⋅ 10 13 cm 2 .We calculate the transit depth governed by the D 2line with δ Na = exp( τ λ NaD2 ). ( 17)

Model Validation
To verify the output of our model, we use Jupiter's moon Io and its Na exosphere.Through gravitational tidal forces from Jupiter and the other Galilean moons, Io is the most volcanic active body in the Solar System.The neutral Na exosphere of Jupiter is emergent from Io's volcanism interacting with the plasma from Jupiter's magnetosphere.Remote and in-situ measurements allow for a validation of our simulation results obtained with the DTFE.
State-of-the-art models that are specifically designed for the simulation of Jupiter's particle environment, including Io and Europa, are more sophisticated than our approach in this work.As our goal is the simulation of exoplanet systems, many facets of modern Galilean satellite models are currently simply unknowns for exoplanet modeling.They depend on parameters that we presently cannot quantify to a satisfying degree around planets & moons outside the Solar System.Development of SERPENS may include improvements and more parameters to adjust in the future.For the time being, we analyze Io in a more crude fashion.Our aim is to confirm the general morphology and density ranges of sodium in the Io torus and apply our methods to the eight candidate systems at hand. Figure 6 depicts the sputtered particle distribution making the characteristic "banana cloud" from low velocity sputtered particles visible (at few km/s).We additionally plot the Delaunay tessellation showing the position of all particles in the simulation.Burger et al. (2001) measured the column density profile of Io's sodium corona at 10 mutual eclipses between the Galilean satellites.Their measurements in column density ranged from 1010 to 10 12 cm 2 .Similar values have been reported in the past (Brown, 1974;Brown & Chaffee, 1974) and by Schneider et al. (1991).We are able to reproduce these values using the DTFE and a sodium mass-loss rate of 10 kg/s.

Results of Toroidal Simulations
Equipped with a novel DTFE algorithm, SERPENS generates gas density distributions driven by the mass-loss rate of an evaporating exomoon.The emergent toroidal exosphere adheres to a profile with scale height (Johnson & Huggins, 2006), where v ej /v orb is the ratio between atom ejection speed and orbital velocity of the gas or venting exomoon.In order for the toroidal exosphere to emerge, gas needs to vent beyond the gravitational influence of the satellite, that is, the ejection velocity needs to surpass the escape velocity of the satellite.A minimal scale-height can therefore be calculated as v ej ≈ v esc (Table 5, 6th column).Increasing the temperature of the system generally leads to higher ejection speeds which in return increase the scale height of an emerging torus.In terms of our simulation, this is captured by the dependency of the sputtering parameter v b on temperature.In contrast to thermal evaporation, the sputtering model includes the ejection of high-energy particles increasing the effective size of the toroidal exosphere.Exomoon phase-curves capture the temporal evolution of particle distributions using the DTFE at each simulation advance once a quasiperiodic steady-state has been established.These phase-curves carry observational utility in that the exomoon gas density evolution may be compounded over multiple exoplanet transits.For Ṁ ϕ ☽ ) ∼ Ṁ0 constant, the resulting phase-curves 10 and 11 are source-limited by an average Ṁ, such that our density simulations can be described analytically via a simple scaling law, where various values of Ṁ are chosen based on the expected thermal or tidal outgassing rates computed by dishoom.In more intricate cases, where radiation pressure and shadow variations are substantial, this scaling reduces to the offset of the phase-curve.Such cases exhibit radiative forces more than twice as strong as gravity from the star, or they involve localized clouds that remain entirely concealed during exomoon occultation.
For a simplified scenario where these effects are negligible, the particle density distribution is expected to follow a cosine-like behavior dependent on the orbital phase angle ϕ ☽ of the satellite, motivated by the geometry and periodicity of illumination (cf.Marconi, 2007;Oza et al., 2018): where Ṁ0 is a constant (or averaged) mass-loss rate, m the mass of a volatile, τ its survival timescale, V tor the volume of a torus and ϕ ☽ the orbital phase angle of the satellite.A phase angle of ϕ ☽ = 90°represents the trailing phase, symbolizing dawn, while a phase angle of ϕ ☽ = 270°signifies the leading phase during elongation, associated with dusk.However, our simulations show that the frequency of the particle density oscillation is not bound to one orbital period, but may be higher depending on the system in question (e.g.HD-209458 b and WASP-69 b), that is, ω tor ≠ ω ☽ .This feature is not unexpected and we find that in the absence of radiation Note.(M, R) x denotes the type of exomoon considered in the cells following to the right.〈⋅〉 ϕ is the exomoon orbit average for a quantity.We calculate these averages for line-of-sight column densities N [cm 2 ] and particle densities n [cm 3 ].For each orbit averaged value, we show the ranges to the maximum and minimum with + .ΔN ☽,occ are the orders of magnitude the column densities drop off by when the exomoon is occulted by the exoplanet.We further table the analytical expression for the minimal torus scale-height H = a ☽ ⋅ v esc / v orb and the wavelength Doppler shift at the sodium D 2 line Δλ NaD2 .Table A1 provides results for simulations without the presence of radiative effects.
pressure, the frequency of the particle density oscillation is correlated with the maximum sputtering velocity of the system (cf.Table 4 and Figures 7-9).Analytical fits of the form A cos ω(ϕ ☽ ϕ ☽,0 ) + 〈log n〉 ϕ incorporate a periodic signal, but not necessarily orbit-periodic, superimposed on an average steady-state background density that is attained after multiple orbital periods.We apply such fits to the density evolution in the case of an exo-Io venting gas into orbit (see Figures 7-9).A large value of A implies significant dusk-dawn variations in the particle density distribution due to stellar gravity decreasing or amplifying acceleration toward the exomoon and exoplanet.
Taking into account radiation pressure, results are no longer governed by a (co-)sinusoidal signal, but by dominant density peaks near noon.This reflects in poor analytical fits of the seen in Figures 10 and 11.
We find that the LOS column density follows a trend governed by the combination of an orbit periodic and orbit double-periodic signal where v esc /v orb is the ratio of escape-to orbit velocity, related to the torus scale height H = a ☽ v esc / v orb (Johnson & Huggins, 2006).C 1 and C 2 are different amplitudes that can be associated with a dawn-dusk asymmetry (C 2 tends to zero for perfect symmetry).Exo-Io column density fits are constructed as accordingly, fixing the periodicity of the signal, but allowing for dawn-dusk asymmetries that are superimposed on a steady-state background.Radiation pressure amplifies this asymmetry by pushing particles toward the planet (during transit this is toward our line of sight) while they move laterally along the orbit.Particle clouds are accelerated close to the exoplanet and expand at the moment when the exoplanet no longer blocks the line of sight to the exomoon.This results in a dense, fast moving cloud (peak) that quickly disperses afterward (dip).The resulting strong noon-dusk variation is therefore not only correlated with an expected uneven thermal evaporation due to surface temperature inhomogeneity, but also with the magnitude of radiation pressure.Sometimes, as in the case of HAT-P-1 b, an Earth-sized exomoon's gravitational pull may be strong enough to withstand this effect.If the fast moving clouds near dusk are visible during a transit event, that is, if the exomoon passes through noon and dusk phases during exoplanet transit, significant radial motion toward the observer will present itself in high-resolution spectra by inducing Doppler shifts.Our simulations are able to explain Doppler shifts of up to Δλ NaD2 ∼ 0.13-0.50Å depending on the system (see last column of Table 5), which correspond to radial motion of ∼10-30 km/s.From the extracted line-of-sight column densities we are able to calculate expected sodium spectral transit depths according to Equation 17for both cloud and torus scenarios (see Table 6).In addition to spectral transit depths associated to a toroidal sodium exosphere, we list geometric transit-and occultation depths of the exoplanet following Equations 14 and 15 in Table 6 (also see the transit curve 17 for WASP-69 b).
Visual representations for the exo-Io simulations are shown as a collection of plots for each system, including Figures 12 and 13 for HD-189733 b and WASP-96 b, respectively, and allow for a clearer picture during interpretation.These plots are arranged in a 4 × 4 grid and show distributions at exomoon phases ϕ ☽ = 0°, 90°, 180°, 270°.Equivalent figures for HAT-P-1 b, HD-209458 b, WASP-49 b, WASP-17 b, WASP-69 b, and XO-2N b can be found in Appendix A2.
Upon studying the exomoon phase curves (Figures 10 and 11) for each exoplanet host system, we find that there is no single parameter that generalizes the simulated density evolution.For instance, some systems display the highest column densities for exomoons with an Earth-like mass and radius, whereas others surprisingly have density maxima for an Enceladus sized body (cf.Tables 5 and 6).We want to quickly mention the most important results for each system.

WASP-69 b I:
An Earth-like exomoon produces a sodium cloud that reaches maximum density approaching noon.Here we have radiation pressure compressing an otherwise elongated particle distribution.WASP-69 b I's sputtering distribution has the lowest maximum velocity from the considered systems which is unable to act against this compression.The lower mass exomoons have a smaller gravitational acceleration to hold onto particles.The initial particle cloud is less confined, such that the noon effect is not as prominent.Here, highest density is achieved at leading phases (dusk), where a denser particle tail is created through the gravitational pull of the star and the Poynting-Roberston effect (note that WASP-69 b is the lowest mass and coldest exoplanet system considered).

WASP-17 b I:
WASP-17 b I has the smallest semi-major axis among the exomoons studied.Many particles fall into the exoplanet and radiative processes at noon amplify this effect.During trailing phases (dawn) particle densities are highest, as particles accumulate during the shadow phases (increased lifetime and no radiation pressure).This way the particles are more free to escape the system than afterward.The semi-major axis of the exomoon resides at the boundary of the Alfvénic sphere of its host (see Table 3).Outflowing particles may therefore experience additional plasma effects that we do not consider here.

HD-189733 b II:
From noon until midnight a half-torus reminiscent of Io's "banana cloud" is formed.After midnight radiative forces act against further fueling of the banana cloud which disrupts the build-up of a complete torus.However, particles remain in line-of-sight, which results in only a minor effect on the LOS column and thus only a small dawn-dusk asymmetry.A visualization can be found in Figure 12.With the exception of an Earthlike exomoon, WASP-69 b I's baseline behavior is similar.upon.This evolution is characterized by a peak density in both particle and LOS density near noon, followed by a decline.
In addition to the detailed orbit-phase dependent density calculations captured in the exomoon phase-curves, we calculate orbit-averaged values for column-and particle densities for all three exomoon mass ranges.Our main results for the photoionization limited case, captured in Table 5, demonstrate that large variations in column density are present, often over multiple orders of magnitude.Furthermore, within the line of sight, column densities may drop multiple magnitudes ΔN ☽,occ , due to occultation of the exomoon by the planet, especially in the sodium cloud scenarios.The phases-curves in Figure 16 show results that incorporate this effect.Care should be taken when interpreting values very close to exomoon occultation ingress and egress.Here, large variances in particle densities are present.Sampling noise will propagate immediately into fluctuations in the reconstructed density field, as the discrete point sample itself will be the source for the density value estimates.
Photoionization lifetimes usually result in very strictly confined particle clouds around the exomoon, effectively an extended exosphere.Corresponding phase-curves show strong variations in contrast to simulations, where lifetimes are comparable to that of sodium within the Jupiter-Io plasma torus.Over these extended timescales, tori can form around exoplanets (cf.Figures 12 and 13, and Appendix A2).Nevertheless, we notice sudden jumps in density at dawn and after noon in multiple systems.

Formation of an Exomoon Torus
We aim to explore the mechanisms by which toroidal formations similar to Jupiter's Io-torus can create fluctuations in the distribution of particles and column densities.In scenarios that result in localized clouds, where the Note.We also include R for the photoionization-limited sodium cloud.(M,R) x denotes the type of exomoon considered in the cells following to the right.δ Na,torus is the orbit-averaged sodium transit depth for the torus scenario (sodium lifetime of 3 hr).
abundances of neutral sodium are limited by rapid photoionization, it is often challenging for exomoons to maintain such a torus.However, if we assume that the particle lifetime is extended to approximately match that of sodium in Io's plasma torus, a toroidal structure emerges.By analyzing a small selection of exoplanets that have been investigated, we can identify several criteria that contribute to the formation of a fragmented neutral torus.
1. Increased lifetime, τ i : Our simulations show that neutral lifetimes of just a few minutes are not sufficient to build up an observable torus of neutrals.Other mechanisms need to be present to prolong their lifetimes.2. Low-mass exomoon or high mass-loss rates, Ṁ: Due to tidal locking, a plasma that is corotating with the exoplanet's magnetic field does not reach the high velocities that we know from Jupiter and Io.The resulting sputtering distribution is centered at smaller velocities and more particles are recaptured by the exomoon or fall into the exoplanet.Therefore, sputtering needs to happen on a body that is less massive or additional impacts of highly energetic particles or micrometeoroids are needed to create enough ejecta in the high-energy tail of the distribution.For each system we have run three simulations that consider different satellite masses and radii: Earth-sized (⊕), exo-Io (Io), and Enceladus-sized (Enc).Vertical dotted lines in the LOS plots mark the occultation of the satellite by the exoplanet.The present panels do not include occultation effects that strongly counteract peaks in the LOS column.During transit of the exoplanet, the exomoon will surpass a limited amount of phase-angles illustrated by a horizontal bar on the plots.
3. Sufficiently large magnetic field of the host planet, B p,crit , shielding the exomoon from an extreme stellar plasma environment."Sufficiently large" depends on the model for the magnetic dipole moment, exoplanet thermal mass-loss, and stellar wind parameters.It is therefore different for each system and no universal threshold can be provided in this work.4. Sustained mass-loss at systems with β > 1: By definition, values of β > 1 will strip particles from the exoplanet-exomoon system, since radiation pressure dominates over gravity.If there is no steady supply of particles, a torus will quickly be depopulated.Figure 16.LOS column densities for the considered systems with the inclusion occultation effects.For each system we have modeled exomoons of different masses and radii according to Table 1, representing a timely evolution of a disintegrating body starting as an Earth-mass and -radius exomoon.We ran two sets of simulations with photoionization-limited lifetimes and lifetimes of 3-hr, similar to the lifetime of Na at Io. Included are indicators for the range of phases an exomoon would pass through during an exoplanet transit.Vertical dotted lines mark the phase-angles at which the line-of-sight to the exomoon gets blocked by the planet.
This initial attempt to model toroidal atmospheres, as theorized by Johnson and Huggins (2006), offers a detailed preliminary analysis, intended to serve as a stepping stone for future research.To gain a more precise understanding of the above mentioned criteria, a broader range of simulations using our open-source model is necessary.Constraining the parameter space for exotori formation would greatly improve our understanding of the formation and evolution of these systems.Such a quantitative analysis would also provide valuable insights into the underlying mechanisms driving the formation and stability of these fragmented neutral tori.For that reason, we also encourage theoretical work in studying the different physical and plasma processes that are well known in our solar system (Lopes et al., 2023).

Gravitational and Geological Considerations
Our presented results invoke symmetries to simplify calculations.Keeping track of resulting simplifications is important to judge the simulations' validity.
The geometry in our exoplanet-exomoon systems is set in such a way that the exomoon lies directly in the orbital plane with the physical justification of strong stellar tides forcing the exomoon into the orbital plane of the exoplanet.In addition, at each integration step we provide a gravitational "kick" to tidally lock the exomoon into a circular orbit around its host.In this way tidal energy is dissipated into the heating of the satellite itself, rather than dissipating into orbital decay.Without this approximation of tidal resonance, many systems move toward a hyperbolic or parabolic orbit for the exomoon.The fate of a non-resonant, rogue exomoon is still an open area of dynamical research.Rather than falling into the exoplanet, it is believed that a rogue exomoon would first orbit the star and eventually impact or be captured by the exoplanet as a revived satellite (Hansen, 2022).Observations of ejected exomoons or exocomets may be present at star systems exhibiting transient signatures of continuum and narrow-band sodium light (Chakraborty et al., 2004;Lecavelier Des Etangs et al., 1996;Mamajek et al., 2012;Martinez et al., 2019).Geophysical models are needed to study the precise tidal heating of the Earth-Io-Enceladus sized exomoons here.However, our understanding of the processes is limited, even in the case of Io, where Na/Si and Na/Mg ratios remain uncertain until sample return missions provide insight (Ogliore et al., 2023).Rheology of the interiors, especially tidally heated melting, will influence the tidal Q of the satellite.Consequently, the tidal heating rate Ė☽ ∝ e 2 ☽ / Q, driven by a forced eccentricity e ☽ , contributes to the surface temperature according to ΔT ∝ Ė1/ 4 ☽ (Oza et al., 2019).This directly influences the height of the exobase and the concurrent gravitational binding of particles in the exosphere.Not only can sputtered particles escape more easily, but thermal evaporation also gains in importance.Together with an increase in temperature comes a more intense mass-loss rate, which raises concerns about exomoon evaporation and associated disintegration timescales.The timescale for evaporative destruction of exo-Ios was found to be above the stellar age for most systems considered here (Oza et al., 2019), with the exception of WASP-69 b.In the event of a catastrophic disintegration, an exoring composed of debris, dust, gas, and plasma may form.Under extreme photodesorption conditions, alkaline signatures can be sustained or even enhanced over an extended period.That said, it is important to note that our current simulations solely incorporate the influence of temperature through the velocity distribution parameter v b .As temperature increases, the average velocity of particles in the simulation also increases, resulting in more uniformly distributed particle clouds across larger scale heights.Refer to Appendix A5 for an exploration of significantly higher surface and exobase temperatures, leading to increased sputtering velocity parameters.
For close-in orbiting transit planets, the Keplerian velocity relative to the stellar winds is non-neglectable.This aberration effect results in the magnetic tail not being directed radially away from the star.Therefore, magnetic tail and shadow are not aligned and thus plasma may detach from a neutral cloud and ENAs.Tracking both fast neutrals and the plasma flow may show a comet-like distribution with two tails.Currently, the physics of chargeexchange and electron-impact ionization at exoplanetary ring systems is too poorly understood to incorporate them into our simulations.Nevertheless, SERPENS is already set up with the support of chemical networks, including electron-density scalable reaction rates and velocity vector kicks.We look forward to the exploration of comet-like signatures from exoplanets due to this effect in the future.Mass-loss rates of pure sodium at ∼10 5 kg/s (Table 4) are magnitudes larger than at Io and would result in disintegration on timescales of a few hundred million to a few billion years, as studied in Oza et al. (2019).A larger progenitor is one option for exomoon survival, however, there are other important considerations to take into account: (a) the exomoon might be captured, and so the disintegration is ongoing; (b) the larger end of loss rates in Table 4, may evolve into more of an exoring system of "moonlets"; (c) a component of the sodium may fall back onto the exomoon; (d) there is a lack of established geophysical models capable of accounting for the sodiumto-rock ratio, an essential factor that remains unknown even in the case of Io.Currently, mass-loss rates are constant throughout the orbital motion, which results in an accumulation of particles while the exomoon is in the shadow of the exoplanet.A simulation with phase-dependent mass-loss rates Ṁ = Ṁ ϕ ☽ ) is needed for the validation of this effect as cooling in the shadow-phases may counteract the accumulation (Schmidt et al., 2023).
If an exoplanet has high obliquity, it typically results in an inclination of its satellite's orbital plane.Consequently, the satellite's orbital motion around the planet would cause the circumplanetary toroidal structure to also be inclined relative to the planet's orbital plane, potentially leading to diffuse scattering from an inclined exoring that could modulate the transit signal, as explored by for example, Sucerquia et al. (2020) and Zuluaga et al. (2022).Additionally, a gaseous torus closely encircling the exoplanet could potentially cause an overestimation of the planet's radius, especially in the case of gas giants (Zuluaga et al., 2015).Such instances could further support the hypothesis attributing high-altitude atmospheric features to exorings or exomoons.However, the stability of an inclined toroidal structure remains uncertain due to the extreme tidal forcing present in the studied systems as mentioned above.

Ultraviolet, Optical, and Infrared Observability of Exomoons
Our simulation results here are intended to be tested for follow-up with highresolution spectroscopy (Cosentino et al., 2012;Pepe et al., 2000) where an absolute measurement of the Doppler shift of the vented exomoon gas can be obtained.This is brightest and most analogous to Io for the alkali metals Na I and K I in the optical, however byproduct ions in the UV, excited auroral lines (Kislyakova et al., 2018), as well as molecules in the infrared are currently accessible in low-resolution with HST and JWST.
Examples of observational anomalies which may indicate the presence of a satellite in exoplanet transmission spectroscopy are night-to-night variability.Significant variations of metallic columns in consecutive exoplanet transits are one indication of an evaporating (or disintegrating) satellite.Spectrographs with resolving powers of R ≳ 10 5 including ESPRESSO and UVES on VLT (Dekker et al., 2000;Pepe et al., 2014), HARPS & HARPS-N (Cosentino et al., 2012;Mayor et al., 2003), and PEPSI on LBT (Strassmeier et al., 2015), are sensitive to roughly ∼10 9 -10 11 Na/cm 2 at absorption depths of ≲1% (Draine, 2011) depending on the collecting power.Next decade high-resolution instruments onboard the Extremely Large Telescope (ELT) and 30 m Telescope (TMT) will be able to constrain extremely tenuous quantities of gas to unprecedented thresholds.
Important for LOS observations is the fact that during exomoon occultation by the exoplanet, column densities often drop my multiple orders of magnitude, pointing toward the localization of the sodium cloud toward the satellite.This is more reminiscent of an extended exomoon exosphere than to a neutral torus around the planet.The number of magnitudes by which the column density drops during occultation is given by ΔN ☽,occ found in Table 5.With an instrument sensitivity of ≲10 9 Na/cm 2 multiple exoplanet systems are sofar observed to be consistent with exomoon simulations presented here.WASP-69 b has the most significant column that is always above this if we disregard the occultation of the exomoon.Many other systems reach the threshold only during certain phases that resemble a peak in particle/column density for these systems.A summary of orbit-averaged column densities for both cloud and torus scenario can be found in Table 7.We compare our values with those that have been analytically calculated by Oza et al. (2019) using the principal mass-loss source based on the dishoom outgassing code.Transmission spectra data of WASP-69 b, as presented by N. Casasayas-Barris et al. ( 2017), exhibit noisy light curves that cannot be accounted for by their combined model fits.The characteristics of the data strongly indicate that sodium might not be associated with the exoplanet itself.Our results suggest variations in the equivalent width of sodium due to the orbit of an outgassing satellite or exoring sustaining alkaline signatures (cf.Section 5.2).A translation of column density to spectroscopic transit depth according to Equation 17 can be found in Figure 17 for the case of WASP-69 b I. Since we do not know the initial phase of the exomoon at first transit, multiple passes are necessary to make an assessment of the exomoon existence plausibility.A sodium column density evolution curve needs to be reconstructed from observations to be compared to our model and predictions.
Our simulations of WASP-17 b I show that a significant amount of material can be deposited in the upper atmosphere of the hot Jupiter due to the strong radiation pressure (β > 5) and close-in orbit at only a ☽ ∼ 1.28 R p .In contrast to our findings for WASP-69 b I, the sodium cloud is strongly localized (compare Figures A3 and A5).This phenomenon could offer a plausible explanation for the observed absence of a sodium signal during Hubble Space Telescope (HST) observations (Alderson et al., 2022), especially if an exomoon was concealed behind the exoplanet during those observations.Furthermore, this localization pattern lends support to the hypothesis that a satellite may account for the presence of silicates on WASP-17 b, as opposed to an exoring or other exogenic sources.The employed mass-loss rate estimate at this exomoon is approximately two orders of magnitude smaller than that at WASP-69 b I, potentially extending the longevity of the putative exomoon, provided its orbital dynamics remain stable.

Follow-Up Investigations
This paper restricts itself to the analysis of sodium for eight exoplanets.We are currently investigating the effects and distributions of K and SO 2 for exoplanet systems with alkali detections using our SERPENS framework.Additionally, great interest exists in simulating the behavior of metals such as Mg and Fe, which have been observed at high altitudes in the exospheres of ultra-hot Jupiters WASP-76 b and WASP-121 b, within a nonhydrostatic environment (Ehrenreich et al., 2020;Hoeijmakers et al., 2023;Sing et al., 2019).Photoionization lifetimes are especially short at ultrahot systems, an orbiting exomoon may help provide a sufficiently high neutral atom rate, at high altitudes.Moreover, we aim to investigate the possibility of disintegrating metal-rich exomoons or debris systems surrounding these exoplanets, as they may offer a plausible explanation for the observed metallic signatures at high altitudes, where an exomoon can orbit.

Conclusion
We have combined an atmospheric evolution and escape model with a three-dimensional test particle Monte Carlo simulation code plus time, to predict volcanic gas column densities (e.g., Na/K) in neutral exospheres and tori created by putative exomoons.We have applied this formalism to a set of eight candidate star systems to predict an exomoon phase curve for the first time.An exomoon phase curve tracks the density evolution of a volatile which may be seen near transit, occultation, or in orbit based on the expected lightcurves of a venting exomoon.Of course, this setup of a volcanic exo-Io system, can also be applied to outgassing Solar System satellites, as well as a star-planet "super-Io" system, considering for instance a closein exoplanet like 55 Cnc e (Demory et al., 2023;Morris et al., 2021;Valdés et al., 2022) forming a torus based on similar physics studied here.Modifying mass and radius of the evaporating object to closer resemble an Earth-or Enceladus-sized object can significantly alternate the distribution of particles along the orbit.Consequently, tori associated with rapidly evaporating or disintegrating satellites experience an equally accelerated evolution in the presence of strong tidal forces and radiative processes.Conversely, the density evolution trend can provide insights into the size and mass of an active satellite.However, further analysis is necessary to establish the precise correlation.Our approach allows for scaling with mass-loss rates as long as the mechanism of evaporation remains consistent.Nonetheless, our simulations demonstrate a wide variety of possible distributions, influenced by the specific environment of the exomoon.It is therefore difficult to predict the evolution of a particular system on the basis of another.Nevertheless, we have identified two primary types of baseline particle distribution evolution, roughly separated by a radiation pressure parameter of β ∼ 2, upon which further variations are built.Below this threshold, systems exhibit an anticipated periodic and double-periodic evolution in particle and column density, respectively.The density of particles evaporated from exomoons in an environment of more significant β becomes more pronounced near noon, followed by a decrease around ϕ ☽ ∼ +60°during dusk (elongation).To further constrain this boundary, a larger set of exoplanets must be studied.It is important to note that the β values employed for the exoplanet systems in this study are derived from the analysis of Lyman-α spectra and are thus not yet finely tuned for sodium gas (see Section 3.4).
The exoplanet WASP-69 b stands out as the most promising target among the considered systems due to its large column density ≥10 9 cm 2 .This high density exists in an environment with comparatively small radiation pressure and reasonably sized Alfvén sphere, featuring an estimated magnetic field strength of B p ∼ 15 G (Narang et al., 2023;Yadav & Thorngren, 2017).These factors contribute to the potential stability of a torus.We look forward to JWST cycle-1 GTO observations and encourage monitoring for any potential exomoon signatures.In our simulations, HD-189733 b emerges as the second most favorable system for exomoon exploration with instruments capable of detecting sodium columns at ≳10 7 cm 2 .Analyzing the system with a plasma code would allow to investigate the possibility of particles being carried away to distances of approximately ∼16 R p , thus establishing a connection between our approach and the results presented in (Ben-

A2. Simulation Visualizations
The following Figures A1-A6 provide visualizations of our results for exomoons that have both mass and radius alike Jupiter's moon Io.More comments on the individual systems can be found in the main body of this work, in particular Section 4. Each figure incorporates 16 subfigures that exhibit the exomoon at orbit phases 0°90°, 180°, 270°, both for sodium lifetimes of 3 hr (characteristic to Io) and photoionization limited lifetimes.Top-down and line-of-sight perspectives are provided with particle and column densities, respectively.Note.(M,R) x denotes the type of exomoon considered in the cells following to the right.〈⋅〉 ϕ is the orbit averaged value.We calculate these averages for line-of-sight column densities N [cm 2 ] and particle densities n [cm 3 ].ΔN shadow are the orders of magnitude the column densities drop off by when the exomoon is occulted by the exoplanet.We further table the analytical expression for the torus scale-height H = a ☽ ⋅ v esc / v orb and the wavelength Doppler shift at the sodium D2 line Δλ NaD2 .which is equal to the average density at the vertices of the simplex S n times its volume.Therefore, in which S ij is one of the simplices of which nucleus i is a vertex.As the complete set S ij constitutes the contiguous Voronoi cell W i of nucleus i, we get This directly leads to Equation A3 for w i ≡ m i .
We can use the Delaunay tessellation as a multidimensional linear interpolation grid for a field f(x).Each Delaunay simplex is a region with constant field gradient ∇f, such that Given the points substituting the vertices of a simplex and the value of the sampled field at each of these locations, we can calculate the gradient ∇f from the inversion relation (here in d = 3) with the quantities Δx i ≡ x i x 0 etc. and Δf i ≡ f i f 0 (i = 1, 2, 3).This linear interpolation allows us to assign a field value f(x) to every location x inside the tessellation.Locations outside the tessellated space are set to have a field value of 0.

A4. Defining the Exoplanet Shadow
As photoionization is the main loss mechanism of neutral sodium, we need to check if particles are moving inside the exoplanet's shadow.Here, stellar radiation is blocked and lifetimes are increased to 3 hr, typical for charge exchange at Jupiter's moon Io.The following derives our shadow-condition for the position r of any particle inside the simulation boundaries.4. Upper right: Increase of v b to 10 km/s and v M to 60 km/ s.A more equalized distribution emerges, effectively resulting in a constant LOS density as long as the line of sight is not blocked by the exoplanet.Note that the upper panels show the particle densities on a logarithmic scale log n [cm 3 ], whereas the phase curves describe the column densities.
D MC code simulates the sputtering & outgassing of an exomoon at 8 candidate systems to provide density maps of the neutral environment • Exomoon phase-curves that is, spatial & temporal variations of orbiting alkali clouds & tori imprinted in transmission spectra are predicted

Figure 1 .
Figure 1.Characteristic distances of the exoplanet system HD-189733 b: a p semi-major axis of the planet, a s ≡ a ☽ semi-major axis of the exomoon, R Roche Roche-limit, and R Hill Hill sphere radius.All lengths are to scale.A zoom-in on exoplanet can be found in Figure 2.

Figure 2 .
Figure 2. HD-189733 b and the depiction of an exomoon in three different sizes residing in the stability region between Roche limit and Hill sphere.(M,R) ⊕ denotes an exomoon of fundamental parameters equal to Earth, whereas (M,R) Io and (M,R) Enceladus are exomoons characteristic to Jupiter's moons Io and Enceladus.Masses and radii are listed in Table1.

Figure 4 .
Figure 4. Different velocity distributions for the case of WASP-49 b I. Equilibrium temperature is set to 1400 K, tidal heating temperature to 1902 K.For the sputtering distribution we take the corresponding values from Table4.The two sputtering cases correspond to a collisional cascade (α = 3) and single elastic collisional ejection (α = 7/3).The sputtering velocity distributions show a characteristic high-energy tail.

Figure 5 .
Figure 5. Delaunay Tessellation Field Estimator (DTFE) triangular grid against a static grid.The adaptive DTFE tessellation is spanned by the positions of particles in the simulation.It removes the ambiguity of bin size encountered in conventional static grids.The sourcing satellite is located at the center of the particle cloud and orbits the planet represented by the central circle.

Figure 6 .
Figure6.Io orbiting Jupiter with the famous sodium "banana cloud"(Wilson et al., 2002) clearly visible.x-and y-axes show distance from Jupiter in Jupiter radii and the coloring refers to the sodium column density in cm 2 on a logarithmic scale.The underlying Delaunay tessellation connects all particles in the simulation and is used for density calculations.The simulation includes the gravitational accelerations from the icy Galilean moons Europa, Ganymede and Callisto.

Figure 8 .
Figure 8. Equivalent to Figure 7, but for the HD-209458 b system.

Figure 9 .
Figure 9. Equivalent to Figure 7, but for the HD-189733 b system.

Figure 10 .
Figure 10.Particle and line-of-sight column density phase curves for the WASP exoplanets studied in this paper in a photoionization-limited cloud scenario.Particles inside the exoplanet shadow have an increased lifetime of 3 hr.For each system we have run three simulations that consider different satellite masses and radii: Earthsized (⊕), exo-Io (Io), and Enceladus-sized (Enc).We have fit analytic models according to Equations 18 and 19 to phase-curves associated with an exo-Io.Vertical dotted lines in the LOS plots mark the occultation of the satellite by the exoplanet.During transit of the exoplanet, the exomoon will surpass a limited amount of phaseangles illustrated by a horizontal bar on the plots.

Figure 11 .
Figure 11.Continuation of Figure 10 for the other exoplanets considered in this study in a photoionization-limited cloud scenario.

Figure 12 .
Figure 12.HD-189733 b II simulation results.Plotted are the distributions at exomoon phases ϕ s ≡ ϕ ☽ ∈ {0°, 90°, 180°, 270°} from left to right.The first row shows a top-down view of the cloud scenario, where sodium lifetimes are limited by photoionization.The second row corresponds to the line-of-sight view on the cloud.The third and forth row are equivalent to the first two but show the torus scenario, where sodium lifetimes are set to 3 hr.Even in the cloud scenario, where a strong localization due to rapid ionization is expected, a partial torus emerges.Increasing the lifetime results in a denser torus which is depopulated mostly by radiative effects.

Figure 14 .
Figure14.Particle-and line-of-sight (LOS) column density phase curves for the WASP exoplanets studied in this paper in a 3-hr-lifetime torus scenario.For each system we have run three simulations that consider different satellite masses and radii: Earth-sized (⊕), exo-Io (Io), and Enceladus-sized (Enc).Vertical dotted lines in the LOS plots mark the occultation of the satellite by the exoplanet.The present panels do not include occultation effects that strongly counteract peaks in the LOS column.During transit of the exoplanet, the exomoon will surpass a limited amount of phase-angles illustrated by a horizontal bar on the plots.

Figure 15 .
Figure 15.Continuation of Figure 14 in a 3-hr-lifetime torus scenario.The behaviors show less variability, with a steady decline from shadow egress to shadow ingress.Inside the shadow, radiation pressure and Poynting-Robertson drag are absent, which results in an increase in density.At HD-209458 b I, the exomoon orbits very close to the planet at a ☽ = 1.93 R p .Due to significant radiative effects in this system (β = 4), the densities experience a sharp increase after noon as particles pushed toward the exoplanet get gravitationally accelerated and slingshot toward the observer.

Figure 17 .
Figure 17.Consecutive transits and occultations for the WASP-69 b system.Both the geometric exoplanet light curve, ΔF ⋆ / F ⋆ ∼ (R p / R ⋆ ) 2 , and the modeled periodic exomoon spectroscopic transit depth δ Na are shown.The initial phase-shift of the exomoon curve is unknown.We assume a conjunction of star-exoplanet-exomoon -observer at the first transit.Following transits are shown for comparison.Note that the two light curves are not additive as they represent two different phenomena.Available sodium spectral data represent the sum of what is plotted as the orange curve and what might be intrinsic to the exoplanet atmosphere.Distinguishing the contributions remains a challenging task.The light curve of the exoplanet is, therefore, more akin to a reference time axis.

Figure A1 .
Figure A1.WASP-49 b I simulation results.Plotted are the distributions at exomoon phases ϕ ☽ ∈ {0°, 90°, 180°, 270°} from left to right.The first row shows a topdown view of the cloud scenario, where sodium lifetimes are limited by photoionization.The second row corresponds to the line-of-sight view on the cloud.The third and forth row are equivalent to the first two but show the torus scenario, where sodium lifetimes are set to be 3 hr.

Figure A2 .
Figure A2.XO-2N b I simulation results.Plotted are the distributions at exomoon phases ϕ s ≡ ϕ ☽ ∈ {0°, 90°, 180°, 270°} from left to right.The first row shows a topdown view of the cloud scenario, where sodium lifetimes are limited by photoionization.The second row corresponds to the line-of-sight view on the cloud.The third and forth row are equivalent to the first two but show the torus scenario, where sodium lifetimes are set to be 3 hr.

Figure A3 .
Figure A3.WASP-69 b I simulation results.Plotted are the distributions at exomoon phases ϕ s ≡ ϕ ☽ ∈ {0°, 90°, 180°, 270°} from left to right.The first row shows a top-down view of the cloud scenario, where sodium lifetimes are limited by photoionization.The second row corresponds to the line-of-sight view on the cloud.The third and forth row are equivalent to the first two but show the torus scenario, where sodium lifetimes are set to be 3 hr.Just like at HD-189733 b II, the exomoon is capable of building up a partial torus even in a photoionization limited environment.

Figure A4 .
Figure A4.HD-209458 b I simulation results.Plotted are the distributions at exomoon phases ϕ s ≡ ϕ ☽ ∈ {0°, 90°, 180°, 270°} from left to right.The first row shows a top-down view of the cloud scenario, where sodium lifetimes are limited by photoionization.The second row corresponds to the line-of-sight view on the cloud.The third and forth row are equivalent to the first two but show the torus scenario, where sodium lifetimes are set to be 3 hr.

Figure A5 .
Figure A5.WASP-17 b I simulation results.Plotted are the distributions at exomoon phases ϕ s ≡ ϕ ☽ ∈ {0°, 90°, 180°, 270°} from left to right.The first row shows a top-down view of the cloud scenario, where sodium lifetimes are limited by photoionization.The second row corresponds to the line-of-sight view on the cloud.The third and forth row are equivalent to the first two but show the torus scenario, where sodium lifetimes are set to be 3 hr.

Figure A6 .
Figure A6.HAT-P-1 b I simulation results.Plotted are the distributions at exomoon phases ϕ s ≡ ϕ ☽ ∈ {0°, 90°, 180°, 270°} from left to right.The first row shows a topdown view of the cloud scenario, where sodium lifetimes are limited by photoionization.The second row corresponds to the line-of-sight view on the cloud.The third and forth row are equivalent to the first two but show the torus scenario, where sodium lifetimes are set to be 3 hr.

Figure A7 .
Figure A7.Effect of increasing the velocity distribution parameters for the HD-189733 b system with the inclusion of exomoon occultation at noon.Upper left: Values taken from Table4.Upper right: Increase of v b to 10 km/s and v M to 60 km/ s.A more equalized distribution emerges, effectively resulting in a constant LOS density as long as the line of sight is not blocked by the exoplanet.Note that the upper panels show the particle densities on a logarithmic scale log n [cm 3 ], whereas the phase curves describe the column densities.

Table 1
Masses and Radii Used as Simulation Parameters

Table 3
shows the systems we analyze.The fundamental parameters necessary for our simulations are summarized in Table4.From average line-of-sight (LOS) column density observations 〈N〉 over the surface of the star (notations 〈⋅〉 and ⋅ are used interchangeably), we can estimate a mass loss rate via

Table 3
Characteristics of the Eight Exoplanet Systems, Sorted by Exoplanet Semi-Major Axis Tabled are mass M p , radius R p , and semi-major axes a p of the exoplanet.τsyncis the tidal locking timescale of the exoplanet.Ṁth is the calculated thermal mass loss rate responsible for the Alfvén sphere equatorial radius R A .The Hill-sphere radius is given by a Hill in units of planetary radii.ip is the tilt of the orbital plane of the exoplanet relative to the observer.Stellar parameters are captured in Table2.

Table 2
Stellar Parameters of the Eight Systems Considered in This Study, Including Spectral Type, Mass, and Radius in Units Normalized to the Sun

Table 4
Input Parameters for the Particle Evolution Simulation for Each Exoplanet System

Journal of Geophysical Research: Planets
MEYER ZU WESTRAM ET AL.

Table 6
Results of Our Simulations.R Is Equal to 1 δ for Transit Depths δ, Either in Transit or Occultation

Table 7
Oza et al. (2019)ium Columns to Calculations From Oza et al. (2019)Note.Our simulations are based on the same mass-loss rates derived by the model dishoom.Our sputtering simulations show smaller values in the case of a photoionization limited sodium cloud.A torus scenario comparable to Io, on the other hand, overshoots the analytical values fromOza et al. (2019).+ denote the largest and smallest achieved column densities and therefore the range of values that could be retrieved during a single exomoon orbit.