# Devolatilization of extrasolar planetesimals by ${}^{60}\text{Fe}$ and
${}^{26}\text{Al}$ heating

Joseph W. Eatson1, Tim Lichtenberg2, Richard J. Parker1 & Taras V. Gerya3
1Department of Physics and Astronomy, The University of Sheffield, Hicks
Building, Hounsfield Road, Sheffield, S3 7RH, UK
2Kapteyn Astronomical Institute, University of Groningen, P.O. Box 800, 9700
AV Groningen, NL
3Institute of Geophysics, Department of Earth Sciences, ETH Zurich,
Sonneggstrasse 5, 8092 Zurich, CH
E-mail<EMAIL_ADDRESS>Society Dorothy Hodgkin Fellow

###### Abstract

Whilst the formation of Solar system planets is constrained by meteoritic
evidence, the geophysical history of low-mass exoplanets is much less clear.
The bulk composition and climate states of rocky exoplanets may vary
significantly based on the composition and properties of the planetesimals
they form from. An important factor influenced by planetesimal composition is
water content, where the desiccation of accreting planetesimals impacts the
final water content of the resultant planets. While the inner planets of the
Solar system are comparatively water-poor, recent observational evidence from
exoplanet bulk densities and planetary formation models suggest that rocky
exoplanets engulfed by substantial layers of high-pressure ices or massive
steam atmospheres could be widespread. Here we quantify variations in
planetesimal desiccation due to potential fractionation of the two short-lived
radioisotopes ${}^{26}\text{Al}$ and ${}^{60}\text{Fe}$ relevant for internal
heating on planetary formation timescales. We focus on how order of magnitude
variations in ${}^{60}\text{Fe}$ can affect the water content of
planetesimals, and how this may alter the formation of extrasolar ocean
worlds. We find that heating by ${}^{26}\text{Al}$ is the dominant cause of
planetesimal heating in any Solar system analogue scenario, thus validating
previous works focussing only on this radioisotope. However,
${}^{60}\text{Fe}$ can become the primary heating source in the case of high
levels of supernova enrichment in massive star-forming regions. These
diverging scenarios can affect the formation pathways, bulk volatile budget,
and climate diversity of low-mass exoplanets.

###### keywords:

exoplanets – planets and satellites: terrestrial planets – planets and
satellites: interiors – planets and satellites: composition – planets and
satellites: formation

††pubyear: 2023††pagerange: Devolatilization of extrasolar planetesimals by
${}^{60}\text{Fe}$ and ${}^{26}\text{Al}$ heating–References

##  1 Introduction

Astronomical observations (Lacedelli et al., 2022; Luque & Pallé, 2022;
Diamond-Lowe et al., 2022; Piaulet et al., 2023) and formation models (Zeng et
al., 2019; Venturini et al., 2020b) have recently opened new frontiers in
characterising low-mass exoplanets, highlighting a potential deviation between
volatile-rich and under-dense versus volatile-poor and mostly rocky worlds
(Dorn & Lichtenberg, 2021; Zilinskas et al., 2023; Piette et al., 2023). In
the early Solar System, short-lived radioisotopes (SLRs) sensitively affected
the thermal evolution, differentiation, and volatile composition of
planetesimals and forming planets (Fu et al., 2017; Monteux et al., 2018;
Lichtenberg et al., 2023; Krijt et al., 2022), which can shape system-to-
system deviations in exoplanet volatile trends (Lichtenberg et al., 2019a;
Lichtenberg & Krijt, 2021).

Substantial internal heating of planetesimals by SLRs, in particular
${}^{26}\text{Al}$, in the Solar System was first theorised by Urey (1955).
Evidence of a substantial ${}^{26}\text{Al}$ enrichment in the early Solar
System in the form of an excess of the ${}^{26}\text{Al}$ decay product
${}^{26}\text{Mg}$ in calcium–aluminium-rich inclusions (CAIs) was discovered
by Gray & Compston (1974) and Lee et al. (1976). With such levels of
enrichment, ${}^{26}\text{Al}$ was shown to be a dominant heating mechanism
for planetesimals and asteroidal precursor bodies (Grimm & McSween, 1989,
1993). The excess levels of ${}^{60}\text{Fe}$ have been less clear from
meteoritic data (Tang & Dauphas, 2012; Telus et al., 2018; Trappitsch et al.,
2018; Cook et al., 2021), with estimates varying over two orders of magnitude
(from ${}^{60}\mathrm{Fe}/^{56}\mathrm{Fe}=10^{-8}$ to $10^{-6}$) in recent
years, with the lower ${}^{60}\text{Fe}$ levels favoured in most recent
measurements. At the lower end limit, the abundances of ${}^{60}\text{Fe}$ are
not expected to be sufficient to cause substantial heating and volatile loss
of planetesimals, however, substantial variations of the
${}^{26}\text{Al}$/${}^{60}\text{Fe}$ ratio are expected across star-forming
regions (Vasileiadis et al., 2013; Boss, 2019, 2022; Fatuzzo & Adams, 2022;
Parker et al., 2023).

Multiple formation mechanisms exist for both ${}^{26}\text{Al}$ and
${}^{60}\text{Fe}$ (Lugaro et al., 2018; Diehl et al., 2021; Laird et al.,
2023). While spallation from stellar cosmic rays has been suggested as a
potential mechanism for enrichment of ${}^{26}\text{Al}$ (Adams, 2021), the
low levels of 10Be (McKeegan et al., 2000; Dunham et al., 2022) and mostly
homogenous distribution of decay products throughout the Solar System (Desch
et al., 2022) suggest that this mechanism is not the most influential (Parker,
2020). Mechanisms external to a stellar system are more probable, such as
through supernovae (Chevalier, 2000; Adams, 2010; Ouellette et al., 2010;
Parker et al., 2014; Pfalzner et al., 2015), stellar winds (Gaidos et al.,
2009), or AGB stars (Lugaro et al., 2018). Supernovae in particular produce
both ${}^{60}\text{Fe}$ and ${}^{26}\text{Al}$ in abundance, with usually
elevated ${}^{60}\text{Fe}$/${}^{26}\text{Al}$ ratios relative to the Solar
System (Kuffmeier et al., 2016; Lichtenberg et al., 2016), while stellar winds
from massive evolved stars (such as Wolf-Rayets) can produce significant
quantities of ${}^{26}\text{Al}$ over their lifetime (Limongi & Chieffi, 2006,
2018). As the distribution of almost all SLRs are approximately homogenous,
the half life times are short compared to the lifetime of the protoplanetary
disks, and ${}^{26}\text{Al}$-devoid CAIs (FUN-CAIs) are rare (Villeneuve et
al., 2009; Desch et al., 2023), it is likely that SLRs are deposited
immediately prior to or during the formation of the first CAIs. As such, we
can infer that the enrichment of SLRs in a stellar system is sensitive to the
local stellar environment immediately leading up to the formation of the
stellar system (Adams, 2010; Portegies Zwart, 2019; Parker, 2020).

Heating from SLRs has a significant impact on the volatile evolution of
initially ice-rich planetesimals (Castillo-Rogez & Young, 2017; Golabek &
Jutzi, 2021; Newcombe et al., 2023), which can have a knock-on effect on the
properties of subsequent planets (Grimm & McSween, 1993; Lichtenberg et al.,
2021). Of importance for this paper, planetesimal internal heating from the
decay of SLRs can result in melting, heating, vaporisation and subsequent
outgassing of volatiles – particularly water and volatile carbon compounds. In
the case of stellar systems with significant SLR enrichment, this can result
in a significant reduction in surface and internal water content, which is a
prime marker for bulk composition and affects the long-term climatic evolution
of rocky planets. Constraining the influence of radiogenic heating for
multiple radioisotopes on the volatile loss from planetesimals during
accretion would allow for improvements in planetary formation models, better
estimations of the ocean world population, and quantify the possible
fractionation between major atmophile elements across planetary populations
(Wang et al., 2019; Spaargaren et al., 2023).

Previous simulations of planetesimal devolatilization primarily focused on the
influence of ${}^{26}\text{Al}$ on planetesimal water and carbon compounds
(Lichtenberg et al., 2019a; Lichtenberg & Krijt, 2021; Lichtenberg & Clement,
2022). This work, in contrast, specifically focusses on the influence of
${}^{60}\text{Fe}$ and thus varying ${}^{60}\text{Fe}$/${}^{26}\text{Al}$
ratios on planetesimal desiccation. A parameter space exploration is performed
to determine how ${}^{60}\text{Fe}$ heating can affect planetesimals of
varying sizes, iron contents and degree of ${}^{60}\text{Fe}$ enrichment
relative to the Solar System. Internal core-mantle differentiation and thus
${}^{60}\text{Fe}$ distribution is simulated by use of two distinct internal
structure models, which are discussed in Section 2. As there are multiple
parameters to be explored, the number of simulations was extremely large, as
such, 2D models with a temperature-based rather than a chemistry-based model
were utilised. Additionally, constraining some of these parameters for
subsequent simulation sets was the basis for a number of simulation subsets,
which are described in Section 3. Finally, the results are discussed and
conclusions are made in Section 4 & 5.

##  2 Methodology

###  2.1 Numerical & heating simulation

(a) A comparison of retained water fraction over time for a number of
simulations of different resolutions. (b) A comparison of final water fraction
over simulation cell count. The final result is insensitive to simulation
resolution.

Figure 1: A comparison of H2O desiccation for a set of simulations with
varying resolution. These simulations were conducted to determine whether
simulation resolution impacted the results of the simulations. We found that
there was no significant change as resolution varies, as such a value of
$250^{2}$ cells was utilised as a good compromise between speed and accuracy.

