#  Cronomoons: origin, dynamics, and light-curve features of ringed exomoons

Mario Sucerquia,1,2 Jaime A. Alvarado-Montes,3,4 Amelia Bayo,1,2 Jorge
Cuadra,5,2 Nicolás Cuello,6 Cristian A. Giuppone,7 Matías Montesinos,8,2 J.
Olofsson,1,2,9 Christian Schwab,3,4 Lee Spitler,3,4,10 and Jorge I. Zuluaga11.
1Instituto de Física y Astronomía, Facultad de Ciencias, Universidad de
Valparaíso, Av. Gran Bretaña 1111, 5030 Casilla, Valparaíso, Chile
2Núcleo Milenio Formación Planetaria - NPF, Universidad de Valparaíso, Av.
Gran Bretaña 1111, Valparaíso, Chile
3Department of Physics & Astronomy, Macquarie University – Sydney, NSW 2109,
Australia.
4Centre for Astronomy, Astrophysics and Astrophotonics, Macquarie University –
Sydney, NSW 2109, Australia.
5Departamento de Ciencias, Facultad de Artes Liberales, Universidad Adolfo
Ibáñez, Avenida Padre Hurtado 750, Viña del Mar, Chile.
6Univ. Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France.
7Universidad Nacional de Córdoba, Observatorio Astronómico - IATE. Laprida
854, 5000 Córdoba, Argentina
8Escuela de Ciencias, Universidad Viña del Mar, Viña del Mar, Chile.
9Max Planck Institute for Astronomy, Königstuhl 17, D-69117 Heidelberg,
Germany
10Australian Astronomical Optics, 105 Delhi Rd, North Ryde, NSW 2113,
Australia
11SEAP research group, Instituto de Física, FCEN, Universidad de Antioquia –
Calle 70 No. 52-21, Medellín, Colombia<EMAIL_ADDRESS>
(Accepted 30 November 2021. Received 03 August 2021; in original form 22 July
2021)

###### Abstract

In recent years, technical and theoretical work to detect moons and rings
around exoplanets has been attempted. The small mass/size ratios between moons
and planets means this is very challenging, having only one exoplanetary
system where spotting an exomoon might be feasible (i.e. Kepler-1625b i). In
this work, we study the dynamical evolution of ringed exomoons, dubbed
cronomoons after their similarity with Cronus (Greek for Saturn), and after
Chronos (the epitome of time), following the Transit Timing Variations (TTV)
and Transit Duration Variation (TDV) that they produce on their host planet.
Cronomoons have extended systems of rings that make them appear bigger than
they actually are when transiting in front of their host star. We explore
different possible scenarios that could lead to the formation of such
circumsatellital rings, and through the study of the dynamical/thermodynamic
stability and lifespan of their dust and ice ring particles, we found that an
isolated cronomoon can survive for time-scales long enough to be detected and
followed up. If these objects exist, cronomoons’ rings will exhibit gaps
similar to Saturn’s Cassini Division and analogous to the asteroid belt’s
Kirkwood gaps, but instead raised due to resonances induced by the host
planet. Finally, we analyse the case of Kepler-1625b i under the scope of this
work, finding that the controversial giant moon could instead be an Earth-mass
cronomoon. From a theoretical perspective, this scenario can contribute to a
better interpretation of the underlying phenomenology in current and future
observations.

###### keywords:

planets and satellites: rings – techniques: photometric – methods: analytical
– techniques: photometric – planets and satellites: dynamical evolution and
stability – planets and satellites: detection.

††pubyear: 2021††pagerange: Cronomoons: origin, dynamics, and light-curve
features of ringed exomoons–Cronomoons: origin, dynamics, and light-curve
features of ringed exomoons

##  1 Introduction

Figure 1: Schematic illustration of the system considered in our model. We
highlight the semi-major axis of the planet ($a_{\rm p}$) and the moon’s
($a_{\rm\mathrm{stop}}$); the mass and radius of the star, the planet, and the
moon ($M_{\rm\star}$, $M_{\rm p}$, and $M_{\rm m}$; $R_{\star}$, $R_{\rm p}$
and $R_{\mathrm{m}}$, respectively). The ring’s tilt, its inner and outer
radii are also shown ($i_{\rm R},R_{\mathrm{m}},r_{\rm i}$ and $r_{\rm o}$),
as well as the transit impact parameter $b$.

The first confirmation of  an exomoon would certainly be  an important step to
refine both observational and theoretical predictions of exoplanetary systems.
Finding an ensemble of moons would enable a new tool, giving us more leverage
to better understand the formation and dynamical evolution of planetary
systems. According to formation models, moons should be common among Jupiter-
like planets discovered within the circumstellar habitable zones of their host
stars (see e.g. Heller & Barnes 2013; Heller et al. 2014; Martínez-Rodríguez
et al. 2019). As this class of exoplanet is known to be numerous, a successful
confirmation (or lack thereof) will have strong implications for models of
migration and habitability of moons.

There were recently six spurious claims for exomoon candidates by Fox &
Wiegert (2021), which have been proven incorrect from both an observational
(Kipping & Teachey, 2020) and a stability perspective (Quarles et al., 2021).
Does this reflect a lack of moons? or is it insufficient sensitivity of
current observations? Direct detection of such satellites could be achievable
for the next generation of telescopes like the European Extremely Large
Telescope (E-ELTS; Peters & Turner, 2013), whereas radial velocity and
astrometric methods are probably incapable of reaching these targets
(Vanderburg et al., 2018). Perhaps the only feasible way is the analysis of
light curve data in systems with multiple transits of planets orbited by large
exomoons. This would involve the detection of one or more second-order
dynamical effects on planetary light curves: Transit Duration Variations
(TDVs), Transit Timing Variations (TTVs; Kipping et al., 2012; Kipping et al.,
2013; Kipping et al., 2014), and Transit Radius Variations (TRVs; Rodenbeck,
Heller & Gizon, 2020).

The aforementioned effects (TDVs, TTVs, TRVs) could be useful to assess the
physical and orbital properties of exoplanets and exomoons. So far, these are
the most promising techniques to reveal the presence of (yet undetected) low-
mass companions such as rocky planets, rings, moons, and even smaller objects
(see e.g. Seager & Mallén-Ornelas 2003; Kipping 2009a; Zuluaga et al. 2015;
Sucerquia et al. 2020a, b, and references therein). TTVs are related to the
changes in the planet’s orbital position at the beginning of its transit which
could reveal the presence of nearby planets and/or satellites; while TDVs are
variations in the duration of transits, emerging as a consequence of changes
in the velocity of the planet during the transit, produced by the presence of
a companion satellite. Finally, TRVs are due to small time-dependent
contributions of moons to the transit depth, which happens during a transit
when moons lie in the line of sight of the observer. Further information about
the presence of rings, cometary tails, or other circumplanetary structures
could be encoded in the photometric shape of the ingress and egress stages of
planetary transits (see e.g, Zuluaga et al. 2015; Rappaport et al. 2018;
Arkhypov, Khodachenko & Hanslmeier 2021).

Kepler-1625b i (Teachey & Kipping, 2018) is currently the most convincing
exomoon candidate spotted by the detection of a mutual transit,  and by re-
drawing a new linear ephemeris through all data points as of yet, Kipping
(2021) has calculated a TTV of $\sim 17$ min for this system. However,
interpretation of this observational result is still being debated (see e.g.
Heller, Rodenbeck & Bruno 2019; Kreidberg, Luger & Bedell 2019; Teachey et al.
2020). The main uncertainty in interpretation lies in the planet–moon physical
and orbital characteristics estimated for this system: the mass-ratio $M_{\rm
m}/M_{\rm p}\sim 10^{-3\,\text{to}\,-2}$ (where $M_{m}$ and $M_{p}$ are the
mass of the moon and the planet, respectively) in particular suggests that
Kepler-1625b i is unlikely to have formed in the Galilean moon formation
framework (Canup & Ward, 2006).

The aim of this work is to study the dynamical stability of a gravitationally
bound system composed of a single-planet being orbited by a ringed moon (see
Fig. 1). We name these objects cronomoons after the planet Saturn (Cronus, in
Greek mythology) and Chronos (the epitome of time) in relation to the set of
unique observational signatures that this kind of system might exhibit, namely
large TRVs but small TTVs and TDVs. However, this situation could also occur,
for example, in exoplanet systems where moons have large, dense, and thick
atmospheres sufficient to obscure a large portion of the stellar disc
(Rodenbeck et al., 2020), but with considerably small bulk masses unable to
cause detectable perturbations on the TTV and TDV signals.

Figure 2: Possible pathways to form cronomoons, according to Section 2: i)
corresponds to circumsatellital disc remnant, ii) to moon-moon collisions,
iii) to exocomet collisions, and iv) to ploonets. Note that the scenarios from
i) to iii) are planetocentric frameworks, whereas iv) is centred on the star.

Using approaches such as pure-dynamical simulations including tidal
interactions, the stability of hierarchical systems of moons (i.e. moons
around moons) has been proposed and tested (see Forgan 2018; Kollmeier &
Raymond 2019; Rosario-Franco et al. 2020). In this regard, to analyse the
stability and lifespan of large sets of particles in disc-like structures
around moons, it is important to construct a framework able to explain
observations displaying baffling signals. For instance, recent analysis of
some light curves have suggested the existence of planets whose densities are
abnormally low, also called super-inflated planets (see. e.g. Akinsanmi et al.
2020; Piro & Vissapragada 2020). However, such planets could also have a
normal density and be surrounded by systems of rings opaque enough to eclipse
a significant amount of the stellar disc (see e.g. Zuluaga et al. 2015 and
references therein). Also, the magnitude and detectability of TTV and TDV
signals for different theoretical populations of exomoons have been
constrained (Sucerquia et al., 2019, 2020a) and delimited (Kipping & Teachey,
2020) in terms of the semi-major axes and masses of systems of moons,
indicating that any TTVs or TDVs larger than these constrained values must
have a different origin.

