# Water Evolution & Inventories of Super-Earths Orbiting Late M-Dwarfs

Keavin Moore  Department of Earth & Planetary Sciences, McGill University,
3450 rue University, Montréal, QC, H3A 0E8, Canada  Benjamin David  Department
of Physics, McGill University, 3600 rue University, Montréal, QC, H3A 2T8,
Canada  Albert Yian Zhang  Department of Physics, McGill University, 3600 rue
University, Montréal, QC, H3A 2T8, Canada  Nicolas B. Cowan  Department of
Earth & Planetary Sciences, McGill University, 3450 rue University, Montréal,
QC, H3A 0E8, Canada  Department of Physics, McGill University, 3600 rue
University, Montréal, QC, H3A 2T8, Canada

###### Abstract

Super-Earths orbiting M-dwarf stars may be the most common habitable planets
in the Universe. However, their habitability is threatened by intense
irradiation from their host stars, which drives the escape of water to space
and can lead to surface desiccation. We present simulation results of a box
model of water cycling between interior and atmosphere and loss to space, for
terrestrial planets of mass 1–8 $M_{\oplus}$ orbiting in the habitable zone of
a late M-dwarf. Energy-limited loss decreases with planetary mass, while
diffusion-limited loss increases with mass. Depending on where it orbits in
the habitable zone, a 1 $M_{\oplus}$ planet that starts with 3–8 Earth Oceans
can end up with an Earth-like surface of oceans and exposed continents; for an
8 $M_{\oplus}$ super-Earth, that range is 3–12 Earth Oceans. Planets
initialized with more water end up as waterworlds with no exposed continents,
while planets that start with less water have desiccated surfaces by 5 Gyr.
Since the mantles of terrestrial planets can hold much more water than is
currently present in Earth’s atmosphere, none of our simulations result in
Dune planets — such planets may be less common than previously thought.
Further, more water becomes sequestered within the mantle for larger planets.
A super-Earth at the inner edge of the habitable zone tends to end up as
either a waterworld or with a desiccated surface; only a narrow range of
initial water inventory yields an Earth-like surface.

planets and satellites: atmospheres – planets and satellites: interiors –
planets and satellites: tectonics – planets and satellites: terrestrial
planets – planets and satellites: oceans – stars: low-mass

##  1 Introduction

###  1.1 Habitability of Earth-like Planets

Rocky planets around M-dwarfs have become prime targets for astrobiology.
These stars are the most abundant in the Galaxy and are known to host many
rocky planets (e.g., Dressing & Charbonneau 2015; Sabotta et al. 2021), but
M-dwarfs are significantly more active than FGK-dwarfs (Scalo et al., 2007;
Shields et al., 2016), emitting abundant X-ray and extreme ultraviolet (EUV)
radiation, collectively known as XUV radiation. As a result, the habitability
of rocky planets orbiting M-dwarfs is at risk due to higher irradiation during
the extended early evolutionary periods of M-dwarfs (for a review, see Shields
et al. 2016).

The significant radiation impinging on the atmospheres of these rocky planets
can lead to rapid and significant loss of water to space, or even surface
desiccation (e.g., Wordsworth & Pierrehumbert 2013, 2014; Ramirez &
Kaltenegger 2014; Luger & Barnes 2015; Tian & Ida 2015). Two previous studies
used box models to investigate the fate of planetary surface water as a proxy
for habitability (Moore & Cowan, 2020; Moore, Cowan, & Boukaré, 2023). These
models include atmospheric loss to space, and track water evolution between
planetary reservoirs during an early magma ocean stage and a plate-tectonics-
driven deep-water cycling phase. These studies concluded that sequestering
water in the solid mantle or a magma ocean can improve the habitability
prospects of Earth-mass planets by reducing loss of water to space.

###  1.2 Habitability of Super-Earths

Rocky planets larger than the Earth are expected to be more common than Earth-
size planets around M-dwarfs (e.g., Dressing & Charbonneau 2015; Mulders et
al. 2015; Hsu et al. 2020; Chachan & Lee 2023) and as such, their water
evolution must be investigated (e.g., Kruijver et al. 2021). These “super-
Earths” have masses in the range $1~{}M_{\oplus}<M_{\mathrm{p}}\lesssim
10~{}M_{\oplus}$, and super-Earths orbiting low-mass stars may be the most
common habitable planets. Regardless of their intrinsic frequency, super-
Earths are easier to discover and characterize than smaller terrestrial
planets. Furthermore, terrestrial planets forming around low-mass stars may
initially contain 10–30% water by mass, depending on the details of their
migration through the protoplanetary disk (Alibert & Benz, 2017).

How does water evolution and partitioning operate on massive terrestrial
planets? The manuscript continues as follows. In Section 2, we summarize the
box model of Moore, Cowan, & Boukaré (2023) before emphasizing mass dependence
and our expectations. Section 3 then presents temporal water evolution and
parameter space results for various terrestrial planets. These results are
then discussed in Section 4 along with the nuances and exclusions of our
super-Earth box model, before we summarize our results and conclude in Section
5.

##  2 Methods

The model used in this study is a modified version of that in Moore, Cowan, &
Boukaré (2023), without a persistent basal magma ocean; we instead test
various magma ocean water partition coefficients across different planetary
masses, and a basal magma ocean would be represented by a higher partition
coefficient. Moore, Cowan, & Boukaré (2023) showed that a short-lived basal
magma ocean could benefit surface habitability, which could help rapidly-
desiccated planets if substantial water becomes sequestered in the mantle;
however, a long-lived basal magma ocean leads to the majority of water
becoming trapped in the mantle. Basal magma oceans are probably more likely
for super-Earths (C.-É Boukaré, priv. comm.), which we discuss in Section 4.

We direct the reader to Moore, Cowan, & Boukaré (2023) for a thorough
explanation of our model of water evolution, but we note key assumptions
below. We adopt an XUV saturation timescale for the host star of
$t_{\mathrm{sat}}=$ 1 Gyr, and use the stellar tracks of Baraffe et al. (2015)
to determine the stellar luminosity at a given age, offsetting the track by 5
Myr to account for the lag between star and planet formation. Significant
changes to the model and the mass dependence in certain equations are outlined
below. The overall planetary evolution pathway remains the same as our
previous study (Moore, Cowan, & Boukaré, 2023): the planet begins with a magma
ocean, extending from surface to core, and which completely solidifies by the
end of the runaway greenhouse phase. The runaway greenhouse phase is
concurrent with the magma ocean phase. When stellar irradiation is greater
than $325$ W/m2 (Turbet et al., 2021), all water exists in a steam atmosphere;
when the irradiation falls below this limit, we set the surface temperature
$T_{\mathrm{surf}}=293.15$ K $=20^{\circ}$ C for simplicity. Following surface
solidification, plate tectonics begins, initiating the deep-water cycle, i.e.,
exchange of water between interior and surface/atmosphere through degassing
(interior to surface) and regassing (surface to interior) mediated by seafloor
pressure and mantle temperature. Our deep-water cycling parameterization is
based on an Earth-twin with plate tectonics (Moore, Cowan, & Boukaré, 2023).

Water loss to space is energy-limited (Watson et al., 1981) when the planet is
in a runaway greenhouse; during this time, all water is readily available in
the atmosphere to be photodissociated and lost to space. Water loss to space
becomes diffusion-limited (Walker, 1977) once the planet exits the runaway
greenhouse and the magma ocean solidifies; at this point, the more modest
surface temperatures mean that water loss is limited by its vertical diffusion
through the atmosphere, towards the exobase, above which it may be lost. Note
that, although the diffusion limit is typically not valid for steam
atmospheres, our atmosphere contains water, hydrogen, and oxygen, and hence
the diffusion limit corresponds to hydrogen diffusing through an oxygen
background.

###  2.1 Expected Mass Dependence

We use the scaling relations of Valencia et al. (2006) to calculate planetary
radius, surface area, and mantle volume as a function of mass. These scalings
do not explicitly account for molten phases, and thus cannot fully represent
the magma ocean stage (Dorn & Lichtenberg, 2021); however, we neglect the
contribution of surface & atmospheric water on planetary radius within our
model, although we do explore this assumption in Section 4. The relation of
Valencia et al. (2006) between planetary radius, $R_{\mathrm{p}}$, and mass,
$M_{\mathrm{p}}$ is:

$R_{\mathrm{p}}=R_{\oplus}\left(\frac{M_{\mathrm{p}}}{M_{\oplus}}\right)^{0.27},$
(1)

where $R_{\oplus}$ and $M_{\oplus}$ are the radius and mass of the Earth,
respectively.

We can substitute this relation into the equations for water loss to space.
The energy-limited escape rate (Watson et al., 1981) depends on the XUV energy
required to overcome the gravitational potential of the planet:

$\dot{M}_{\mathrm{energy-limited}}=\frac{\epsilon_{\mathrm{XUV}}\pi
F_{\mathrm{XUV}}R_{\mathrm{p}}R_{\mathrm{XUV}}^{2}}{GM_{\mathrm{p}}K_{\mathrm{tide}}}.$
(2)

Here, $F_{\mathrm{XUV}}$ is the XUV flux at the orbital distance of the
planet, and we hold the orbital distance fixed throughout each simulation for
simplicity; we set $K_{\mathrm{p}}=1$, meaning no tidal contribution. The XUV
absorption efficiency, $\epsilon_{\mathrm{XUV}}$, is typically 1–10% (e.g.,
Ercolano & Clarke 2010; Barnes et al. 2016; Lopez 2017), and we hold this
value fixed at $\epsilon_{\mathrm{XUV}}=0.1$; we previously discussed the
impact of this assumption in Moore, Cowan, & Boukaré (2023). Making the
planetary mass dependence explicit, the energy-limited escape equation
becomes:

$\dot{M}_{\mathrm{energy-limited}}\propto
R_{\mathrm{p}}^{2}\frac{R_{\mathrm{p}}}{M_{\mathrm{p}}}=\frac{R_{\mathrm{p}}^{3}}{M_{\mathrm{p}}}=M_{\mathrm{p}}^{-0.19}.$
(3)

Here, we have approximated the surface area at which the XUV energy is
deposited, $R_{\mathrm{XUV}}^{2}$, as $R_{\mathrm{p}}^{2}$, while the
additional $R_{\mathrm{p}}/M_{\mathrm{p}}$ corresponds to the gravitational
potential; see Section 4 for a sensitivity analysis of our XUV deposition
radius, which may be much greater for hot, extended atmospheres, especially
during the magma ocean/runaway greenhouse phase. This equation implies that
the energy-limited escape rate will decrease with increasing planetary mass.

The diffusion-limited escape rate (Walker, 1977) can be determined using the
following equation:

$\dot{M}_{\mathrm{diffusion-limited}}=m_{\mathrm{H}}\pi
R_{\mathrm{p}}^{2}\frac{bg(m_{\mathrm{O}}-m_{\mathrm{H}})}{k_{\mathrm{B}}T_{\mathrm{therm}}(1+X_{\mathrm{O}}/X_{\mathrm{H}})}.$
(4)

Within this equation, $m_{\mathrm{H}}$ and $m_{\mathrm{O}}$ are the atomic
masses of hydrogen and oxygen, respectively, the binary diffusion coefficient
between hydrogen and oxygen is $b=4.8\times
10^{17}(T_{\mathrm{therm}})^{0.75}$ cm-1 s-1, $g$ is the planetary surface
gravity, and $X_{\mathrm{O}}/X_{\mathrm{H}}=1/2$ arises since
photodissociation of water releases two hydrogen atoms for every one oxygen
atom. We further assume that oxygen either escapes with hydrogen, or reacts
away with the surface, such that it does not build up in the atmosphere. We
adopt an average thermospheric temperature of $T_{\mathrm{therm}}=400$ K
(Hunten et al., 1987), following Luger & Barnes (2015) and Moore, Cowan, &
Boukaré (2023). We discuss the sensitivity of our results to a constant
$T_{\mathrm{therm}}$ in Section 4.