Numerical simulations of planetesimals were conducted using the I2ELVIS
thermo-mechanical geodynamic
code111https://github.com/FormingWorlds/i2elvis_planet, which uses a
combination of the marker-in-cell and conservative finite difference schemes
to simulate planetesimal evolution (Gerya & Yuen, 2003, 2007), with a similar
setup as in Lichtenberg et al. (2019a) and Lichtenberg et al. (2021).
Simulations were conducted with a fixed grid resolution of $250\times 250$
cells. Higher resolutions were tested, but found to produce similar results at
the expense of significantly increased processing time (Fig. 1). Simulations
evolve using a variable timestep, which is constrained by the following
parameters:

  * •

maxxystep: The maximum displacement distance for material in a cell ($\partial
x/\partial t\cdot\Delta t$).

  * •

maxtkstep: The maximum temperature change per step.

  * •

maxtmstep: The maximum total timestep.

Additionally, the maximum timestep is set to $2500$ years throughout the
simulation, in order to prevent numerical instability.

Upon initialisation I2ELVIS assigns properties to cells in order to make up
the planetesimal and its surrounding environment. These properties are
analogous to geological compositions within the planetesimal:

  * •

“Sticky air”, the surrounding medium that ensures thermal stability, a proxy
for the medium surrounding the planetesimal (Crameri et al., 2012).

  * •

Hydrous silicates which contain water and the simulation’s supply of
${}^{26}\text{Al}$.

  * •

Iron, which contains the simulation’s supply of ${}^{60}\text{Fe}$.

Hydrous silicates have three important states in a simulation that determine
the water content of the planetesimal. Initially water in the wet silicates
are assumed to be stored as ice, which can undergo phase changes based on
temperature. Beyond a temperature of $273\text{\,}\mathrm{K}$
($T_{\text{melt}}$) the ice melts, which delays further temperature increase
as the enthalpy of fusion needs to be overcome. This phase change is
reversible, and re-freezing can occur if the cell falls below the threshold
temperature. Beyond a temperature of $1223\text{\,}\mathrm{K}$, the upper
limit of the amphibolite stability field, the entire water content of the cell
becomes totally and irreversibly vaporised, and is assumed to leave the
planetesimal through out-gassing processes (though these are not simulated);
this temperature is referred to as $T_{\text{vap}}$. Water retention across
the planetesimal is stored as a fraction of the form:

$\mathcal{H}=\frac{f_{\text{H}_{2}\text{O},t}}{f_{\text{H}_{2}\text{O},i}}\equiv\frac{N_{\text{wet}}}{N_{\text{dry}}},$
(1)

where $f_{\text{H}_{2}\text{O},t}$ is the volume of hydrous silicates that has
never been heated above $T_{\text{vap}}$, $f_{\text{H}_{2}\text{O},i}$ is the
initial volume of hydrous silicates, $N_{\text{wet}}$ is the number of hydrous
silicate cells still containing water in the liquid or ice phases, and
$N_{\text{dry}}$ is the number of hydrous silicate cells that have been
removed of water. This ratio is calculated and stored for every timestep.

An important variable discussed is the final retained water fraction,
$\mathcal{H}_{f}$, which is the value of $\mathcal{H}$ at the end of the
simulation
($\mathcal{H}_{f}=f_{\text{H}_{2}\text{O},f}/f_{\text{H}_{2}\text{O},i}$).
Measuring the exact volume of remaining water introduces a significant number
of additional parameters, and as such is beyond the scope of this project.
This method can be performed for a number of other chemicals, such as CH4,
CO2, or CO (Lichtenberg & Krijt, 2021), but for simplicity we focus here on
water as the main volatile.

Heating of the planetesimal is provided by the decay of ${}^{26}\text{Al}$ and
${}^{60}\text{Fe}$ isotopes. The specific heating of an SLR is given by the
formula:

$H_{\text{SLR}}(t)=\frac{N_{\text{A}}E_{\text{SLR}}\lambda}{m_{\text{SLR}}}e^{-\lambda
t},$ (2)

where $N_{\text{A}}$ is Avogadro’s
constant222$6.022\,140\,76\text{\times}{10}^{23}$, $E_{\text{SLR}}$ is the SLR
decay energy, $\lambda$ is the decay constant of the SLR
($\lambda=\ln(2)/\tau_{\text{1/2}}$) and $m_{\text{SLR}}$ is the molar mass of
the SLR. The specific heating rate for SLRs embedded within a rocky body can
be calculated with the equation:

$Q_{\text{SLR}}(t)=f_{\text{E,CI}}Z_{\text{SLR}}H_{\text{SLR}}(t),$ (3)

where $f_{\text{E}}$ is the chrondritic elemental abundance and
$Z_{\text{SLR}}$ is the radioisotopic enrichment (Castillo-Rogez et al.,
2009). Enrichment is measured relative to the early Solar System estimates of
${}^{60}\text{Fe}$ and ${}^{26}\text{Al}$ enrichment; the Solar System
enrichment values for ${}^{26}\text{Al}$ and ${}^{60}\text{Fe}$ are detailed
in Table 1. The total heating rate can then be calculated through the formula:

$Q_{\text{T}}(t)=Q_{\text{26Al}}(t)+Q_{\text{60Fe}}(t),$ (4)

where $Q_{\text{26Al}}$ and $Q_{\text{60Fe}}$ are the individual heating
contributions per unit mass of ${}^{26}\text{Al}$ and ${}^{60}\text{Fe}$.

Fig. 2 shows decay heating from ${}^{26}\text{Al}$ and ${}^{60}\text{Fe}$
sources over time, as well as specific heating rates for the SLRs with Solar
System abundances. As can be seen, ${}^{26}\text{Al}$ specific heating is
significantly higher for the first 2 million years after CAI formation,
however this is supplanted by ${}^{60}\text{Fe}$ after this point. In the case
of Solar System abundances, however, ${}^{26}\text{Al}$ heating dominates due
to the lower quantity of ${}^{60}\text{Fe}$ over the entire period of
planetesimal formation and cooling. In other star systems that had
interactions with supernovae during their formation we may find a greater
abundance of ${}^{60}\text{Fe}$, which would impact radiogenic heating and the
water budget of subsequently forming planets.

Throughout this paper the SLR enrichment is described relative to the
estimated isotopic enrichment of the Solar System after CAI formation. This is
parametrised as the SLR enrichment ratio, $\Lambda$.

$\Lambda_{\text{SLR}}=\frac{Z_{\text{SLR},\star}}{Z_{\text{SLR},\odot}},$ (5)

where $Z_{\text{SLR},\star}$ is the post-CAI formation SLR enrichment of the
system described in the simulation and $Z_{\text{SLR},\odot}$ is the
equivalent for the Solar System. Table 1 notes the enrichment values for each
for Solar System estimates. Simulations are run for a total of
$10\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$, we assumed that the planetesimal
being simulated is formed at $1\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$ after
formation of the first CAIs. As such
$1\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$ of radioactive decay is simulated
prior to the start of hydrodynamical simulations to reflect this.

(a) Specific heating, based on Equation 3. (b) Planetary body specific heating
rate, based on Equation 4.

Figure 2: A comparison of the specific heating, $H_{\text{SLR}}$, and specific heating rates, $Q_{\text{SLR}}$, of the SLRs ${}^{26}\text{Al}$ and ${}^{60}\text{Fe}$ at Solar System abundances and enrichments. Whilst the specific heating rate of ${}^{60}\text{Fe}$ remains higher for longer, Solar System abundances result in a minimal impact on the total heating rate. The tables are produced using Eq. 2 and Eq. 3 using parameters detailed in Table 1. Isotope | $f$ | $\tau$ | $E$ | $Z_{\text{SLR},\odot}$
---|---|---|---|---
${}^{26}\text{Al}$ | 0.0085 | $0.717\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$ | $3.210\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$ | $5.250\text{\times}{10}^{-5}$
${}^{60}\text{Fe}$ | 0.1828 | $2.600\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$ | $2.712\text{\,}\mathrm{M}\mathrm{e}\mathrm{V}$ | $1.150\text{\times}{10}^{-8}$
Table 1: SLR properties, with the associated Solar System SLR-to-stable
isotope ratios (Kita et al., 2013; Tang & Dauphas, 2012).

###  2.2 Assumptions & limitations

Our simulations are based on static planetesimals that do not accrete further,
meaning that ongoing accretion of, for example, varying composition of pebbles
that grow planetesimals is not taken into account (Lichtenberg et al., 2021;
Johansen et al., 2023). The addition of accreting pebble layers on top of
planetesimals has been shown to alter their thermal evolution (Sturtz et al.,
2022). Furthermore, we do not model geochemical reactions operating in the
host rock, such as the oxidation of iron metal by water flow. This is an
important aspect, as the diversity of FeO and hydrogen-bearing phases in inner
Solar System asteroids bear evidence of rapid water-rock reactions during
planetary formation (Sarafian et al., 2017; Monteux et al., 2018; McCubbin &
Barnes, 2019; Grant et al., 2023), and may influence the final volatile
fractionation on planetesimals. Following up from this work, we will thus
investigate the more detailed multi-phase dynamics of water-rock aggregates
(Travis & Schubert, 2005; Bland & Travis, 2017). This is a substantial
simplification as mobilization and transport of volatile ices, fluids, and
gases strongly influences their redistribution between forming metal core,
planetary mantle, and outgassed reservoir (Suer et al., 2023). Additionally,
the physical dislocation in multi-phase fluid treatments can differ
substantially if the initial aggregate consists of a mixture that becomes
comparable in ice and rock composition (Gerya, 2019). Zhang (2023) recently
suggested sublimation of volatile ice phases as an additional devolatilization
mechanism for ice-rich planetesimals that cross to the inside of the water
snow line (due to the snow line migrating outwards during the Class I stage of
the solar protoplanetary disk, Lichtenberg et al., 2021), further highlighting
the importance of multiphasic treatments of planetesimals. In a similar vein,
ongoing metal-silicate differentiation is not treated directly, but bracketed
by our assumptions on Fe distribution, which we discuss in the next
subsection.