Given the unusual characteristics of the exomoon candidate Kepler-1625b i, is
it possible that we may be observing a different type of system? May this moon
have a ringed structure making it look bigger than it really is? What are the
detectable imprints of such a system on the light curve? To answer these
questions, we explain in Section 2 some possible scenarios for the formation
of hypothetical ringed exomoons, while we test their long-term dynamical and
thermodynamic stability in Section 3. The detectability of such moons using
current and future instrumentation is studied in Section 4, and we analyse in
Section 5 the system  Kepler-1625b i under the scope of this work. Finally, we
discuss in Section 6 the main findings of this research and its possible
consequences on the search and characterisation of exomoons and exorings.

##  2 Formation scenarios of cronomoons

Rings seem to be an exclusive feature of the giant planets in the Solar
System. Each of them has a set of circumplanetary rings with quite dissimilar
characteristics regarding composition, surface density, and optical
properties. These differences suggest the existence of different formation and
evolutionary pathways of these ring systems (see, e.g. Hyodo & Charnoz 2017).
Moreover, it is already known that some minor Solar System bodies, such as the
centaurs Chariklo and Chiron as well as the trans-Neptunian object Haumea,
bear ring system that have been detected by stellar occultations. The
formation mechanisms of these rings are still unknown, although it is presumed
that they could be the result of a collision between small bodies (see e.g
Braga-Ribas et al. 2014; Ortiz et al. 2015, 2017 and citing articles).

The variety of ring species known so far, the fact that small-mass bodies may
possess them, plus the stability studies developed in Section 3 lead us to
think that moons may also have a ring system around them. We discuss below on
what these formation routes might be and their possible evolutionary paths.

  1. 1.

Circumsatellital disc remnant: One of the main but still debated scenarios to
explain the origin of rings-systems around Saturn draws an analogy to
circumstellar discs during the process of planetary formation, proposing that
Saturn’s rings are the remnants of Saturn’s sub-nebula disc (Pollack, Summers
& Baldwin, 1973; Pollack, 1975). However, because the average chemical
composition of these rings is different from that of Saturn’s satellites, it
has been argued that the formation of the rings occurred at a later stage, as
explained in detail by Charnoz et al. (2009).

Still, it is worth noting that growing planets are surrounded by remnant discs
which can form regular satellites, and during their formation some parts of
the discs can thrive around the satellites (see case i) in Fig. 2). For
instance, a dusty circumplanetary disc with a bulk mass of $\sim
10^{-3}\;M_{\oplus}$ has been recently observed around PDS70 b, and there is
evidence of another one surrounding a the planet c, which have been
interpreted as a moon-forming region (Isella et al., 2019). Therefore, planets
and regular satellites share an analogous formation channel where discs
composed of dust (and sometimes ice) can survive the coalescence of material
and last for long time-scales that would allow us to detect them either via
their emission properties (as PDS70 b ,c), or the imprints they would produce
on the light curve of their host planet/satellite.

  2. 2.

Moon-moon collisions: The rate at which moons migrate due to tidal angular
momentum exchange depends mainly on their masses and the properties of the
planet. Supposing a steady state for planetary migration, giant planets
harbouring moons with dissimilar compositions and masses (e.g. one dense and
rocky, the other light and composed of lighter material) will induce different
rates of orbital migration through tidal migration as they are gravitationally
perturbed. The worst-case scenario of this combination of gravitational
effects could lead the system to chaotic orbital configurations that end in
collisions with each other, or between the moons and the planet, leading to
the formation of rings of particles around them.

In the Solar System, giant planets have a large number of satellites around
them. Similarly, we expect massive exoplanets (at large circumstellar
distances) to harbour large numbers of satellites, rings, and exomoons.
However, different evolutionary models suggest that short-period (Alvarado-
Montes & García-Carmona, 2019) and ultra-short-period planets (Alvarado-Montes
et al., 2021) have undergone a slow but progressive migration to reach the
currently observed positions. In such processes, different planetary
properties will be affected, thereby affecting their surroundings. For
example, the planetary Hill radius is compressed (Spalding, Batygin & Adams,
2016), allowing a stronger gravitational interaction between its satellites
and other external bodies such as the host star. Also, the angular momentum
exchange between planets and their satellites is affected by migration as
tidal interaction become stronger, which will modify the planetary rotation
rate and the orbital mean motion of moons.

At present, the problem of the orbital stability and fate of multiple moons
migrating in conjunction with the planet has not yet been studied in depth,
but processes such as those described above could move these moons away from
their possible protective resonances and introduce chaos into the system (see
e.g. Quarles et al. 2021, and references therein), which may result in
collisions and mergers between them. The by-product of such interactions could
be the depletion of these moons, or the survival of a few of them (or a single
massive moon) which transiently possess surrounding ring-like particle systems
after the collision (see case ii in Fig. 2), and are expected to survive for
longer than hundreds of Myrs (Hesselbrock & Minton, 2017; Crida et al., 2019).
However, the moon-moon collisions giving rise to such rings are stochastic
processes that are not yet constrained by observations, with a lack of
stringent theoretical computations from population synthesis studies. Whether
the rings resulting from these collisions are made of dust or ice depends on
the composition, density and their reservoirs of icy material, which will be
given by the distance of the host planet from the star.

  3. 3.

Exo-cometary collisions: Currently, cometary activity has been detected in
about thirty stars. These objects have been identified, on the one hand, due
to their asymmetric light curves that reveal the presence of their tails, and
on the other hand, because of the chemical fluctuations from starlight passing
through their tails and transiently modifying the stellar spectrum (for a
complete summary of these systems see Strøm et al. 2020).

In the Solar System, long-period comets emerged after the orbital exchange of
Neptune and Uranus described in the Nice model (Tsiganis et al., 2005).
Exoplanetary systems are also prone to processes where systems of low-mass icy
or dusty objects have high-eccentricity orbits around their star. Under the
presence of short-period planets retaining their moons (Namouni, 2010;
Spalding et al., 2016), it is likely that exocomets will collide with such
planets and moons (this scenario will be the subject of a future work, led by
the same authors of this paper). Also, after a long-lasting gravitational
interaction with the planet, such exocomets may end up located in short-period
orbits which increase the probability of collision to form tilted rings around
planets or moons. Icy exocomets colliding against moons and planets motivated
the calculations done in this work using icy particles as possible
constituents of close-in exorings (see case iii in Fig. 2), whose life time-
scales will be computed in Section 3.2.

  4. 4.

Ploonets: Ploonets are moons that have been moved from circumplanetary to
circumstellar orbits by the effect of, for example, tides (Sucerquia et al.
2019). These ploonets manage to settle transiently into orbits close to their
host planet, but they may have eccentric orbits. These orbital characteristics
make them very prone to collide either with their host planet or other nearby
planets and moons, as depicted in the middle panel of Fig. 2, which is a
scenario that we will study in a future paper. Depending on the orbital
inclinations of ploonets, collisions and close approaches could form rings-
systems around planets and their moons.

  5. 5.

Other formation channels: 1\. A common feature between the giant planets of
the Solar System is the presence of irregular satellites whose origin is
explained via gravitational capture of wanderer objects in the planetary
neighbourhood. If rings around bodies such as Chariklo, Chiron, and Haumea
survive to the capture process of giant planets, this might form a cronomoon
in the Solar System. 2. Under the assumption of the existence of submoons
(Kollmeier & Raymond, 2019), in analogy to satellital migration around short-
period planets, a submoon could undergo inward migration until reaching its
Roche limit with respect to the host moon, spreading its material and creating
a system of rings (Canup, 2010).

##  3 Ring’s lifespan, shape, and size

Figure 3: Upper panel: asymptotic moon’s semi-major axis
($a_{\rm\mathrm{stop}}$) as a function of the planetary mass (blue solid
lines). For reference, the principal and secondary (stability limit) Hill
radii are drawn in red and orange lines, respectively. The intersection
between dotted lines are the values used in our simulations (Section 3.2.2).
Lower panel: final position $a_{\rm\mathrm{stop}}$ for planet–moon systems
having different moon-to-planet mass ratios. We assume here the physical
properties of a Jupiter-like planet and a Titan-like moon.

In this section, we test the dynamical stability of a system composed of a
Solar-like host star and a close-in giant planet orbited by a massive ringed
moon (see Fig. 1 and the nomenclature therein used throughout this work).
Planet-moon mass ratios are assumed to be within the framework of Galilean’s
moons formation, where $M_{\rm m}/M_{\rm p}\sim 10^{-4}$. Also, moon’s orbital
positions (i.e, moon’s semi-major axes) will be constrained by adopting the
satellite migration model around close-in exoplanets described by Alvarado-
Montes, Zuluaga & Sucerquia (2017) and tested in Sucerquia et al. (2020a) (see
Section 3.1), allowing us to reduce the number of unknown parameters.

Figure 4: Left-hand panel: Sublimation time-scales for water–ice spheres (with
bond albedos of $0.95$) of different sizes and placed at average distances
$\bar{a}$ from the star. Right-hand panel: Poynting-Robertson decay time-
scales for different tilts of the rings, placed at $\bar{a}$ from the star.