Again using the Valencia et al. (2006) scaling relations, the mass dependence
of diffusion-limited escape can be clarified:

$\dot{M}_{\mathrm{diffusion-limited}}\propto R_{\mathrm{p}}^{2}g\propto
R_{\mathrm{p}}^{2}\left(\frac{M_{\mathrm{p}}}{R_{\mathrm{p}}^{2}}\right)=M_{\mathrm{p}},$
(5)

so diffusion-limited escape will increase with planetary mass. The first
factor of $R_{\mathrm{p}}^{2}$ again corresponds to planetary surface area. We
calculate the expected water loss over 5 Gyr for each planetary mass in Table
1.

Planetary | Location Within | Orbital | Orbital | Expected Water Loss [Earth Oceans]
---|---|---|---|---
Mass | Habitable Zone (HZ) | Distance [AU] | Period [days] | EL during RG | DL rest of sim | Total | Maximum
1 $M_{\oplus}$ | Hot, Inner Edge | 0.025 | 4.81 | 7.72 | 2.82 | 10.5 | 19.7
| Middle (Mid) | 0.037 | 8.67 | 2.18 | 2.95 | 5.13 | 8.92
| Cold, Outer Edge | 0.050 | 13.6 | 0.83 | 2.99 | 3.83 | 5.66
2 $M_{\oplus}$ | Hot, Inner Edge | 0.025 | 4.81 | 6.77 | 5.64 | 12.4 | 17.8
| Middle (Mid) | 0.037 | 8.67 | 1.91 | 5.91 | 7.82 | 9.71
| Cold, Outer Edge | 0.050 | 13.6 | 0.73 | 5.98 | 6.71 | 7.42
4 $M_{\oplus}$ | Hot, Inner Edge | 0.025 | 4.81 | 5.93 | 11.3 | 17.2 | 19.4
| Middle (Mid) | 0.037 | 8.67 | 1.67 | 11.8 | 13.5 | 13.9
| Cold, Outer Edge | 0.050 | 13.6 | 0.64 | 12.0 | 12.6 | 12.8
8 $M_{\oplus}$ | Hot, Inner Edge | 0.025 | 4.81 | 5.20 | 22.5 | 27.8 | 27.8
| Middle (Mid) | 0.037 | 8.67 | 1.47 | 23.6 | 25.1 | 25.1
| Cold, Outer Edge | 0.050 | 13.6 | 0.56 | 23.9 | 24.5 | 24.5
Table 1: Expected amount of water that will be lost as a function of planetary
mass, $M_{\mathrm{p}}$, and orbital distance, neglecting a magma ocean and
deep-water cycle. Each HZ orbital distance fixed throughout a given
simulation, and determined using the calculator of Kopparapu et al. (2013) at
4.5 Gyr, comparable to the age of the Solar System; the Inner Edge corresponds
to the runaway greenhouse limit, the Outer Edge corresponds to freeze-out of
carbon dioxide, and the Middle is taken to be the average orbital distance
between these two limits. Here, “EL” refers to energy-limited loss, which is
expected to decrease with $M_{\mathrm{p}}$; “RG” is the runaway greenhouse;
and “DL” is diffusion-limited loss, expected to increase with
$M_{\mathrm{p}}$. The total expected loss is the sum of the previous two
columns, while the maximum expected loss is calculated using the maximum of EL
and DL loss at each timestep. In practice, these planets lose much less water
due to the cycling of water between interior, surface, and atmosphere.

###  2.2 Model Inputs & Constraints

Table 2 shows the parameter space explored to investigate the effect of
planetary mass on water retention and water partitioning for the first 5 Gyr
of the lifetime of super-Earths orbiting a late M-dwarf; we choose 5 Gyr as it
is roughly the age of the Solar System, allowing direct comparison with the
habitable Earth, but also to capture the early, intensely-irradiated lifetime
of the planet when loss to space would be energy-limited and presumably most
substantial, if not for a magma ocean, since the pre-main sequence of an
M-dwarf may persist for Gyrs. Four planetary masses are investigated:
$M_{\mathrm{p}}=1~{}M_{\oplus}$, $2~{}M_{\oplus}$, $4~{}M_{\oplus}$, and
$8~{}M_{\oplus}$. Four initial water mass fractions are explored, ranging from
$4.7\times 10^{-4}$ to $1.9\times 10^{-3}$; these represent the end-member
scenario that water delivery is a consequence of accretion, and hence, the
initial water inventory will scale approximately with planetary mass.
Terrestrial planets around low-mass stars may form with more than 10% water by
mass (Alibert & Benz, 2017), but we find that for water mass fractions above
$2\times 10^{-3}$, all our simulated planets end as waterworlds.

We also run a set of simulations that begin with the same initial water
inventory independent of mass, 2, 4, 6, and 8 Earth Oceans; since planets
forming in-situ in the hot, volatile-depleted habitable zone around M-dwarf
stars form dry (e.g., Schneeberger et al. 2024), planetary water inventory is
instead delivered later by e.g., comets, and hence planets will begin with
roughly the same water inventory regardless of planetary mass. Three fixed
orbital distances within the habitable zone (HZ) are simulated — Inner HZ, Mid
HZ, and Outer HZ — calculated using the HZ relations of Kopparapu et al.
(2013) for an M8 host star at 4.5 Gyr, i.e., approximately the age of the
Earth. The water saturation limit of the magma ocean, which controls the
timing of atmospheric degassing, is held constant at $C_{\mathrm{sat}}=0.01$
due to its thorough exploration in Moore, Cowan, & Boukaré (2023); this means
that an atmosphere is only degassed once the magma ocean becomes
supersaturated with water.

Name | Parameter | Values Tested
---|---|---
Planetary mass |  $M_{\mathrm{p}}$ [$M_{\oplus}$] | 1, 2, 4, 8
Initial water mass fraction | $(M_{\mathrm{init}}/M_{\mathrm{p}})$ |  $4.7\times 10^{-4}$, $9.4\times 10^{-4}$, $1.4\times 10^{-3}$, $1.9\times 10^{-3}$
Initial water inventory |  $M_{\mathrm{init}}$ [Earth Oceans] | 2, 4, 6, 8
Magma ocean water partition coefficient | $D$ | 0.001, 0.01, 0.1
Orbital distance/location within HZ |  $a_{\mathrm{orb}}$ [AU] | 0.025 (Inner HZ), 0.037 (Mid HZ), 0.050 (Outer HZ)
Host stellar type | — | M8
Table 2: Parameter space explored in this study. Planetary mass is expressed
in units of Earth masses, where $M_{\oplus}\approx 5.97\times 10^{24}$ kg.
Initial water mass fraction is based on [2, 4, 6, 8] Earth Oceans on a 1
$M_{\oplus}$ planet, where 1 Earth Ocean $\approx 1.4\times 10^{21}$ kg. We
run two sets of simulations: in the first, planets begin with the same initial
water mass fraction, while the second initializes planets with the same mass
of water, measured in Earth Oceans. These roughly represent end-members of
possible water delivery scenarios.

The mantle water capacity will increase as planetary mass increases. As such,
we adopt the 12 Earth Ocean limit on mantle water capacity for a 1
$M_{\oplus}$ planet from Moore & Cowan (2020) and assume the mantle water
capacity scales as $(M_{\mathrm{p}}/M_{\oplus})\times 12$ Earth Oceans,
consistent with Cowan & Abbot (2014); this is a conservative assumption given
the high water capacity of post-perovskite (Townsend et al., 2016). Further,
the different surface water regimes — waterworld, Earth-like, Dune planet, and
desiccated — will also vary with planetary mass, as shown in Fig. 1. A
waterworld will have a completely inundated surface with no exposed
continents, while an Earth-like planet has both exposed continents & surface
oceans. A Dune planet would have up to ${\sim}$1% of an Earth Ocean on its
surface, while a desiccated planet has been stripped of all of its surface
water, either through loss to space or regassing to the interior. These
regimes can be delineated by the mass of water on the planet’s surface,
$W_{\mathrm{s}}$. A planet is a waterworld if $W_{\mathrm{s}}\geq
3.7(M_{\mathrm{p}}/M_{\oplus})^{0.08}$ Earth Oceans (Cowan & Abbot, 2014),
Earth-like if
$0.01(M_{\mathrm{p}}/M_{\oplus})^{0.54}<W_{\mathrm{s}}<3.7(M_{\mathrm{p}}/M_{\oplus})^{0.08}$
Earth Oceans, a Dune planet if $10^{-5}(M_{\mathrm{p}}/M_{\oplus})^{0.54}\leq
W_{\mathrm{s}}\leq 0.01(M_{\mathrm{p}}/M_{\oplus})^{0.54}$ (Abe et al., 2011),
and desiccated if $W_{\mathrm{s}}<10^{-5}(M_{\mathrm{p}}/M_{\oplus})^{0.54}$
Earth Oceans, where $10^{-5}$ Earth Oceans is the amount of water currently in
the Earth’s atmosphere (Gleick, 1993). An Earth-like planet would maintain
habitable surface conditions through the silicate weathering thermostat which
regulates the climate over geological timescales (e.g., Walker et al. 1981;
Sleep & Zahnle 2001). A Dune planet or a waterworld, while lacking a silicate
weathering thermostat, may also be habitable, which we discuss further in
Sections 4 and 5.

Figure 1: Surface water regimes for terrestrial planets ranging from
$0.1~{}M_{\oplus}$ to $10~{}M_{\oplus}$. Each limit between regimes is
explained in the text. Planets in either the waterworld, Earth-like, or Dune
planet regime may be habitable, while desiccated surfaces are uninhabitable.
The transition between surface water regimes depends weakly on planetary mass
(cf. Cowan & Abbot 2014).

##  3 Results

###  3.1 Same Initial Water Mass Fraction

We begin with simulations initiated with the same water mass fraction (see
Table 2), which is relevant if water is delivered as part of the planetary
building blocks: a $1~{}M_{\oplus}$ planet begins with [2, 4, 6, 8] Earth
Oceans, while an $8~{}M_{\oplus}$ planet begins with [16, 32, 48, 64] Earth
Oceans. Two such simulations are compared in Fig. 2 for the same initial water
mass fraction, for planets orbiting at the Inner HZ; the magma ocean/runaway
greenhouse phase is shaded grey (left), while the longer deep-water cycling
period is shown on the right. In Fig. 2, although more water becomes
sequestered in the mantle for the super-Earth, it is a waterworld by 5 Gyr,
while the $1~{}M_{\oplus}$ planet is desiccated and uninhabitable.

We also compare the impact of magma ocean-solid mantle water partition
coefficient $D$ on the water cycling between planets. As expected, $D$
controls water partitioning during the magma ocean only, by definition, and
does not change the outcome after 5 Gyr. We discuss this further in Section 4.

Figure 2: Comparison of evolution of water inventories for terrestrial planets
of mass $1~{}M_{\oplus}$ and $8~{}M_{\oplus}$. Each simulation corresponds to
a planet orbiting at the inner edge of the habitable zone (Inner HZ)
initialized with a water mass fraction of $4.7\times 10^{-4}$, and assumes a
magma ocean water partition coefficient of $D=0.001$. Although the
$1~{}M_{\oplus}$ planet becomes desiccated, the $8~{}M_{\oplus}$ planet
retains significant surface water by 5 Gyr. Although the amount of water lost
to space increases with planetary mass (see Table 1), so too does the initial
water endowment and hence the final surface water inventory. More massive
planets have a higher surface gravity, so more water becomes sequestered (and
trapped) within the mantle for larger super-Earths.