###  2.3 ${}^{60}\text{Fe}$ distribution models

Shortly after the formation of the stellar system, subsequent rapid formation
of planetesimals, and initial radogenic heating, the hot, undifferentiated
liquidus phase material of the planetesimal begins to undergo differentiation
(Nimmo & Kleine, 2015). Segregation of silicates and iron occur on an
approximately similar timescale, which occur over a period of
$0.4\,$\mathrm{M}\mathrm{y}\mathrm{r}$$ to
$10\,$\mathrm{M}\mathrm{y}\mathrm{r}$$ (Lichtenberg et al., 2019b), with the
upper limit being the approximate timescale of our simulation. The segregation
time of the planetesimal material is dependent on the initial formation time,
degree of radiogenic heating (Dodds et al., 2021), and mode of core formation
(Walte et al., 2020). As ${}^{26}\text{Al}$ heating drives melting of
silicates, a strongly convective mantle forms, which can inhibit the growth of
the core (Neumann et al., 2012). As the planetesimal cools this convective
flow decreases and eventually the mantle solidifies, this occurs on a
timescale of $10\,$\mathrm{M}\mathrm{y}\mathrm{r}$$ to
$50\,$\mathrm{M}\mathrm{y}\mathrm{r}$$ (Neumann et al., 2018). There are two
qualitative types of core formation in planetesimals, which differ in
timescale: in percolative core formation, S-rich metals may form an
interconnected vain network, which can remove mantle metals quickyl via
gravitational drainage before the appearance of the first silicate melts
(Yoshino et al., 2003; Ghanbarzadeh et al., 2017). However, laboratory
experiments find that typically some metals remain stranded in the silicate
matrix (Bagdassarov et al., 2009; Cerantola et al., 2015; Walte et al., 2023),
which are then removed from the mantle once the silicates reach the
rheological transition. As the range of core formation times across these two
major modes intersects the epoch our simulations occur in, and a full model of
core formation and material segregation are beyond the scope of this work, two
distinct planetesimal structures are simulated in this paper:

  1. 1.

The _core model_ : Iron is contained in a large metal core surrounded by
silicates. This represents an idealised end-state for a pre-differentiated
planetesimal, in a scenario where metal-silicate separation operates rapidly
and early.

  2. 2.

The _grain model_ : Iron is randomly distributed throughout the planetesimal.
This represents the end-member scenario of metal core formation completing
slowly, after the main stage of internal devolatilization.

Fig. 3 illustrates the physical differences of these models, as well as their
model-dependent parameters. An additional model with a mantle containing iron
grains as well as a small, still forming core was considered, as it
represented the mid-point between the two models. However, simulations with
these models were not conducted as the results from the first two models were
very similar, as discussed in Section 3.2. In addition to the parameters
dictating iron abundance in each model, $\Psi$ and $\Phi$, the planetesimal
radius and isotopic enrichment are varied, in order to explore a broad
parameter space.

Figure 3: Comparison of the initial conditions of each model with an
approximately equivalent iron fraction of $N_{\text{Fe}}/N_{\text{T}}=0.25$.
The _core model_ distributes all iron (yellow) into a central core surrounded
by silicates (grey), with a planetesimal radius ratio of $\Psi$ (Eq. 6). The
_grain model_ randomly distributes iron throughout the planetesimal with a
cell fraction of $\Phi$ (Eq. 7).

####  2.3.1 Core model

The first model considered for this research consists of an iron core
surrounded by hydrous silicates. This model is comparatively easy to
implement, and required no modification of the underlying I2ELVIS code.
However, this model is not a realistic case for planetesimal formation while
still retaining water, as iron core formation by liquid iron alloy percolation
requires heating of a planetesimal above the dehydration temperature of
$1223\text{\,}\mathrm{K}$ to $1700\text{\,}\mathrm{K}$ (Neumann et al., 2012).
Whilst this case therefore does not represent and accurate account of how
metal and volatile parts of a planetesimal would behave, it can approximate
SLR repartitioning during core formation, which occurs in tandem with
dehydration. As such, it represents an idealised end-member case, and can be
compared with the successor grain model to determine how different
distributions of iron and its related ${}^{60}\text{Fe}$ heating throughout
the planetesimal body may affect the final water fraction.

The size of the core is controlled by the core radius ratio parameter, given
by the equation:

$\Psi=\frac{r_{\text{c}}}{r_{\text{pl}}},$ (6)

where $r_{\text{c}}$ is the core radius and $r_{\text{pl}}$ is the total
planetesimal radius. Throughout these simulations this parameter is varied
from 0.0 up to 0.99; values beyond this were considered redundant as this is
already an extreme value, and that there would be no hydrous silicate cells
left to measure desiccation from.

Whilst simple to execute, this model is less physically accurate than the
_grain model_ also described in this paper. Primarily, the degree of
differentiation for even large planetesimals is comparatively slow, and occurs
over a timescale of $\sim 1\,$\mathrm{M}\mathrm{y}\mathrm{r}$$ to $\sim
10\,$\mathrm{M}\mathrm{y}\mathrm{r}$$ after CAI formation. As such, a
differentiated body is unlikely for a large planetesimal so early after its
formation. However, with more extreme degrees of radiogenic heating from
${}^{60}\text{Fe}$ the process of iron melt segregation may be accelerated due
to faster melting of the body.

####  2.3.2 Grain model

The _grain model_ provides a more realistic simulation of the planetesimal
shortly after formation. Upon initialisation a percentage of the cells
containing hydrous silicates are converted at random to have an iron
composition instead. The model therefore produces planetesimals where core
formation has not occurred. A pseudorandom number generator is used to
determine which cells are converted, as true random numbers were not deemed
necessary.

As the percentage of cells affected by the change increases the amount of iron
in the planetesimal increases, we define the Fe grain volume ratio, $\Phi$, as
the ratio of the number of cells changed at initialisation to the number of
grains that were unaffected:

$\Phi=\frac{N_{\text{Fe}}}{N_{\text{S}}},$ (7)

where $N_{\text{Fe}}$ is the number of cells with the iron marker and
$N_{\text{S}}$ is the number of cells with the hydrous silicate marker upon
initialisation. This parameter is used as a stand-in for the iron content in
the planetesimal; similarly to the _core model_ , ${}^{60}\text{Fe}$
enrichment is a separate parameter.

This model offers significant improvements over the iron _core model_ in terms
of physical accuracy, as it offers a closer analogue to the post-formation
conditions and properties of a typical planetesimal. Additionally,
planetesimals under these conditions should evenly heat faster than
planetesimals with an iron core, due to the even distribution of
${}^{60}\text{Fe}$ throughout the rocky body of the planetesimal. However,
whilst more accurate than the iron _core model_ , the _grain model_ still has
some physical inaccuracies. In particular, iron and silicates cannot segregate
due to limitations in the one-phase fluid approximation, though below the iron
alloy melting temperature this is not a concern (Keller & Suckale, 2019; Zhang
et al., 2021).

###  2.4 Data recording & processing

2D outputs are performed every 50 timesteps, and contain the primitive
variables (density, velocity & pressure), temperature and radiogenic heating
rate. These outputs can then be used to calculate other values and can be
averaged for specific regions of the planetesimal. Additionally, hydrous
fraction, ice fraction and water fraction are dumped for the entire
planetesimal at every time step.

##  3 Results

For simulations where desiccation occurs we found that the evolution of the
planetesimals was broadly the same. As the planetesimal begins to heat up over
the first $10^{5}\,$\mathrm{y}\mathrm{r}$$ of the simulation some cells exceed
the melting temperature, before subsequently exceeding the vaporisation
temperature. Desiccation occurs rapidly at this point, before the planetesimal
begins to thermalise, then cool, preventing any further desiccation from
occurring. This progression can be seen in Fig. 4, where this radpid
desiccation can be observed.

Figure 4: _Core model_ comparison of water fraction, $\mathcal{H}$, as a
function of simulation time for planetesimals of radius
$50\text{\,}\mathrm{k}\mathrm{m}$ and an un-enriched iron core with a radius
ratio of $0.25$. Simulations follow a similar structure, with most desiccation
occurring over a very short period, before reaching a minimum as radiogenic
heating reduces.

The parameter space exploration of this paper was broken into a series of
simulations. Firstly, the _core model_ was utilised in order to get a baseline
of how ${}^{60}\text{Fe}$ heating affects water content, and its relative
heating impact compared to ${}^{26}\text{Al}$. This model was also used to
constrain the ideal size of planetesimal for subsequent simulations, and also
to determine how iron content can affect the efficiency of ${}^{26}\text{Al}$
as a heating element. Subsequently, the _grain model_ was used for a separate
set of simulations, in order to determine how iron distribution affects
heating, in order to determine if heating is the only major variable in
desiccation.

###  3.1 Core model

####  3.1.1 Core size ratio