Seeking a more significant contribution of cronomoons in planetary light
curves, we set up our simulations so that planets have close-in circular
orbits from their host star, with distances ranging from 0.3 to 1.0 au. At
such orbital positions the system is exposed to large amounts of stellar
radiation, which makes volatile material highly unstable. For reasons
explained in Sections 3.2, 4, and 6; and in spite of the high irradiation
received by the system, we will not restrict our study to rings composed of
refractory material (the so called warm-exorings described by Schlichting &
Chang 2011); that is, we will also include those composed of ices. In this
work we have dubbed them dusty and icy rings when composed of silicates and
ices, respectively.

###  3.1 Tidal migration of host moons

To study the dynamical stability of cronomoons we need to constrain the
orbital and physical features of host moons. The semi-major axis of the moon
will be assumed at a position also known as Satellite’s Tidal Orbital Parking
($a_{\rm\mathrm{stop}}$, Sucerquia et al. 2020a), reached due to orbital
migration triggered by tidal interactions with the parent planet. According to
the model in Alvarado-Montes et al. (2017), which we adopt in this work prior
to the evolution of cronomoons, both the evolution of physical size (i.e. the
planetary radius $R_{\rm p}$) and interior dissipative properties of the
planet modify the final fate of migrating moons. This can have important
consequences for the survival of exomoons in compact systems and even for the
orbital evolution of short-period (Alvarado-Montes & García-Carmona, 2019) and
ultra-short-period planets (Alvarado-Montes et al., 2021).

For this work, as some important criteria for orbital migration of moons from
Alvarado-Montes et al. (2017) share some of its attributes with previous
classical models (e.g. Barnes & O’Brien 2002 and Sasaki, Barnes & O’Brien
2012), the evolution of the planetary rotational rate ($\Omega_{\rm p}$) and
the planet and moon’s mean motion ($n_{\rm p}$ and $n_{\rm m}$, respectively)
follow:

$\begin{split}\frac{\mathrm{d}\Omega_{\rm
p}}{\mathrm{d}t}=-\frac{3}{2}\frac{k_{\rm 2}}{Q}\frac{R_{\rm
p}^{3}}{\kappa^{2}G}\left[\frac{(GM_{\rm\star})^{2}}{M_{\rm p}a_{\rm
p}^{6}}\,\mathrm{sgn}\,(\Omega_{\rm p}-n_{\rm p})+\right.\\\
\left.\frac{M_{\rm m}^{2}}{M_{\rm p}^{3}}n_{\rm
m}^{4}\,\mathrm{sgn}\,(\Omega_{\rm p}-n_{\rm m})\right].\end{split}$ (1)
$\frac{\mathrm{d}n_{\rm m}}{\mathrm{d}t}=-\frac{9}{2}\frac{k_{\rm
2}}{Q}\frac{M_{\rm m}R_{\rm p}^{5}}{G^{5/3}M_{\rm p}^{8/3}}n_{\rm
m}^{16/3}\mathrm{sgn}\,(\Omega_{\rm p}-n_{\rm m}).$ (2)

Where the moon’s semi-major axis ($a_{\rm m}$) can be found via Kepler’s third
law as $a_{\rm m}\approx(GM_{\rm p}/n_{\rm m}^{2})^{1/3}$ and $\kappa$ is the
planetary gyration radius. Also, the dissipative properties of the planet are
represented by the Love number $k_{\rm 2}$, a coefficient standing for its
deformation capabilities; and the tidal quality factor $Q$, which represents
the tidal energy dissipated per rotational period.  For both $k_{\rm 2}$ and
$Q$ we use the formalism of Alvarado-Montes et al. (2017) which follows
Ogilvie (2013) and Guenel et al. (2014). Within this formalism we assume
initial $k_{\rm 2}$ and $Q$ as those of a Jupiter-like planet (Mathis, 2015)
and make both quantities evolve during the orbital migration of moons. This
locates moons at a final stable semi-major axis (Sucerquia et al., 2020a). On
the contrary, if both of these quantities are assumed static, moons would
migrate inwards and eventually cross their Roche limit as in Barnes & O’Brien
(2002) and Sasaki et al. (2012).

Tidal models state that moons can migrate outwards until they reach an orbital
position where synchronisation occurs, i.e. $\Omega_{\rm p}\approx n_{\rm m}$.
Then, inward migration follows since the torques on the tidal bulge of the
planet are reverted. However, in Alvarado-Montes et al. (2017) inward
migration is suppressed because $\Omega_{\rm p}$ never reaches $n_{\rm m}$,
and this might have notable consequences for the moon’s final fate, including
but not limited to: detachment from the host planet (Sucerquia et al., 2019),
or orbital parking on an asymptotic semi-major axis (Sucerquia et al., 2020a).
We calculate $a_{\rm\mathrm{stop}}$ under the model by Alvarado-Montes et al.
(2017), and assume a departing moon’s semi-major axis of $a_{\rm m}=3\,R_{\rm
p}$ and a planetary rotational period of 13 h. $a_{\rm\mathrm{stop}}$ is the
value we use to set up the numerical simulations (Section 3) and to assess
their detectability (Section 4). The main outcomes of this section are
presented in Fig. 3: the upper panel presents the evolution of
$a_{\rm\mathrm{stop}}$, the planetary Hill radius $R_{\rm Hill}$, and the
stability limit of $0.48\,R_{\rm Hill}$ (Domingos, Winter & Yokoyama, 2006)
for planets at $0.3$, $0.5$, $0.8$ and $1.0$ au; while the lower panel shows
the relationship between $a_{\rm p}$ and $a_{\rm m}$ for $M_{\rm m}/M_{\rm
p}=10^{-5},10^{-4}$ and $10^{-3}$, from blue to yellow, respectively.

###  3.2 Thermodynamic stability: rings’ lifespan

In general, isolated planetary rings are gravitationally stable and their
spatial domain is externally constrained by the Roche Limit111This is defined
as the orbital limit inside which a body’s self-attraction (particles) is
exceeded by the tidal forces of a primary body (host moon)., $a_{\rm Roche}$,
to prevent any rings’ particles to form larger bodies (Roche, 1849; Crida &
Charnoz, 2012),

$a_{\rm
Roche}=r_{\mathrm{m}}\left(2\frac{\rho_{\mathrm{m}}}{\rho_{\mathrm{part}}}\right)^{\frac{1}{3}}$
(3)

where $r_{\mathrm{m}}$ and $\rho_{\mathrm{m}}$ are the radius and density of
the satellite, respectively; and $\rho_{\mathrm{part}}$ is the density of the
ring’s particles. Still, at short distances from the star, it is expected that
icy rings are promptly sublimated leaving behind rings composed by dusty
particles (‘dusty rings’). Such lifespan will be limited by the Poynting-
Robertson effect. Below, we present the framework to calculate both effects in
the context of this work.

####  3.2.1 Sublimation of icy rings

Given that rings around moons can be composed of particles with icy envelopes,
we study the evolution of spherical blocks of ice adopting the model proposed
in Sucerquia et al. (2019). This will allow us to estimate the typical time-
scales of icy rings’ orbital evolution around close-in giant planets, while
calculating the sublimation of icy envelopes to constrain the typical lifespan
of this type of rings systems. We assume these particles to be rapidly
rotating spheres (rotation period much shorter than orbital period) located in
circular orbits. The equilibrium temperature of these particles is given by:

$T_{\rm eq}=\left[\frac{L_{\rm\star}\,(1-A)}{16\,\pi\,\sigma_{\rm
SB}\,a^{2}}\right]^{1/4},$ (4)

where $A$ is the bond albedo, $L_{\rm\star}$ the stellar luminosity, $a$ the
orbital radius, and $\sigma_{\rm SB}$ the Stefan–Boltzmann constant. In
addition, water from the icy particles sublimates at temperature $T$ at a rate
$Z_{\mathrm{r}}$ (in kg m-2 s-1) given by (Estermann, 1955; Bohren & Albrecht,
1998),

Figure 5: Comparison between icy (left-hand panels) and dusty (right-hand
panels) rings for a moon orbiting a planet following a circular orbit at
$a_{\rm p}=0.3$ au. The rings have an initial inclination of $i=0$° and
circular orbits (i.e. $e=0$). The lower panels show $\delta e$ and $\delta a$
of the rings’ particles, which are the maximum semi-amplitude of their
eccentricity and the semi-major axis, respectively. The state showing the
surviving (white), collided (red), and eccentric (blue) particles correspond
to a system evolution of 100 yr.

$Z_{\mathrm{r}}=e_{\mathrm{sat}}(T)\left(\frac{m_{\mathrm{w}}}{2\pi
RT}\right)^{1/2}\exp\left[\frac{2\,m_{w}\,\sigma_{i}\,(T)}{\rho_{i}(T)\,r\,RT}\right]\,\,\,.$
(5)

with $e_{\mathrm{sat}}(T)$, $\sigma_{\mathrm{i}}(T)$, $\rho_{\mathrm{i}}(T)$
the saturation vapour pressure, the surface ice tension, and the surface ice
density, respectively. $m_{\mathrm{w}}$ is the molecular weight of water and
$R$ the universal constant of ideal gas. All these quantities are adopted from
standard thermodynamic expressions by Andreas (2007), which have widely been
applied to study sublimation processes in bodies like the Moon, Saturn’s rings
(Hedman, 2015), and comets (Vincent et al., 2016). With this in mind, the
total water mass loss rate, $\dot{M}_{\mathrm{s}}$, can be calculated as
follows:

$\dot{M}_{\mathrm{s}}=-Z_{\mathrm{r}}S_{\rm
pl}(T)=-Z_{\mathrm{r}}\left[\frac{6\,\sqrt{\pi}\,m}{\rho_{\mathrm{i}}(T)}\right]^{2/3},$
(6)

where the sphere’s surface area $S_{\rm pl}(T)$ has been written in terms of
its mass $m$ and $\rho_{\rm i}(T)$ (in kg m-3), which is a function of
temperature:

$\rho_{\rm i}(T)=916.7-0.175\,(T-273.15)-5.0\times 10^{-4}(T-273.15)^{2}$ (7)

and agreeing within about 1% of the tabulated values for ice density for
temperatures from 13 to 273.15 K (Hobbs, 1974). Then, the temporal evolution
of the ice block’s mass reads:

$m^{1/3}-m_{\rm o}^{1/3}=Z_{\mathrm{r}}\left(\frac{1}{\rho_{\rm
i}(T)}\right)^{2/3}\left(\frac{4\pi}{3}\right)^{1/3}\;(t_{\rm o}-t),$ (8)

where the initial mass $m_{\rm o}$ corresponds to an initial time $t_{\rm o}$.
Thus, the time necessary to fully sublimate the particles’ ice envelope is:

$t_{\rm ev}=\frac{1}{Z_{\rm r}}\left(\frac{4\pi\;m_{\rm o}}{3\;\rho_{\rm
i}(T)^{2}}\right)^{1/3}.$ (9)

The lifespan of icy rings around exomoons is presented in the left-hand panel
of Fig. 4. There, the constituent particles of the ring are assumed to be
spheres with radii ranging from $1$ to $10$ m (parallel to the particles in
Saturn’s rings, see Brilliantov et al. 2015) and placed at semi-major axes
$\bar{a}=0.5-1.0$ au. Note that distances $\bar{a}<0.6$ au are harmful to icy
rings, with corresponding lifespans of less than a year.  However, less
luminous stars might allow cronomoons with icy rings and orbiting close-in
planets to exist for longer time-scales. Also, icy rings located at larger
distances can survive up to time-scales of the order of Myr. In fact, the
calculations presented in Fig. 4 correspond to an underestimation of the
actual values for ice sublimation (i.e. we are not taking into account the
actual shape of the particles): irregular or asymmetric surfaces would produce
regions where the incidence and absorption of stellar radiation is not
homogeneous, affecting thus the sublimation of the rings’ particles.
Furthermore, when rings’ tilts are either low or co-planar to the moon’s
orbit, the shadow projected on the rings’ particles from the planet and the
moon will shield the rings from stellar radiation.

####  3.2.2 Poynting-Robertson decay

We assume that once icy rings are depleted of volatile material, only dusty
particles are left in the ring. Then, the Poynting-Robertson effect acts on
these dusty particles and triggers their orbital decay and lost. The time-
scales for the Poynting-Robertson particle decay ($t_{\mathrm{PR}}$) for
‘warm-exorings’, as derived by Schlichting & Chang (2011), reads:

$t_{\mathrm{PR}}\sim\frac{\pi c^{2}\Sigma}{\sin{i}\>(L_{\rm\star}/4\pi a_{\rm
p}^{2})\;Q_{\rm PR}(5+\cos{i}^{2})}$ (10)

where $\Sigma=109$ g cm-2 is the mass surface density of the ring, assumed
similar to Saturn’s B-ring (Hedman & Nicholson, 2016), $c$ the speed of light,
$i$ the ring’s tilt (i.e. normal angle with respect to the planet’s orbital
plane), $a_{\rm p}$ the planetary semi-major axis, and $Q_{\rm PR}$ the
radiation pressure efficiency factor which was assumed equal to 0.5. The
right-hand panel of Fig. 4 shows the value of $t_{\mathrm{PR}}$ for several
combinations of averaged star–moon distances $\bar{a}$ and ring’s tilt $i$.
According to Equation 10 and both panels of Fig. 4, ringed moons with high
inclinations can have lifespans of the order of Myr.

Figure 6: Numerical simulation outcomes after $\sim 750.000$ periods (in
average) for icy rings having a tilt of $0^{\circ}$. $\delta a$ and $\delta e$
is the maximum semi-amplitude of the excitation of the ring particles’ semi-
major axis and eccentricity, respectively. Points in red are those particles
that collided with the satellite.

Figure 7: Same as Fig. 6 but for a rings’ tilt of $30^{\circ}$.

###  3.3 Dynamical stability: rings’ size

After orbital migration of the moon due to planet-moon tidal interactions has
already occurred, we adopt the model in Alvarado-Montes et al. (2017) to find
the corresponding $a_{\rm\mathrm{stop}}$ (Sucerquia et al., 2020a) for a given
set of $a_{\rm p}$, as well as planet and moon’s masses, $M_{\rm p}$ and
$M_{\rm m}$ (see Section 3.1 and Fig. 3). Once the parent moon is fixed at the
aforementioned orbital position, the dynamical stability of its rings’
particles is studied by integrating the equations of motion of a set of test-
particles evolving under the gravitational influence of the moon, the planet,
and the star. This allows us to constrain the inner and outer rings’ radii,
rings’ stability, and lifespan.

To set up the initial conditions of the rings’ particles (i.e. positions and
velocities) we used $r_{\mathrm{Roche}}$ as the maximum outer limit of the
rings, assuming that both the planet and the satellite are completely
spherical and their orbits are co-planar. Under these assumptions, the ring-
systems maintain their shape as long as a significant depletion of particles
via thermodynamic processes is averted (see previous Section). Also, the
initial positions of the rings’ particles were generated randomly by
implementing blue-noise sampling algorithms222We use the fibpy algorithm
(https://github.com/matt77hias/fibpy), which is based on the work of Vogel
(1979)., in order to avoid oversampling problems.

We prepared two sets of simulations: 1) icy rings having particles with
densities similar to water ice, $\rho=0.92$ g cm-3 (left-hand panel of Fig.
5); and 2) dusty rings composed of refractory particles with densities similar
to dust, $\rho=1.50$ g cm-3 (right-panel of Fig. 5). The host is a Titan-sized
moon with a density of $\rho_{\mathrm{m}}=3.00$ g cm-3, which orbits a Saturn-
like planet. Both sets of rings were divided into two subsets according to the
rings’ inclination, either $0^{\circ}$ or $39^{\circ}$. These angles were
arbitrarily selected to assess the effect of the inclination on the orbital
evolution of rings’ particles, while avoiding values greater than 39.2º where
Kozai–Lidov resonances appear (Kozai, 1959; Lidov, 1962) and produce evident
effects in orbiting ringed structures (see e.g. Sucerquia et al. 2017). For
each subset, we set the initial position of the planet at $a_{\rm p}=$ 0.3,
0.5, 0.8, and 1.0 au.

The motivation to include icy rings rises from the particles’ lifespan
calculations presented in Fig. 4 and the possible formation mechanisms that
were explained in Section 2. The tangible difference between the simulations
of icy and dusty rings is their initial size which is calculated via equation
3. As it can be seen in Fig. 5, icy rings are much richer dynamically and
present more regions of resonance, while dusty rings are more stable and
present a lower scattering of the constituent particles. The systems were
allowed to evolve over $\sim 750.000$ periods of the rings’ particles, using
the N-body integrator REBOUND with a leapfrog approach (Rein & Liu, 2012) and
a time-step of 10-5 yr. Figs 6 and 7 show the outcomes of the simulations for
inclinations of $0$° and $30$°, respectively.

The upper panels of Fig. 5 show a snapshot of the simulations after 100 yr of
evolution, where particles that collided with the satellite are marked with
red and eccentric particles are depicted in blue. The lower panels of Fig. 5
(as well as Figs 6 and 7) present the semi-amplitude of the excitation in the
eccentricity ($\delta e$) and semi-major axis ($\delta a$) experienced by the
rings’ particles due to the dynamical perturbations produced by the planet and
the star. For the case illustrated in Fig. 5 we used 10.000 particles, while
the systems presented in Figs 6 and 7 were composed of 1000 test particles.

Under the previously mentioned conditions, we found that ringed moons are
stable structures with regions that are strongly affected by the periodic
perturbations of the planet and the star. Increasing the eccentricity or
modifying the semi-major axes of the particles favoured the creation of gaps
in the rings.  For example, in the system with $a_{\rm p}$ = 0.3 au there are
at least three notorious regions of mean motion resonances (MMR) with the
planetary semi-major axis (i.e. 1:15, 1:20 and probably 2:25 MMR in the lower
panels of Fig. 5).  In Figs 6 and 7 we identify the location of the most
important $N/1$ MMR between particles and the planet’s position that roughly
coincides with some spikes. We think that the additional spikes observed in
$\delta_{a}$ and $\delta_{e}$ could be a consequence of higher order MMR (e.g,
$N/2$, $N/3$, etc) and secular resonances that could be present within this
region (e.g. the $\nu 6$ resonance due to the star and mutual inclinations
between the bodies, see Murray & Dermott 2000). However, a thoughtful analysis
of the resonances is beyond the scope of this paper and we leave it for a
forthcoming work.

Particles that were affected by resonant perturbations increased their
eccentricity, becoming prone to collide with surrounding particles of the
system. This, in addition, could increase the particles’ orbital decay rate
and thus the size of resonance gaps. Nevertheless, we did not include the
effect of such collisions in our simulations, so the gaps did not increase
their size in a prominent fashion. Additionally, we have shown here that
particles remained fairly well-confined inside the Roche boundary of the
system and, therefore, it is reasonable to think that their assumed initial
size was an optimal outer radius. This assumption is used in the next section
to establish the detectability of ringed moons.