The results of the parameter space exploration indicate that there may be
significant water present — and potentially trapped — within the mantle below
a desiccated surface. Predictably, larger planets initiated with more water
(i.e., higher water mass fraction) lead to more water trapped in the mantle
and inundated, waterworld surfaces.

Overall, we find that at fixed water mass fraction, although more water is
lost to space for larger super-Earths, fewer planets become desiccated and
more water becomes sequestered in the mantle for larger planets. Super-Earths
are more likely to become waterworlds due to their comparatively large initial
water endowment with respect to the surface water regimes shown in Fig. 1.

###  3.2 Same Initial Water Inventory

We now show results for simulations that begin with the same initial water
inventory regardless of mass — [2, 4, 6, 8] Earth Oceans — representing water
delivery after planetary formation is complete, possibly through comet
delivery.

Fig. 3 is the same as Fig. 2, but now each planet begins the simulation with
the same water inventory of 4 Earth Oceans and orbits in the middle of the
habitable zone (Mid HZ). The $8~{}M_{\oplus}$ planet becomes desiccated
shortly after 3 Gyr; although energy-limited loss during magma ocean is lower
and larger planetary mantles sequester more water early in the deep-water
cycling phase, this water is degassed to the surface and rapidly lost to
space. Conversely, the $1~{}M_{\oplus}$ planet ends up in a habitable, Earth-
like surface water regime.

Figure 3: Same as Fig. 2, but for planets initiated with the same mass of
water (4 Earth Oceans) and orbiting in the middle of the habitable zone (Mid
HZ). Although the $1~{}M_{\oplus}$ planet (top) ends in an Earth-like surface
water regime, the $8~{}M_{\oplus}$ planet (bottom) has a desiccated,
uninhabitable surface by 5 Gyr. Energy-limited loss during magma ocean
decreases with planetary mass while diffusion-limited loss during deep-water
cycling increases. The effect of surface gravity on mantle water inventory is
visible near the beginning of the deep-water cycling period: more water
becomes sequestered in the mantle of the super-Earth. Due to enhanced
diffusion-limited loss, the surface of the $8~{}M_{\oplus}$ planet becomes
desiccated although some water remains trapped in the solid mantle.

The parameter space results for fixed initial water inventory are plotted as
Initial Water Inventory vs. orbital distance in Fig. 4. This figure also
includes purple boxes around surface water “survivors”. These represent
planets that began the simulations with less water than they would naively be
expected to lose, and so would become desiccated without a deep-water cycle
(see Table 1), but which remain habitable in our simulations. For larger
planets, more water is lost to space, and more simulations result in a
desiccated surface, as indicated by a black ‘X’. However, over half of the
simulated planets in Fig. 4 are in a habitable, Earth-like surface water
regime by 5 Gyr. Water loss is more substantial for super-Earths due to the
increase in diffusion-limited loss. Furthermore, larger planets sequester more
water in their mantles due to their enhanced gravity. If water inventory is
independent of planetary mass, then higher-mass terrestrial planets are more
likely to be desiccated. Indeed, by 5 Gyr, there are many more desiccated
planets and barely any waterworlds compared to simulations initiated with
fixed water mass fraction.

Figure 4: Final water inventories for planets initiated with the same water
inventory (i.e., initial water inventory independent of mass, reasonable for
the dry formation environment around M-dwarfs; see, e.g., Schneeberger et al.
2024). The initial water inventory is indicated by a dashed orange circle, and
the final water inventory by a thick black circle, coloured based on the
surface water regimes of Fig. 1. Additionally, the final mantle water
inventory is plotted — again scaled to the initial inventory — as a filled
brown circle. The purple boxes highlight planets that were expected to lose
all of their water yet remain habitable (see text for details). Desiccated
surfaces are more common on super-Earths by 5 Gyr, with fewer habitable
outcomes for larger planets when initiated with fixed water mass.

We plot the results of all simulations for both formation scenarios as initial
water mass fraction vs. planetary mass in Fig. 5; here, different symbols
represent the three HZ orbital positions, coloured based on their final
surface water regime, and desiccated planets are represented by a black ‘X’.
Note that we run additional simulations to fill some gaps in this figure:
initial water masses of [3, 12] Earth Oceans for each planetary mass. We
discuss the implications of Fig. 5 in the following section.

Figure 5: Final surface water inventory as a function of initial water mass
fraction, planetary mass, and habitable zone orbital distances (slightly
offset for visualization). To first-order, the outcome of a simulation depends
on the planet’s initial water mass, denoted by dashed grey lines; this result
is consistent with Fig. 1. To second-order, we see that at the inner edge of
the HZ, a super-Earth is less habitable than an Earth-sized planet: the larger
planet tends to absorb water in its mantle, in addition to losing it to space,
usually resulting in a desiccated surface.

##  4 Discussion

We have extended the $1~{}M_{\oplus}$ models of Moore, Cowan, & Boukaré (2023)
to investigate water evolution on super-Earths up to $8~{}M_{\oplus}$. We
briefly describe some additional mass-dependent parameters in our model below.
Many aspects, although important, are not feasible for implementation within a
box model, and as such, should be explored in future studies.

Energy-limited loss was expected to decrease with increasing planetary mass,
and diffusion-limited loss was expected to increase with increasing mass.
These behaviours are especially apparent in Fig. 3, when initiating planets
with the same water endowment regardless of their mass. Because the duration
of diffusion-limited loss is much longer than that of energy-limited loss, the
amount of water lost should increase with increasing $M_{\mathrm{p}}$, which
we see in Fig. 4; for our tested planetary masses, the decrease in energy-
limited loss with mass is also consistent with the atmospheric escape
calculations of Chin et al. (2024). While it has been previously predicted
that super-Earth atmospheres are likely stable against XUV-driven thermal
escape (Tian, 2009), subsequent studies predicted immense water loss and
substantial atmospheric oxygen build-up for super-Earths (Luger & Barnes,
2015); however, in all three studies, planetary interior reservoirs were
neglected. Furthermore, diffusion-limited loss depends on the thermospheric
temperature $T_{\mathrm{therm}}$. We adopt $T_{\mathrm{therm}}=400$ K based on
previous studies (e.g., Luger & Barnes 2015), but we note the linear
dependence on $T_{\mathrm{therm}}$ in Eqn. 4, such that increasing
$T_{\mathrm{therm}}$ to ${\sim}1000$ K would result in a decrease of
diffusion-limited loss by a factor of 2.5.

Our energy-limited loss parameterization presumes that the XUV deposition
radius $R_{\mathrm{XUV}}$ is simply the planetary radius $R_{\mathrm{p}}$.
However, $R_{\mathrm{XUV}}$ could be greater than $R_{\mathrm{p}}$, especially
for a puffed-up hot, wet atmosphere above a magma ocean. We therefore test
three additional XUV deposition radii to elucidate its effect on our overall
results: $R_{\mathrm{XUV}}=$ [$R_{\mathrm{p}}+500$ km, $R_{\mathrm{p}}+1000$
km, $R_{\mathrm{p}}+2000$ km]. We find that the 1 $M_{\oplus}$ planet
initiated with 4 Earth Oceans loses more water when $R_{\mathrm{XUV}}$ is
increased, becoming desiccated by 5 Gyr, while this planet survived the
original $R_{\mathrm{XUV}}=R_{\mathrm{p}}$ simulations as a habitable planet.
We see this same behaviour for the 2 $M_{\oplus}$ planet initiated with 6
Earth Oceans; however, these are the only two simulations that are
significantly changed. Regardless, our loss rates during the magma
ocean/runaway greenhouse period may still be too high, as a recent study
including H2O photochemistry during the stellar pre-main sequence found that
water loss is suppressed (Kawamura et al., 2024).

Super-Earths form with greater accretional energy than their Earth-mass
counterparts, and likely have deeper magma oceans (Stixrude et al., 2020).
Further, the water capacity of the early magma ocean, within which water is
protected from loss to space, should increase with increasing mass (cf. Cowan
& Abbot 2014). A more massive planet will cool slower since volume is a
stronger function of mass than surface area (e.g., Seales & Lenardic 2021). A
slowly-solidifying super-Earth may have a longer magma ocean phase than its
corresponding runaway greenhouse phase. This could prolong the magma ocean
phase, and maintain a hotter solid mantle during deep-water cycling.

We note that the higher surface gravity of super-Earths may lead to the
formation of high-pressure ices, inhibiting the deep-water cycling and
seafloor weathering, and hence the silicate weathering thermostat; however,
even our highest tested initial water inventory of 64 Earth Oceans is far
below the conservative high-pressure ice limit of ${\sim}$100 Earth Oceans
(Alibert, 2014; Nakayama et al., 2019).

We initiate every deep-water cycling simulation with the same solid mantle
temperature of $T_{\mathrm{m,0}}=3000$ K. This value should probably scale
with planetary mass, and a hotter mantle would presumably lead to more intense
early degassing during deep-water cycling and, if the mantle remains hydrated,
an extended period of plate tectonics (if plate tectonics begins at all on
super-Earths; see, e.g., Valencia et al. 2007; Valencia & O’Connell 2009;
Korenaga 2010; Foley et al. 2012). We can estimate $T_{\mathrm{m,0}}=4880$ K
for an $8~{}M_{\oplus}$ planet by extrapolating Table 5 of Schaefer & Sasselov
(2015). Using this as the initial solid mantle temperature for deep-water
cycling does lead to more early degassing in our model, which slightly offsets
the amount of water loss to space for a period; however, the final water
inventories and partitioning remain unchanged by 5 Gyr regardless of
$T_{\mathrm{m,0}}$. Indeed, as noted by McGovern & Schubert (1989) and
Schaefer & Sasselov (2015), the effect of initial mantle temperature is
minimal after a few hundred Myrs.

In this study we have adopted the HZ limits of Kopparapu et al. (2013), but
habitable planets orbit very close to their M-dwarf hosts and are likely
tidally-locked (e.g., Kasting et al. 1993). GCM simulations suggest that the
inner edge moves inwards to 2$\times$ the stellar flux at the runaway
greenhouse limit (Yang et al., 2013, 2014). We therefore perform a sensitivity
analysis by instead taking the inner edge to be at 2$\times$ the RG limit of
Kopparapu et al. (2013), i.e.,
$a_{\mathrm{InnerHZ,tl}}=\sqrt{(a_{\mathrm{InnerHZ}}^{2})/2}=0.017$ AU. For
this closer Inner HZ limit, our model predicts that the planet is permanently
in a runaway greenhouse and never becomes habitable (cf. Turbet et al. 2023).
In principle, a wider habitable zone incorporating a closer Inner HZ edge
could benefit the super-Earths in our results (see Fig. 5) if the planet is
able to exit the runaway greenhouse phase while retaining a water inventory.

The magma ocean water partition coefficient $D$ controls the amount of water
in the solid mantle as the magma ocean solidifies. As the magma ocean
solidifies from the core-mantle boundary towards the surface, the magma ocean
water partition coefficient should change since it is a function of pressure
(Papale, 1997, 1999); indeed, internal pressures of super-Earths will also be
higher than for Earth-mass planets. We hence test $D=0.001$, 0.01, and 0.1,
since we cannot fully represent this pressure-dependence in our 0-D box model
as we do not resolve the 1-D internal structure. As expected, increasing $D$
leads to more substantial water sequestered in the solid mantle during the
magma ocean phase.