First, we performed a set of simulations in order to explore how core size –
and therefore iron abundance – affects desiccation of planetesimals. The core-
to-radius ratio, $\Psi$, was varied between 0.05 and 0.95 in steps of 0.05,
with additional simulations of 0.01 and 0.99 for completeness. Simulations for
$\Psi=0$ and $\Psi=1$ were not conducted as these were determined to be
redundant (Section 2.3.1). The common parameters for these simulations were a
radius of $100\text{\,}\mathrm{k}\mathrm{m}$ with ${}^{60}\text{Fe}$
enrichment a factor of $10^{3}$ greater than Solar System estimates. No
${}^{26}\text{Al}$ enrichment was included in these simulations, in order to
focus solely on desiccation due to ${}^{60}\text{Fe}$ heating. The results of
this set of simulations were then used to constrain the core size parameter
space, in order to reduce the number of required simulations.

Fig. 5 shows that initially desiccation is fairly limited, but rapidly
increases as the core-to-planetesimal ratio increases. Eventually a point is
reached where increasing the core size does not lead to any appreciable
increases in desiccation amount. Desiccation increases rapidly when the
planetesimal is $\sim 1\%$ iron by volume, plateaus at $\sim 10\%$, and
decreases again at $\gtrsim 50\%$ iron by volume. As a result, two values of
$\Psi$ were chosen for further simulations:

  * •

$\Psi=0.25$, which is in the range of the initial increase in desiccation, and
corresponds to a Fe volume fraction of $1.56\%$.

  * •

$\Psi=0.50$, which is in the trough of maximum desiccation, and corresponds to
a Fe volume fraction of $12.5\%$.

Figure 5: _Core model_ simulation of final retained H2O fraction,
$f_{\text{H}_{2}\text{O},f}/f_{\text{H}_{2}\text{O},i}$, peak mean
planetesimal temperature and peak mean mantle temperature for simulations
varying core size, $\Psi$, with ${}^{60}\text{Fe}$ enrichment $10^{3}$ times
greater than Solar System estimates. Desiccation becomes significant after the
core radius, $r_{\text{c}}$, is greater than 20% of the radius of the
planetesimal, $r_{\text{p}}$, quickly reaches a maxima and decreases rapidly
beyond 90% due to only crustal deposits remaining.

####  3.1.2 Radius comparison

After determining how desiccation amount varies through modification of the
core size, we then progressed to determining how desiccation amount varies
through the size of the planetesimal itself. Lichtenberg et al. (2019a) notes
that there is a strong dependence on planetesimal radius and desiccation, with
larger planetesimals undergoing greater desiccation due to the significantly
greater mass of ${}^{26}\text{Al}$. Planetesimals were varied from
$1\text{\,}\mathrm{k}\mathrm{m}$ to $100\text{\,}\mathrm{k}\mathrm{m}$ with 1
dex spacing. The core-to-planetesimal radius ratio was maintained at
$\Psi=0.5$ and ${}^{60}\text{Fe}$ enrichment was varied between
$1\leq\Lambda_{\text{60Fe}}\leq 10^{4}$ with $1/3$ dex spacing between each
simulation. ${}^{26}\text{Al}$ enrichment was not included, in order to focus
on the influence of ${}^{60}\text{Fe}$ on water content.

Fig. 6 shows that desiccation is strongly dependent on planetesimal radius.
Small planetesimals $<$10\text{\,}\mathrm{k}\mathrm{m}$$ need an extremely
high degree of ${}^{60}\text{Fe}$ enrichment in order to undergo significant
desiccation, while above $10\text{\,}\mathrm{k}\mathrm{m}$ there is a rapid
increase in desiccation amount, which tapers out by
$100\text{\,}\mathrm{k}\mathrm{m}$ as water becomes increasingly rarefied.
Even in the case of large planetesimals high ${}^{60}\text{Fe}$ enrichment is
still required for desiccation; this is a recurring theme in later simulation
subsets. Lines for the $\Lambda_{\text{60Fe}}=1$ and
$\Lambda_{\text{60Fe}}=10$ simulation subsets were not included, as no
desiccation was observed in any of these simulations.

In order to constrain the parameter space of our simulations, we use a common
planetesimal sizes of $50\text{\,}\mathrm{k}\mathrm{m}$ and
$100\text{\,}\mathrm{k}\mathrm{m}$ for all subsequent simulations, which is
approximately at the peak of the birth size distribution of planetesimals
generated by the streaming instability mechanism (Li et al., 2019; Simon et
al., 2022).

Figure 6: Final retained H2O fraction, $\mathcal{H}_{f}$ for simulations with
varying planetesimal radius, $r_{\text{pl}}$. The _core model_ is utilised
with a core-to-planetesimal radius ratio of $\Psi=0.5$. There is a strong
inverse correlation between $r_{\text{pl}}$ and $\mathcal{H}_{f}$, similar to
what was observed in Lichtenberg et al. (2019a) with water fraction, radius
and ${}^{26}\text{Al}$ enrichment, though with significantly less pronounced
desiccation.

####  3.1.3 Isotopic enrichment

Once planetesimal geometry had been constrained, we conducted an additional
set of simulations in order to detail the influence of isotopic enrichment on
desiccation. Fig. 7 shows the results of simulations where only
${}^{60}\text{Fe}$ enrichment occurs, $\Lambda_{\text{60Fe}}$ is varied from 1
to $10^{4}$ and $\Psi$ is varied between 0.25 and 0.5. The results are
consistent with the previous simulations, and shows that desiccation becomes
significant above $200$ times the lower limit for Solar System
${}^{60}\text{Fe}$ enrichment (recall that the lower limit is $\Lambda_{\rm
60Fe}=1.15\times 10^{-8}$), with a larger core resulting in more desiccation.

Figure 7: Final retained H2O fraction for _core model_ simulations undergoing
radiogenic heating solely from ${}^{60}\text{Fe}$. ${}^{60}\text{Fe}$ must be
enriched to greater than approximately $3\times 10^{2}$ times the lower limit
of the Solar System abundance for desiccation to occur.

To compare the influence of ${}^{60}\text{Fe}$ and ${}^{26}\text{Al}$ on water
content, a number of simulations with both SLRs at varying degrees of
enrichment were conducted. Similar to the previous simulations, two subsets
were run, one where $\Psi=0.25$ and another where $\Psi=0.50$,
${}^{60}\text{Fe}$ abundance was varied from $1$ to $10^{3}$ times Solar
System enrichment in a sequence of the form $\\{10^{k}\\}$, planetesimal
radius was not varied, and set as $50\text{\,}\mathrm{k}\mathrm{m}$.
Meanwhile, ${}^{26}\text{Al}$ enrichment was restricted to 0, 1, 1.78, 5.62
and 10 times Solar System enrichment, a sequence of the form $\\{10^{k/4}\\}$.

Fig. 8 is a contour plot of the final retained water fraction,
$\mathcal{H}_{f}$, for the $\Psi=0.25$ simulation set. These plots use a
continuous gradient interpolated between these simulations using Delauney
triangulation, which can be used to generate contours between sparse, non-
linearly spaced data. We find that ${}^{26}\text{Al}$ is very effective at
reducing the final water content of a planetesimal, with near total
desiccation happening in planetesimals with only marginally greater than Solar
System ${}^{26}\text{Al}$ abundance. Conversely, ${}^{60}\text{Fe}$ does have
an impact, but only for simulation with little to no ${}^{26}\text{Al}$. In
the case of a simulation with total ${}^{26}\text{Al}$ depletion and
$\Lambda_{\text{60Fe}}=1000$ a final water fraction of $\mathcal{H}_{f}=0.796$
is found, which is higher than the Solar System enrichment estimation
simulation, which reports a value of $\mathcal{H}_{f}=0.509$. Fig. 9 details
the temperature change due to ${}^{26}\text{Al}$ and ${}^{60}\text{Fe}$
enrichment. An identical change is observed, with ${}^{60}\text{Fe}$
enrichment having a minimum impact on the peak mean planetesimal temperature,
while ${}^{26}\text{Al}$ enrichment rapidly pushes the planetesimal
temperature through the vaporisation point.

Fig. 10 is similar to Fig. 8 for the $\Psi=0.5$ simulation set, a similar lack
of influence due to ${}^{60}\text{Fe}$ enrichment is observed. However,
low-${}^{26}\text{Al}$-enrichment simulations do show markedly less
desiccation, primarily due to having $\sim 30\%$ of the hydrous silicate mass
of the $\Psi=0.25$ planetesimals. Conversely, as there is significantly more
${}^{60}\text{Fe}$, we observe a final water abundance fraction for the
$Z_{\text{26Al}}=0$, $\Lambda_{\text{60Fe}}=1000$ simulation of
$\mathcal{H}_{f}=0.153$, far below the maximum highest ${}^{60}\text{Fe}$-only
desiccation amount of the $\Psi=0.25$ simulations.

To summarise,${}^{26}\text{Al}$ is typically the primary radiogenic heating
source for planetesimals with a clearly differentiated core, requiring
${}^{60}\text{Fe}$ enrichment far greater than the lower Solar System
estimates to result in significant H2O outgassing. However, in planetary
systems that are highly enriched in ${}^{60}\text{Fe}$, for example by
supernovae enrichment, ${}^{60}\text{Fe}$ can start to significantly
contribute to planetesimal devolatilization.