##  4 Effects in light curves

Cronomoons have effects on the light curves of transiting exoplanets that are
similar to those produced by single moons, but adding special and unique
features. For example, they can modify the relationship between depth,
duration, and timing of the transit by miscalculations of some physical
parameters of the satellite, the planet, and the star.

Moons orbiting around planets produce TTV and TDV effects on planetary light
curves (Kipping, Fossey & Campanella, 2009). In the case of massive moons,
these would push their parent planet further away from their mutual centre of
mass than less massive moons do, while also changing more significantly the
planet’s position and velocity. However, if the actual size of the moon is
overestimated (see Zuluaga et al. 2015), TTVs and TDVs could be imperceptible.
An example of this is a low-mass moon with a ring system that makes it look
much larger than it really is (see Section 5).

###  4.1 Transit depth and TRV of ringed exomoons

Owing to their extended rings-system, some transiting exoplanets can exhibit
apparent sizes that exceed their actual size. This is due to the fact that the
transit technique only cares about the portion of the stellar disc that is
obscured during a transit; thus, larger areas would immediately mean bigger
planets. However, ringed planets can be smaller than what is predicted from
their transit properties (see Zuluaga et al. 2015). In those cases, the
fractional change in stellar flux or transit depth, given by
$\delta=(F_{\star}-F)/F_{\star}$, mainly depends on the tilt or orbital
inclination of the rings towards the observer ($i_{\mathrm{R}}$), the inner
and outer radii of the rings ($f_{\mathrm{i}}$ and $f_{\mathrm{o}}$
respectively, in units of $R_{\rm p}$); and the opacity of the rings ($\tau$)
defined in terms of the light that they absorb (Barnes & Fortney, 2004).

Assuming a spherical moon and a constant optical opacity for its rings-system,
the transit depth can be estimated as follows (Zuluaga et al., 2015),

$\delta=\frac{A_{\mathrm{Rp}}}{A_{\star}},$ (11)

where $A_{\star}$ is the area of the stellar disc and

$A_{\mathrm{Rp}}\;=\;\pi R_{\rm
p}^{2}+\pi\left[r^{2}(f_{\mathrm{o}})-r^{2}(f_{\mathrm{i}})\right]$ (12)

is the effective area projected by the planet and its rings. $r^{2}(f)$ stands
for the effective rings’ radius projected towards the star, and it is a
function of the rings’ size and opacity. For more details on $r^{2}(f)$ see
Eq. 3 in (Zuluaga et al., 2015).

Figure 8: Upper panel: Fraction between the observed and true transit depth
for a ringed exomoon as a function of the projected rings’ tilt ($\cos i_{\rm
R}$) over the observer line of sight, and different values of the rings’
normal opacity ($\tau$) for icy and dusty rings. Lower panel: rms amplitudes
of TTV (orange) and TDV (blue) signals for the system studied in section
3.2.2. The coloured dots are the true values obtained from the moon’s mass and
size used in the aforementioned simulation, while the shadows are the expected
values for massive moons.  For reference, we include the observational limits
for TESS. That is, the Near Faint Detection Limit (black-dashed line) which
corresponds to a faint star with a small planet, and the Optimal Detection
Limit (black-dotted line) related to a star 10x brighter with a planet 2.5x
larger than the previous case (Deeg & Tingley, 2017).

By applying Equations 11 and 12 to the rings’ sizes estimated in the previous
section; and for several orientations and compositions of the system (see the
upper panel of Fig. 8), we found that ringed moons crossing the stellar disc
during the flat part of the planet’s transit will induce flux changes from
$\sim 1.5$ to $3.2$ times larger than non-ringed moons. For a system composed
of a solar-type star and a Saturn-like planet with an Titan-like moon, the
planet and moon would induce a drop in stellar flux of $\delta
F_{\mathrm{p}}\sim 7000$ ppm and $\delta F_{\mathrm{m}}\sim 14$ ppm,
respectively. Over different transit measurements, this could potentially lead
to apparent Transit Radius Variations (Rodenbeck et al., 2020) of about $0.21$
per cent. In the scenario of a ringed moon, the stellar flux variation
associated solely to the moon will be up to $\delta F_{\mathrm{m-ringed}}\sim
33$ and $\sim 47$ ppm for dusty and icy rings, respectively. Such a system
would have a variation in the transit radius of $\sim 0.47$ (dusty) and $\sim
0.77$ (icy) per cent over consecutive periods. Also, moons with highly
inclined orbits could induce variations on the inferred transit’s impact
parameter over consecutive eclipses.

###  4.2 TTV and TDV of ringed exomoons

A planet with its massive satellite will orbit a common centre of mass, and
the mean orbital advance of the planet will be given by the orbital advance of
the centre of mass plus/minus the planet’s velocity relative to the centre of
mass. A distant observer can notice that the planet may either delay or
anticipate its transit due to changes in the position of the planet’s orbit
caused by its moon, producing the so-called Transit Timing Variation (or TTV).
Also, during a transit event any change in the velocity of the planet can
induce a shortening or widening of the transit duration during different
eclipses, something also known as Transit Duration Variation (or TDV). The
optimal TTVs occur when orbits are circular and co-planar for both the planet
and the moon. For edge-on orbits and from the line of view of an observer on
Earth, the barycentric root-mean-squared amplitude for TTVs
($\delta_{\mathrm{TTV}}$) is given as follows,

$\delta_{\mathrm{TTV}}=\frac{1}{2\pi}\frac{a_{\rm\mathrm{stop}}\,M_{\rm
m}}{a_{\rm p}\,M_{\rm p}}P_{\rm p}\sqrt{\frac{\Phi_{\mathrm{TTV}}}{2\pi}},$
(13)

where $\Phi_{\mathrm{TTV}}\rightarrow\pi$ is the scaling factor for co-aligned
circular orbits, as detailed in Kipping (2009b). In a similar way, the optimal
rms velocity-induced $\delta_{\mathrm{TDV}}$ for the aforementioned
configuration reads

$\delta_{\mathrm{TDV}}=\bar{\tau}\frac{a_{\rm\mathrm{stop}}\,M_{\rm m}P_{\rm
p}}{a_{\rm p}\,M_{\rm p}\;P_{\rm m}}\sqrt{\frac{\Phi_{\mathrm{TDV}}}{2\pi}},$
(14)

with $P_{\rm m}$ the orbital period of the moon and $\bar{\tau}$ the averaged
transit duration given by:

${\color[rgb]{0,0,0}\bar{\tau}=\frac{P_{\rm
p}}{\pi}\arcsin{\left[\sqrt{\frac{(R_{*}+R_{\rm p})^{2}-(b\;R_{*})^{2}}{a_{\rm
p}^{2}-b^{2}R_{*}^{2}}}\right]}},$ (15)

where $R_{*}$ is the stellar radius. We have assumed here the most favourable
case for the impact parameter (i.e. $b=0$). Also, it is worth noting that
Equations 13 and 14 do not depend on the moon’s actual or apparent size;
hence, the rings-system (whose mass is much smaller than that of the moon) has
a negligible effect on the amplitude of TTVs and TDVs  for a given semi-major
axis. However, a deep transit (such as that generated by a ringed moon) is
expected to be associated with a massive moon and to generate large
$\delta_{\mathrm{TTV}}$ and $\delta_{\mathrm{TDV}}$ amplitudes which are
directly proportional to the moon’s mass.

Although there exists a degeneracy as evidenced from equations 13 and 14, the
ratio of TDV/TTV signals (or TVs for short) can allow us to recover the mass
of moons. To predict TVs after a single observation it is also possible to
estimate the mass of moons through a simple mass-radius scaling-law assuming a
bulk composition and density for the moon. By following Zeng et al. (2019),
ringed moons with radii inferred from the transit light curves introduced in
Section 4.1 and fully made of refractory material (MgSiO3) will have estimated
masses of $4.17$ and $8.08$ times the Earth’s mass (for dusty and icy rings,
respectively), producing significant effects on TVs as shown by the middle and
lower panels of Fig. 8. There, expected TVs induced by massive moons (denoted
as exp. in the plots) are of the order of minutes and are depicted with light
green/blue for TTVs/TDVs. The lower limit of each region corresponds to rings
made of dust (brown dashed lines), whereas the upper limit represents rings
made of ices (blue dotted lines). For the parameters adopted in our
simulations (see Section 3.2.2), the lower panel of Fig. 8 shows the TVs
labelled as ‘true’ (dot markers).  Finally, in the same figure we include the
Near Faint Detection Limit (NFDL; 292 s) for the Transiting Exoplanet Survey
Satellite (TESS), which corresponds to a faint a star with a small planet, as
well as the Optimal Detection Limit (ODL; 12 s) for a star 10x brighter with a
planet 2.5x larger. For a full description of the parameters involved in each
scenario the reader can refer to table 1 in Deeg & Tingley (2017).

To summarise this section, cronomoons can produce deep imprints on the stellar
flux (i.e large TRVs) but at the same time  small/undetectable transit
secondary effects such as TTVs and TDVs. This is a result of having icy or
dusty rings-systems around them.

##  5 The system Kepler-1625b i