However, $D$ only controls partitioning during the relatively short-lived
magma ocean phase; less of the total water inventory may be lost to space than
expected, but water sequestered in the mantle does not necessarily enhance
surface habitability. All else remaining the same, the amount of water present
within the solid mantle by 5 Gyr is nearly identical regardless of $D$ due to
the short lifetime of the magma ocean compared to the much longer deep-water
cycling phase and the significant degassing from the hot mantle following
surface solidification, so varying $D$ has nearly no effect on the final
planetary water inventories, their surface water regime, nor their
habitability prospects. A Gyr-lived basal magma ocean — posited by Moore,
Cowan, & Boukaré (2023) for an Earth-mass planet, and probably more likely to
exist in super-Earths (C. É. Boukaré, priv. comm.) — would presumably cause
even more water to become trapped in the solid mantle than our current
simulations indicate. This reinforces our conclusion that substantial water
becomes locked in the solid mantle for super-Earths, in agreement with
Schaefer & Sasselov (2015).

##  5 Summary & Conclusions

The surface of a super-Earth should accommodate slightly more water before
becoming inundated. For the same initial water mass fraction, more massive
planets are likely to avoid surface desiccation and become waterworlds with
substantial water inventories trapped in their mantle; indeed, only 13/48 such
simulations end in a habitable Earth-like surface water regime, while a
staggering 30/48 end up as waterworlds. This is mainly due to the very large
water inventories for super-Earths and the weak dependence of surface water
regime on mass shown in Fig. 1. For this formation scenario, i.e., water as a
consequence of accretion, although smaller planets are more likely to become
Earth-like, super-Earths are less prone to desiccation and retain very wet
mantles. Hence, super-Earths are possibly more likely to be habitable if
waterworlds can support the origin and evolution of life. Indeed, recent
studies support the long-term habitability of waterworlds in the absence of
silicate weathering through CO2-driven buffering (Ramirez & Levi, 2018; Kite &
Ford, 2018), which can even operate in the absence of geochemical cycling;
however, these results do not apply to an M8 host star, as tested throughout
the current study.

Fig. 4, instead initializing all planets with 2–8 Earth Oceans, reveals that
surface desiccation is more likely for super-Earths due to the enhanced
diffusion-limited loss during deep-water cycling, but also because more water
becomes sequestered in the solid mantle. In this dry, in-situ formation
scenario (e.g., Schneeberger et al. 2024), super-Earths are more likely to
become desiccated and less likely to end up in an Earth-like surface water
regime.

Although the wide swath in the log-plotted Fig. 1 makes Dune planets seem like
a very likely outcome of water evolution, our coarse parameter space search
does not produce any. Since loss exceeds degassing late in our simulations,
Dune planets may be very rare around late M-dwarfs. Although Dune planets lack
a silicate weathering thermostat, the low planetary water content and dry,
unsaturated atmosphere would both reduce water loss to space and widen the
predicted habitable zone (Abe et al., 2011); since we do not find Dune planets
in our parameter search, however, we leave this detail to future model
iterations. Regardless, we can still estimate rough desiccation boundaries for
a given planetary mass from Fig. 5: 1 $M_{\mathrm{\oplus}}$ becomes desiccated
between ${\sim}$2–3 Earth Oceans, while this boundary falls between 2–8 Earth
Oceans for an 8 $M_{\mathrm{\oplus}}$ super-Earth depending on orbital
distance, with super-Earths at the Inner HZ most likely to become desiccated.
Waterworlds arise for 1 $M_{\mathrm{\oplus}}$ above about 8 Earth Oceans, and
for 8 $M_{\mathrm{\oplus}}$ above about 12 Earth Oceans.

Plate tectonics requires both water dissolved in the mantle and moderate
surface temperatures. Our model is restricted to the plate-tectonic mode,
although it is possible that super-Earths are instead in a stagnant lid regime
similar to modern-day Venus, especially if surface oceans are not formed
through magma ocean degassing (O’Neill & Lenardic, 2007; Miyazaki & Korenaga,
2022). Insofar as stagnant lid planets degas water more effectively than they
regas it, this tectonic state could help super-Earths release water that would
be otherwise become trapped in their mantle by plate tectonics. However, we
leave investigation of alternate tectonic modes to future modeling efforts.

## Acknowledgments

KM and NBC thank Sean Raymond and Charles-Édouard Boukaré for useful
exchanges, and the anonymous referees for insightful comments and suggestions.
NBC acknowledges support from an NSERC Discovery Grant, a Tier 2 Canada
Research Chair, and an Arthur B. McDonald Fellowship. The authors also thank
the Trottier Space Institute (TSI) and l’Institut de recherche sur les
exoplanètes (iREx) for their financial support and dynamic intellectual
environment. This research was partially funded by the Heising-Simons
Foundation.

## Data Availability

The model and simulation results presented within this manuscript are
available in the following Zenodo repository: 10.5281/zenodo.12548229 (Moore,
2024).

## References

  * Abe et al. (2011) Abe, Y., Abe-Ouchi, A., Sleep, N. H., & Zahnle, K. J. 2011, Astrobiology, 11, 443
  * Alibert (2014) Alibert, Y. 2014, A&A, 561, A41. doi:10.1051/0004-6361/201322293
  * Alibert & Benz (2017) Alibert, Y. & Benz, W. 2017, A&A, 598, L5. doi:10.1051/0004-6361/201629671
  * Baraffe et al. (2015) Baraffe, I., Homeier, D., Allard, F., & Chabrier, G. 2015, A&A, 577, A42
  * Barnes et al. (2016) Barnes, R., Deitrick, R., Luger, R., Driscoll, P. E., et al. 2016, arXiv e-prints, 1608.06919v2
  * Chachan & Lee (2023) Chachan, Y. & Lee, E. J. 2023, arXiv:2305.00803. doi:10.48550/arXiv.2305.00803
  * Checlair et al. (2019) Checlair, J. H., Salazar, A. M., Paradise, A., et al. 2019, ApJ, 887, L3. doi:10.3847/2041-8213/ab5957
  * Chin et al. (2024) Chin, L., Dong, C., & Lingam, M. 2024, arXiv:2401.16211. doi:10.48550/arXiv.2401.16211
  * Cowan & Abbot (2014) Cowan, N. B., & Abbot, D. S. 2014, ApJ, 781, 27
  * Dorn & Lichtenberg (2021) Dorn, C., & Lichtenberg, T. 2021, ApJ, 922, L4. doi:10.3847/2041-8213/ac33af
  * Dressing & Charbonneau (2015) Dressing, C. D., & Charbonneau, D. 2015, ApJ, 807, 45
  * Ercolano & Clarke (2010) Ercolano, B. & Clarke, C. J. 2010, MNRAS, 402, 2735. doi:10.1111/j.1365-2966.2009.16094.x
  * Foley et al. (2012) Foley, B. J., Bercovici, D., & Landuyt, W. 2012, Earth and Planetary Science Letters, 331, 281. doi:10.1016/j.epsl.2012.03.028
  * Gleick (1993) Gleick, P. H. 1993, Water in Crisis: A Guide to the World’s Fresh Water Resources (New York: Oxford University Press)
  * Hsu et al. (2020) Hsu, D. C., Ford, E. B., & Terrien, R. 2020, MNRAS, 498, 2249
  * Hunten et al. (1987) Hunten, D. M., Pepin, R. O., & Walker, J. C. G. 1987, Icarus, 69, 532. doi:10.1016/0019-1035(87)90022-4
  * Kasting et al. (1993) Kasting, J. F., Whitmire, D. P., & Reynolds, R. T. 1993, Icarus, 101, 108
  * Katz et al. (2003) Katz, R. F., Spiegelman, M., & Langmuir, C. H. 2003, GGG, 4, 1073
  * Kawamura et al. (2024) Kawamura, Y., Yoshida, T., Terada, N., et al. 2024, ApJ, 967, 95. doi:10.3847/1538-4357/ad3e7e
  * Kite & Ford (2018) Kite, E. S. & Ford, E. B. 2018, ApJ, 864, 75. doi:10.3847/1538-4357/aad6e0
  * Komacek & Abbot (2016) Komacek, T. D., & Abbot, D. S. 2016, ApJ, 832, 54
  * Kopparapu et al. (2013) Kopparapu, R. K., Ramirez, R., Kasting, J. F., Eymet, V., et al. 2013, ApJ, 765, 131
  * Korenaga (2010) Korenaga, J. 2010, ApJ, 725, L43
  * Krissansen-Totton et al. (2021a) Krissansen-Totton, J., Galloway, M. L., Wogan, N., Dhaliwal, J. K., & Fortney, J. J. 2021, ApJ, 913, 107
  * Krissansen-Totton et al. (2021b) Krissansen-Totton, J., Fortney, J. J., & Nimmo, F. 2021, , 2, 216. doi:10.3847/PSJ/ac2580
  * Krissansen-Totton et al. (2022) Krissansen-Totton, J., & Fortney, J. J. 2022, ApJ, 933, 115
  * Kruijver et al. (2021) Kruijver, A., Höning, D., & van Westrenen, W. 2021, The Planetary Science Journal, 2, 208. doi:10.3847/PSJ/ac24aa
  * Li et al. (2022) Li, H.-F., Oganov, A. R., Cui, H., et al. 2022, Phys. Rev. Lett., 128, 035703. doi:10.1103/PhysRevLett.128.035703
  * Lopez (2017) Lopez, E. D. 2017, MNRAS, 472, 245. doi:10.1093/mnras/stx1558
  * Luger & Barnes (2015) Luger, R., & Barnes, R. 2015, Astrobiology, 15, 119
  * McGovern & Schubert (1989) McGovern, P. J., & Schubert, G. 1989, Earth and Planetary Science Letters, 96, 27
  * Miyazaki & Korenaga (2022) Miyazaki, Y. & Korenaga, J. 2022, Astrobiology, 22, 713. doi:10.1089/ast.2021.0126
  * Moore & Cowan (2020) Moore, K., & Cowan, N. B. 2020, MNRAS, 496, 3786
  * Moore, Cowan, & Boukaré (2023) Moore, K., Cowan, N. B., & Boukaré, C.-É. 2023, MNRAS, 526, 6235. doi:10.1093/mnras/stad3138
  * Moore (2024) , Model & data associated with “Water Evolution & Inventories of Super-Earths Orbiting Late M-Dwarfs”, v1.0, Zenodo, doi:10.5281/zenodo.12548229
  * Mulders et al. (2015) Mulders, G. D., Pascucci, I., & Apai, D. 2015, ApJ, 798, 112
  * Nakayama et al. (2019) Nakayama, A., Kodama, T., Ikoma, M., & Abe, Y. 2019, MNRAS, 488, 1580
  * O’Neill & Lenardic (2007) O’Neill, C., & Lenardic, A. 2007, Geophysical Research Letters, 34, L19204
  * Papale (1997) Papale, P. 1997, Contrib. Mineral Petrol, 126, 237
  * Papale (1999) Papale, P. 1999, American Mineralogist, 84, 477
  * Ramirez & Kaltenegger (2014) Ramirez, R. M. & Kaltenegger, L. 2014, ApJ, 797, L25. doi:10.1088/2041-8205/797/2/L25
  * Ramirez & Levi (2018) Ramirez, R. M. & Levi, A. 2018, MNRAS, 477, 4627. doi:10.1093/mnras/sty761
  * Sabotta et al. (2021) Sabotta, S., Schlecker, M., Chaturvedi, P., Guenther, E. W., et al. 2021, A&A, 653, A114. doi:10.1051/0004-6361/2021409
  * Scalo et al. (2007) Scalo, J., Kaltenegger, L., Segura, A. G., Fridlund, M. Ribas, I., Kulikov, Y. N., Grenfell, J. L., Rauer, H., Odert, P., Leitzinger, M., Selsis, F., Khodachenko, M. L., Eiroa, C., Kasting, J., & Lammer, H. 2007, Astrobiology, 7, 85
  * Schaefer & Sasselov (2015) Schaefer, L., & Sasselov, D. 2015, ApJ, 801, 40
  * Schneeberger et al. (2024) Schneeberger, A., Mousis, O., Deleuil, M., et al. 2024, arXiv:2401.15981. doi:10.48550/arXiv.2401.15981
  * Seales & Lenardic (2021) Seales, J. & Lenardic, A. 2021, Icarus, 367, 114560. doi:10.1016/j.icarus.2021.114560
  * Shields et al. (2016) Shields, A. L., Ballard, S., & Johnson, J. A. 2016, Phys. Rep., 663, 1. doi:10.1016/j.physrep.2016.10.003
  * Sleep & Zahnle (2001) Sleep, N. H., & Zahnle, K. 2001, Journal of Geophysical Research, 106, 1373
  * Stamenković et al. (2012) Stamenković, V., Noack, L., Breuer, D., & Spohn, T. 2012, ApJ, 748, 41
  * Stixrude (2014) Stixrude, L. 2014, Phil. Trans. R .Soc. A 372:20130076. http://dx.doi.org/10.1098/rsta.2013.0076
  * Stixrude et al. (2020) Stixrude, L., Scipioni, R., & Desjarlais, M. P. 2020, Nature Communications, 11, 935. doi:10.1038/s41467-020-14773-4
  * Tian (2009) Tian, F. 2009, ApJ, 703, 905. doi:10.1088/0004-637X/703/1/905
  * Tian & Ida (2015) Tian, F. & Ida, S. 2015, Nature Geoscience, 8, 177. doi:10.1038/ngeo2372
  * Townsend et al. (2016) Townsend, J. P., Tsuchiya, J., Bina, C. R., et al. 2016, Earth and Planetary Science Letters, 454, 20. doi:10.1016/j.epsl.2016.08.009
  * Turbet et al. (2021) Turbet, M., Bolmont, E., Chaverot, G., et al. 2021, Nature, 598, 276. doi:10.1038/s41586-021-03873-w
  * Turbet et al. (2023) Turbet, M., Fauchez, T. J., Leconte, J., et al. 2023, A&A, 679, A126. doi:10.1051/0004-6361/202347539
  * Valencia et al. (2006) Valencia, D., O’Connell, R. J., & Sasselov, D. 2006, Icarus, 181, 545
  * Valencia et al. (2007) Valencia, D., O’Connell, R. J., & Sasselov, D. D. 2007, ApJ, 670, L45. doi:10.1086/524012
  * Valencia & O’Connell (2009) Valencia, D. & O’Connell, R. J. 2009, Earth and Planetary Science Letters, 286, 492. doi:10.1016/j.epsl.2009.07.015
  * Walker (1977) Walker, J. C. G. 1977, Evolution of the Atmosphere (New York: MacMillan)
  * Walker et al. (1981) Walker, J. C. G., Hays, P. B., & Kasting, J. F. 1981, Journal of Geophysical Research, 86, 9776
  * Watson et al. (1981) Watson, A. J., Donahue, T. M., & Walker, J. C. G. 1981, Icarus, 48, 150
  * Wordsworth & Pierrehumbert (2013) Wordsworth, R. D., & Pierrehumbert, R. T. 2013, ApJ, 778, 154
  * Wordsworth & Pierrehumbert (2014) Wordsworth, R., & Pierrehumbert, R. 2014, ApJ, 785, L20
  * Yang et al. (2013) Yang, J., Cowan, N. B., & Abbot, D. S. 2013, ApJ, 771, L45. doi:10.1088/2041-8205/771/2/L45
  * Yang et al. (2014) Yang, J., Boué, G., Fabrycky, D. C., et al. 2014, ApJ, 787, L2. doi:10.1088/2041-8205/787/1/L2