Figure 8: Comparison of final H2O fraction, $\mathcal{H}_{f}$, for simulations
with a core ratio, $\Psi$, of 0.25 in the _core model_. Desiccation increases
dramatically as ${}^{26}\text{Al}$ enrichment increases; however
${}^{60}\text{Fe}$ enrichment only becomes significant with greater than
approximately $3\times 10^{2}\,Z_{\text{60Fe,ss,low}}$ or higher. Figure 9:
Comparison of peak planetesimal mean temperature, $\bar{T}_{\text{peak,pl}}$,
for simulations where $\Psi=0.25$ in the _core model_. Similarly to Fig. 8 we
see that there is a strong dependence on planetesimal temperature and
${}^{26}\text{Al}$ enrichment but an extremely weak dependence on
${}^{60}\text{Fe}$ enrichment and planetesimal temperature. This shows that in
the case of the _core model_ ${}^{60}\text{Fe}$ has a weak dependence. Figure
10: Comparison of $f_{\text{H}_{2}\text{O},f}/f_{\text{H}_{2}\text{O},i}$ for
_core model_ planetesimals where $\Psi=0.50$, as there is a significantly
lower mass of Al, there is a lower dependence on ${}^{26}\text{Al}$ enrichment
and a greater ${}^{60}\text{Fe}$ enrichment dependence compared to Fig. 8.
However, it is still not significant, as ${}^{60}\text{Fe}$ enrichment must be
2 orders of magnitude higher than Solar System estimates to have any
reasonable impact.

###  3.2 Grain model

For the second set of simulations we implemented the _grain model_ instead of
the _core model_. Two major parameters are varied over these simulations:
${}^{60}\text{Fe}$ enrichment $\Lambda_{\text{60Fe}}$ and iron volume fraction
$\Phi$. ${}^{26}\text{Al}$ enrichment $\Lambda_{\text{26Al}}$ is also varied,
though only between 0, 1 and 10. The parameter space of ${}^{26}\text{Al}$
enrichment was significantly compressed compared to the _core model_
simulation set, in order to accommodate the expanded ${}^{60}\text{Fe}$
enrichment space, and introduce the $\Phi$ parameter without a drastic
increase in simulation count. $\Lambda_{\text{60Fe}}$ was varied from 1 to
$10^{4}$ – with 1 dex steps – while $\Phi$ was varied between 0.01 and 0.9.
Planetesimal radius was increased from the _core model_ simulation set to
$100\text{\,}\mathrm{k}\mathrm{m}$. The results of these simulations are very
similar to the results of the _core model_ , ${}^{60}\text{Fe}$ enrichment has
to be significantly enriched compared to the Solar System in order for
desiccation to occur in the case of ${}^{26}\text{Al}$ depletion. Meanwhile,
${}^{26}\text{Al}$ enrichment relative to canonical Solar System levels
results in a greater deal of desiccation, and enrichment beyond that results
in near-total desiccation. A divergence from the previous simulation set
arises from varying $\Phi$, where high values impede desiccation from
${}^{26}\text{Al}$, as there is significantly less ${}^{26}\text{Al}$
available to heat the planetesimal. This is similar to the results observed
for the _core model_ $\Psi=0.5$ subset as shown in Fig. 10.

Fig. 11 shows $\mathcal{H}_{f}$ for varying values of enrichment parameters
$\Lambda_{\text{26Al}}$, $\Lambda_{\text{60Fe}}$ and grain volume fraction,
$\Phi$. Data is interpolated between simulations in the same manner as Section
3.1.3. Similarly to the _core model_ simulations, ${}^{60}\text{Fe}$
enrichment has a significantly lower impact on final water content than
${}^{26}\text{Al}$ enrichment. In fact, the link between desiccation and
${}^{60}\text{Fe}$ enrichment is less pronounced than with the _core model_.
In the case of fully ${}^{26}\text{Al}$-depleted simulations, we find that no
appreciable desiccation occurs until $\Lambda_{\text{60Fe}}>1000$ even for
simulations where $\Phi=0.9$. Fig. 12 shows the associated temperature of
these simulations, here we can see similar dependencies as with Fig. 11, and
that $\bar{T}_{\text{peak}}$ is not exceeded outside of simulations with
${}^{26}\text{Al}$ enrichment. In summary, the _grain model_ produces similar
results to the _core model_ , as such we can infer that bulk heating is more
important than heating specific areas of the planetesimal. Furthermore, this
lends further evidence that ${}^{26}\text{Al}$ enrichment is significantly
more influential with ${}^{60}\text{Fe}$/${}^{26}\text{Al}$ fractionation
factors close to the Solar System, as ${}^{60}\text{Fe}$ enrichment must be
multiple orders of magnitude higher to produce any significant effect, while a
similar effect occurs with a single-digit factor increase in
${}^{26}\text{Al}$ enrichment.

Figure 11: _Grain model_ comparison of $\mathcal{H}_{f}$ for simulations
varying ${}^{26}\text{Al}$ enrichment, ${}^{60}\text{Fe}$ enrichment and Fe
grain volume ratio, $\Phi$. Desiccation is similarly dependent on
${}^{26}\text{Al}$ enrichment as in the case of the _core model_ , while
having a similar weak dependence on ${}^{60}\text{Fe}$ enrichment. In order
for desiccation to begin the planetesimal requires a large amount of iron
grains by volume and an extremely large amount of enrichment ($>10^{3}\cdot
Z_{\text{Fe,ss}}$). This suggests that even with a more physically accurate
distribution of iron in a planetesimal it is still unlikely that
${}^{60}\text{Fe}$ is significant as a radiogenic heating source in planetary
formation. Figure 12: _Grain model_ comparison of peak mean planetesimal
temperature for simulations varying ${}^{26}\text{Al}$ enrichment,
${}^{60}\text{Fe}$ enrichment and Fe grain volume ratio, $\Phi$. We note that
cooling decreases in the case of ${}^{26}\text{Al}$ enriched simulations with
large quantities of low enriched Fe grains can actually lower planetesimal
temperature and hinder heating and hence desiccation, as there is a lower mass
of ${}^{26}\text{Al}$ to heat the planetesimal.

##  4 Discussion

###  4.1 ${}^{60}\text{Fe}$ versus ${}^{26}\text{Al}$

Our simulation results suggest that radiogenic heating from ${}^{26}\text{Al}$
is significantly more important compared to ${}^{60}\text{Fe}$ for
${}^{60}\text{Fe}$/${}^{26}\text{Al}$ fractions close to the Solar System
value. In certain cases, with ${}^{60}\text{Fe}$ enrichment more than $10^{2}$
to $10^{3}$ relative to the Solar System, ${}^{60}\text{Fe}$ heating can
result in significant desiccation of a planetesimal, though these values may
only be realised in massive star-forming regions, or with distinct supernovae
enrichment in individual systems. Future considerations for determining which
SLRs are important for the process of planetary formation would need to narrow
down a typical galactic enrichment range, and whether such ${}^{60}\text{Fe}$
enrichment levels are common on a planetary system level. Additional
simulations can narrow down this parameter space further by incorporating more
comprehensive models of stellar feedback on local scales (Nicholson et al.,
2019; Parker et al., 2023; Patel et al., 2023).

####  4.1.1 Other SLRs

Whilst other SLRs formed by stellar nucleosynthesis could be considered for
simulation, such as ${}^{41}\text{Ca}$ and ${}^{53}\text{Mn}$ (Russell et al.,
2001), these would have a lower abundance and enrichment to
${}^{60}\text{Fe}$, resulting in an even lower influence on desiccation. As
such, these were not included in our simulations.

###  4.2 ${}^{60}\text{Fe}$ as a temperature sustainer?

Figure 13: A comparison of mean power density, $P$, as a function of
simulation time for various simulations. The _grain model_ is utilised where
$\Phi=0.25$, while $\Lambda_{\text{60Fe}}$ and $\Lambda_{\text{26Al}}$ are
varied. In all cases specific power drops significantly between 2 and 10 Myr
after the start of the simulation as ${}^{26}\text{Al}$ decays, though
simulations with a higher content of ${}^{60}\text{Fe}$ still have a
significant power density.

Whilst our initial results show that there is a minimal impact of
${}^{60}\text{Fe}$ on planetesimal heating outside highly enriched systems,
the longer half-life of ${}^{60}\text{Fe}$ could potentially sustain internal
temperatures if the planetesimal was sufficiently enriched. This can be seen
in Fig. 2 where the specific heating of ${}^{60}\text{Fe}$ is greater than the
specific heating rate of ${}^{26}\text{Al}$
$\sim$2\text{\,}\mathrm{M}\mathrm{y}\mathrm{r}$$ from CAI formation. In the
case of Solar System abundances, however, this does not produce a significant
amount of heating, and is still almost 2 orders of magnitude less than the
specific heating rate of ${}^{26}\text{Al}$ at the specific heating crossing
point. For more ${}^{60}\text{Fe}$-enriched planetesimals this could result in
sustained heating. A continued heating source would prevent the planetesimal
from cooling and continue the vaporisation and out-gassing processes. In order
to infer the influence of ${}^{60}\text{Fe}$ over longer time scales from our
simulations, we calculated the mean averaged radiogenic heating power density:

$\overline{P}=\frac{\sum{Q_{\text{T,cell}}(t)/\rho_{\text{cell}}}}{N_{\text{cells}}},$
(8)

over the planetesimal. Fig. 13 shows a comparison of $\overline{P}$ over time
between simulations utilising the _grain model_ where $\Phi=0.25$ with varying
isotopic enrichment. After $2\,$Myr the average power density falls in the
case of all simulations, becoming similar for simulations with matching
${}^{60}\text{Fe}$ enrichment amounts. It is clear that in all cases with
significant ${}^{60}\text{Fe}$ enrichment that ${}^{60}\text{Fe}$ heating
becomes the strongest heating mechanism near the end of the simulations, and
should retain adequate heating for some time after the simulations have
concluded. Whilst this enduring heating would not change the results with our
desiccation model, more complex thermochemical models with out-gassing and
without a discrete evaporation temperature would lead to greater desiccation
rates.