We studied the system Kepler-1625b i under the physical scenario proposed in
this work.  Despite Kepler-1625b i is the most accepted exomoon candidate
whose dip was first reported by Teachey & Kipping (2018), there is still no
consensus to date on the existence of such an exomoon nor its physical/orbital
properties. For instance, Kreidberg et al. (2019) disagrees with the existence
of the exomoon’s dip, although their work has been shown to be afflicted by
higher systematic noise issues (Teachey et al., 2020). Still, there exists
some convergence: Rodenbeck et al. (2020) independently recovered the dip and
the three teams have analysed the same data and agree with the TTV evidence.
Kepler-1625b i is thought to be a massive Neptune-sized moon (Teachey, Kipping
& Schmitt, 2018), but this moon might well be a less massive satellite
surrounded by an optically thick rings-system whose dynamical features could
be inferred by following the formalism developed in the previous Sections.

If we assume that a moon has actually been detected through a second stellar
flux drop in the light curve of Kepler-1625b i, its depth would allow us to
deduce the moon’s size $R_{\mathrm{m,o}}$ and consequently infer the moon’s
mass $M_{\rm m,o}$ through scaling laws (e.g. that of Chen & Kipping 2017). In
addition, the moon’s semi-major axis can also be constrained from the duration
of the dip produced by the moon. Table 1 summarises some attempts to calculate
the physical and orbital properties of the planet and the moon. Then,
complementary methods to transit photometry can be used to fully characterise
the system, and after consecutive transits the signals associated with
secondary effects should conform to those predicted by Equations 13 and 14. If
that is not the case the predicted physical parameters must be reviewed. In
particular, if large TRVs are found along with small TTVs and TDVs, we could
be in the case described in the previous sections. What would these signals be
if Kepler-1625b i were a cronomoon?

$M_{\mathrm{star}}$ | $R_{\mathrm{star}}$ | $M_{\mathrm{p}}$ | $R_{\mathrm{p}}$ | $a_{\mathrm{p}}$ | $M_{\rm m,o}$ | $R_{\mathrm{m,o}}$ | $a_{\rm m,o}$
---|---|---|---|---|---|---|---
$M_{\odot}$ | $R_{\odot}$ | $M_{\mathrm{Jup}}$ | $R_{\mathrm{Jup}}$ | au | $M_{\mathrm{Nep}}$ | $R_{\mathrm{Nep}}$ | $R_{\mathrm{p}}$
$1.079$ | $1.793$ | $10.0$ | $1.18$ | $0.87$ | $1.0$ | $1.0$ | $40.0$
This work, with $i_{\rm R}\sim 90$° |  |  |
$1.079$ | $1.793$ | $10.0$ | $1.18$ | $0.87$ | $0.18$ | $0.38$ | $40.0$
Table 1: Physical and orbital properties of the system Kepler-1625 as given by
Morton et al. (2016), Mathur et al. (2017), Teachey et al. (2018) and Teachey
& Kipping (2018). Properties compiled by Moraes & Vieira Neto (2020) as input
for their simulations.

Our first consideration is that the observed moon’s size
$R_{\mathrm{m,o}}$333Hereafter we will use the notation observed (subindex o),
true (subindex t), and expected (subindex e), to refer to the observed
parameters, the true parameters (i.e. those deduced after applying the method
presented in this paper), and the expected parameters from extra measurements
like TVs (see Section 4) obtained with the observed values, respectively. does
not correspond to the actual size of the moon but to its Roche limit (or a
projection of it onto the observer’s line of view), which is the outer radius
of the rings-system. Also, we have assumed an edge-on projection of the rings
towards the star (i.e. $i_{\rm R}\sim 90$°, see Section 4.1), so from Equation
3 the moon’s true size ($R_{\mathrm{m,t}}$) reads,

$R_{\mathrm{m,t}}=R_{\mathrm{m,o}}\left(\frac{\rho_{\mathrm{part}}}{2\;\rho_{\mathrm{m}}}\right)^{\frac{1}{3}}.$
(16)

Note that for a non edge-on transit, the true size of the moon must be
obtained by inverting Equation 12 instead of Equation 3.

Equation 16 will vary according to the composition ($\rho_{\mathrm{m}}$)
chosen for the moon. However, if regular satellites have cumulative masses
($\mathcal{O}[10^{-4}]$, Canup & Ward 2006) these will range between $\sim
0.31$ and $3.14$ the Earth’s mass, corresponding to a moon’s density between
the Earth’s density and that of newfound Super-Earths (i.e. from 5 to 8 g
cm-3). Therefore, according to Equation 16 the minimum satellite radius is
about $0.38$ times the radius of Neptune which equates to a moon’s true mass
of $\sim 3.14$ times the Earth’s mass, dubbed here as $\pi$-Earth. That being
said, within the framework of this work the exomoon candidate Kepler-1625b i
could be a $\pi$-Earth accompanied by an optically thick rings-system made of
ices and extending to distances similar in size to Neptune.

The rms values of transit secondary effects ($\delta_{\mathrm{TTV}}$ and
$\delta_{\mathrm{TTV}}$) can be estimated using Equations 13 and 14; and
compared with the values  in Kipping (2021) who computed a TTV rms of 17.08
min. When using the values mentioned in the paragraph above, we found that
$\delta_{\mathrm{TTV,e}}\sim 9.03$ min and $\delta_{\mathrm{TTV,t}}\sim 1.67$
min. Moreover, we obtained rms TDVs equal to $\delta_{\mathrm{TDV,e}}\sim
2.70$ min and $\delta_{\mathrm{TDV,t}}\sim 0.50$ min. It is worth mentioning
that no TDVs have been measured yet, so the difference between expected and
obtained values supports the hypothesis of an unseen planetary companion in a
resonant orbital configuration (see Kreidberg et al. 2019; Teachey et al.
2020).

Figure 9: Simulation results for the system Kepler-1625b i assumed as a
cronomoon with icy rings in face–on transit. $R_{\mathrm{obs}}$ is reported as
the observed size of the moon.

Figure 10:  TTVs and $\Delta$TTVs (upper panel), and TDVs and $\Delta$TDVs
(lower panel) of the system Kepler-1625b i obtained via N-body numerical
simulations. Cyan dots represent the values of a massive giant moon, whereas
red dots correspond to a cronomoon.  $\Delta$TTVs and $\Delta$TDVs stand for
the elapsed time between consecutive transits and these are compared to the
NFDL and ODL for TESS. Values inside the shaded regions are undetectable.

Adopting the parameters found under the hypothesis of a cronomoon listed on
Table 1, we performed numerical simulations using the N-body code REBOUND
(Rein & Liu, 2012) and a 15th order Gauss-Radau integrator (IAS15; Rein &
Spiegel 2015). In the absence of any other perturbing body, unlike the cases
studied in Section 3.2.2, we found that the orbits of the rings’ particles
were very regular (lower panel in Fig. 9) and did not present any orbital
resonances or gaps (upper panel in Fig. 9). Also, we extracted TTVs and TDVs
for Kepler-1625b i when orbited separately by a Neptune-like moon and a
cronomoon, represented by the blue and orange dots in Fig. 10, respectively.
Also, the difference in TTVs and TDVs between consecutive transits are shown
in Fig. 10, whose values are compared to the NFDL and ODL for TESS.
$\Delta$TTVs and $\Delta$TDVs outside these limits could be detected by TESS.
We found that both sets of TTV amplitudes are in agreement with those
calculated via equation 13, contrary to TDV signals whose values are much
larger than expected. The origin of this difference between the measured and
expected TDV values lies in the long orbital period of the planet (with its
moon) around the common centre of mass, as compared to the transit duration.
We computed the period for both the Neptune-like moon and the cronomoon as
12.6 d, with an expected transit duration of 13.5 h.

It is important to note that the parameters compiled in Table 1 were obtained
assuming the existence of a moon with no rings. Therefore, a new analysis of
the data by assuming a moon with rings could lead to very different
conclusions regarding the estimates of mass and semi-major axis of the moon,
while reproducing too the deep drop in stellar brightness originally
associated with a giant moon.

Additionally, if it is assumed in such calculations the 17-min TTV reported by
Kipping (2021) which differs from the 9-min one obtained with the parameters
of the previously mentioned table, a less massive moon (or a cronomoon) than
that reported should be located at a much more distant position than
previously thought. Although the analysis of the Kepler-1625b i system under
the framework of this research including the tidal migration and orbital
stability of the satellite will be explored in a future work, Fig. 11 shows
the orbital separations that reproduce a TTV of 17 min for the range of lunar
masses compatible with dust and icy rings. In these circumstances, the
satellite should be located between 4 and 5 times the initially reported
position of 40 times the planet’s size (i.e. the white square in Fig. 11).
These new positions are inside the planet’s Hill limit ($R_{\mathrm{Hill}}$),
but beyond its secondary Hill limit (i.e. 0.48 $R_{\mathrm{Hill}}$), so that
orbital stability would not be fully guaranteed.

Figure 11:  TTV map (equation 13) for Kepler-1625b i for a range of masses and
planet-moon separations. The white circle represents the
$\delta_{\mathrm{TTV}}$ found with data from Table 1, while the dashed yellow
line presents the mass-distance combinations that account for the TTV obtained
by Kipping (2021). The white square is the region where a cronomoon is
compatible with such a measurement. Vertical red lines are the satellite
stability limits.

##  6 Summary and discussion

We studied the likelihood for ringed satellites around close-in exoplanets to
exist and the feasilbity for them to be detectable. We have dubbed these
objects as cronomoons. These systems could explain some light curves with
atypical behaviour (give res to works with atypical behavior). We explore
proportionality between the inferred size of cronomoons and their dynamical
effects on their host planet in the context of Transit Timing Variation (TTV)
and Transit Duration Variation (TDV). We predict the TTV and TDV signals would
be much smaller than variations in the measured moon’s radius (TRV).