# Deep Learning-based Analysis of Basins of Attraction

David Valle<EMAIL_ADDRESS>Nonlinear Dynamics, Chaos and Complex Systems
Group, Departamento de Física, Universidad Rey Juan Carlos, Tulipán s/n, 28933
Móstoles, Madrid, Spain
Alexandre Wagemakers<EMAIL_ADDRESS>Nonlinear Dynamics, Chaos
and Complex Systems Group, Departamento de Física, Universidad Rey Juan
Carlos, Tulipán s/n, 28933 Móstoles, Madrid, Spain
Miguel A.F. Sanjuán<EMAIL_ADDRESS>Nonlinear Dynamics, Chaos and
Complex Systems Group, Departamento de Física, Universidad Rey Juan Carlos,
Tulipán s/n, 28933 Móstoles, Madrid, Spain

###### Abstract

This study showcases the effectiveness of convolutional neural networks (CNNs)
in characterizing the complexity and unpredictability of basins of attraction
for diverse dynamical systems. This novel method is optimal for exploring
different parameters of dynamical systems since the conventional methods are
computationally expensive for characterizing multiple basins of attraction.
Additionally, our research includes a comparison of different CNN
architectures for this task showing the superiority of our proposed
characterization method over the conventional methods, even with obsolete
architectures.

###### pacs:

$05.45.-a$, $05.45.Df$, $07.05.Mh$

> The application of machine learning algorithms to complex dynamical systems
> has opened up new possibilities for predicting the behavior of these systems
> under different initial conditions. One important aspect of studying
> dynamical systems is understanding their basins of attraction which are the
> regions in the phase space where initial conditions lead towards a
> particular attractor. Since we are interested in the asymptotic behavior of
> trajectories, the basins of attraction provide a visual tool to have a clear
> image of which set of initial conditions bring their associated trajectories
> to specific attractors in phase space. In chaotic systems, the boundary sets
> between these basins can be fractal, what can make the long-term prediction
> of initial conditions challenging. Basins of attraction can be characterized
> using various metrics, such as the fractal dimension, the basin entropy, the
> boundary basin entropy, and the Wada property. These measurements provide a
> quantitative understanding of the behavior of dynamical systems. Moreover,
> their prediction involves complex nonlinear operations that are well adapted
> for machine learning algorithms. Additionally, convolutional neural networks
> are particularly suitable for this task due to the representation of basins
> as images, where each pixel corresponds to an initial condition and each
> color represents a different attractor. The use of machine learning
> algorithms in this context constitutes an important step towards predicting
> the sensitivity of the system to initial conditions, and improving our
> understanding of the behavior of complex dynamical systems.

##  I Introduction

Complex dynamical systems are typically modeled through differential or
difference equations and therefore follow deterministic rules. However, even
though they are deterministic, it does not necessarily mean that these systems
are easy to predict. As well known, in the case of chaotic dynamics, a small
variation in the initial conditions may result that trajectories would go to a
different attractor after a sufficiently long run. In particular, in an
experimental scenario a slight variation in the initial conditions can
completely alter the long-term evolution when the system happens to be
chaotic. This inherent property of chaos is closely related to fractal
geometry, so it is not surprising that both disciplines are so closely
connected Peitgen, Jürgens, and Saupe (1992).

Fractal geometry and basins of attraction are related concepts in the study of
dynamical systems. Fractal geometry is a mathematical framework for describing
objects that display self-similar patterns at different scales. In the context
of dynamical systems, fractal sets such as the Mandelbrot set are often used
to illustrate the behavior of iterated functions Peitgen, Jürgens, and Saupe
(1992). Basins of attraction refer to the regions in the phase space of a
dynamical system where initial conditions lead towards a particular attractor,
and boundary sets refer to the separation between the different basins Nusse
and Yorke (1996a); Aguirre, Viana, and Sanjuán (2009). The relationship
between fractal geometry and basins of attraction lies in the fact that the
boundary sets can be fractal, providing a manifestation of the complex and
intricate structure of the phase space. In some multistable dynamical systems,
the basin boundaries display complicated fractal structures at every scale,
hindering the long-term prediction of initial conditions starting in the
vicinity of such boundaries.

The analysis of basins of attraction with different metrics brings practical
and theoretical information on the dynamical system. The asymptotic behavior
of the system is embodied into a single structure that matches an initial
condition to an attractor. By tracking the properties of this structure as a
parameter of the system vary, we gain insight into how the system evolves with
that parameter modification. It enhances the understanding of the
predictability and stability of the system and also eases the identification
of the bifurcation points. Obviously, the applications of these ideas can be
useful for the analysis of different problems in science and technology that
are modeled using dynamical systems.

The unpredictability of a basin can be characterized using various metrics
such as the fractal dimension, the basin entropy, the boundary basin entropy,
and the Wada property. The fractal dimension can be used as a metric to
quantify the long term uncertainty through the irregularity of the boundary
between different basins of attraction McDonald  _et al._ (1985). The basin
entropy is a scalar measure of the unpredictability of phase space inspired by
the classical Gibbs entropy in statistical mechanics. The entropy of a small
discretized region of the phase space will depend on the number of attractors
and their proportion in the sampled region. Basin entropy is an average of
this local entropy over all the phase space, providing an overall value for
the unpredictability of the basins of attraction Daza  _et al._ (2016). The
boundary basin entropy is calculated using the same procedure as the basin
entropy, but the discretization occurs only in regions where a boundary is
present, therefore this metric reflects the uncertainty associated to the
boundaries of the basins of attraction Daza  _et al._ (2016).

Finally, the Wada property is a topological property associated to dynamical
systems with three or more basins of attraction, where a boundary point
belongs at the same time to all the boundaries. In other words, every point in
the boundary is in the boundary of all basins which is unique. The interest in
the Wada basin boundaries lies on the fact that they are the most entangled
boundaries one can imagine, since they separate all the basins at the same
time. Furthermore, the existence of this property implies a higher degree of
unpredictability for the dynamical system. As a matter of fact, small
perturbations near the boundaries can lead to any of the different attractors
of the system. Recently, an intermediate situation has been analyzed in Daza,
Wagemakers, and Sanjuán (2018); DazaClasifying2022 ; Wagemakers, Daza, and
Sanjuán (2020), where some regions have the Wada property but others do not,
so that the basins are called partial Wada basins.

Basins of attraction are well suited for the application of machine learning
techniques to predict the metrics previously described. The prediction of
these measurements involves complex nonlinear operations. Since machine
learning algorithms infer nonlinear relationships from data, they fit
perfectly for the task. From a more specific point of view, basins can be
represented as images where each pixel corresponds to an initial condition,
and the color represents the attractor. Therefore, we can clearly distinguish
patterns in basins, such as filaments or boundaries meandering repeatedly. Due
to the presence of these patterns, machine learning algorithms perform
remarkably well to analyze these objects.

The use of machine learning and deep learning techniques opens a wide field of
applications for the classification of basins of attraction. Indeed, machine
learning has already been used to help solve problems in nonlinear dynamics.
Among them, we can mention improving the prediction time of chaotic
trajectories Fan  _et al._ (2020); Pathak  _et al._ (2018), reconstructing the
long-term dynamics of a system Lu, Hunt, and Ott (2018), as well as
classifying data with associated fractal properties Kirichenko, Radivilova,
and Bulakh (2019); Shi (2018), and many other applications described in a
recent focus issue of this journal Tang  _et al._ (2020).

An innovative use of neural network capabilities was recently published in
Valle  _et al._ (2022). The estimation of the fractal dimension of the basins
of attraction of the Duffing oscillator was compared with the classical method
proposed by McDonald et al. McDonald  _et al._ (1985). The new method showed
similar robustness but required less computation time. Our main objective here
is to efficiently quantify the uncertainty of the phase space of common
classes of dynamical systems. To achieve this, we will quantify the
aforementioned metrics on the basins of attractions that provides information
on the predictability of a dynamical system. The result is the quantification
of previously computed basins of attraction using trained convolutional neural
networks.