The sustenance of high temperatures is relevant for late-formed planetesimals
and cometesimals in extrasolar systems. In the Solar System, there is evidence
for prolonged planetesimal formation, essentially until the very end of the
disk phase (Kleine et al., 2020; Simon et al., 2022; Lichtenberg et al.,
2023). Recent evidence from JWST for water-enriched inner disks in both low-
mass and high-mass star-forming regions (Perotti et al., 2023; Ramirez-Tannus
et al., 2023) illustrates that desiccation of planetesimals until very late
stages of the disk could impact the final composition of rocky and terrestrial
exoplanets. How effective the sustenance by ${}^{60}\text{Fe}$ is will depend
on the late-stage pebble flux and dislocation of vapour from ice-rock
mixtures. Our simulations show that a very high ${}^{60}\text{Fe}$ enrichment
relative to the solar system produces temperatures amenable for
devolatilization, which would be sustained until the very late stages of disk
evolution. During the late stages of the disk the varied composition of comet-
like objects can impact the thermal evolution of still forming bodies within
the disk (see Section 2.2) (Golabek & Jutzi, 2021; Arakawa & Wakita, 2023).
This effect is not taken into account in the simulations discussed in this
paper. As such, further investigation of the conditions where
${}^{60}\text{Fe}$ can still maintain devolatilization in planetary systems
should be performed. Such investigations should ideally have a focus on
smaller objects, high ice fractions, and varying composition (such as a varied
C/O ratio) (Davidsson, 2021; Lichtenberg & Krijt, 2021).

###  4.3 Possible influence on exoplanet populations

Gaining information on the SLR distribution across planet-forming systems has
been a continuous challenge, which complicates assessing their influence on
exoplanet formation (Lugaro et al., 2018; Parker, 2020). However, it has been
suggested that planetary debris in polluted white dwarf systems may serve as
additional constraint on the SLR distribution across planetary systems if the
effects of accretion energy and SLR heating can be distinguished (Jura et al.,
2013; Jura & Young, 2014). With increasing measurement precision and numerical
modelling, research has shown that polluted white dwarfs indeed seem to
conserve the footprint of planetesimal differentiation across a significant
fraction of planetary systems (Bonsor et al., 2020, 2023b; Curry et al.,
2022).

From a star formation perspective, for a typical Initial Mass Function (e.g.
Maschberger, 2013), massive stars that would produce 60Fe are produced in
significant numbers (i.e. $>$5) in regions that also form more than 1000 low-
mass stars. Notable examples of such regions in our Galaxy are Cyg OB2 (e.g.
Wright et al., 2015), Westerlund 1 (e.g. Clark et al., 2005) and the Arches
cluster (Figer et al., 2002). However, stars likely form stochastically, and
in this scenario the only minimum limit on the most massive stars that can
form in a region is the total mass of the region itself (Elmegreen, 2006),
though see Weidner & Kroupa (2006). However, in any region that produces
massive stars, significant photoevaporation of the gaseous component of
protoplanetary discs occurs (e.g. Scally & Clarke, 2001; Adams et al., 2004;
Nicholson et al., 2019; Concha-Ramírez et al., 2019; Winter & Haworth, 2022),
meaning that the planetary systems that could be heavily enriched in 60Fe
would likely be devoid of gas giants (Patel et al, submitted).

An alternative enrichment scenario, which avoids the destructive
photoevaporation from massive stars, is from Asymptotic Giant Branch stars
(e.g. Karakas & Lugaro, 2016; Lugaro et al., 2018; Parker & Schoettler, 2023).
AGB stars are a stellar evolutionary phase that all 1 – 8 M⊙ stars undergo
(Herwig, 2005), and these stars produce significant yields of
${}^{60}\text{Fe}$ and ${}^{26}\text{Al}$ via their winds (Ventura et al.,
2018). AGB stars were previously discounted as a viable source of SLR
enrichment due to the supposed low probability of an evolved star encountering
a young star as it forms its planetary system (Kastner & Myers, 1994).
However, Parker & Schoettler (2023) report the serendipitous Gaia DR3
discovery of an AGB star interloping through the young star-forming region NGC
2264, and furthermore, show that under reasonable assumptions for the initial
conditions of the star-forming region, an AGB star could enrich the Solar
System in both ${}^{60}\text{Fe}$ and ${}^{26}\text{Al}$, with some stars
attaining higher ${}^{60}\text{Fe}$ ratios than in the Solar system (e.g.
$\Lambda_{\text{60Fe}}\sim 10^{-4}$).

In addition, studies of planet-forming disks indicate that planetesimal
formation is rapid, typically operating on a timescale of $10^{5}$ yr
(Drazkowska et al., 2022; Manara et al., 2022). Combined, this suggests that a
significant fraction of exoplanetary systems are expected to undergo
significant SLR-driven heating during the formation of their planets, with
substantial fractionation possible between ${}^{60}\text{Fe}$ and
${}^{26}\text{Al}$.

Exoplanet surveys have just started to reach into the low-mass regime, where
it is possible to distinguish bulk-level under- and over-densities in the size
regime of super-Earths and below (Wordsworth & Kreidberg, 2022; Piette et al.,
2023). Further investigation of these population trends will constrain the
effects of accretion environment and geophysical evolution on bulk
composition. While typically astronomical studies tend to focus on atmospheric
escape and disk migration as the main effects on volatile bulk composition
(Venturini et al., 2020a; Bean et al., 2021), we have shown here that
variations in the main SLRs can drive substantial devolatilization of the
planetesimal building blocks. From a Solar System perspective, meteoritic
evidence of initial incorporation of water, and other volatiles and rapid
devolatilization of the building blocks of Earth and its planetesimal
precursors is accumulating (Sarafian et al., 2017; McCubbin & Barnes, 2019;
Lichtenberg et al., 2021; Hirschmann et al., 2021; Grewal et al., 2021;
Peterson et al., 2023; Lewis et al., 2022; Newcombe et al., 2023; Stephant et
al., 2023; Grant et al., 2023). This presents a challenge for a wide spread of
potentially Earth-like rocky planets with a similar geodynamic an climatic
regime, as the mass fraction of surface water to allow for the coexistence of
oceans and exposed land with a geodynamic regime as the Earth is limited to
within a relatively narrow range (Cowan & Abbot, 2014; Schaefer & Sasselov,
2015; Noack et al., 2017). Observational insights on the water content and
devolatilization trends of planetary debris around polluted white dwarfs
(Farihi et al., 2013; Curry et al., 2022) and in debris disks (Lichtenberg &
Krijt, 2021; Marino, 2022; Bonsor et al., 2023a) will be crucial to narrow
down the potential spread in these effects.

Ultimately, increased understanding of the influence of planetesimal degassing
is required to understand the information received from exoplanet population
analyses, planetary debris, and detailed characterization of individual low-
mass exoplanets.

##  5 Conclusions

Overall we conclude that the influence of radiogenic heating on water content
in protoplanetary systems due to ${}^{60}\text{Fe}$ is weaker compared to
heating from the typically more abundant and energetic ${}^{26}\text{Al}$.
However, the contribution of ${}^{60}\text{Fe}$ to the devolatilization of
early-formed planetesimals becomes substantial for enrichment levels of
${}^{60}\mathrm{Fe}/^{56}\mathrm{Fe}\gtrsim 10^{-6}$, which is the upper end
range of enrichment values that has previously been inferred for the Solar
System. In exoplanet systems that formed in high-mass star-forming regions
with substantial supernovae feedback, ${}^{60}\text{Fe}$ may thus be an
important contributor to SLR-driven heating. Our parameter space explored
${}^{26}\text{Al}$ enrichment over an order of magnitude, while
${}^{60}\text{Fe}$ was explored over four orders of magnitude, we found that
${}^{60}\text{Fe}$ did not become dominant unless the system was
${}^{26}\text{Al}$ depleted, ${}^{60}\text{Fe}$ enriched, or possessed a high
iron content. Each of these characteristics may be fullfilled in exoplanetary
systems forming in diverse star-forming regions, however, we anticipate
${}^{26}\text{Al}$ to be the main driving factor of planetesimal
differentiation and devolatilization.

Fractionation in ${}^{60}\text{Fe}$ and ${}^{26}\text{Al}$ as a function of
the birth star-forming region is expected due to intrinsic stochasticity of
star formation (where variations in the IMF lead to significantly different
numbers of massive stars, and therefore enrichment, e.g. Nicholson & Parker,
2017), the stochasticity of dynamical evolution of young star-forming regions,
which affects the amount of SLR ejecta that may be captured by an individual
planetary system (Patel et al., submitted) and the short half-lives of the
main SLRs. Our results thus suggest that SLR-driven internal desiccation of
planetesimals may contribute to compositional scatter across exoplanetary
systems, following the initial SLR enrichment of exoplanetary systems at their
birth. Future work shall explore the geochemical effect of water flow and iron
oxidation inside volatile-rich planetesimals, and the compositional effects of
ongoing accretion on the final composition and atmospheric diversity of low-
mass exoplanets.

## Acknowledgements

JWE and RJP acknowledge support from the Royal Society from a Dorothy Hodgkin
Fellowship and an Enhancement Award. TL was suppported by a grant from the
Branco Weiss Foundation.

## Data availability

The data underlying this article is stored at The University of Sheffield ORDA
repository. The data and scripts to reproduce the figures included in this
paper are stored as well as the initial conditions for each simulation, such
that the simulations can be performed with a suitable workstation. The
repository is located at https://doi.org/10.15131/shef.data.23815599.v1.