To analyse the dynamical stability of cronomoons we used satellite migration
models to determine their asymptotic semi-major axis, a position where
cronomoons will remain for most of their lifespan. This allowed us to
accurately assess the dynamics of their rings by fixing the gravitational
effect of the planet and the star on the rings’ particles. Additionally, this
provides information that can be used to estimate the effects that moons
produce on TTVs and TDVs, while helping reduce the degeneracy between these
two kinds of signals. We also studied the life expectancy of the constituent
particles of rings when they are subject to sublimation processes and orbital
decay, finding that such particles have a sufficient life expectancy that make
their discovery feasible by recognising morphological changes on planetary
light curves (after successive observations and depending upon mean
circumstellar distances).

By using N-body numerical simulations we found that periodic perturbations of
both the star and the planet are not significant enough to disrupt cronomoons.
Still, such perturbations can modify their rings and enrich their morphology:
the most evident effect is the opening of gaps which arise due to resonant
perturbations with the planet. As expected, those scenarios where host moons
and planets are located in close-in orbits may potentially be more harmful to
rings-systems because of the strength of gravitational interactions. However,
if cronomoons orbit more distant planets and evolve under the migration model
used in this work (Alvarado-Montes et al., 2017), they will have more stable
orbits so long as they keep within the dynamical stability limit of the
system. In contrast, scenarios where the physical parameters of the planet are
fixed can result with moons located at any distance (Dobos et al., 2021),
including short-period orbits where rings can rapidly be dissipated by
perturbing mechanisms.