Convolutional Neural Networks (CNN) are a type of artificial neural networks
designed specifically for image classification tasks. They were developed by
LeCun et al. Lecun  _et al._ (1998) in 1998 as a class of deep feed‐forward
artificial neural networks. They are composed of multiple layers, including
convolutional layers, pooling layers, and fully connected layers. The
convolutional layers apply a set of filters to the input image producing a
feature map. These filters can be learned during the training process and are
responsible for detecting specific features in the image, such as edges,
shapes, and textures. The pooling layers perform down-sampling operations on
the feature maps, reducing their spatial dimensions while retaining the most
important information. This helps to reduce the computation required by the
network and also provides some degree of translation invariance, allowing the
network to recognize features regardless of their position in the image
Goodfellow, Bengio, and Courville (2016); Szeliski (2011).

The fully connected layers receive the feature maps from the pooling layers
and produce the final output of the network. These layers perform a nonlinear
mapping from the feature maps to the output, using multiple activation
functions and weight matrices. The weights of the network are learned during
the training process, using labeled training data, and are updated to minimize
the error between the predicted and actual outputs Goodfellow, Bengio, and
Courville (2016); Szeliski (2011). This supervised learning approach is
applied to the set of pre-labeled basins that act as targets. Labels are
simply the metrics of interest computed previously with an standard algorithm
from the literature. With this setup, we obtain a trained CNN that can
accurately predict the desired measurements for at least all of the systems
analyzed in this paper.

The rest of the paper is organized as follows. In Sec. II, we provide insight
about basins of attraction and their characterization metrics, followed by
information about the computation of these metrics and the dynamical systems
used. Section III presents the architectures of several CNNs used in our
research, the image pre-processing, the image dataset distribution, and the
hyperparameter settings for the training of the CNNs. Section IV presents a
comparison of the performance of the proposed CNNs between themselves and
reference methods for calculating the metrics, demonstrating the superior
efficiency of our method.

##  II Complex dynamical systems and their characterization metrics

We have computed basins of attraction from the Duffing oscillator, the forced
damped pendulum, the Newton fractal, the magnetic pendulum, and finally the
exit basins from the Hénon-Heiles system (which in our work fulfill the same
purpose as the basins of attraction). Examples of basins from these systems
can be seen in Fig. 1. Each basin is made of $110889$ trajectories, therefore
obtaining basins of $333\times 333$ pixels for different choices of parameters
of such systems. Additional information about the computation performed here
can be found in Appendix A.

Figure 1: Examples of basins from the dynamical systems used in our work. (a)
This figure corresponds to a basin of attraction from the Duffing oscillator,
(b) This one corresponds to a basin from the Newton fractal, (c) corresponds
to a basin of attraction of the forced damped pendulum, (d) is an exit basin
from the Hénon-Heiles Hamiltonian system, and finally (e) corresponds to an
example of a basin of attraction of the magnetic pendulum with $4$ magnets.
Further details about these dynamical systems can be found in Appendix A.

We summarize the main features about the computation of the metrics selected
for the study of the basins, and provide some details:

  * •

The fractal dimension is a positive real number $FDim\in[D-1,D]$ where D is
the dimension of the phase space. For $D=2$, a fractal basin is characterized
by $FDim<2$, meaning that the structures have a dimension somewhere between
the line and the plane. To compute this number, we have implemented a Monte
Carlo adaptation of the the box-counting method proposed by McDonald et al.
McDonald  _et al._ (1985). Briefly, the method consists of detecting if a box
of side length $\varepsilon$ located at random contains a piece of the
boundary. If it does, the box is labeled as uncertain. The ratio of uncertain
boxes over the total number of boxes $f(\varepsilon)$ depends on the box size.
The fractal dimension is computed with a linear fit of the curve $\log
f(\varepsilon)$ against $\log\varepsilon$. In our implementation, the fractal
dimension was computed by randomly throwing boxes of side length ranging from
$3$ to $33$ pixels on the basin. The box length size was increased by $3$
pixels every $350000$ boxes.

  * •

The basin entropy is associated to the degree of uncertainty of the basin
represented by a positive number which satisfies $Sb\in[0,log(N_{a})]$ where
$N_{a}$ refers to the number of attractors in the basin. For basins with
$Sb=0$ there is only one attractor, whereas for basins with $Sb=log(N_{a})$
there are $N_{a}$ basins with intermingled boundaries DazaClasifying2022 . To
compute this number, we have implemented a Monte Carlo adaptation of the
method proposed by Daza et al. in Daza  _et al._ (2016) where we measure the
entropy of a box located at random in the basin, assuming that the
trajectories inside that box are statistically independent. This measure is
then repeated over a fixed number of random boxes, and the final basin entropy
is the average entropy of all the boxes. In our implementation, the basin
entropy was computed by randomly throwing $350000$ boxes of side length $15$
pixels.

  * •

The boundary basin entropy $S_{bb}$ is associated to the degree of uncertainty
of the boundaries of the basin DazaClasifying2022 . It is calculated using the
same procedure as the basin entropy taking into account only the regions where
a boundary is present. This metric provides a sufficient condition to assess
easily that the boundaries are fractal, being $S_{bb}>log(2)$. This is a
sufficient but not necessary condition, since there might be fractal
boundaries with $Sbb<log2$. Nevertheless, this threshold can be very useful to
assess quickly the fractality of some boundaries Daza  _et al._ (2016).

  * •

Wada boundaries separate three or more basins at a time. From the dynamical
point of view, the most interesting feature of Wada basins is the fact that an
arbitrarily small perturbation of a system with initial conditions lying in a
Wada boundary can drive the trajectory to any of the at least three possible
attractors. It implies a special kind of unpredictability. In our work the
presence of the Wada property was verified using the merging method proposed
by Daza et al. Daza, Wagemakers, and Sanjuán (2018). This method rely on the
surprising fact that merging two of the Wada basins leaves unchanged the
structure of the boundary. To asses the Wada property, it suffices to check
that the boundary remains unaltered under the operation of merging any two
basins. In our implementation of the merging method, the pixels of the merged
boundaries and the original boundary are compared with some margin of error
defined as the fattening parameter $r=5$.

We have computed each of these metrics for all our basins. The first three
metrics were generated $10$ times for each basin since their calculations were
done as a Monte Carlo method. By computing the same metric $10$ times for a
basin, we have defined the true value of the metric as the mean value of those
$10$ measurements. The robustness and deviation of the measurements with
respect to their mean value are discussed later in Sec. V.

Once we characterized all the basins used on our work, we got a total of
$84842$ basins from the Duffing oscillator, $14276$ basins from the Newton
fractal, $18926$ basins from the forced damped pendulum and $1230$ basins from
the Hénon-Heiles system, and $14384$ basins from the magnetic pendulum basins.
All these basins will be later used by the CNN to train and test for the
accuracy of the predictions on each of the characterization metrics.

##  III Convolutional neural networks and their training

We compare different CNN architectures for characterizing basins of attraction
using various metrics. The architectures compared in our work are: AlexNet
Krizhevsky, Sutskever, and Hinton (2012), VGG16 Simonyan and Zisserman (2015),
VGG19 Simonyan and Zisserman (2015), GoogLeNet Szegedy  _et al._ (2015), and
ResNet50 He  _et al._ (2016). Each architecture is schematized in Fig. 2. This
choice of architectures is due to historical reasons, since they were all
breakthroughs in their respective fields, introducing novel techniques that
improved the accuracy of image classification tasks and allowed for deeper and
more complex neural networks to be trained effectively.

To calculate the different characterization metrics, we have trained separate
CNNs for each architecture mentioned above. In particular, we trained four
CNNs of the same architecture, with each CNN returning a single metric. We use
the following terminology: FDim for the fractal dimension, Sb for the basin
entropy, Sbb for the boundary basin entropy, and Wada for the presence of the
Wada property.

The different architectures were implemented using the Python programming
language libraries Tensorflow (2.7.0) and Keras (2.7.0). These libraries
provide the building blocks necessary to develop the desired architectures,
and all these were built following the same approach as explained by the
authors of their corresponding papers Krizhevsky, Sutskever, and Hinton
(2012); Simonyan and Zisserman (2015); Szegedy  _et al._ (2015); He  _et al._
(2016), with the exception of AlexNet and VGG architectures. In the case of
AlexNet, the top dense layers were modified adding dropout Srivastava  _et
al._ (2014) to prevent overfitting and reducing the amount of neurons in half
due to insufficient RAM memory. Whereas for the VGG architectures, the top
dense layers were removed due to their poor results.

Regarding the image pre-processing, after the transformation of the basins to
images, we started with matrices sized $333$ $\times$ $333$ which values were
the classified attractors corresponding to each initial condition. We then
reshaped each basin to $333$ $\times$ $333$ $\times$ $1$ as an image with same
width and height as the number of initial conditions on each axis of the
basin. The extra dimension corresponds to a grayscale color channel, where
each pixel intensity is associated to the previously classified attractor.

Due to the number of images used in the training process, loading the whole
image set overflows the RAM. Hence, we have developed a custom generator that
has allowed us to randomly select groups of $16$ images, eventually granting
the CNN to train with all the images in the set. Additionally, during the
random grouping of images for the training, the images are subjected to two
possible data augmentation techniques. They correspond to a horizontal and a
vertical flip, each of them occurring with a probability of 50$\%$. So even by
randomly sampling the training set, in case that the same image is repeated
for the training, we still have a 50$\%$ chance of getting flipped
horizontally, and a 50$\%$ of getting flipped vertically, or a 25$\%$ chance
getting flipped both ways.

In terms of the distribution of images used to form the training, validation,
and test sets, the basins from the Duffing oscillator and the Newton fractal
compose the training image set. The basins from the forced damped pendulum and
the Hénon-Heiles Hamiltonian system are used as the validation image set,
which leaves the basins from the magnetic pendulum to compose the test image
set. The distribution of systems in each set was adjusted based on the number
of images from each system, since not all systems have the same number of
images and a balanced distribution of images across the training, validation,
and test sets is necessary.

With this distribution, the training set for the neural network amounts to
$74.16\%$ of the total images with $99118$ images. The validation set was kept
for improving the results through the validation step and amounts $15.08\%$ of
the total $20156$ images. And finally, the test set checked the prediction
accuracy of the network with the remaining $10.76\%$ consisting of $14384$
basins.

The choice of parameters and systems influences the distribution of the
metrics values. However, as can be observed in Fig. 3, in spite of these
uneven metrics, the distributions of the three sets are quite similar.
Therefore the neural network is able to train accordingly and shows precise
results.

With regards to the training of the different neural networks, an Nvidia
Quadro RTX 5000 GPU was used. The weights and biases of the networks were
randomly initialized and adjusted during the training through the optimizer
provided by the tensorflow library Adam Kingma and Ba (2015) (Adaptive Moment
Estimation Method). This optimizer was used with the default parameters:
learning rate $=0.001$, $\beta_{1}=0.9$, $\beta_{2}=0.999$,
$\epsilon=10^{-8}$, $\text{decay}=0$, $amsgrad$ = true.

We trained each CNN for $100$ epochs while setting a batch size of $16$
images. Beyond this number of epochs we have noticed that the CNNs started to
overfit. During this process we have measured the loss function of the
validation image set as the mean square error (mse), when estimating the
fractal dimension, the basin entropy, and the boundary basin entropy. This
loss function measures the square of the difference between the predictions
and the true values, and hence tends to be minimized. The mathematical
expression is given by:

$\text{mse}=\frac{1}{N}\sum_{t=1}^{N}\left(\hat{x}_{t}-x_{t}\right)^{2},$ (1)

where $N$ is the total number of observations, $\hat{x_{t}}$ is the predicted
value, and $x_{t}$ the expected value.