## References

  * Adams (2010) Adams F. C., 2010, ARA&A, 48, 47
  * Adams (2021) Adams F. C., 2021, ApJ, 919, 10
  * Adams et al. (2004) Adams F. C., Hollenbach D., Laughlin G., Gorti U., 2004, ApJ, 611, 360
  * Arakawa & Wakita (2023) Arakawa S., Wakita S., 2023, arXiv e-prints, p. arXiv:2401.00231
  * Bagdassarov et al. (2009) Bagdassarov N., Solferino G., Golabek G. J., Schmidt M. W., 2009, Earth and Planetary Science Letters, 288, 84
  * Bean et al. (2021) Bean J. L., Raymond S. N., Owen J. E., 2021, Journal of Geophysical Research (Planets), 126, e06639
  * Bland & Travis (2017) Bland P. A., Travis B. J., 2017, Science Advances, 3, e1602514
  * Bonsor et al. (2020) Bonsor A., Carter P. J., Hollands M., Gänsicke B. T., Leinhardt Z., Harrison J. H. D., 2020, MNRAS, 492, 2683
  * Bonsor et al. (2023a) Bonsor A., Wyatt M. C., Marino S., Davidsson B. J. R., Kral Q., 2023a, arXiv e-prints, p. arXiv:2307.01574
  * Bonsor et al. (2023b) Bonsor A., Lichtenberg T., Drążkowska J., Buchan A. M., 2023b, Nature Astronomy, 7, 39
  * Boss (2019) Boss A. P., 2019, ApJ, 870, 3
  * Boss (2022) Boss A. P., 2022, ApJ, 933, 1
  * Castillo-Rogez & Young (2017) Castillo-Rogez J., Young E. D., 2017, in Elkins-Tanton L. T., Weiss B. P., eds, , Planetesimals: Early Differentiation and Consequences for Planets. Cambridge University Press, pp 92–114, doi:10.1017/9781316339794.005
  * Castillo-Rogez et al. (2009) Castillo-Rogez J., Johnson T. V., Lee M. H., Turner N. J., Matson D. L., Lunine J., 2009, Icarus, 204, 658
  * Cerantola et al. (2015) Cerantola V., Walte N. P., Rubie D. C., 2015, Earth and Planetary Science Letters, 417, 67
  * Chevalier (2000) Chevalier R. A., 2000, The Astrophysical Journal, 538, L151
  * Clark et al. (2005) Clark J. S., Negueruela I., Crowther P. A., Goodwin S. P., 2005, A&A, 434, 949
  * Concha-Ramírez et al. (2019) Concha-Ramírez F., Wilhelm M. J. C., Portegies Zwart S., Haworth T. J., 2019, MNRAS, 490, 5678
  * Cook et al. (2021) Cook D. L., Meyer B. S., Schönbächler M., 2021, ApJ, 917, 59
  * Cowan & Abbot (2014) Cowan N. B., Abbot D. S., 2014, ApJ, 781, 27
  * Crameri et al. (2012) Crameri F., et al., 2012, Geophysical Journal International, 189, 38
  * Curry et al. (2022) Curry A., Bonsor A., Lichtenberg T., Shorttle O., 2022, MNRAS, 515, 395
  * Davidsson (2021) Davidsson B. J. R., 2021, MNRAS, 505, 5654
  * Desch et al. (2022) Desch S. J., Young E. D., Dunham E. T., Fujimoto Y., Dunlap D. R., 2022, arXiv e-prints, p. arXiv:2203.11169
  * Desch et al. (2023) Desch S. J., Dunlap D. R., Williams C. D., Mane P., Dunham E. T., 2023, Icarus, 402, 115611
  * Diamond-Lowe et al. (2022) Diamond-Lowe H., et al., 2022, AJ, 164, 172
  * Diehl et al. (2021) Diehl R., et al., 2021, Publ. Astron. Soc. Australia, 38, e062
  * Dodds et al. (2021) Dodds K. H., Bryson J. F. J., Neufeld J. A., Harrison R. J., 2021, Journal of Geophysical Research: Planets, 126, e2020JE006704
  * Dorn & Lichtenberg (2021) Dorn C., Lichtenberg T., 2021, ApJ, 922, L4
  * Drazkowska et al. (2022) Drazkowska J., et al., 2022, arXiv e-prints, p. arXiv:2203.09759
  * Dunham et al. (2022) Dunham E. T., et al., 2022, Geochimica Cosmochimica Acta, 324, 194
  * Elmegreen (2006) Elmegreen B. G., 2006, ApJ, 648, 572
  * Farihi et al. (2013) Farihi J., Gänsicke B. T., Koester D., 2013, Science, 342, 218
  * Fatuzzo & Adams (2022) Fatuzzo M., Adams F. C., 2022, ApJ, 925, 56
  * Figer et al. (2002) Figer D. F., et al., 2002, ApJ, 581, 258
  * Fu et al. (2017) Fu R. R., Young E. D., Greenwood R. C., Elkins-Tanton L. T., 2017, in Elkins-Tanton L. T., Weiss B. P., eds, , Planetesimals: Early Differentiation and Consequences for Planets. Cambridge University Press, pp 115–135, doi:10.1017/9781316339794.006
  * Gaidos et al. (2009) Gaidos E., Krot A. N., Williams J. P., Raymond S. N., 2009, ApJ, 696, 1854
  * Gerya (2019) Gerya T., 2019, Introduction to numerical geodynamic modelling. Cambridge University Press
  * Gerya & Yuen (2003) Gerya T. V., Yuen D. A., 2003, Physics of the Earth and Planetary Interiors, 140, 293
  * Gerya & Yuen (2007) Gerya T. V., Yuen D. A., 2007, Physics of the Earth and Planetary Interiors, 163, 83
  * Ghanbarzadeh et al. (2017) Ghanbarzadeh S., Hesse M. A., Prodanović M., 2017, Proceedings of the National Academy of Science, 114, 13406
  * Golabek & Jutzi (2021) Golabek G. J., Jutzi M., 2021, Icarus, 363, 114437
  * Grant et al. (2023) Grant H., Tartèse R., Jones R., Piani L., Marrocchi Y., King A., Rigaudier T., 2023, Meteoritics & Planetary Science
  * Gray & Compston (1974) Gray C. M., Compston W., 1974, Nature, 251, 495
  * Grewal et al. (2021) Grewal D. S., Dasgupta R., Marty B., 2021, Nature Astronomy, 5, 356
  * Grimm & McSween (1989) Grimm R. E., McSween Harry Y. J., 1989, Icarus, 82, 244
  * Grimm & McSween (1993) Grimm R. E., McSween H. Y., 1993, Science, 259, 653
  * Herwig (2005) Herwig F., 2005, ARA&A, 43, 435
  * Hirschmann et al. (2021) Hirschmann M. M., Bergin E. A., Blake G. A., Ciesla F. J., Li J., 2021, Proceedings of the National Academy of Science, 118, e2026779118
  * Johansen et al. (2023) Johansen A., Ronnet T., Schiller M., Deng Z., Bizzarro M., 2023, A&A, 671, A74
  * Jura & Young (2014) Jura M., Young E. D., 2014, Annual Review of Earth and Planetary Sciences, 42, 45
  * Jura et al. (2013) Jura M., Xu S., Young E. D., 2013, ApJ, 775, L41
  * Karakas & Lugaro (2016) Karakas A. I., Lugaro M., 2016, ApJ, 825, 26
  * Kastner & Myers (1994) Kastner J. H., Myers P. C., 1994, ApJ, 421, 605
  * Keller & Suckale (2019) Keller T., Suckale J., 2019, Geophysical Journal International, 219, 185
  * Kita et al. (2013) Kita N. T., et al., 2013, Meteoritics & Planetary Science, 48, 1383
  * Kleine et al. (2020) Kleine T., Budde G., Burkhardt C., Kruijer T. S., Worsham E. A., Morbidelli A., Nimmo F., 2020, Space Sci. Rev., 216, 55
  * Krijt et al. (2022) Krijt S., Kama M., McClure M., Teske J., Bergin E. A., Shorttle O., Walsh K. J., Raymond S. N., 2022, arXiv e-prints, p. arXiv:2203.10056
  * Kuffmeier et al. (2016) Kuffmeier M., Frostholm Mogensen T., Haugbølle T., Bizzarro M., Nordlund Å., 2016, ApJ, 826, 22
  * Lacedelli et al. (2022) Lacedelli G., et al., 2022, MNRAS, 511, 4551
  * Laird et al. (2023) Laird A. M., et al., 2023, Journal of Physics G Nuclear Physics, 50, 033002
  * Lee et al. (1976) Lee T., Papanastassiou D. A., Wasserburg G. J., 1976, Geophys. Res. Lett., 3, 41
  * Lewis et al. (2022) Lewis J. A., Jones R. H., Brearley A. J., 2022, Geochimica Cosmochimica Acta, 316, 201
  * Li et al. (2019) Li R., Youdin A. N., Simon J. B., 2019, ApJ, 885, 69
  * Lichtenberg & Clement (2022) Lichtenberg T., Clement M. S., 2022, ApJ, 938, L3
  * Lichtenberg & Krijt (2021) Lichtenberg T., Krijt S., 2021, ApJ, 913, L20
  * Lichtenberg et al. (2016) Lichtenberg T., Parker R. J., Meyer M. R., 2016, Monthly Notices of the Royal Astronomical Society, 462, 3979
  * Lichtenberg et al. (2019a) Lichtenberg T., Golabek G. J., Burn R., Meyer M. R., Alibert Y., Gerya T. V., Mordasini C., 2019a, Nature Astronomy, 3, 307
  * Lichtenberg et al. (2019b) Lichtenberg T., Keller T., Katz R. F., Golabek G. J., Gerya T. V., 2019b, Earth and Planetary Science Letters, 507, 154
  * Lichtenberg et al. (2021) Lichtenberg T., Drążkowska J., Schönbächler M., Golabek G. J., Hands T. O., 2021, Science, 371, 365
  * Lichtenberg et al. (2023) Lichtenberg T., Schaefer L. K., Nakajima M., Fischer R. A., 2023, in Inutsuka S., Aikawa Y., Muto T., Tomida K., Tamura M., eds, Astronomical Society of the Pacific Conference Series Vol. 534, Protostars and Planets VII. p. 907 (arXiv:2203.10023), doi:10.48550/arXiv.2203.10023
  * Limongi & Chieffi (2006) Limongi M., Chieffi A., 2006, The Astrophysical Journal, 647, 483
  * Limongi & Chieffi (2018) Limongi M., Chieffi A., 2018, The Astrophysical Journal Supplement Series, 237, 13
  * Lugaro et al. (2018) Lugaro M., Ott U., Kereszturi Á., 2018, Progress in Particle and Nuclear Physics, 102, 1
  * Luque & Pallé (2022) Luque R., Pallé E., 2022, Science, 377, 1211
  * Manara et al. (2022) Manara C. F., Ansdell M., Rosotti G. P., Hughes A. M., Armitage P. J., Lodato G., Williams J. P., 2022, arXiv e-prints, p. arXiv:2203.09930
  * Marino (2022) Marino S., 2022, arXiv e-prints, p. arXiv:2202.03053
  * Maschberger (2013) Maschberger T., 2013, MNRAS, 429, 1725
  * McCubbin & Barnes (2019) McCubbin F. M., Barnes J. J., 2019, Earth and Planetary Science Letters, 526, 115771
  * McKeegan et al. (2000) McKeegan K. D., Chaussidon M., Robert F., 2000, Science, 289, 1334
  * Monteux et al. (2018) Monteux J., Golabek G. J., Rubie D. C., Tobie G., Young E. D., 2018, Space Sci. Rev., 214, 39
  * Neumann et al. (2012) Neumann W., Breuer D., Spohn T., 2012, Astronomy & Astrophysics, 543, A141
  * Neumann et al. (2018) Neumann W., Kruijer T. S., Breuer D., Kleine T., 2018, Journal of Geophysical Research: Planets, 123, 421
  * Newcombe et al. (2023) Newcombe M. E., et al., 2023, Nature, 615, 854
  * Nicholson & Parker (2017) Nicholson R. B., Parker R. J., 2017, MNRAS, 464, 4318
  * Nicholson et al. (2019) Nicholson R. B., Parker R. J., Church R. P., Davies M. B., Fearon N. M., Walton S. R. J., 2019, MNRAS, 485, 4893
  * Nimmo & Kleine (2015) Nimmo F., Kleine T., 2015, Washington DC American Geophysical Union Geophysical Monograph Series, 212, 83
  * Noack et al. (2017) Noack L., Snellen I., Rauer H., 2017, Space Sci. Rev., 212, 877
  * Ouellette et al. (2010) Ouellette N., Desch S. J., Hester J. J., 2010, ApJ, 711, 597
  * Parker (2020) Parker R. J., 2020, Royal Society Open Science, 7, 201271
  * Parker & Schoettler (2023) Parker R. J., Schoettler C., 2023, The Astrophysical Journal Letters, 952, L16
  * Parker et al. (2014) Parker R. J., Church R. P., Davies M. B., Meyer M. R., 2014, MNRAS, 437, 946
  * Parker et al. (2023) Parker R. J., Lichtenberg T., Patel M., Polius C. K. M., Ridsdill-Smith M., 2023, MNRAS, 521, 4838
  * Patel et al. (2023) Patel M., Polius C. K. M., Ridsdill-Smith M., Lichtenberg T., Parker R. J., 2023, MNRAS, 525, 2399
  * Perotti et al. (2023) Perotti G., et al., 2023, Nature, 620, 516
  * Peterson et al. (2023) Peterson L. D., Newcombe M. E., Alexander C. M. O., Wang J., Sarafian A. R., Bischoff A., Nielsen S. G., 2023, Geochimica Cosmochimica Acta, 340, 141
  * Pfalzner et al. (2015) Pfalzner S., et al., 2015, Phys. Scr., 90, 068001
  * Piaulet et al. (2023) Piaulet C., et al., 2023, Nature Astronomy, 7, 206
  * Piette et al. (2023) Piette A. A. A., Gao P., Brugman K., Shahar A., Lichtenberg T., Miozzi F., Driscoll P., 2023, arXiv e-prints, p. arXiv:2306.10100
  * Portegies Zwart (2019) Portegies Zwart S., 2019, A&A, 622, A69
  * Ramirez-Tannus et al. (2023) Ramirez-Tannus M. C., et al., 2023, arXiv e-prints, p. arXiv:2310.11074
  * Russell et al. (2001) Russell S. S., Gounelle M., Hutchison R., 2001, Philosophical Transactions: Mathematical, Physical and Engineering Sciences, 359, 1991
  * Sarafian et al. (2017) Sarafian A. R., et al., 2017, Philosophical Transactions of the Royal Society of London Series A, 375, 20160209
  * Scally & Clarke (2001) Scally A., Clarke C., 2001, MNRAS, 325, 449
  * Schaefer & Sasselov (2015) Schaefer L., Sasselov D., 2015, ApJ, 801, 40
  * Simon et al. (2022) Simon J. B., Blum J., Birnstiel T., Nesvorný D., 2022, arXiv e-prints, p. arXiv:2212.04509
  * Spaargaren et al. (2023) Spaargaren R. J., Wang H. S., Mojzsis S. J., Ballmer M. D., Tackley P. J., 2023, ApJ, 948, 53
  * Stephant et al. (2023) Stephant A., Zhao X., Anand M., Davidson J., Carli C., Cuppone T., Pratesi G., Franchi I. A., 2023, Earth and Planetary Science Letters, 615, 118202
  * Sturtz et al. (2022) Sturtz C., Limare A., Chaussidon M., Kaminski É., 2022, Icarus, 385, 115100
  * Suer et al. (2023) Suer T.-A., Jackson C., Grewal D. S., Dalou C., Lichtenberg T., 2023, Front. Earth. Sci., 11, 1159412
  * Tang & Dauphas (2012) Tang H., Dauphas N., 2012, Earth and Planetary Science Letters, 359, 248
  * Telus et al. (2018) Telus M., Huss G. R., Nagashima K., Ogliore R. C., Tachibana S., 2018, Geochimica Cosmochimica Acta, 221, 342
  * Trappitsch et al. (2018) Trappitsch R., et al., 2018, ApJ, 857, L15
  * Travis & Schubert (2005) Travis B. J., Schubert G., 2005, Earth and Planetary Science Letters, 240, 234
  * Urey (1955) Urey H. C., 1955, PNAS, 41, 127
  * Vasileiadis et al. (2013) Vasileiadis A., Nordlund Å., Bizzarro M., 2013, ApJ, 769, L8
  * Ventura et al. (2018) Ventura P., Karakas A., Dell’Agli F., García-Hernández D. A., Guzman-Ramirez L., 2018, MNRAS, 475, 2282
  * Venturini et al. (2020a) Venturini J., Ronco M. P., Guilera O. M., 2020a, Space Sci. Rev., 216, 86
  * Venturini et al. (2020b) Venturini J., Guilera O. M., Haldemann J., Ronco M. P., Mordasini C., 2020b, A&A, 643, L1
  * Villeneuve et al. (2009) Villeneuve J., Chaussidon M., Libourel G., 2009, Science, 325, 985
  * Walte et al. (2020) Walte N. P., Solferino G. F. D., Golabek G. J., Silva Souza D., Bouvier A., 2020, Earth and Planetary Science Letters, 546, 116419
  * Walte et al. (2023) Walte N. P., Howard C. M., Golabek G. J., 2023, Earth and Planetary Science Letters, 617, 118247
  * Wang et al. (2019) Wang H. S., Liu F., Ireland T. R., Brasser R., Yong D., Lineweaver C. H., 2019, MNRAS, 482, 2222
  * Weidner & Kroupa (2006) Weidner C., Kroupa P., 2006, MNRAS, 365, 1333
  * Winter & Haworth (2022) Winter A. J., Haworth T. J., 2022, European Physical Journal Plus, 137, 1132
  * Wordsworth & Kreidberg (2022) Wordsworth R., Kreidberg L., 2022, ARA&A, 60, 159
  * Wright et al. (2015) Wright N. J., Drew J. E., Mohr-Smith M., 2015, MNRAS, 449, 741
  * Yoshino et al. (2003) Yoshino T., Walter M. J., Katsura T., 2003, Nature, 422, 154
  * Zeng et al. (2019) Zeng L., et al., 2019, Proceedings of the National Academy of Science, 116, 9723
  * Zhang (2023) Zhang Z., 2023, ApJ, 956, L25
  * Zhang et al. (2021) Zhang Z., Bercovici D., Jordan J. S., 2021, Journal of Geophysical Research (Planets), 126, e06754
  * Zilinskas et al. (2023) Zilinskas M., Miguel Y., van Buchem C. P. A., Snellen I. A. G., 2023, arXiv e-prints, p. arXiv:2301.05190