The viability of cronomoons as a physical system strongly depends on the
environment where they develop and dwell in, as described in the different
formation scenarios of Section 2. For example, mergers between sibling
satellites or close encounters with exocomets could populate exomoons with
rings, and even result in different types of circumplanetary and
circumsatellital structures that can be classified according to their
composition and morphology. Collisions between moons could generate co-planar
dusty rings, whereas encounters or flybys with exocomets could generate icy
rings in random orbital inclinations (taking as a proxy the orbital properties
of the comets in the Solar System444See https://ssd.jpl.nasa.gov/?sb_elem.).
Additionally, any deviation from coplanarity between the orbits of moons and
their rings could be interpreted as a signpost of a recent tidal disruption
event that delivered fresh material around the moon to form the rings, given
that dynamical evolution tends to a fast coplanarization of this type of
systems (Sucerquia et al., 2017). Interestingly, studying transit secondary
effects of systems with tilted configurations would allow us to follow them up
and monitor them over very long baselines.

Beyond their dynamical stability, cronomoons can significantly alter first-
and second-order effects on planetary light curves. We showed that they can
exhibit transits with large TRV but dimmed TTV and TDV signals (contrary to
what is expected for massive systems). In addition, the transit geometry of
cronomoons at their ingress and egress stages, as well as their rings’ gaps
may reveal the presence of these objects, but it was also found that large
transit depths generated by circumsatellital rings can be related to moderate
masses that do not necessarily perturb the dynamics of host planets.

Among the ten current claims for exomoon candidates (Ben-Jaffel & Ballester,
2014; Bennett et al., 2014; Teachey & Kipping, 2018; Oza et al., 2019;
Kipping, 2020; Quarles et al., 2020; Fox & Wiegert, 2021), we chose to study
the system Kepler-1625b i (Section 5) – the most convincing exomooon candidate
so far – and compare its properties to those of a stable cronomoon. Despite an
accurate fitting model for the light curve of this system was not used, we
tested its dynamical evolution and detectable features by assuming a face-on
transit: this implies that the sum of the rings’ tilt and the inclination of
the system is $\sim 90$°, helping us illustrate another phenomenological
explanation of the observations by predicting small values for TTVs and TDVs.

In this work we mainly studied the case of rings around moons which orbit
close-in planets, as this kind of systems could be easier to detect. However,
giant planets located in farther circumstellar orbits could also bear ringed
moons whose existence is possible via different mechanisms like capture of
moons or destructive orbital resonances. Events such as the Late Heavy
Bombardment (LHB) as detailed by the Nice Model of the Solar System, or the
recent impact of the Shoemaker–Levy 9 comet on Jupiter (see e.g. Zahnle & Mac
Low 1994) could give rise to fortuitous collisions that are, in fact, frequent
in the mature stages of planetary systems. Also, despite impacts of comets
against the Earth during the K/T mass extinction event (Siraj & Loeb, 2021)
are still debated due to a lack of well-established evidence in the Solar
System (Desch et al., 2021), this scenario might be somewhat frequent
especially during the early stages of planetary evolution and might lead to
the formation of rings around moons.

The fact that no cronomoons have been detected in the Solar System could mean
that their lifespan was shorted due to the dispersion of particles in their
orbital plane, produced by the presence of other planets and moons (a scenario
we will study in a future work). In addition, the small density of most moons
in the Solar System could have produced small cross sectional areas for
cronomoon formation after the passage of an outsider object (e.g. during the
LHB event), so that capturing material to form a cronomoon might have been
feasible but also unlikely. On the other hand, cronomoons formed after moon-
moon collision events might be more likely and frequent to form in the early
stages of the Solar System, but these objects might have already been depleted
via the aforementioned processes.

Adding to the scenarios mentioned in Section 2, a feasible avenue for the
formation of cronomoons in the current Solar System could be the capture of
objects already bearing rings. For instance, cronomoons might be formed if
giant planets trapped ring-bearing bodies such as the Centaurus Chariklo and
Chiron, or the trans-Neptunian object Haumea given that their rings-systems
survive to the gravitational stresses produced by the capture process. Future
geological studies applied on surfaces of moons searching for ancient impacts
may detect a disparity between their natural composition and the particles
that were deposited on their surfaces by older rings-systems. These may
possibly form ridge-like structures in a similar way to the process of ring’s
particles falling on a planet’s atmosphere, evidenced by the Cassini’s Grand
Finale phase as described by Waite et al. (2018) and references therein.

To summarise our work, we studied a feasible scenario whereby rings around
moons might flourish and thrive. We dubbed these objects cronomoons after
their similarity with Cronus (the Greek name for Saturn), and after Chronos
(the epitome of time), following the effects on light curves in terms of
transit timing and duration. Cronomoons, in addition to endowing us with a new
interpretation of light curves, give us a glimpse about the exotic
extraterrestrial landscapes that could take place around exoplanets, and even,
in our own Solar System.

## Acknowledgements

The authors thank the referee, David Kipping, for the helpful feedback and
positive remarks, which allowed us to improve the quality of this research.
Mario Sucerquia acknowledges support by Agencia Nacional de Investigación y
Desarrollo (ANID) through FONDECYT postdoctoral 3210605. Mario Sucerquia,
Amelia Bayo, Johan Olofsson, Jorge Cuadra, and Matías Montesinos thank ANID –
Millennium Science Initiative Program – NCN19_171. Jaime A. Alvarado-Montes
acknowledges funding support from Macquarie University through the
International Macquarie University Research Excellence Scholarship (‘iMQRES’).
This project has received funding from the European Union’s Horizon 2020
research and innovation programme under the Marie Skłodowska-Curie grant
agreement Nº 210021. Jorge I. Zuluaga is funded by Vicerrectoría de Docencia,
UdeA. Cristian Giuppone acknowledge the computational resources from CCAD –
UNC (https://ccad.unc.edu.ar), which are part of SNCAD – MinCyT, República
Argentina. The authors of this paper call for respect of fundamental rights in
order to achieve equity, social justice, and peace needed in Colombia.
This work is dedicated to Martín.

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Akinsanmi et al. (2020) Akinsanmi B., Santos N. C., Faria J. P., Oshagh M., Barros S. C. C., Santerne A., Charnoz S., 2020, A&A, 635, L8
  * Alvarado-Montes & García-Carmona (2019) Alvarado-Montes J. A., García-Carmona C., 2019, MNRAS, 486, 3963
  * Alvarado-Montes et al. (2017) Alvarado-Montes J. A., Zuluaga J. I., Sucerquia M., 2017, MNRAS, 471, 3019
  * Alvarado-Montes et al. (2021) Alvarado-Montes J. A., Sucerquia M., García-Carmona C., Zuluaga J. I., Spitler L., Schwab C., 2021, MNRAS, 506, 2247
  * Andreas (2007) Andreas E. L., 2007, Icarus, 186, 24
  * Arkhypov et al. (2021) Arkhypov O. V., Khodachenko M. L., Hanslmeier A., 2021, A&A, 646, A136
  * Barnes & Fortney (2004) Barnes J. W., Fortney J. J., 2004, ApJ, 616, 1193
  * Barnes & O’Brien (2002) Barnes J. W., O’Brien D. P., 2002, ApJ, 575, 1087
  * Ben-Jaffel & Ballester (2014) Ben-Jaffel L., Ballester G. E., 2014, ApJ, 785, L30
  * Bennett et al. (2014) Bennett D. P., et al., 2014, ApJ, 785, 155
  * Bohren & Albrecht (1998) Bohren C. F., Albrecht B. A., 1998, Atmospheric thermodynamics
  * Braga-Ribas et al. (2014) Braga-Ribas F., et al., 2014, Nature, 508, 72
  * Brilliantov et al. (2015) Brilliantov N., Krapivsky P. L., Bodrova A., Spahn F., Hayakawa H., Stadnichuk V., Schmidt J., 2015, Proceedings of the National Academy of Sciences, 112, 9536
  * Canup (2010) Canup R. M., 2010, Nature, 468, 943
  * Canup & Ward (2006) Canup R. M., Ward W. R., 2006, Nature, 441, 834
  * Charnoz et al. (2009) Charnoz S., Morbidelli A., Dones L., Salmon J., 2009, Icarus, 199, 413
  * Chen & Kipping (2017) Chen J., Kipping D., 2017, ApJ, 834, 17
  * Crida & Charnoz (2012) Crida A., Charnoz S., 2012, Science, 338, 1196
  * Crida et al. (2019) Crida A., Charnoz S., Hsu H.-W., Dones L., 2019, Nature Astronomy, 3, 967
  * Deeg & Tingley (2017) Deeg H. J., Tingley B., 2017, A&A, 599, A93
  * Desch et al. (2021) Desch S., Jackson A., Noviello J., Anbar A., 2021, Astronomy and Geophysics, 62, 3.34
  * Dobos et al. (2021) Dobos V., Charnoz S., Pál A., Roque-Bernard A., Szabó G. M., 2021, arXiv e-prints, p. arXiv:2105.12040
  * Domingos et al. (2006) Domingos R. C., Winter O. C., Yokoyama T., 2006, MNRAS, 373, 1227
  * Estermann (1955) Estermann I., 1955, Gases at low densities. In Rossini, F.D., Thermodynamics and Physics of Matter. Princeton University Press, Princeton NJ
  * Forgan (2018) Forgan D., 2018, Research Notes of the American Astronomical Society, 2, 191
  * Fox & Wiegert (2021) Fox C., Wiegert P., 2021, MNRAS, 501, 2378
  * Guenel et al. (2014) Guenel M., Mathis S., Remus F., 2014, A&A, 566, L9
  * Hedman (2015) Hedman M. M., 2015, ApJ, 801, L33
  * Hedman & Nicholson (2016) Hedman M. M., Nicholson P. D., 2016, Icarus, 279, 109
  * Heller & Barnes (2013) Heller R., Barnes R., 2013, Astrobiology, 13, 18
  * Heller et al. (2014) Heller R., et al., 2014, Astrobiology, 14, 798
  * Heller et al. (2019) Heller R., Rodenbeck K., Bruno G., 2019, A&A, 624, A95
  * Hesselbrock & Minton (2017) Hesselbrock A. J., Minton D. A., 2017, Nature Geoscience, 10, 266
  * Hobbs (1974) Hobbs P. V., 1974, Ice physics
  * Hyodo & Charnoz (2017) Hyodo R., Charnoz S., 2017, AJ, 154, 34
  * Isella et al. (2019) Isella A., Benisty M., Teague R., Bae J., Keppler M., Facchini S., Pérez L., 2019, ApJ, 879, L25
  * Kipping (2009a) Kipping D. M., 2009a, MNRAS, 392, 181
  * Kipping (2009b) Kipping D. M., 2009b, MNRAS, 396, 1797
  * Kipping (2020) Kipping D., 2020, ApJ, 900, L44
  * Kipping (2021) Kipping D., 2021, MNRAS, 500, 1851
  * Kipping & Teachey (2020) Kipping D., Teachey A., 2020, Serbian Astronomical Journal, 201, 25
  * Kipping et al. (2009) Kipping D. M., Fossey S. J., Campanella G., 2009, MNRAS, 400, 398
  * Kipping et al. (2012) Kipping D. M., Dunn W. R., Jasinski J. M., Manthri V. P., 2012, MNRAS, 421, 1166
  * Kipping et al. (2013) Kipping D. M., Forgan D., Hartman J., Nesvorný D., Bakos G. Á., Schmitt A., Buchhave L., 2013, ApJ, 777, 134
  * Kipping et al. (2014) Kipping D. M., Nesvorný D., Buchhave L. A., Hartman J., Bakos G. Á., Schmitt A. R., 2014, ApJ, 784, 28
  * Kollmeier & Raymond (2019) Kollmeier J. A., Raymond S. N., 2019, MNRAS, 483, L80
  * Kozai (1959) Kozai Y., 1959, AJ, 64, 367
  * Kreidberg et al. (2019) Kreidberg L., Luger R., Bedell M., 2019, ApJ, 877, L15
  * Lidov (1962) Lidov M. L., 1962, Planet. Space Sci., 9, 719
  * Martínez-Rodríguez et al. (2019) Martínez-Rodríguez H., Caballero J. A., Cifuentes C., Piro A. L., Barnes R., 2019, ApJ, 887, 261
  * Mathis (2015) Mathis S., 2015, in Martins F., Boissier S., Buat V., Cambrésy L., Petit P., eds, SF2A-2015: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics. pp 283–288 (arXiv:1510.05639)
  * Mathur et al. (2017) Mathur S., et al., 2017, ApJS, 229, 30
  * Moraes & Vieira Neto (2020) Moraes R. A., Vieira Neto E., 2020, MNRAS, 495, 3763
  * Morton et al. (2016) Morton T. D., Bryson S. T., Coughlin J. L., Rowe J. F., Ravichandran G., Petigura E. A., Haas M. R., Batalha N. M., 2016, ApJ, 822, 86
  * Murray & Dermott (2000) Murray C. D., Dermott S. F., 2000, Solar System Dynamics. Cambridge
  * Namouni (2010) Namouni F., 2010, ApJ, 719, L145
  * Ogilvie (2013) Ogilvie G. I., 2013, MNRAS, 429, 613
  * Ortiz et al. (2015) Ortiz J. L., et al., 2015, A&A, 576, A18
  * Ortiz et al. (2017) Ortiz J. L., et al., 2017, Nature, 550, 219
  * Oza et al. (2019) Oza A. V., et al., 2019, ApJ, 885, 168
  * Peters & Turner (2013) Peters M. A., Turner E. L., 2013, ApJ, 769, 98
  * Piro & Vissapragada (2020) Piro A. L., Vissapragada S., 2020, AJ, 159, 131
  * Pollack (1975) Pollack J. B., 1975, Space Sci. Rev., 18, 3
  * Pollack et al. (1973) Pollack J. B., Summers A., Baldwin B., 1973, Icarus, 20, 263
  * Quarles et al. (2020) Quarles B., Li G., Rosario-Franco M., 2020, ApJ, 902, L20
  * Quarles et al. (2021) Quarles B., Eggl S., Rosario-Franco M., Li G., 2021, AJ, 162, 58
  * Rappaport et al. (2018) Rappaport S., et al., 2018, MNRAS, 474, 1453
  * Rein & Liu (2012) Rein H., Liu S. F., 2012, A&A, 537, A128
  * Rein & Spiegel (2015) Rein H., Spiegel D. S., 2015, MNRAS, 446, 1424
  * Roche (1849) Roche E., 1849, Acad. Montpellier, 1, 243
  * Rodenbeck et al. (2020) Rodenbeck K., Heller R., Gizon L., 2020, A&A, 638, A43
  * Rosario-Franco et al. (2020) Rosario-Franco M., Quarles B., Musielak Z. E., Cuntz M., 2020, AJ, 159, 260
  * Sasaki et al. (2012) Sasaki T., Barnes J. W., O’Brien D. P., 2012, ApJ, 754, 51
  * Schlichting & Chang (2011) Schlichting H. E., Chang P., 2011, ApJ, 734, 117
  * Seager & Mallén-Ornelas (2003) Seager S., Mallén-Ornelas G., 2003, ApJ, 585, 1038
  * Siraj & Loeb (2021) Siraj A., Loeb A., 2021, Scientific Reports, 11, 3803
  * Spalding et al. (2016) Spalding C., Batygin K., Adams F. C., 2016, ApJ, 817, 18
  * Strøm et al. (2020) Strøm P. A., et al., 2020, PASP, 132, 101001
  * Sucerquia et al. (2017) Sucerquia M., Alvarado-Montes J. A., Ramírez V., Zuluaga J. I., 2017, MNRAS, 472, L120
  * Sucerquia et al. (2019) Sucerquia M., Alvarado-Montes J. A., Zuluaga J. I., Cuello N., Giuppone C., 2019, MNRAS, 489, 2313
  * Sucerquia et al. (2020a) Sucerquia M., Ramírez V., Alvarado-Montes J. A., Zuluaga J. I., 2020a, MNRAS, 492, 3499
  * Sucerquia et al. (2020b) Sucerquia M., Alvarado-Montes J. A., Zuluaga J. I., Montesinos M., Bayo A., 2020b, MNRAS, 496, L85
  * Teachey & Kipping (2018) Teachey A., Kipping D. M., 2018, Science Advances, 4, eaav1784
  * Teachey et al. (2018) Teachey A., Kipping D. M., Schmitt A. R., 2018, AJ, 155, 36
  * Teachey et al. (2020) Teachey A., Kipping D., Burke C. J., Angus R., Howard A. W., 2020, AJ, 159, 142
  * Tsiganis et al. (2005) Tsiganis K., Gomes R., Morbidelli A., Levison H. F., 2005, Nature, 435, 459
  * Vanderburg et al. (2018) Vanderburg A., Rappaport S. A., Mayo A. W., 2018, AJ, 156, 184
  * Vincent et al. (2016) Vincent J. B., et al., 2016, A&A, 587, A14
  * Vogel (1979) Vogel H., 1979, Mathematical Biosciences, 44, 179
  * Waite et al. (2018) Waite J. H., et al., 2018, Science, 362, aat2382
  * Zahnle & Mac Low (1994) Zahnle K., Mac Low M.-M., 1994, Icarus, 108, 1
  * Zeng et al. (2019) Zeng L., et al., 2019, Proceedings of the National Academy of Science, 116, 9723
  * Zuluaga et al. (2015) Zuluaga J. I., Kipping D. M., Sucerquia M., Alvarado J. A., 2015, ApJ, 803, L14