Figure 2: Schematic representation of the CNNs architectures used in this
article. Each architecture is built as the authors of the corresponding papers
explain Krizhevsky, Sutskever, and Hinton (2012); Simonyan and Zisserman
(2015); Szegedy  _et al._ (2015); He  _et al._ (2016), except for the fully
connected layers. Since the authors proposed 4096 neuron layers but given our
image resolution it would vastly exceed memory usage, so these dense layers
were reduced to 2048 neuron layers and removed for performance reasons in the
case of the VGG architectures. The output layers also denote fully connected
layers, being 1-neuron layers for the estimation of the fractal dimension, the
basin entropy, and boundary basin entropy, whereas for the verification of the
Wada property we use an output of 2 neurons. The inception blocks,
convolutional blocks, and identity blocks that appear on the GoogLeNet and
ResNet architectures denote blocks of convolutions that occur on those
architectures, and the exact process is shown at the bottom of the image. It
is also worth mentioning that the main layers of each CNN appear in the
representation, however the normalization layers, and activation layers with
the rectified linear unit activation function are not shown in the schematic
representation. Nevertheless, these layers always appear in that order after
each convolutional layer. Another layer that also does not appear in the
representation are the dropout layers after each fully connected layer other
than the output one, each dropout layer has been established with a value of
$p=0.5$. Finally, in the case of the GoogLeNet architecture, 3 outputs appear
at different depths, and therefore each basin is predicted 3 values over a
metric. In our work, the best performing branch of the architecture is the
deepest one for all the studied metrics, hence all the results shown later
from this architecture belong to the output of the deepest branch.

Figure 3: Distribution of values of each metric on all image sets. The blue
histogram denotes the training set, the orange denotes the validation image
set, and the green denotes the test image set. It can be seen that the
histograms do not follow a uniform distribution, which indirectly introduces a
bias during training. However, by having the three sets with a similar
distribution, the bias is minimized in our results.

When checking the Wada property, the mean square error cannot be used as a
loss function, since basins are labeled as either Wada or not Wada. A more
useful loss function for this case is the categorical cross-entropy (cce).
This loss function measures the difference between the predicted probability
distribution and the true probability distribution of the target class, given
by

$\text{cce}=-\frac{1}{N}\sum_{t=1}^{N}\sum_{c=1}^{C}y_{tc}\log\left(\hat{y}_{tc}\right).$
(2)

From this expression, $N$ is the total number of observations, $C$ the total
number of classes on the set (Wada or not Wada), $y_{tc}$ is the observation
belonging to class $c$, and finally $\hat{y}_{tc}$ is the predicted
observation $t$ belonging to class $c$.

Once the training is finished, the final set of weights and biases that
performed best for the validation set were chosen as the final result of the
training. All the code for the training and use of the CNNs can be found in:
Github repository: https://github.com/RedLynx96/Deep-Learning-based-Analysis-
of-Basins-of-Attraction.

##  IV Results and Performance

The most precise results of the measured metrics in each difference
architectures are shown in Table 1. They show the mean square error for the
fractal dimension (FDim), the basin entropy (Sb), and the boundary basin
entropy (Sbb), while for the Wada property we show the failure probability,
which indicates the percentage of basins wrongly labeled. All these
predictions are made on images of the validation set which, as specified in
Sec. III, consist of basins from the forced damped pendulum and the Hénon-
Heiles Hamiltonian system.

Architecture | FDim $(10^{-4})$ | Sb $(10^{-4})$ | Sbb $(10^{-4})$ | Wada $(\%)$
---|---|---|---|---
AlexNet | 1.756 | 4.258 | 2.103 | 2.44
VGG16 | 3.668 | 0.863 | 3.980 | 1.95
VGG19 | 3.692 | 1.408 | 3.406 | 1.33
GoogLeNet | 1.015 | 1.024 | 5.924 | 2.39
ResNet50 | $\mathbf{0.464}$ | $\mathbf{0.240}$ | $\mathbf{0.463}$ | $\mathbf{0.84}$
Table 1: Performance of each architecture for estimating the desired metrics
on the validation dataset. The metric shown in the table for the fractal
dimension (FDim), the basin entropy (Sb), and the boundary basin entropy (Sbb)
is the mean square error value at the end of the training for the validation
set, whereas for the Wada property we show the failure probability. It can be
seen from the table that the ResNet50 architecture is the most accurate one
for all the metrics, as it has the least associated error.

As shown in the table, the ResNet50 architecture exhibits the lowest error. It
can also be inferred in the cases of the estimation of the mean squared error,
that the error of the neural networks increases with deeper network
architecture except in the cases of GoogLeNet, and ResNet50. In the case of
GoogLeNet, note that this architecture has three outputs and therefore each
basin is predicted three times over a metric. Overall, the best performing
branch of the architecture

Figure 4: Performance of the ResNet50 architecture at predicting the
$\mathbf{14384}$ basins of attraction from the test set. The linear
regressions show the comparison between the true values and the predicted
values for each metric studied in our work. Each regression is followed in its
right side with a histogram where the difference of the true and predicted
values is shown. Such histogram approximates to a Gaussian distribution with
mean $\mu$ and standard deviation $\sigma$. Since the estimation of the Wada
property is a classification task, a more practical representation is the
confusion matrix between the predictions and the true values. Such confusion
matrix is shown in the bottom right of the image, where it is shown how the
ResNet50 is able to correctly classify $99.75\%$ of basins, and therefore
proving the efficiency of the ResNet50 architecture.

is the deepest one, and therefore all the results shown in this table belong
to the output of the deepest branch. In Ref. He  _et al._ (2016), the authors
argued how deep convolutional neural networks make good predictions up to a
certain point, where deeper networks involve more error in the predictions.
This is the key point of ResNet, where a convolutional neural network with
residual terms avoids this problem.

Regarding the performance of our CNNs on the testing set, we used $14384$
images from the basins of attraction of the magnetic pendulum. We ran this
image set through all our CNNs to assess their performance. To achieve this,
we have computed all the metrics of each basin using all of our CNNs and
compared the predicted values with the true values established by the
reference methods. We have obtained the error in the prediction of each metric
for a given basin by subtracting the predicted and true values, and then
calculated the average error and its standard deviation for the $14384$
images. In Fig. 4, we show the predictions of our yet to prove best performing
architecture: the ResNet50. We compare the predicted and true values of all
the metrics in the test image basins, where the ideal case is to obtain a
straight line which would indicate that the predicted and true values are the
same. These linear regressions are followed by a histogram showing the
distribution error of each prediction. With regards to the Wada property, a
more useful approach is shown in the bottom right of Fig. 4. The confusion
matrix reflects how each basin of attraction has been classified and how it
should have been classified, being the elements of the diagonal correctly
classified and those outside wrongly classified.

Finally, speed and performance tests for each CNN are shown in Table 2, where
we show every architecture tested alongside the mean value error in the
predictions $(\mu)$, the standard deviation of the errors $(\sigma)$, and the
time $t$ required to characterize each basin from the testing set $(t)$. These
results are also compared to the reference methods used to label our ground
truths. As discussed in Sec. II, each metric is computed $10$ times for a
single basin from where the mean value and standard deviation are calculated.

The results in Table 2 show that the architecture with the fastest computation
time is AlexNet, which makes sense since it is the simplest architecture among
those compared. The next fastest architecture is ResNet50, and even though it
is not the fastest architecture, it is the one with the minimal error in its
predictions. Therefore, we consider it to be the optimal architecture for
performing this task, as the relationship between computation time and
precision is the most advantageous among all the proposed architectures.

---
Table 2: Performance comparison among CNN architectures. Each architecture
predicts the fractal dimension (FDim), the basin entropy (Sb), the boundary
basin entropy (Sbb), and the presence of the Wada property (Wada) for the test
set of $14384$ images of basins of attraction from the magnetic pendulum. The
predicted values are compared to the ground truth values established by the
reference methods explained in Sec. II, giving the mean ($\mu$) and standard
deviation ($\sigma$) error for each architecture. The time ($t$) each
architecture takes to predict all images is also included. It can be seen how
ResNet50 has the lowest prediction error and is faster than all architectures
except for AlexNet. However, ResNet50 still outperforms AlexNet in terms of
precision, making it the best architecture for this task. Therefore, it can be
concluded that ResNet50 is the best performing architecture for this task.

Although the ResNet50 has a slightly higher error compared to the reference
methods, it provides a significant improvement in computation time. The
ResNet50 characterizes basins much faster, taking only 7 minutes and 25
seconds, compared to a total time of 42 hours 49 minutes 38 seconds using the
reference methods. This makes the ResNet50 approximately $\mathbf{342}$ times
faster for the calculation of the $14384$ images from basins of attraction of
the magnetic pendulum. Among the implemented solutions, this is the optimal
architecture for a massive characterization of basins of attraction, such as
exploring the evolution of a dynamical system.

One way to improve the accuracy and robustness of these results would be to
improve the training conditions and labeling of the basins. As shown in Fig.
3, although the histograms of the training, validation, and test sets have the
same shape in each of the metrics, the distribution should be uniform to
produce optimal results in the output of the network. To minimize the error in
the predictions, it would be necessary to compute more basins from different
dynamical systems so that the proportion remains balanced. However, achieving
uniformity in one metric does not necessarily lead to uniformity in all of
them, and to achieve this a greater number of basins and dynamical systems
would be necessary.

Another problem with the labeling is that the fractal dimension, the basin
entropy, and the boundary basin entropy are metrics calculated via the Monte
Carlo method. Hence, part of the error associated with these labels is due to
the random sampling conditions. One way to minimize this error has been done
by measuring the same magnitude several times and taking the mean as the true
value. Finally, regarding the Wada property, as we have explained before, we
believe that the main source of error in the estimation of the metric comes
from our ground truth. So a way to improve the utility of our work is labeling
these basins with a different method such as the Nusse-Yorke Nusse and Yorke
(1996b) or the saddle-straddle Wagemakers, Daza, and Sanjuán (2020), which
despite being more computationally expensive classifies basins in a more
robust way.

##  V Conclusions

We have trained a variety of CNNs using basins of attraction of the Duffing
oscillator, the Newton fractal, the forced damped pendulum, and the Hénon-
Heiles Hamiltonian system. These CNNs proved to be capable of characterizing
basins from complex dynamical systems quickly and accurately. Additionally, we
have compared the different architectures of CNNs in terms of speed,
precision, and robustness between them and against the reference
characterization methods. In doing so, we found that the ResNet50 architecture
was the most efficient for this task.

In Table 2, we can see the results of testing all the studied architectures
and the reference methods, which prove that the ResNet50 architecture is the
best for characterizing the $14384$ images from the basins of attraction of
the magnetic pendulum. We can observe that the ResNet50 error is slightly
higher than the reference algorithms error for the studied metrics.

We finally conclude that from all the CNNs architectures used in this work,
the optimal architecture for this problem is the ResNet50 due to its precision
and speed. Even though all the methods proposed here need a previously
computed basin of attraction, the heavy computational time required for a
massive characterization of basins when exploring one or more dynamical
systems, can be easily minimized by the deep learning tool. The trade-off
between computing speed and loss of precision is worth it for this task.

##  Appendix A Dynamical systems parameters

We discuss the algorithms and the choice of parameters used for the
computation of the basins from the different dynamical systems that we had
discussed. All the basins and metrics were computed using an Intel Core
i7-10700 CPU 2.90 GHz. The basins of attraction for the Duffing oscillator and
the forced damped pendulum were computed numerically on the Julia programming
language using a standard RK45 integrator and the method proposed by Datseris
and Wagemakers in Datseris and Wagemakers (2022). On the other hand, basins
for the Hénon-Heiles Hamiltonian system, the magnetic pendulum, and the Newton
fractal were computed numerically on C++ using a standard 4th order Runge-
Kutta method. These basins of attraction were initially computed with a
resolution of $1000\times 1000$ pixels. However, to improve the size of the
image set, each basin was split into 9 basins of $333\times 333$ pixels, and
also the original basin was downsampled to that resolution, thus obtaining a
total of $10$ basins of $333\times 333$ pixels from an original $1000\times
1000$ basin.

As for the Duffing image set, we have computed $84842$ images corresponding to
basins of the Duffing oscillator defined by

$\ddot{x}+0.15\dot{x}-x+x^{3}=\gamma\cos\omega t.$ (3)

We have considered the ranges $0.1\leq\gamma\leq 0.5$ and $0.2\leq\omega\leq
2.5$. This choice of parameters in the differential equation describes the
motion of a unit mass particle in a double-well potential of the form
$V(x)=x^{4}/4-x^{2}/2$ , where the particle also experiences a linear drag
with the environment $(0.15)$ and a external periodic forcing
($\gamma\cos\omega t$).

Regarding the forced damped pendulum, $18926$ images of basins of attraction
were obtained from different parameters of the dynamical system. The system is
defined by the following equation

$\ddot{\theta}+0.2\dot{\theta}+\sin\theta=F\sin\omega t.$ (4)

Here, we have considered $0.5\leq\omega\leq 2.0$ and $0.5\leq F\leq 2.0$. This
choice of parameters in the differential equation covers a variety of basins
of attraction, such as simple non-fractal basins, riddled basins, or Wada
basins.

For the Hénon-Heiles Hamiltonian system, we have computed $1230$ exit basins,
which are analogous to the basins of attraction for open Hamiltonian systems
Nieto2019 . The Hamiltonian associated with the Hénon-Heiles system is defined
by:

$H=\frac{1}{2}\left(p_{x}^{2}+p_{y}^{2}\right)+\frac{1}{2}\left(x^{2}+y^{2}\right)+x^{2}y-\frac{y^{3}}{3}.$
(5)

This system describes the motion associated with a star around a galactic
center Hénon and Heiles (1964), the parameters $x$, $y$, $p_{x}$ , $p_{y}$ are
associated with the positions and momenta of a particle in the system, and
$H=E$ is a constant that represents the total energy of the system. For the
computation of the exit basins, we have considered values of the energy
$0.5\leq E\leq 2.0$, where the initial launch position of the particle has
been modified $(x,y)$ and the associated moments $(p_{x},p_{y})$ have been
computed using the tangential shooting as in Aguirre, Vallejo, and Sanjuán
(2001).

Regarding the Newton’s fractal, $14276$ basins of attraction have been
obtained. The polynomial used to calculate the basins of attraction was a
fifth degree polynomial

$p(z)=a_{0}+a_{1}z+a_{2}z^{2}+a_{3}z^{3}+a_{4}z^{4}+a_{5}z^{5},$ (6)

where $z$ is an initial condition in the complex plane and $a_{n}$ is a real
random number that satisfies $0.0\leq a_{n}\leq 1.0$. Each iteration of the
initial condition $z$ is performed following the rule

$z_{n+1}=z_{n}-b\frac{p\left(z_{n}\right)}{p^{\prime}\left(z_{n}\right)}.$ (7)

The basins of attraction have been computed for different initial conditions
of $z$, given that $-2.5\leq z\leq 2.5$, for both real and imaginary
parameters. In addition, the parameter $b$ has also been considered as a
complex number, such that $0.0\leq b\leq 1.0$ for both its real and imaginary
parts. Thus, the basins of attraction have been obtained by computing
different fifth-degree polynomials with random values for both $a_{n}$ and
$b$.

Finally, with regards to the magnetic pendulum $14384$ basins of attraction of
this system have been used. The system consists of a pendulum oscillating
above a plane with magnets that is also influenced by the air friction. The
differential equations that govern the movement of the pendulum in the system
are defined by

$\begin{split}&\ddot{x}+b\dot{x}+0.2x-\\\
&\sum_{i=1}^{n}\frac{x_{i}-x}{\left(\left(x_{i}-x\right)^{2}+\left(y_{i}-y\right)^{2}+0.2^{2}\right)^{3/2}}=0\\\
\\\ &\ddot{y}+b\dot{y}+0.2y-\\\
&\sum_{i=1}^{n}\frac{y_{i}-y}{\left(\left(x_{i}-x\right)^{2}+\left(y_{i}-y\right)^{2}+0.2^{2}\right)^{3/2}}=0.\end{split}$
(8)

These differential equations describe the position of the pendulum $(x,y)$ at
a distance $0.2$ above the plane in which $n$ magnets are located, the
position of these magnets is given by the coordinates $(x_{i},y_{i})$ as they
are placed having all the same separation between them and at the same
distance $a$ from the center of the plane, which attracts the pendulum due to
gravity and being represented with the scalar $0.2$. We also assume that the
pendulum has linear drag $b$, and that the magnetic force between a magnet and
the pendulum is proportional to $1/r^{3}$.

To compute the basins of attraction of this system we have considered
$-0.01\leq b\leq 0.25$, and the distance from the center to the magnets
$1.5\leq a\leq 3.0$. We have computed basins with this variation of parameters
for the case in which $2,3$ and $4$ magnets appear in the system.

## ACKNOWLEDGMENTS

This work has been financially supported by the Spanish State Research Agency
(AEI) and the European Regional Development Fund (ERDF, EU) under Project No.
PID2019-105554GB-I00 (MCIN/AEI/10.13039/501100011033).

## References

  * Aguirre, Vallejo, and Sanjuán (2001) Aguirre, J., Vallejo, J. C., and Sanjuán, M. A. F., “Wada basins and chaotic invariant sets in the Hénon-Heiles system,” Phys. Rev. E 64 (2001).
  * Aguirre, Viana, and Sanjuán (2009) Aguirre, J., Viana, R. L., and Sanjuán, M. A. F., “Fractal structures in nonlinear dynamics.” Rev. Mod. Phys. 81, 333 (2009).
  * Datseris and Wagemakers (2022) Datseris, G. and Wagemakers, A., “Effortless estimation of basins of attraction,” Chaos 32, 023104 (2022).
  * Daza  _et al._ (2016) Daza, A., Wagemakers, A., Georgeot, B., Guéry-Odelin, D., and Sanjuán, M. A. F., “Basin entropy: a new tool to analyze uncertainty in dynamical systems,” Sci. Rep. 6, 31416 (2016).
  * Daza, Wagemakers, and Sanjuán (2018) Daza, A., Wagemakers, A., and Sanjuán, M. A. F., “Ascertaining when a basin is wada: the merging method,” Sci. Rep. 8, 9954 (2018).
  * Fan  _et al._ (2020) Fan, H., Jiang, J., Zhang, C., Wang, X., and Lai, Y.-C., “Long-term prediction of chaotic systems with machine learning,” Phys. Rev. Research 2, 012080 (2020).
  * Goodfellow, Bengio, and Courville (2016) Goodfellow, I., Bengio, Y., and Courville, A.,  _Deep Learning_ (MIT Press, 2016).
  * He  _et al._ (2016) He, K., Zhang, X., Ren, S., and Sun, J., “Deep residual learning for image recognition,” in  _Proceedings of the IEEE conference on computer vision and Pattern Recognit._ (2016) pp. 770–778.
  * Hénon and Heiles (1964) Hénon, M. and Heiles, C., “The applicability of the third integral of motion: Some numerical experiments,” The Astron. J. 69, 73 (1964).
  * Kingma and Ba (2015) Kingma, D. P. and Ba, J., “Adam: A method for stochastic optimization,” in  _3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings_ , edited by Y. Bengio and Y. LeCun (2015).
  * Kirichenko, Radivilova, and Bulakh (2019) Kirichenko, L., Radivilova, T., and Bulakh, V., “Machine learning in classification time series with fractal properties,” Data 4, 5 (2019).
  * Krizhevsky, Sutskever, and Hinton (2012) Krizhevsky, A., Sutskever, I., and Hinton, G. E., “Imagenet classification with deep convolutional neural networks,” in  _Advances in Neural Information Processing Systems 25_ , edited by F. Pereira, C. J. C. Burges, L. Bottou, and K. Q. Weinberger (Curran Associates, Inc., 2012) pp. 1097–1105.
  * Lecun  _et al._ (1998) Lecun, Y., Bottou, L., Bengio, Y., and Haffner, P., “Gradient-based learning applied to document recognition,” Proceedings of the IEEE 86, 2278–2324 (1998).
  * Lu, Hunt, and Ott (2018) Lu, Z., Hunt, B. R., and Ott, E., “Attractor reconstruction by machine learning,” Chaos 28, 061104 (2018).
  * McDonald  _et al._ (1985) McDonald, S. W., Grebogi, C., Ott, E., and Yorke, J. A., “Fractal basin boundaries,” Physica D 17, 125–153 (1985).
  * Nusse and Yorke (1996a) Nusse, H. and Yorke, J., “Basins of attraction,” Science 271, 1376–1380 (1996a).
  * Nusse and Yorke (1996b) Nusse, H. E. and Yorke, J. A., “Wada basin boundaries and basin cells,” Physica D 90, 242–261 (1996b).
  * Pathak  _et al._ (2018) Pathak, J., Wikner, A., Fussell, R., Chandra, S., Hunt, B. R., Girvan, M., and Ott, E., “Hybrid forecasting of chaotic processes: Using machine learning in conjunction with a knowledge-based model,” Chaos 28, 041101 (2018).
  * Peitgen, Jürgens, and Saupe (1992) Peitgen, H.-O., Jürgens, H., and Saupe, D.,  _Chaos and Fractals: New Frontiers of Science_ (Springer, New York, NY, 1992).
  * Shi (2018) Shi, C.-T., “Signal pattern recognit. based on fractal features and machine learning,” Applied Sciences 8, 1327 (2018).
  * Simonyan and Zisserman (2015) Simonyan, K. and Zisserman, A., “Very deep convolutional networks for large-scale image recognition,” in  _3rd International Conference on Learning Representations_ (2015) pp. 1–14.
  * Srivastava  _et al._ (2014) Srivastava, N., Hinton, G., Krizhevsky, A., Sutskever, I., and Salakhutdinov, R., “Dropout: A simple way to prevent neural networks from overfitting,” Journal of Machine Learning Research 15, 1929–1958 (2014).
  * Szegedy  _et al._ (2015) Szegedy, C., Liu, W., Jia, Y., Sermanet, P., Reed, S., Anguelov, D., Erhan, D., Vanhoucke, V., and Rabinovich, A., “Going deeper with convolutions,” in  _Proceedings of the IEEE conference on computer vision and Pattern Recognit._ (2015) pp. 1–9.
  * Szeliski (2011) Szeliski, R.,  _Computer vision algorithms and applications_ (Springer, New York, NY, 2011).
  * Tang  _et al._ (2020) Tang, Y., Kurths, J., Lin, W., Ott, E., and Kocarev, L., “Introduction to focus issue: When machine learning meets complex systems: Networks, chaos, and nonlinear dynamics.” Chaos 30, 063151 (2020).
  * Valle  _et al._ (2022) Valle, D., Wagemakers, A., Daza, A., and Sanjuán, M. A. F., “Characterization of fractal basins using deep convolutional neural networks,” Int J Bifurcat Chaos 32, 2250200 (2022).
  * Wagemakers, Daza, and Sanjuán (2020) Wagemakers, A., Daza, A., and Sanjuán, M. A. F., “The saddle-straddle method to test for wada basins,” Commun. Nonlinear Sci. Numer. Simul. 84, 105167 (2020).
  * (28) Daza, A., Wagemakers, A., and Sanjuán, M. A. F., “Classifying basins of attraction using the basin entropy”. Chaos, Solitons and Fractals, 159, 112112 (2022).
  * (29) Nieto, A.R., Zotos, E.E., Seoane, J.M., and Sanjuán, M. A. F., “Measuring the transition between nonhyperbolic and hyperbolic regimes in open Hamiltonian systems.” Nonlinear Dyn 99, 3029–3039 (2020).