# Serpentinization in the thermal evolution of icy Kuiper belt objects in the
early Solar system

Anikó Farkas-Takács Konkoly Observatory, Research Centre for Astronomy and
Earth Sciences, Konkoly Thege 15-17, H-1121 Budapest, Hungary  Eötvös Loránd
University, Faculty of Science, Pázmány P. st. 1/A, 1171 Budapest, Hungary
Csaba Kiss Konkoly Observatory, Research Centre for Astronomy and Earth
Sciences, Konkoly Thege 15-17, H-1121 Budapest, Hungary  ELTE Eötvös Loránd
University, Institute of Physics, Pázmány P. st. 1/A, 1171 Budapest, Hungary
Sándor Góbi Eötvös Loránd University, Faculty of Science, Pázmány P. st. 1/A,
1171 Budapest, Hungary  MTA-ELTE Lendület Laboratory Astrochemistry Research
Group, Institute of Chemistry, ELTE Eötvös Loránd University, H-1518 Budapest,
Hungary  Ákos Kereszturi Konkoly Observatory, Research Centre for Astronomy
and Earth Sciences, Konkoly Thege 15-17, H-1121 Budapest, Hungary

###### Abstract

Here we present an improved algorithm to model the serpentinization process in
planetesimals in the early Solar system. Although it is hypothesized that
serpentinization-like reactions played an important role in the thermal
evolution of planetesimals, few and restricted models are available in this
topic. These process may be important as the materials involved were abundant
in these objects. Our model is based on the model by Góbi & Kereszturi (2017),
and contains improvements in the consideration of heat capacities and
lithospheric pressure, and in the calculation of the amount of interfacial
water. Comparison of our results with previous calculations show that there
are significant differences in the e.g. the serpentinization time – the time
necessary to consume most of the reactants at specific initial conditions – or
the amount of heat produced by this process. In a simple application we show
that in icy bodies, under some realistic conditions, below the melting point
of water ice, serpentinization reaction using interfacial water may be able to
proceed and eventually push the local temperature above the melting point to
start a ’runaway’ serpentinization. According to our calculations in objects
with radii R $\gtrsim$ 200 km serpentinization might have quickly reformed
nearly the whole interior of these bodies in the early Solar system.

planetary systems: planets and satellites: composition, interiors

##  1 Introduction

We only have indirect information on the internal material properties of the
small bodies in the Solar system, coming chiefly from the analysis of
meteorites. Some samples show that, in addition to heat from accretion and
radioactive decay, there may have been chemical processes that at some point
in the early history of the Solar system significantly altered the
mineralogical and lithological characteristics of the objects and contributed
to internal heat production. One of these chemical processes is the hydration
of silicates. As it is an exothermic reaction, it may have contributed to the
heating of planetesimals in the early Solar system (Gail et al., 2014),
especially in bodies with water ice content. The aqueous alteration found in
carbonaceous chondrites is partly explained by the formation of Mg–serpentine
via the serpentinization and we use here serpentinization as a model approach,
as it is a moderately simple reaction and uses abundant reactants.
Serpentinite consists of one or more serpentine group minerals resulting from
the hydration of silicates. The presence of serpentine has been detected in
meteorites, especially in CM chondrites with high carbon content whose primary
minerals are members of the serpentine group and account for 55–58% of the
meteorite by volume (Scott et al., 1988).

According to some views, serpentinization is a very rare process, and
serpentine found inside meteorites were formed already in the planetary nebula
or during accretion (Lunine, 2006). There may have also been other changes
inside the planetesimals that could alter serpentine. For instance, members of
the serpentine family may have been dehydrated and/or altered by heat.
Serpentinization itself can produce significant heat, which can cause
serpentine to become amorphous, as it has already been found in some
meteorites (Zega et al., 2003). This is probably the reason why no well-
crystallized serpentine minerals have been observed.

There are several reactions to form serpentinite from olivine. In these
reactions the rock absorbs a large amount of water and consequently destroys
the structure of the original minerals while it increases its volume and
decreases its density. In one of the main serpentinization processes Mg-
pyroxenes (enstatite, MgSiO3) is also required in addition to the Mg - rich
olivine (forsterite, Mg2SiO4); the end product of the reaction is purely
serpentinite (Mg3Si2O5(OH)4). The following reaction shows the stoichiometric
equation of the formation of serpentinite:

$\ce{Mg2SiO4+MgSiO3+2H2O->Mg3Si2O5(OH)4}$ (1)

This reaction takes place only in the presence of liquid water at a
temperature- and pressure-dependent reaction rate (Wegner & Ernst, 1983). In
this case, the enthalpy is $69\,kJ\,mol^{-1}$ (Robie & Waldbaum, 1968), with a
weak dependence on the temperature (Fyfe, 1974).

In another reaction brucite (Mg(OH)2) is produced in addition to serpentine
(Martin & Fyfe, 1970). This has not been found in large amounts in meteorites,
however, it might have been thermally decomposed after its formation.

$\ce{2Mg2SiO4+3H2O->Mg3Si2O5(OH)4+Mg(OH)2}$ (2)

The actual rate of serpentinization depends on the composition of the rock and
the ability of the liquid to transport magnesium and other elements during the
process.

Based on meteorite samples, it is generally estimated that these aqueous
changes occur on a short time scale of about 100 years, at low temperatures,
but above the melting temperature of H2O, supported by previous models
(Dufresne & Anders, 1962; Grimm & Mcsween, 1989; Zolensky et al., 1989). Góbi
& Kereszturi (2017) pointed out the importance of interfacial water. This type
of water can exist at temperatures below the melting point of bulk ”classical”
water as a microscopic liquid film along mineral–H2O interfaces at subzero
temperatures, mainly due to van der Waals forces, enabling the
serpentinization process to proceed. Based on the simulations they estimated
that the reaction duration was comparable to time scales from previous
studies.

When constructing a serpentinization heat model we have to take into account
the amount of heat released, the heat consumption required to melt the ice,
and the heat transfer as well. There are several models for estimating heat
production. The simplest is the heat-balance model (Lowell & Rona, 2002) which
assumes the rock-water system is stationary and completely ignores thermal
conduction. The temperature dependence of reaction enthalpies and heat
capacities has already been taken into account in the improved model by Allen
& Seyfried (2004). With these methods, only the increase in the temperature of
the surrounding material can be estimated. Heat loss during ice melting is
only taken into account in dynamic models (Cohen & Coker, 2000). In heat-
balance models, convection of liquid water is considered, while in dynamic
models, conduction in solid components and diffusion of molecules are
considered to be the main modes of heat transfer.

The heat transfer significantly affects the total heat balance and can be
influenced by a number of parameters, such as the porosity of the
planetesimal. The value of porosity in carbon-containing chondrites is around
20% (Consolmagno et al., 2008) while for smaller trans-Neptunian objects
(TNOs) it can reach up to 60% if they formed after the decay of most 26Al
(Bierson & Nimmo, 2019). The porosity can affect melting and heat transfer due
to the lack of continuity of the solid material as it reduces the conductive
heat transfer but promotes convection by allowing fluid to migrate in the
interconnected gaps. Higher porosity can increase the rate of the reaction, as
the effective surface area is larger, thus the contact surface of olivine and
H2O molecules also increases. In this way, the amount of interfacial liquid
water can also increase at temperatures below the melting point of water.
Furthermore, the interconnected porosity allows the liquid water to reach
places where it has already run out, thus allowing the reaction to continue.

The known meteorite samples originated mainly from the main asteroid belt
where the average densities of objects are much higher than in the outer Solar
System ($\rho$ $\geq$ 2 g cm-3 in contrast with the transneptunian region
where even densities of $\rho$ $\leq$ 1 g cm-3 are common). In large Kuiper
belt objects current estimates point to a common primordial (bulk) density of
$\sim$1.8 g cm-3 (Barr & Schwamb, 2016; Grundy et al., 2019), indicating a
rock-to-ice ratio of $\sim$42:58 in volume, and $\sim$70:30 in mass. This is
notably higher ice content than in the main belt, where e.g. the $\rho$ = 2.16
g cm-3 density of (1) Ceres (Park et al., 2016) indicates of rock-to-ice mass
ratio of $\sim$80:20, and the high $\rho$ = 3.5 g cm-3 density of (4) Vesta
(Russell et al., 2012) suggests a very low water ice content. If the
primordial density in the Kuiper belt was really close to the value obtained
from the bulk density of large objects ($\sim$1.8 g cm-3), small bodies in the
D $<$ 500 km range should have a notable macroporosity to have the observed
bulk densities below $\sim$1 g cm-3 (see e.g. Grundy et al., 2019, for a
recent evaluation). The abundance of ice and the level of porosity indicate
that aqueous alteration processes may have played an important role in the
evolution of these objects, at least for a short time, early in the evolution
of the Kuiper belt, when radiogenic decay provided enough heat for these
reactions to start.

In this paper, we present a revised model of the serpentinization process
inside planetesimals. This model is incorporated into a general thermal
evolution model. The main outline of the model is presented in Sect. 2. We
apply it to explore the role of this chemical process in the early evolution
of planetesimals in the Kuiper belt in Sect. 3. A detailed list of the
equations used is given in the Appendix.

##  2 Thermal evolution model of planetesimals

Our serpentinization model is based on Góbi & Kereszturi (2017) which is an
improved serpentinization model from heat balance (Lowell & Rona, 2002; Allen
& Seyfried, 2004) and dynamic models (Cohen & Coker, 2000). Góbi & Kereszturi
(2017) used a kinetic approach to estimate the rate of serpentinization and
follow its evolution. Their model implemented interfacial water but included
some simplifications and neglected several effects that did not cause a
significant difference in the results for the smaller object ($\sim$15 km
radius), low water ice content and low porosity they studied. In our cases,
however, these neglected effects are important as we aim to consider larger
planetesimals and icy bodies as well. In the description of the
serpentinization model (in Sect. 2.1. and 2.2.), we present which effects were
taken into account in the current model as an improvement, compared with the
previous model in which these parameters were simplified.

To create a complete heat development model, several other effects must be
considered which are important in the early solar system: the decay of
radionuclides (dominated by the short-lived 26Al), the accretion heat and the
time of accretion with respect to the onset of radiogenic decay. The heat-
generating capabilities of these heat sources were incorporated into the
improved serpentinization model, and the heat loss/cooling of the
planetesimals and the heat conduction were taken into account (in Sects. 2.3.
and 2.4.).

###  2.1 Modelling the serpentinization process

In this part of the paper we focus solely on the algorithm which is able to
follow the serpentinization process itself. We assume that the reaction
described by Eq. 1 is the dominant process for serpentinite production in the
planetesimals and we take only this reaction into account when we calculate
the heat released, following Góbi & Kereszturi (2017). The planetesimal used
in the model is made of seven components: silicate rocks are forsterite (Mg-
rich olivine), enstatite (Mg pyroxen), hydrated rock (serpentinite), non-
reactive solid material and H2O in the three-phase state: liquid, solid and
void space filled with H2O vapor. The melting point of H2O (268 K) is obtained
from the properties of a saturated solution of MgSO4 (Kargel, 1998), which was
used in model in earlier studies, too (Cohen & Coker, 2000; Góbi & Kereszturi,
2017). While MgSO4 is quite common in some chondrites (Zolensky et al., 1999)
in the case of TNOs other solutes could also be considered. However, the
effect of any salt or salt combination on the melting temperature of water is
probably small.

To be comparable with the Góbi & Kereszturi (2017) model, the serpentinization
algorithm itself is tested in a selected layer of a sphere, with the
homogeneous and isotropic distribution of the components within that layer.
Due to the simple forward Euler scheme, the model remains sensitive to the
choice of the time step. We explore this in Sect. 2.2 to find the optimal
(largest allowable) time step for our calculations. Apart from the
improvements discussed below, this is the same scheme as was used in Góbi &
Kereszturi (2017) (see the model scheme in Fig 1.).

Figure 1: The serpentinization model scheme.

#### Heat capacity

Heat properties of materials and the average values of the physical properties
of planetesimal depend on the temperature and composition (see in Sect. A.1).
We calculated the heat capacity of H2O in a different, more accurate way than
in the base model. After fixing it, the final heat production was higher by a
few Kelvin. However, the difference remained within the uncertainty range,
even at longer timescales.

#### Lithospheric pressure

Several quantities depend on the total pressure which is calculated as the sum
of lithospheric and vapor pressures (see in Sect. A.3).

Góbi & Kereszturi (2017) used a simplification to determine the lithospheric
pressure due to the small size (R = 15 km) of their main test objects, as the
contribution of the lithospheric pressure to total pressure is insignificant
for small objects. This was replaced by Eq. A35 which gives reliable results
for larger objects, too. Using this method the reaction rate becomes lower
which makes a more significant difference for larger objects. We compare our
results with the Góbi & Kereszturi (2017) litospheric pressure values in Fig.
2.

Figure 2: Lithospheric pressure in the center of the test object versus its
radius. Orange diamonds represent the simplified calculation of lithosphere
pressure (Góbi & Kereszturi, 2017) while the black curve shows the results of
the calculations using our more accurate model.

#### Interfacial liquid water

To obtain the heat gain we need to know the exact amount of liquid water and
the latent heat of water vaporization (Sect. A.4). The latent heat was
calculated in a way different from the base model and it caused a few tenths
of Kelvin difference in the final temperature, with a very small amount of
vapor formed under typical conditions. In those cases when the temperature is
lower than the melting point of H2O, ice is present instead of liquid water
which would not allow serpentinization to take place. However, Góbi &
Kereszturi (2017) pointed out that the presence of microscopic scale
interfacial water is possible on the surface of olivine grains at low
temperatures making serpentinization possible, and it may gradually melt the
icy surrounding of silicate particles due to the heat produced in this
exothermic reaction. The model of Góbi & Kereszturi (2017) allowed the
formation of more interfacial water than it would have actually been possible
at the given temperature. They did not take into account the rate of
serpentinization, although this effect is significant in the temperature range
below the melting point of bulk ice. In our model we consider four values for
the amount of water: i) the actual amount of ice in the layer examined, ii)
the maximum value of interfacial water that can be formed, iii) the amount of
water that is able to react in a specific step, and iv) the amount of ice that
can be melted by the serpentinization heat (see Eq. A44). The minimum of these
four values determines the actual amount of interfacial water which can be
formed.

In the Góbi & Kereszturi (2017) paper the improper handling of the interfacial
water also caused an issue in selecting the proper simulation time step: if
larger time steps for the same timespan were used, the calculated initial
temperatures required for the same serpentinization level were lower and the
final temperatures higher, causing a greater overall heat production in the
simulations. The reason behind this phenomenon was that in the case of
decreased time steps the serpentinization reaction started later in time. By
incorporating the new condition in the melting ice calculation, the
simulations did not produce more interfacial water any more than it was
possible at the given temperature, and there was no difference in the amount
of heat produced as a function of time scale at the same time.

###  2.2 Serpentinization model results compared with previous results

To demonstrate the capabilities of our improved algorithm to follow the
serpentizination process we used the same setup as Góbi & Kereszturi (2017),
but with the corrections listed above. We consider the deepest 100 m radius,
without taking other heat sources and heat/material transfer into account. In
each case the amount of non-reacting material of the test object was constant,
14% of the volume of planetesimals. We investigated how the different
parameters of the model can influence the outcome compared with the Góbi &
Kereszturi (2017) results.

#### Initial temperature

The output of any evolutionary model that includes chemical reactions like
serpentinization is very sensitive to the choice of initial temperature. We
examined the temperature increase ($\Delta$T) during the reaction as well as
the time needed to consume 90% of the reagent material during the
serpentinization reaction ($t_{90}$) at different initial temperature (Tini)
values (Fig. 3). The temperature dependence of chemical reactions is known in
general: the higher the initial temperature, the faster the reaction. The
dependence of the temperature increase on the initial temperature is a result
of the temperature dependence of the heat capacity of the constituent minerals
(equations for calculating the heat capacities can be found in Sect. A.1).

Figure 3: Serpentinization time (t90) and temperature increase ($\Delta$T)
during the serpentinization process. Orange diamonds represent the results by
Góbi & Kereszturi (2017) for an object with 15 km radius and the black dots
represents our improved model for a same object, using a time step of 0.1
year. Both test objects have 16% porosity and 1:2 olivine-to-water ratio. Here
we examined the innermost 100 m radius of the object, to be comparable with
the previous results.

The serpentinization time t90 is significantly longer than in the previous
studies (Fig. 3), and the serpentinization rate is higher towards higher
initial temperatures. Mainly due to the accurate calculation of the
lithospheric pressure and the consideration of the maximum value of sub-
freezing interfacial water t90 increased. In the case of a planetesimal with a
radius of 15 km, the process can be up to 3-4 times longer than in the
previous calculations, depending on the initial temperature.

This is due to the fact that only a small amount of interfacial liquid water
is present in the system below the freezing point and its production depends
on the amount of water reacted, thus the process is slow. When the temperature
reaches the melting point the heat production is entirely used to melt the
ice. Above the melting point, the reaction speeds up due to the accessibility
of large amounts of liquid water. Starting the reaction below the melting
temperature, most of the heat is consumed by melting the ice and this results
in an overall smaller temperature increase.

#### Importance of microscopic liquid water

Interfacial liquid water is an important component as it promotes the progress
of the reaction in the early stage when the temperature is below the melting
point of ice. Compared with previous results, it can be seen that in the sub-
freezing initial temperature range, the process is slower and the heat
production rate is also lower (Fig. 3). The reaction rate increases after the
melting point.

#### Size of planetesimals

The reaction rate increases with increasing pressure and it depends both on
the size of the test object through the lithospheric pressure and the vapor
pressure (Sect. A.3). The pressure dependence of serpentinization was studied
previously (Martin & Fyfe, 1970; Wegner & Ernst, 1983; Jones & Brearley, 2006;
Cohen & Coker, 2000), and it was found to be nearly linear in the range of
1-200 MPa which corresponds to the pressure expected in the size range in our
investigation.

Figure 4: The initial temperature (T90) and temperature increase ($\Delta$T)
for different sizes of planetesimals (robj) in $t_{90}$, using 1:2 olivine to
water ratio. The filled circles and the orange diamonds mark the results from
this present work and from the previous study by (Góbi & Kereszturi, 2017),
respectively.

The initial temperature required for serpentinization to consume 90% of the
materials in 10,000 years (T90) reaches values similar to those determined in
the previous work (Góbi & Kereszturi, 2017), the deviation is only a few K
(see Fig. 4). The largest deviation occurs at the size/pressure value where
the initial temperature falls below the freezing point.

There is an increasing trend in temperature change ($\Delta$T) with increasing
size (see Fig. 4) due to heat capacities which are lower due to the lower
initial temperature (T90). When the required initial temperature falls below
the freezing point the heat production rate decreases significantly as some of
the heat produced is used to melt ice. From this point on, the $\Delta$T again
shows an increasing trend as a function of object size.

#### Olivine to H2O ratio and initial porosity

The component ratio of olivine-to-water is one of the most important initial
parameters of the serpentinization reaction which significantly influences
both the course and the outcome of the reaction. In the case of higher
olivine-to-water ratios, the reaction can be faster than in the cases of
higher water content (see e.g. Fig. 1. in Góbi & Kereszturi, 2017). For each
olivine-to-water ratio value, a higher T90 is needed (see e.g. in Fig. 4), but
the difference is only a few K. The rate of heat production is closely related
to heat capacity which increases with water content as the heat capacity of
water is higher than that of rocky components. This slows down the reaction
towards higher water content.

#### Time step

As we use a simple forward-Euler scheme in our model, the results are expected
to be sensitive to the time step applied. Choosing a large timestep ($\Delta$t
= 10–100 yr, or larger) results in considerable instability in the
calculations, i.e. the final results (e.g. T90 or $\Delta$T we investigated
above) do not show a consistent trend, and depend strongly on the time step
chosen. This annoying effect disappears when the timestep is decreased, and
the calculations become stable for $\Delta$t $\lesssim$ 1 yr (see Fig. 5). To
avoid this problem we used a timestep of $\Delta$t = 0.5 yr. Reducing the
timestep further did not cause a considerable change in the final results.

Figure 5: Upper panel: temperature increase due to serpentinization as a
function of the time required (when one of the reagent materials is 100%
consumed). Different colors show different model time steps. Bottom panel, the
time required for chemical reaction as a function of the model time steps.

###  2.3 Decay of radionuclides

The main source of heat within planetesimals is the radioactive decay of both
short- and long-lived radionuclides. In Fig. 6 we plot the total energy output
as a function of time for some short-lived radionuclides. Because
serpentinization is a rapid process, only a few 10,000 years, we investigated
the early stages of the heat development of planetesimals. During this period,
as shown in Fig. 6, the most significant portion of the heat production from
radioactive decay was provided by 26Al (Lugaro et al., 2018), and only this
isotope was considered in our thermal evolution model as a heat source, with
the following parameters: half-life is t${}_{1/2}=0.717$ Myr, total energy is
TE$=5.07\times 10^{6}$ J/kg, and we assume that ${}^{26}Al$ was homogeneously
distributed in the early solar system (Lugaro et al., 2018).

Figure 6: Total energy output as a function of time for some short-lived
radionuclides. Time scale starting 50,000 years after the time when the
radionuclides start to decay.

###  2.4 Heat transfer

There are two major ways of heat transfer inside a planetary body: thermal
conduction and convection. Both of them were considered, and the thermal
radiation from the surface was also included in the thermal evolution model.

#### Thermal radiation

The power output was calculated from the Stefan-Boltzmann law using a constant
emissivity factor ($\epsilon=0.9$), and assuming an ambient temperature of
Tamb = 50 K which is a typical surface temperature in the outer Solar System
(Eq. A50).

#### Thermal conduction

The heat conduction rate was estimated with the Eq. A51, in all boundaries of
all layers following Hussmann et al. (2006).

#### Thermal convection

If the temperatures within the layers are sufficiently high and water is
present in a liquid state convective motions can occur in the porous media
(Hewitt et al., 2014). The requirement for starting convection is that the
Rayleigh number (Ra) exceeds a specific value. Ra is calculated in a porous
material as:

$\displaystyle\begin{split}Ra&=\frac{\rho\,\beta_{0}\,\Delta
T\,k_{0}\,g\,dr}{\eta\,\alpha}\\\ \end{split}$ (3)

where $\rho$ is the density in the layer, $\beta_{0}$ is the thermal expansion
coefficient and we used $\beta_{0}\,=\,10^{-3}$ as a safe upper limit for any
of the possible constituents, $\Delta$T is the temperature difference across
distance $dr$, which is the thickness of the layer, $k_{0}$ is the
permeability, and we used two values: $10^{-12}\,m^{2}$ following Cohen &
Coker (2000) and $10^{-7}\,m^{2}$ following Jacob (1972). The results were the
same because the Ra did not reach the critical value required for the start of
convection with any of the permeability values. $g$ is the local gravitational
acceleration, $\alpha=k/(\rho\cdot c_{p})$ is the thermal diffusivity ($k$ is
thermal conductivity and $c_{p}$ is specific heat capacity) and $\eta$ is the
dynamic viscosity of the fluid:

$\displaystyle\eta(T)=\eta_{0}\,exp(25\times(T_{melt}/T-1))$ (4)

where $\eta_{0}=10^{13}$ Pa s is the melting point viscosity (Hussmann et al.,
2006) and Tmelt is the melting point. A critical Rayleigh number of 1000 was
used for the onset of convection following Hussmann et al. (2006). However
this critical value can be as low as 40 (Nield & Bejan, 1999) in a porous
medium. In our simulations these Ra values have not been reached at T
$\lesssim$ Tmelt.

The main purpose of our present work is to follow the evolution of the
serpentinization process, and to determine the conditions and onset timescales
this process can operate at. As we show it below, before and during the active
serpentinization stage the temperature in all of our simulations is close to
or below the melting point of ice, and the heat transfer is governed by these
’icy’ conditions. During serpentinization the melted water is quickly consumed
by the chemical reaction, preventing convection (which would be caused by the
otherwise significantly decreased viscosity of liquid water). After the end of
the serpentinization stage radiogenic decay can still provide extra heat to
further melt the ice. In these conditions the low viscosity of liquid water
may lead to convection, generating a faster heat transfer within the inner
layers and may also lead to rearrangements in the structure of these layers.
The consideration of this process, as well as e.g that of the diffusion of
water into solid grains, are beyond the scope of this paper.

##  3 Application to transneptunian objects

###  3.1 Main model outline

Serpentine is a primary candidate to explain the reflectance spectra of
transneptunian objects (e.g. Protopapa et al., 2009), and it has also been
considered as primordial material in impact simulations aimed to explain the
formation of transneptunian binary system (e.g. the Pluto-Charon system;
Canup, 2005). Therefore it is an intriguing question how serpentine can form
under the conditions in the transneptunian region where objects have typically
high ice contents and possibly high porosity, and how far it can contribute to
the heat budget and thermal evolution of these planetesimals.

We consider a spherical body with spherical symmetry (all variables depend on
the radius only), with a size-dependent number of layers. We assumed a
homogeneous and isotropic composition at the start, and a common Tini in the
whole planetesimal. The effect of deviations from this latter assumption is
investigated in Sect. 3.3. Heat production by radiogenic decay and heat
transfer is considered as described in Sects. 2.3 and 2.4.

###  3.2 Initial parameters

#### Composition

The low density of TNOs is due partly to their composition, as the H2O content
of these distant objects is significantly higher than those in the inner Solar
system, but it is also due to their higher porosity. In the case of smaller
TNOs, when the radius does not exceed $\sim$150 km, the porosity can reach 60%
(Bierson & Nimmo, 2019). Such a high porosity is only possible if they formed
after the decay of 26Al to maintain their high porosity as there is no
internal transformation by heat at this stage. In the absence of accurate
knowledge of the internal composition, a simplified composition was used,
which has been used in previous serpentinization studies (See in Góbi &
Kereszturi, 2017; Cohen & Coker, 2000, and in the A).

#### Porosity

Porosity is a very important factor in this model, the initial temperature
(Tini) and the serpentinization time also strongly dependent on it. To
estimate the porosity, we used the calculations from Yasui and Arakawa (2009)
where they determine a size/pressure-dependent porosity for small icy bodies.
The equation defined for the range Regime 3 (P $>$ 2 MPa) was applied due to
the pressure conditions in our objects:

$\displaystyle\phi=a_{3}P^{b_{3}}$ (5)

where P is the lithospheric pressure (in MPa) in a middle layer in the
measured bodies and $a_{3}$ and $b_{3}$ are constants. We used the approximate
value of $a_{3}=0.5$ and $b_{3}=-0.2$ (see Yasui and Arakawa, 2009). With
these calculations, the porosity values are between 13 and 34% in the examined
size range for the 42:58 rock/water ratio specified earlier (see Fig. 7)

#### Temperature

In the early solar system the initial temperature of the objects was
determined by the size, composition, accretion heat, and the heat produced by
the decay of short-lived radioactive nuclei. The heating from solar
irradiation may be important closer to the Sun but it is negligible in the
outer regions. In our model, the initial temperature was taken to be the heat
from the accretion (Hanks & Anderson, 1969):

$\displaystyle T_{acc}=\frac{3}{5}\frac{GM}{RC_{p}}$ (6)

where $G$ is the gravitational constant, $C_{p}$ is the average heat capacity,
$M$ is the mass and $R$ is the radius of the planetesimal (see Fig. 7.). At t
= 0 the planetesimal has a homogeneous temperature distribution with Tacc.
Note that as $C_{p}$ is temperature-dependent Eq. 6 is an implicit equation
for Tacc, and it is solved iteratively; it also determines the shape of the
Tacc curve presented in Fig. 7.

Figure 7: The accretion heat product (Tacc) and porosity versus the size of
the planetesimal (robj), for the case of 42:58 rock/water ratio with size
dependent porosity.

In a next step we tested the sensitivity of the radioactive decay heat product
on the time step and layer thickness. In these tests we used the rock/water
ratio of 42:58 that corresponds to an olivine to water ratio of 0.12.
Radiogenic heat production was not very sensitive to the variation to time
step and layer thickness.

In all simulation configurations, we used a non-reactive rock content of 14%
with a density of 3630 kg m-3 like the earlier studies (Góbi & Kereszturi,
2017). In each case, the olivine water ratio was taken as 0.12, which
corresponds to the assumed rock-water ratio of 42:58 in TNOs. The calculations
are performed with the size-dependent porosity (Eq. 5.), the initial
temperature is calculated from the accretion heat (Eq. 6.), and we choose a
time step of 0.5 yr and a layer thickness of 20 km. The starting time t = 0
was considered to be the end of the accretion.

###  3.3 The effect of serpentinization in the thermal history of trans-
Neptunian objects

Objects with a radius below 150 km are expected to have high porosity (up to
60%) and the initial (accretion) temperature would remain very low ($<$ 42 K),
as indicated by our calculations. These objects may have formed after the
depletion of 26Al (as discussed in Bierson & Nimmo, 2019), and the
serpentinization reaction could likely not produce a significant amount of
heat. Therefore we did not consider these objects in our further calculations.
Very large objects (R $>$ 1000 km) were also excluded from the study size
range because they are likely formed by multiple accretion events and undergo
fast chemical differentiation early in their evolution.

Figure 8: Heat production by serpentinization without radioactive decay during
the first 100,000 years, for different object sizes, which corresponds to
different porosities and initial temperatures as discussed in Sect3.2. The
black symbols represent the initial temperature (due to accretion heat) and
the red symbols show the final temperature in the center of the planetesimals.

We examined the conditions under which serpentinization can produce
significant heat over a few tens of thousands of years in the size range of
$150\leq R\leq 1000$ km. We found that for objects smaller than 600 km, the
initial temperature from the accretion will be so low that the reaction cannot
start or will be very slow, even in the core, in the absence of radiogenic
decay. This appears in Fig. 8 as a ’jump’ in temperature at this specific
size. In larger objects, the process is already able to produce a notable
amount of heat without the contribution of radioactive elements. Above an
initial temperature of $\sim$150 K, assuming the previously determined
material composition and pressure-dependent porosity (Fig. 7) the reaction
rate increases significantly (Fig. 8). As the size increases, both the initial
temperature and the internal pressure increase, resulting in an increase in
the rate of serpentinization. This results in the extension of the central
region where the serpentinization process can produce significant heat in a
few tens of thousand years. The upper panel in Fig. 9 shows the ratio of the
internal volume in which 90% olivine has been consumed in the first 100,000
years for a planetesimal with a specific size. For radii R $\geq$ 650 km the
serpentinized zone extends almost to the surface while it remains in the core
below R = 600 km, in agreement with the slow reaction rate indicated by the
small temperature increase seen in Fig.8. The bottom panel in Fig. 9 shows an
example of how serpentinization proceeds inside an object of R = 620 km in the
first 100,000 years. As shown also in the upper panel, a maximum of 480 km
radius of the serpentinized region is reached corresponding to a maximum
volume ratio of 46%. This R $\approx$ 600 km critical size limit for the
efficient progress of the serpentinization process was obtained without the
consideration of radiogenic decay, and it already indicates that
serpentinization can be an efficient process for the large objects (the dwarf
planets of the trans-Neptunian region) even if they formed after 107 yr or
later, when the heat produced by radiogenic decay is significantly lower (see
Fig. 6). This also suggests that the size limit of efficient serpentinization
is smaller when an additional heat source – radogenic decay – is considered.

Figure 9: Upper panel: ratio of internal volume in which 90% olivine has been
consumed in the first 100,000 years versus the objects size. Bottom panel:
Progress of the serpentinization reaction, presented with the olivine
consumption as a function of depth at different times, for a planetesimal of R
= 620 km. The different colors mark 35,000 (green); 50,000 (blue); 75,000
(blue) and 100,000 (black) years.

Figure 10: Temperature evolution inside planetesimals with radii R = 500, 340
and 200 km (top-left, top-right and bottom-left, respectively), at different
depth in their interior (20 and 80 km, and at their core, as indicated in the
figures). The green curves correspond to a thermal evolution without
serpentinization at the core of the objects (this is representative for most
layers due to the slow heat transport). In each subfigure three points are
marked: t0: the start of the serpentinization process; t1: the time when that
temperature reaches the melting point of water ice in the core, considering
serpentinization; t2: the time when that temperature reaches the melting point
of water ice without serpentinization. The vertical lines represents the end
of the chemical reaction. In the bottom-right subfigure we present the results
for an R = 200 km object, but assuming that the evolution starts 0.7 Myr
later, when the decay of 26Al produces only half of the heat compared with the
object in the bottom-left subfigure.

Figure 11: Upper panel: Contribution of accretion heat (red), radiogenic decay
(green) and serpentinization (blue) to the final temperature of the test
object at the time when the melting temperature of the water ice is reached in
the center; Middle panel: t90 and t0 versus object radius; Bottom panel: The
ratio of the times necessary to reach the melting temperature of water ice,
without and with serpentinization (t2/t1), as a function of the radius of the
object.

To investigate this, we modeled the thermal evolution of test objects with
radii between 150 $<$ R $<$ 500 km considering all heat sources (accretion
heat, radiogenic decay) in addition to serpentinization itself. In Fig. 10 we
present the thermal evolution of three objects with radii 200, 340 and 500 km.
The serpentinization process starts to produce a significant amount of heat
when the temperature of the layer exceeds a critical value of $\sim$180 K, at
time t0 after the start of the simulation. This t0 depends strongly on the
size of the object (see Fig. 11, middle panel) and decreases with size, in a
similar way as the end time of the reaction (t90), represented by vertical
lines in Fig. 10. In these simulations serpentinization completes fully in all
the studied layers even at the melting point, and therefore it does not
contribute to the future thermal evolution. We note that all our simulations
have been run beyond the time of completion of the reaction. In our
simulations serpentinization reduces the time needed to reach the melting
temperature of water by a factor of 1.6-2.6, depending on the size of the
object (Fig. 11, bottom panel). The contribution of the different processes to
the temperature of the object at the end of serpentinization (t90) is
presented in Fig. 11 (upper panel). At smaller sizes radiogenic decay
contributes the most, and the importance of accretion heat increases notably
with size. The heat obtained from serpentinization (50-80 K) becomes more
important with growing size, and exceeds the contribution from radiogenic
decay for R $\approx$ 500 km. In large enough objects (where the critical
reaction-starting temperature of $\sim$180 K is reached) serpentinization
proceeds quickly, and the whole process in finished in $\sim$104 yr.

The results of Wakita and Sekiya (2011) shows that the formation time is very
important because 2.4 Myr after Ca-Al-rich inclusions formation, the icy
planetesimals may not reach the melting temperature of ice. We also examined
the effect of formation time on the thermal evolution, considering the same
setup as above.

In Fig. 10, bottom-right panel we demonstrate the effect of a delayed start
(by 0.7 Myr, the half-life of 26Al) and hence reduced radiogenic heat. In this
simulation the $t_{0}$, $t_{1}$ and $t_{2}$ timescales are notably, by a
factor of $\sim$2 longer.

Despite these longer timescales, serpentinization can still fully proceed, and
the temperature reaches the melting point of ice. As it is expected, the
contribution from radiogenic decay to the final temperature is smaller (11 K
in this specific case) and the relative contribution of the temperature
increase due to serpentinization is higher.

Figure 12: t0 (solid curve), t1 (dashed) and t90 (dash-dotted) times of the
simulations of the thermal evolution of an R = 200 km object versus the
starting time relative to the original starting time when 26Al abundance was
at its maximum. Formation time (the start time of the simulation) is presented
as multiples of the half-life of 26Al (t1/2). The last data point at
7$\times$0.717 Myr $\approx$ 5 Myr corresponds to the expected maximum
lifetime of the protoplanetary disk. Green, red and blue colours correspond to
20 and 80 km depths, and the core of the planetesimal, respectively.

In Fig. 12 we show the effect of a late formation of the planetesimals, and
consequently a late start of the serpentinization process, for different
starting times. The late formation results in a reduced amount of heat from
radiogenic decay, due to the depletion of 26Al. These simulations were
performed for R = 200 km-sized objects, assuming that the formation of the
planetesimal happened at a t = [0,1,…7]$\times$t1/2 after the onset of
isotopic decay, when the abundance of 26Al was at its maximum, as considered
in our previous simulations. Serpentinization needs a much longer time to
start, either in the core or in layers closer to the surface, for later
formation times. While this t0 is in the order of a few thousand years for an
early formation, it is several million years for a formation at t $\approx$5
Myr. The t90 timescale, i.e. the time when 90% of the potentially serpentine
forming material in the core of the planetesimal is consumed, is also notably
longer, while the ratio of t0 to t90 remains roughly the same. Altogether the
lower temperatures due to the smaller radiogenic heat slow down
serpentinization considerably.

We examined what changes occur when the initial temperature distribution is
inhomogeneous and the surface is warmer (Fig. 13). In this test we used an
object with a radius of 240 km and the initial temperature distribution was
set in a way that the surface temperature was 5/3 times the central
temperature (Hanks & Anderson, 1969). Despite the lower starting temperature
the serpentinization process is the fastest in the core due to the higher
lithospheric pressure. The reaction takes longer in the outer layers, and as
the initial temperature is higher, the final temperature will also be higher,
but the temperature difference is smaller than in the homogeneous case.

Figure 13: Thermal evolution of a planetesimal of R = 240 km with an outward
temperature gradient at the start, represented by curves with ’normal’
colours. The ’pale’ colours correspond to the same object/layer, but assuming
a homogeneous temperature distribution at the start of the calculations.

We also examined the case when the initial temperature distribution was
inhomogeneous and the formation occured in two steps (Fig. 14). First, an
object formed with a radius of 160 km, and $\tau_{f}$ = 10,000 or 100,000 yr
later the formation was completed and the object reached a final radius of 320
km. In the time between the two formation/accretion events the object may have
warmed up both from radiogenic decay and serpentinization. We assumed two
values for the time of the first formation event: ts = 0, the maximum 26Al
heat production date, and ts = 3 t1/2, i.e. three 26Al half-life later. For
$\tau_{s}$ = 100,000 yr and ts = 0 the serpentinization could go along all the
way before the second accretion event in the centre and left the core at a
high temperature. Due to the lack of further reactants, the core was heated by
radiogenic decay only from this point on (almost straight blue curve in Fig.
14, bottom left). The final object is built on this warm core in the second
accretion event. Also in our other three cases, the results show that
serpentinization is faster in the core independently of the initial
conditions. The difference between the core and the outer layers is more
significant for ts = 0, and, as expected, everything occurs significantly
later for ts = 3 t1/2 yr.

Figure 14: In this figure we present the temperature evolution obtained in
models with a two-phase formation scenario. In the top row the first
formation/accretion event is followed by another event in $\tau_{f}$ = 10,000
yr, while in the bottom row in $\tau_{f}$ = 100,000 yr. On the left the first
formation event occurs early, at the maximum of the radiogenic heat production
of 26Al; on the right the first formation event occurs three 26Al half-life
later. t = 0 corresponds to the second formation event.

##  4 Conclusions

We presented an improved algorithm to model the serpentinization reaction and
its role in the thermal evolution of planetesimals in the early Solar system,
based on the model by Góbi & Kereszturi (2017). In our model we incorporated
several previously overlooked or neglected effects: (i) the calculation of
heat capacity of water is more accurate in this work; (ii) we have taken into
account the depth dependence of lithospheric pressure; (iii) the latent heat
of vaporization of water was calculated in a more exact way, and (iv) we
improved the method to calculate the amount of microscopic liquid interfacial
water at subzero temperatures which eventually allowed a smaller amount of
water to react compared with previous model. Our model is able to follow the
local chemical evolution of the serpentinization process and the algorithm was
inserted into a more complex internal heat evolution model which considered
radioactive decay and heat transfer. We demonstrated that the improvements in
our models lead to results different from those of earlier models, both in
serpentinization reaction timescales and in final heat production. The
presence of interfacial water at temperatures below the melting point of bulk
water ice may be able to start the reaction, though at a lower rate, and may
eventually be able to produce a notable temperature increase.

Our results suggest that there is a size limit of R $\approx$ 600 km (assuming
Kuiper belt compositions) above which the serpentinization process becomes
efficient, mainly due to the higher accretion heat of these larger objects.
This is even true in later times when radiogenic decay cannot significantly
contribute to the thermal budget because of the depletion of 26Al, several
million years after the onset of radiogenic decay. Serpentinization may
proceed in smaller objects after the onset of radiogenic decay (down to
$\sim$150 km) in the presence of notable heat from radiogenic decay. The
overall importance of serpentinization in the chemical evolution of smaller (R
$\leq$ 500 km) objects in the outer Solar system depends strongly on the time
the process starts at, as accretion heat alone cannot start the process for
these smaller planetesimals. The lifetime of the planetesimal forming disks
around solar-like young stars is $\sim$3-5 Myr, and very few disks survive 10
Myr (see e.g. Russell et al., 2006, for a summary). As we showed above, the
reaction timescales and the overall efficacy depend strongly on the heat
provided by radiogenic decay, and it is reduced strongly if planetesimal
formation is delayed. However, serpentinization is still able to proceed,
although with a notably reduced speed, even if the object is formed a few
million years later, or in multiple accretion events, at the smallest sizes we
investigated (R = 200 km).

The bulk serpentinization efficiency – the ratio of objects in which
serpentinization reformed the interior and those where it could potentially do
it so – also depends strongly on the collisional evolution of planetesimals in
the young trans-Neptunian region. While serpentization is a fast process in
the case of an early formation, it can be considerably slower if the formation
process is delayed. A simple estimate based on Wyatt (2008) shows that for
objects with R $\geq$ 100 km destructive collisions occur on the million year
timescales for a wide range of possible disk parameters (disk mass, disk
extension, mean orbital eccentricity, material strength, etc.), and this
timescale gets longer with the disk dispersal. In this sense the
serpentinization timescale is expected to be much shorter than the collisional
timescale at any time while the disk exists. In our Solar system the trans-
Neptunian population of objects with radii R $>$ 30 km are expected to be
primordial, as they ’decoupled’ early from the collisional evolution, at the
end of the runaway growth, and before the onset of destructive collisions of
smaller objects (Schlichting et al., 2013). This suggests that our R = 150-500
km objects, once formed, are likely not destroyed by collisions, and
serpentinization could take place in their interior.

## Acknowledgements

This research has been supported by the K-125015, K-138962 and K-138594 grants
of the Hungarian Research, Development, and Innovation Office (NKFIH). We are
indebted to our reviewers for their comments and suggestions which have helped
us to notably improve this paper.

## References

  * Allen & Seyfried (2004) Allen, D.E., Seyfried, W.E., Jr 2004, Geochim. Cosmochim. Acta, 68, 1347
  * Anderson et al. (1973) Anderson, D.M., Tice, A.R., McKim,H.L. 1973, in Heider C.M., ed., 2nd Int. Conf., North American Contribution: Permafrost. Printing and Publishing Office, National Academy of Sciences, Washington D.C., p. 289
  * Barr & Schwamb (2016) Barr, Amy C.; Schwamb, Megan E. 2016, MNRAS, 460, 1542
  * Bear (1972) Bear, Jacob 1972, Dynamics of Fluids in Porous Media, Dover
  * Bierson & Nimmo (2019) Bierson, C.J., Nimmo, F. 2019, Icarus, 326, 10
  * Canup (2005) Canup, R.M. 2005, Science, 307, 546
  * Cohen & Coker (2000) Cohen, B.A. and Coker, R.F. 2000, Icarus, 145, 369
  * Consolmagno et al. (2008) Consolmagno, G.J., Britt, D.T., Macke, R.J. 2008 Chem. Erde, 68, 1
  * Desch et al. (2009) Desch, S.J., Cook, J.C., Doggett, T.C., Porter, S.B. 2009, Icarus, 202, 694
  * Dufresne & Anders (1962) Dufresne, E.R., Anders, E. 1962, Geochim. Cosmochim. Acta, 26, 1085
  * Fyfe (1974) Fyfe, W.S. 1974, Geophys. J. R. Astron. Soc., 37, 213
  * Gail et al. (2014) Gail, H.-P., Trieloff, M., Breuer, D., Spohn, T., Early Thermal Evolution of Planetesimals and Its Impact on Processing and Dating of Meteoritic Material, Protostars and Planets VI, Henrik Beuther, Ralf S. Klessen, Cornelis P. Dullemond, and Thomas Henning (eds.), University of Arizona Press, Tucson, 914 pp., p.571-593
  * Góbi & Kereszturi (2017) Góbi, S. and Kereszturi, Á. 2017, MNRAS, 466, 2099
  * Grimm & Mcsween (1989) Grimm, R.E., and Mcsween, H.M., JR. 1989, Icarus, 82, 244
  * Grundy et al. (2019) Grundy, W.M., Noll, K.S., Buie, M.W., et al. 2019, Icarus, 334, 30
  * Hanks & Anderson (1969) Hanks, T.C. and Anderson, D.L. 1969, PEPI, 2, 19
  * Hewitt et al. (2014) Hewitt, D.R., Neufeld, J.A. and Lister, J.L. 2014, Journal of Fluid Mechanics, 748, 879
  * Hussmann et al. (2006) Hussmann, H., Sohl, F., Spohn, T. 2006, Icarus, 185, 258
  * Jacob (1972) Bear, Jacob, 1972. Dynamics of Fluids in Porous Media, Dover
  * Jones & Brearley (2006) Jones, C.L., Brearley, A.J. 2006, Geochim. Cosmochim. Acta, 70, 1040
  * Kargel (1998) Kargel, J.S. 1998, ASSL, 227, 3
  * Lowell & Rona (2002) Lowell, R.P., Rona, P.A. 2002, Geophys. Res. Lett., 29, 26
  * Lugaro et al. (2018) Lugaro, M., Ott, U., Kereszturi, Á. 2018, PrPNP, 102, 1L
  * Lunine (2006) Lunine, J.I. 2006, Meteorites and the Early Solar System II. Univ. Arizona Press, Tucson, AZ, p. 309
  * Malamud & Prialnik (2013) Malamud, U., Prialnik, D., 2013, Icarus, 225, 763
  * Malamud et al. (2017) Malamud, U., Perets, H.B., Schubert, G. 2017, MNRAS, 468, 1056
  * Martin & Fyfe (1970) Martin, B., Fyfe, W.S. 1970, Chem. Geol., 6, 185
  * Müller et al. (2009) Müller, T.G., Lellouch, E., Böhnhardt, H., és mtsai 2009, EM&P, 105, 209
  * Nield & Bejan (1999) Nield, A. and Bejan, A. 1999, Convection in Porous Media, Springer, New York
  * Park et al. (2016) Park, R.S., Konopliv, A.S., Bills, B.G., et al. 2016, Nature, 537, 515
  * Protopapa et al. (2009) Protopapa, S., Alvarez-Candal, A., Barucci, M.A., et al. 2009, A&A, 501, 375
  * Robertson (1988) Eugene C. Robertson, 1988, Thermal properties of Rocks, U.S. Geological Survey
  * Robie & Waldbaum (1968) Robie, R.A., Waldbaum, D.R. 1968, U.S. Geol. Surv. Bull., 1259, 1
  * Russell et al. (2006) Russell, S.S., Hartmann, L., Cuzzi, J., et al. 2006, Meteorites and the Early Solar System II, D. S. Lauretta and H. Y. McSween Jr. (eds.), University of Arizona Press, Tucson, 943 pp., p.233-251
  * Russell et al. (2012) Russell, C.T., Raymond, C. A., Coradini, A., et al. 2012, Science, 336, 684
  * Scott et al. (1988) Scott, E.R.D., Barber, D.J., Alexander, C.M., et al. 1988, in Meteorites and the Early Solar System, eds. J.F. Kerridge & M.S. Matthews (Tucson: Univ. Arizona Press), 718
  * Schlichting et al. (2013) Schlichting, H.E., Fuentes, C.I., Trilling, D.E. 2013, AJ, 146, 36
  * Teiser et al. (2011) Teiser, J., Engelhardt, I., Wurm, G. 2011, ApJ, 742, 5
  * Wakita and Sekiya (2011) Shigeru Wakita and Minoru Sekiya 2011, EP & S,63, 1193
  * Wegner & Ernst (1983) Wegner, W.W., Ernst, W.G. 1983, Am. J. Sci. A, 283, 151
  * Wyatt (2008) Wyatt, M.C., 2008, Annu. Rev. Astron. Astrophys. 46, 339-383
  * Yasui and Arakawa (2009) Yasui, M., Arakawa, M. 2009, JGRE, 114, 9004Y
  * Zega et al. (2003) Zega, T.J.; Garvie, L.A. J. and Buseck, P.R. 2003, AmMin, 88, 1169
  * Zolensky et al. (1989) Zolensky, M.E., Bourcier W.L., and Gooding, J.L. 1989, Icarus, 78, 41
  * Zolensky et al. (1999) Zolensky, M.E., Bodnar, R.J., Gibson Jr., E.K., et al. 1999, Science, 285, 1377

##  Appendix A Serpentinization model

In this model a specific object was considered to be made of seven separate
components: olivine, enstatite, non-reactive solid, serpentinite, and water in
three-phase state as liquid, solid, and vapor in the pore space. These
materials are marked in the subscripts as: $oli$, $ens$, $nre$, $ser$, $wat$,
$ice$ and $vap$. Table 2 presents the variables and constants used in the
serpentinization model. In those cases when the initial values weren’t
specified otherwise, we used the initial values from Table 2.

Table 1: List of variables and constants 1. $C_{p}$ | Heat capacity ($Jkg^{-1}K^{-1}$)
---|---
$g$ | Gravitational acceleration ($m/s^{2}$)
$\textit{k}_{r}$ | Observable serpentinization rate (mol/year)
m | Mass of components ($kg$)
n | Amount of substance (mol)
$\rho$ | Density ($kg/m^{3}$)
$\phi$ | Porosity
$r_{lt}$ | Top of the layer ($m$)
$r_{lb}$ | Bottom of the layer ($m$)
$P_{lit}$ | Lithostatic pressure (Pa)
$P_{vap}$ | Pressure of water vapor (Pa)
$T$ | Temperature ($K$)
$t$ | Time ($year$)
V | Volume of components ($m^{3}$)
$\Delta H_{v}$ | Latent heat of vaporization of water ($J/kg$)
$\Delta H_{r}$ | reaction enthalpy ($J/mol$)
$\Delta h_{vap}$ | vaporization heat ($J$)
$\Delta h_{ser}$ | reaction heat ($J$)
$\Delta n_{ser}$ | Serpentine produced ($mol$)
$\Delta T_{ser}$ | Temperature rise from the serpentinization
$w$ | Content of microscopic liquid water ($g/100\,gsoil$)
$\lambda$ | decay constant ($sec^{-1}$)
$q_{rad}$ | radiogenic heat production rate ($W/kg$)
$Q_{rad}$ | radiogenic heating rate ($J$)
$k$ | thermal konductivity ($W/m/K$)
Table 2: List of variables and constants 2. |  | constant | initial values
---|---|---|---
$G$ | Gravitational constant ($m^{3}kg^{-1}s^{-2}$) | $6.67408\times 10^{-11}$ |
$\varepsilon$ | Emissivity | 0.9 |
$\sigma$ | Stefan-Boltzmann Constant ($W/m^{2}/K^{4}$) | $5.6697\times 10^{-8}$ |
$L$ | Latent heat of fusion of ice ($Jkg^{-1}$) | $3.3\times 10^{5}$ |
$R_{g}$ | Gas constant ($Jmol^{-1}K^{-1}$) | $8.314$ |
$\rho_{oli}$ | Density of olivine ($kgm^{-3}$) | 3210 |
$\rho_{ens}$ | Density of enstatite ($kgm^{-3}$) | 3190 |
$\rho_{nre}$ | Density of non-reactive solid ($kgm^{-3}$) | 3630 |
$\rho_{ser}$ | Density of serpentinite ($kgm^{-3}$) | 2470 |
$S$ | Specific surface area ($m^{2}g^{-1}$) | 100 |
$T_{melt}$ | Melting point of water ($K$) | 268 K |
$t_{1/2}$ | Half-life of ${}^{26}Al$ ($Myr$) | 0.717 |
$TE$ | Total energy of ${}^{26}Al$ ($J/kg$) | $5.07\times 10^{6}$ |
$k_{oli}$ | thermal conductivity of olivine ($W/m/K$) | 5.155 |
$k_{ens}$ | thermal conductivity of enstatite ($W/m/K$) | 5.155 |
$k_{ser}$ | thermal conductivity of serpentinite ($W/m/K$) | 2.95 |
$k_{nre}$ | thermal conductivity of non-reactive solid ($W/m/K$) | 2.8 |
$T_{amb}$ | Ambient temperature ($K$) | 50 |
$n_{oli}/n_{wat}$ | Olivine to water ratio |  | 1:2, 0.12
$n_{nre}$ | Non-reactive material ($\%$) |  | 14
$R$ | Radius of planetesimal ($km$) |  | 15 - 800
$\Delta t$ | Time step ($year$) |  | 0.5

###  A.1 Material properties: heat capacity and density of components

Heat capacity ($c_{p}$) and density ($\rho$) calculations were obtained from
Cohen & Coker (2000). These physical properties of all components are
temperature ($T$) dependent.

$\displaystyle log\,C_{p}\,(olivine)=-11.32+13.58\,x-4.25\,x^{2}+0.44\,x^{3}$
(A1) $\displaystyle
log\,C_{p}\,(enstatite)=-8.620+10.39\,x-3.00\,x^{2}+0.28\,x^{3}$ (A2)
$\displaystyle log\,C_{p}\,(water)=8.250-4.180\,x+1.12\,x^{2}-0.076\,x^{3}$
(A3)

where $x=log\,T$. These expressions are valid from approximately 50 to 500 K.

$\begin{split}C_{p}\,(serpentinite)&=1145+0.048\,T-2.65\times
10^{7}\,T^{-2}\quad(T\,>\,273\,K)\\\
log\,C_{p}\,(serpentinite)&=-0.59-1.51\,x+(2.82\,x^{2})-(0.66\,x^{3})\\\
&(T\,<\,273\,K)\end{split}$ (A4) $C_{p}\,(vapor)=1730.54+0.45\,T\\\ $ (A5)
$\begin{split}C_{p}\,(ice)&=152.46+7.12\,T\quad(T\,>\,150\,K)\\\
&=126.89+7.50\,T\quad(150\,K\,>\,T\,>\,95\,K)\\\
&=-49.97+9.5\,T\quad(95\,K\,>\,T\,>\,50\,K)\end{split}$ (A6)

For a specific layer of the object, we calculated the average heat capacity
from the components’ heat capacities and masses ($m$):

$\begin{split}C_{p}\,(layer)m_{lay}&=C_{p}\,(oli)\,m_{oli}+C_{p}\,(ens)\,m_{ens}+C_{p}\,(ice)\,m_{ice}\\\
&+C_{p}\,(wat)\,m_{wat}+C_{p}\,(vap)\,m_{vap}\\\
&+C_{p}\,(ser)\,m_{ser}+C_{p}\,(oli)\,m_{nre}\end{split}$ (A7)

In the cases of olivine, enstatite, serpentinite and non-reactive material,
constant values of density were used (see in Table 2). For H2O we used:

$\rho_{wat}=-221+13.1\,T-0.0507\,T^{2}+8.49\times 10^{-5}\,T^{3}-5.48\times
10^{-8}\,T^{4}$ (A8)
$\begin{split}\rho_{ice}&=-46.9\,x+1032.71\quad(T\,>\,137\,K)\\\
&=-1.32\,x+935.32\quad(T\,<\,137\,K)\end{split}$ (A9)
$\rho_{vap}=P_{vap}\,0.018\,T^{-1}\,R_{g}^{-1}$ (A10)

where $P_{vap}$ is the vapor pressure of the vapour and $R_{g}$ is the gas
constant. For the examined layer, we calculated the average density from the
components’ masses and volumes ($V$):

$\rho_{lay}=\frac{m_{oli}+m_{ens}+m_{ice}+m_{wat}+m_{vap}+m_{ser}+m_{nre}}{V_{oli}+V_{ens}+V_{ice}+V_{wat}+V_{ser}+V_{voi}+V_{nre}}$
(A11)

where $V_{voi}$ is the void space which is filled with H2O vapor.

###  A.2 Mass and volume of components

This part of the algorithm considers a specific layer in a body, at a specific
depth. The volume of the examined layer is as follows:

$V_{lay}=\bigg{(}\frac{4}{3}\,\pi\,r_{lt}^{3}\bigg{)}-\bigg{(}\frac{4}{3}\,\pi\,r_{lb}^{3}\bigg{)}$
(A12)

where $V_{lay}$ is the volume of the examined layer, $r_{lt}$ is the top of
the layer and $r_{lb}$ is the bottom of the layer.

The initial values of volumes and masses of components (marked with $(0)$)
were calculated from the porosity ($\phi$), olivine-to-water ratio
($n_{oli}/n_{wat}$) and the fraction of non-reactive material ($n_{nre}$):

$\displaystyle V_{voi}(0)=V_{lay}\,\phi$ (A13) $\displaystyle
V_{nre}=V_{lay}\,n_{nre}$ (A14) $\displaystyle m_{nre}=V_{nre}\,\rho_{nre}$
(A15) $\displaystyle m_{oli}/m_{wat}(0)=n_{oli}/n_{wat}\,140.71/18$ (A16)
$\displaystyle m_{lay}=V_{lay}\,\rho_{lay}$ (A17)
$\displaystyle\begin{split}V_{ens}(0)&=(V_{lay}-(V_{nre}+V_{voi}))\frac{m_{oli}/m_{wat}/1.3/\rho_{ens}}{\frac{m_{oli}/m_{wat}}{\rho_{oli}}+\frac{1}{\rho_{ice}}+\frac{m_{oli}/m_{wat}/1.3}{\rho_{ens}}}\quad(T\,<\,T_{melt})\\\
&=(V_{lay}-(V_{nre}+{V_{voi}})\frac{m_{oli}/m_{wat}/1.3/\rho_{ens}}{\frac{m_{oli}/m_{wat}}{\rho_{oli}}+\frac{1}{\rho_{wat}}+\frac{m_{oli}/m_{wat}/1.3}{\rho_{ens}}}\quad(T\,>\,T_{melt})\end{split}$
(A18)
$\displaystyle\begin{split}V_{ice}(0)&=(V_{lay}-(V_{nre}+V_{voi}))\frac{1/\rho_{ice}}{\frac{m_{oli}/m_{wat}}{\rho_{oli}}+\frac{1}{\rho_{ice}}+\frac{m_{oli}/m_{wat}/1.3}{\rho_{ens}}}\quad(T\,<\,T_{melt})\\\
&=0\quad(T\,>\,T_{melt})\\\ \end{split}$ (A19)
$\displaystyle\begin{split}V_{wat}(0)&=0\quad(T\,<\,T_{melt})\\\
&=(V_{lay}-(V_{nre}+V_{voi}))\frac{1/\rho_{wat}}{\frac{m_{oli}/m_{wat}}{\rho_{oli}}+\frac{1}{\rho_{wat}}+\frac{m_{oli}/m_{wat}/1.3}{\rho_{ens}}}\quad(T\,>\,T_{melt})\end{split}$
(A20)
$\displaystyle\begin{split}V_{oli}(0)&=(V_{lay}-(V_{nre}+V_{voi}))\frac{m_{oli}/m_{wat}/\rho_{oli}}{\frac{m_{oli}/m_{wat}}{\rho_{oli}}+\frac{1}{\rho_{ice}}+\frac{m_{oli}/m_{wat}/1.3}{\rho_{ens}}}\quad(T\,<\,T_{melt})\\\
&=(V_{lay}-(V_{nre}+V_{voi}))\frac{m_{oli}/m_{wat}/\rho_{oli}}{\frac{m_{oli}/m_{wat}}{\rho_{oli}}+\frac{1}{\rho_{wat}}+\frac{m_{oli}/m_{wat}/1.3}{\rho_{ens}}}\quad(T\,>\,T_{melt})\end{split}$
(A21)

The melting point of ice ($T_{melt}$ = 268 K) is obtained considering a
saturated solution of MgSO4.

The masses of components (enstatite, ice, water, olivine) are calculated in
the zeroth time step in the following way:

$\displaystyle m_{j}(0)=V_{j}(0)\,\rho_{j}$ (A22)

where $j$ subscript refers to the materials. In all other time steps these
masses are obtained as:

$\displaystyle m_{ser}=m_{ser}^{i-1}+\frac{277.1\,\Delta n_{ser}^{i-1}}{1000}$
(A23) $\displaystyle m_{oli}=m_{oli}^{i-1}-(\Delta n_{ser}^{i-1}\times
0.14071)$ (A24) $\displaystyle m_{ens}=m_{ens}^{i-1}-(\Delta
n_{ser}^{i-1}0\times 0.10039)$ (A25) $\displaystyle
m_{ice}=max(0,m_{ice}^{i-1}-\Delta m_{ice}^{i-1})$ (A26) $\displaystyle
m_{wat}=m_{wat}^{i-1}+\Delta m_{ice}^{i-1}-2\,\Delta n_{ser}^{i-1}\times
0.018)$ (A27) $\displaystyle m_{vap}=\frac{0.018\,P_{vap}\,V_{voi}}{R_{g}\,T}$
(A28)

where the $i-1$ superscript refers to the value of the variable in the
previous time step.

where $\Delta n_{ser}$ is the serpentine produced at the previous moment,
$\Delta m_{ice}$ is the changing mass of ice due to the melting. We calculated
the volumes of components from the masses and densities:

$\displaystyle V=m/\rho$ (A30) $\displaystyle
V_{voi}=V_{voi}^{i-1}+\frac{3}{10^{6}}\times\Delta n_{ser}^{i-1}$ (A31)

Amount of substance of water ($n_{wat}$) and olivine ($n_{oli}$):

$\displaystyle n_{wat}=1000\times m_{wat}/18$ (A33) $\displaystyle
n_{oli}=1000\times m_{oli}/140.71$ (A34)

###  A.3 Amount of serpentinite

To obtain the total pressure, the sum of the lithospheric ($P_{lit}$) and
vapor pressures ($P_{vap}$) we first calculate the gravitational acceleration
($g(l)$) and then $P_{lit}$ in that specific layer:

$\displaystyle\begin{split}g(l)&=G\frac{m_{r}}{r^{2}}\\\
P_{lit}&=\sum_{l=r}^{R}\rho_{lay}(l)g(l)dr\end{split}$ (A35)

where $G$ is the gravitational constant, $r$ is the radius of the specific
layer, $m_{r}$ is mass within the radius $r$ , $R$ is the radius of the
planetesimal and $dr$ is the thickness of the layer.

From the Clausius-Clapeyron relation, the vapor pressure has the approximate
form

$\displaystyle P_{vap}=P_{0}\,e^{\frac{T_{0}}{T}}$ (A36)

where $P_{0}=3.58\times 10^{12}\,Pa$ and $T_{0}=-6140\,K$ in the presence of
ice, while it is changed to $P_{0}=4.7\times 10^{10}\,Pa$ and $T_{0}=-4960\,K$
in the presence of water (Grimm & Mcsween, 1989). In the case of mixed phases
of H2O, the mass-weighted average of $P_{vap}$ is used Góbi & Kereszturi
(2017):

$\displaystyle\begin{split}P_{vap}(0)&=\frac{m_{ice}\,P_{0}\,e^{\frac{T_{0}}{T}}+m_{wat}\,P_{0}\,e^{\frac{T_{0}}{T}}}{m_{ice}+m_{wat}}\\\
P_{vap}&=\frac{m_{ice}^{i-1}\,P_{0}\,e^{\frac{T_{0}}{T}}+m_{wat}^{i-1}\,P_{0}\,e^{\frac{T_{0}}{T}}}{m_{ice}^{i-1}+m_{wat}^{i-1}}\end{split}$
(A37)

The reaction rate depends linearly on pressure and exponentially on
temperature. By knowing the serpentinization rate ($\textit{k}_{r}$), the
number of serpentine produced $\Delta n_{ser}$ can then be calculated in
moles:

$\displaystyle\textit{k}_{r}=4383\,\frac{P_{lit}+P_{vap}}{10^{8}}\,e^{\frac{-3463}{T}}$
(A38) $\displaystyle\Delta n_{ser}=min(n_{oli}(1-e^{-kr\times\Delta
t}),\frac{n_{wat}}{2})$ (A39)

When olivine is in excess, then $n_{wat}/2$ is used as being the limiting
reagent.

###  A.4 Heat budget and temperature increase

To determine the extent of heat production it is necessary to know the heat of
reaction:

$\displaystyle\Delta h_{ser}=\Delta H_{r}\,\Delta n_{ser}$ (A40)

where $\Delta H_{r}=69\,kJ/mol$ is the reaction enthalpy Robie & Waldbaum
(1968)

The vaporation heat is calculated as follows:

$\displaystyle\Delta
H_{v}=3713997.2-7822.6569\,T+17.613373\,T^{2}-0.019018061\,T^{3}$ (A41)
$\displaystyle\Delta h_{vap}=\Delta H_{v}\,\Delta m_{vap}$ (A42)

where $\Delta H_{v}$ is the latent heat of vaporization of water $T$ is capped
at 400 K and $\Delta m_{vap}=m_{vap}-m_{vap}^{i-1}$.

In those cases when the initial temperature is lower than the melting point of
the ice, a certain part of ice can transform into interfacial water Anderson
et al. (1973).

$\displaystyle ln\,w=0.2618+0.5519\,ln\,S-1.449\,S^{-0.264}\,ln(T_{melt}-T)$
(A43)

where $w$ is the content of microscopic liquid water (in g/100 g soil) and $S$
is specific surface area (100 $m^{2}\,g^{-1}$ based on Anderson et al.
(1973)).

In the next step, we calculated the amount of ice that is converted to water
($\Delta m_{ice}$) considering the actual amount of ice in the layer, the
possible amount of interfacial water that can be formed, the amount of water
that is able to react in this specific step, and the amount of ice that can be
melted by serpentinization.

$\displaystyle\begin{split}\Delta
m_{ice}(0)&=min(m_{ice},w\times\frac{m_{oli}+m_{ens}+m_{nre}}{100\times\rho_{wat}})\quad(T\,<\,T_{melt})\\\
&=\frac{\Delta h_{ser}-\Delta h_{vap}}{L}\quad(T\,>\,T_{melt}\quad
and\quad\Delta h_{ser}\,>\,0)\\\ &=0\quad(T\,>\,T_{melt}\quad and\quad\Delta
h_{ser}\,<\,0)\\\ \Delta
m_{ice}&=min(m_{ice},w\times\frac{m_{oli}+m_{ens}+m_{nre}+m_{ser}}{100\,\rho_{wat}},\\\
&2\,\Delta n_{ser}\times 0.018),\frac{\Delta h_{ser}-\Delta
h_{vap}}{L})\quad(T\,<\,T_{melt})\\\ &=\frac{\Delta h_{ser}-\Delta
h_{vap}}{L}\quad(T\,>\,T_{melt}\quad and\quad\Delta h_{ser}\,>\,0)\\\
&=0\quad(T\,>\,T_{melt}\quad and\quad\Delta h_{ser}\,<\,0)\end{split}$ (A44)

where $L$ is the latent heat of the ice-water phase transition

$\displaystyle\Delta h_{ice/wat}=L\Delta m_{ice}$ (A45)

The temperature increase ($\Delta T_{ser}$) due to serpentinization is
calculated as follows:

$\displaystyle\begin{split}\Delta T_{ser}(0)&=\frac{\Delta
h_{ser}}{C_{p}\,(layer)\times m_{lay}}\\\ \Delta T_{ser}&=\frac{\Delta
h_{ser}-\Delta h_{ice/wat}-\Delta h_{vap}}{C_{p}\,(layer)\times
m_{lay}}\end{split}$ (A46)

###  A.5 Decay of radionuclides

We considered solely the ${}^{26}Al$ isotope when calculating the heat from
the radiogenic decay. The radiogenic heat production rate ($Q_{rad}$) is
obtained by the following equations (Desch et al., 2009):

$\displaystyle\begin{split}\lambda&=\frac{ln(2)}{t_{1/2}}\\\
q_{rad}&=TE\times\lambda\\\
Q_{rad}&=q_{rad}\,e^{(-t\times\lambda)}\,m_{rock}\,dt\end{split}$ (A47)

where $\lambda$ is the decay constant, $TE$ is the total energy and $q_{rad}$
is the radiogenic heat production rate (see Table 2).

###  A.6 Heat transfer

#### Thermal conduction

The thermal conductivity of the rocky material components was considered to be
independent of the temperature (Cohen & Coker, 2000; Robertson, 1988). We
calculated the thermal conductivity of the different phases of H2O in the
following way:

$\displaystyle\begin{split}k_{ice}&=9.828\,exp(-0.0057\,T)\\\
k_{vap}&=-0.0143+1.02\times 10^{-4}T\\\ k_{wat}&=-0.581+6.34\times
10^{-3}T-7.93\times 10^{-6}T^{2}\quad(T\,<\,410\,K)\\\
k_{wat}&=0.9721(-0.142+4.12\times 10^{-3}T-5.01\times
10^{-6}T^{2})\quad(T\,>\,410\,K)\end{split}$ (A48)

When determining the average thermal conductivity of the investigated layer, a
composite rock with a homogeneous material distribution was considered. We
calculated the parallel bulk rock conductivity ($k_{p}$, its grains arranged
in a parallel orientation to the direction of heat flow) and the series
conductivity ($k_{s}$, its grains arranged in a layered sequence perpendicular
to the heat flow direction). The mean values of $k_{p}$ and $k_{s}$ fit well
with the observed values, especially for more porous rocks (Robertson, 1988):

$\displaystyle\begin{split}k_{p}&=nV_{1}k_{1}+nV_{2}k_{2}+nV_{3}k_{3}+...\\\
\frac{1}{k_{s}}&=\frac{nV_{1}}{k_{1}}+\frac{nV_{2}}{k_{2}}+\frac{nV_{3}}{k_{3}}+...\\\
k_{lay}&=\frac{k_{p}+k_{s}}{2}\end{split}$ (A49)

where $nV_{1}$, $nV_{2}$, $nV_{3}$ … are fractional volumes of components
($nV_{component}=V_{component}/V_{lay}$).

#### Heat radiation

On the surface of the body the radiated heat is calculated as:

$\displaystyle\begin{split}Q_{th}&=\epsilon\sigma\,(T_{surf}^{4}-T_{amb}^{4})\,A\,dt\end{split}$
(A50)

where $\epsilon$ is the emissivity factor ($\epsilon=0.9$), $T_{surf}$ the
temperature of the surface, $T_{amb}$ is the ambient temperature, which is
assumed to be 50 K, a typical surface temperature of airless bodies due to
solar irradiation in the outer Solar System and $A$ is the radiating surface
area.

#### Thermal conduction

The heat transferred by thermal conduction is calculated as:

$\displaystyle\begin{split}Q_{cond}&=A\,k_{lay}\,\frac{\Delta
T}{dr}\,dt\end{split}$ (A51)

where $\Delta T$ is the temperature difference between the adjacent layers.

#### Thermal convection

If the Rayleigh number exceeded a critical value we also calculated the
Nusselt number and considered the convective heat flow in our calculations:

$\displaystyle\begin{split}Nu&=\left(\frac{Ra}{Ra_{crit}}\right)^{\beta}\\\
Q_{conv}&=\frac{k_{lay}}{dr}Nu\,A\,(T(l+1)-T(l))\,dt\end{split}$ (A52)

where $\beta$ is a dimensionless value that can take values between 0.25 and
1/3 depending on the geometry and boundary conditions. As a reference value,
we use $\beta$ = 0.3 after Hussmann et al. (2006). We used Racrit = 1000 as
critical Rayleigh number.

The final heat equation is:

$\displaystyle Q_{sum}=\Delta h_{ser}-\Delta h_{ice/wat}-\Delta
h_{vap}+Q_{rad}-Q_{cool}+Q_{cond}+Q_{conv}$ (A53)

where $Q_{cool}=Q_{th}$ in the surface and in the inner layers
$Q_{cool}=Q_{cond}$ at the upper boundary of the examined layer.


* Völkel  _et al._ [2022b] S. H. Völkel, N. Franchini, E. Barausse, and E. Berti, Constraining modifications of black hole perturbation potentials near the light ring with quasinormal modes, Phys. Rev. D 106, 124036 (2022b), arXiv:2209.10564 [gr-qc] .
  * Franchini and Völkel [2023b] N. Franchini and S. H. Völkel, Parametrized quasinormal mode framework for non-Schwarzschild metrics, Phys. Rev. D 107, 124063 (2023b), arXiv:2210.14020 [gr-qc] .
  * Ansorg and Panosso Macedo [2016] M. Ansorg and R. Panosso Macedo, Spectral decomposition of black-hole perturbations on hyperboloidal slices, Phys. Rev. D 93, 124016 (2016), arXiv:1604.02261 [gr-qc] .
  * Panosso Macedo  _et al._ [2018] R. Panosso Macedo, J. L. Jaramillo, and M. Ansorg, Hyperboloidal slicing approach to quasi-normal mode expansions: the Reissner-Nordström case, Phys. Rev. D 98, 124005 (2018), arXiv:1809.02837 [gr-qc] .
  * Panosso Macedo [2020a] R. Panosso Macedo, Hyperboloidal framework for the Kerr spacetime, Class. Quant. Grav. 37, 065019 (2020a).
  * Zenginoglu [2008] A. Zenginoglu, Hyperboloidal foliations and scri-fixing, Class. Quant. Grav. 25, 145002 (2008), arXiv:0712.4333 [gr-qc] .
  * Zenginoglu [2011] A. Zenginoglu, A Geometric framework for black hole perturbations, Phys. Rev. D83, 127502 (2011), arXiv:1102.2451 [gr-qc] .
  * Destounis  _et al._ [2020b] K. Destounis, R. D. B. Fontana, and F. C. Mena, Accelerating black holes: quasinormal modes and late-time tails, Phys. Rev. D 102, 044005 (2020b), arXiv:2005.03028 [gr-qc] .
  * Ferrari and Mashhoon [1984] V. Ferrari and B. Mashhoon, New approach to the quasinormal modes of a black hole, Phys. Rev. D 30, 295 (1984).
  * Beyer [1999] H. R. Beyer, On the completeness of the quasinormal modes of the Poschl-Teller potential, Commun. Math. Phys. 204, 397 (1999), arXiv:gr-qc/9803034 .
  * Bizoń  _et al._ [2020] P. Bizoń, T. Chmaj, and P. Mach, A toy model of hyperboloidal approach to quasinormal modes, Acta Phys. Polon. B 51, 1007 (2020), arXiv:2002.01770 [gr-qc] .
  * Price [1972a] R. H. Price, Nonspherical perturbations of relativistic gravitational collapse. 1. Scalar and gravitational perturbations, Phys. Rev. D 5, 2419 (1972a).
  * Price [1972b] R. H. Price, Nonspherical Perturbations of Relativistic Gravitational Collapse. II. Integer-Spin, Zero-Rest-Mass Fields, Phys. Rev. D 5, 2439 (1972b).
  * Gundlach  _et al._ [1994b] C. Gundlach, R. H. Price, and J. Pullin, Late time behavior of stellar collapse and explosions: 2. Nonlinear evolution, Phys. Rev. D 49, 890 (1994b), arXiv:gr-qc/9307010 .
  * Hawking and Ellis [2023] S. W. Hawking and G. F. R. Ellis,  _The Large Scale Structure of Space-Time_, Cambridge Monographs on Mathematical Physics (Cambridge University Press, 2023).
  * Cardoso  _et al._ [2018a] V. Cardoso, J. a. L. Costa, K. Destounis, P. Hintz, and A. Jansen, Quasinormal modes and Strong Cosmic Censorship, Phys. Rev. Lett. 120, 031103 (2018a), arXiv:1711.10502 [gr-qc] .
  * Kim  _et al._ [2013] Y.-W. Kim, Y. S. Myung, and Y.-J. Park, Quasinormal modes and hidden conformal symmetry in the Reissner-Nordström black hole, Eur. Phys. J. C 73, 2440 (2013), arXiv:1205.3701 [hep-th] .
  * Zimmerman  _et al._ [2015] A. Zimmerman, H. Yang, Z. Mark, Y. Chen, and L. Lehner, Quasinormal Modes Beyond Kerr, Astrophys. Space Sci. Proc. 40, 217 (2015), arXiv:1406.4206 [gr-qc] .
  * Richartz [2016] M. Richartz, Quasinormal modes of extremal black holes, Phys. Rev. D 93, 064062 (2016), arXiv:1509.04260 [gr-qc] .
  * Cardoso  _et al._ [2018b] V. Cardoso, J. L. Costa, K. Destounis, P. Hintz, and A. Jansen, Strong cosmic censorship in charged black-hole spacetimes: still subtle, Phys. Rev. D 98, 104007 (2018b), arXiv:1808.03631 [gr-qc] .
  * Destounis [2019a] K. Destounis, Charged Fermions and Strong Cosmic Censorship, Phys. Lett. B 795, 211 (2019a), arXiv:1811.10629 [gr-qc] .
  * Liu  _et al._ [2021] H. Liu, W.-L. Qian, Y. Liu, J.-P. Wu, B. Wang, and R.-H. Yue, Alternative mechanism for black hole echoes, Phys. Rev. D 104, 044012 (2021), arXiv:2104.11912 [gr-qc] .
  * Aretakis [2011a] S. Aretakis, Stability and Instability of Extreme Reissner-Nordström Black Hole Spacetimes for Linear Scalar Perturbations I, Commun. Math. Phys. 307, 17 (2011a), arXiv:1110.2007 [gr-qc] .
  * Aretakis [2011b] S. Aretakis, Stability and Instability of Extreme Reissner-Nordström Black Hole Spacetimes for Linear Scalar Perturbations II, Annales Henri Poincare 12, 1491 (2011b), arXiv:1110.2009 [gr-qc] .
  * Angelopoulos  _et al._ [2018] Y. Angelopoulos, S. Aretakis, and D. Gajic, Horizon hair of extremal black holes and measurements at null infinity, Phys. Rev. Lett. 121, 131102 (2018), arXiv:1809.10037 [gr-qc] .
  * Onozawa  _et al._ [1997] H. Onozawa, T. Okamura, T. Mishima, and H. Ishihara, Perturbing supersymmetric black hole, Phys. Rev. D 55, 4529 (1997), arXiv:gr-qc/9606086 .
  * Okamura [1997] T. Okamura, On scattering off the extreme Reissner-Nordström black hole in N=2 supergravity, Phys. Rev. D 56, 4927 (1997), arXiv:gr-qc/9703050 .
  * Kallosh  _et al._ [1998] R. Kallosh, J. Rahmfeld, and W. K. Wong, One loop supergravity corrections to the black hole entropy and residual supersymmetry, Phys. Rev. D 57, 1063 (1998), arXiv:hep-th/9706048 .
  * Berti [2004] E. Berti, Black hole quasinormal modes: Hints of quantum gravity?, Conf. Proc. C 0405132, 145 (2004), arXiv:gr-qc/0411025 .
  * Moncrief [1974a] V. Moncrief, Odd-parity stability of a Reissner-Nordström black hole, Phys. Rev. D 9, 2707 (1974a).
  * Moncrief [1974b] V. Moncrief, Stability of Reissner-Nordström black holes, Phys. Rev. D 10, 1057 (1974b).
  * Moncrief [1975] V. Moncrief, Gauge-invariant perturbations of Reissner-Nordström black holes, Phys. Rev. D 12, 1526 (1975).
  * Liebling and Palenzuela [2023] S. L. Liebling and C. Palenzuela, Dynamical boson stars, Living Rev. Rel. 26, 1 (2023), arXiv:1202.5809 [gr-qc] .
  * Vaglio  _et al._ [2022] M. Vaglio, C. Pacilio, A. Maselli, and P. Pani, Multipolar structure of rotating boson stars, Phys. Rev. D 105, 124020 (2022), arXiv:2203.07442 [gr-qc] .
  * Destounis  _et al._ [2023b] K. Destounis, F. Angeloni, M. Vaglio, and P. Pani, Extreme-mass-ratio inspirals into rotating boson stars: nonintegrability, chaos, and transient resonances, (2023b), arXiv:2305.05691 [gr-qc] .
  * Mazur and Mottola [2023] P. O. Mazur and E. Mottola, Gravitational Condensate Stars: An Alternative to Black Holes, Universe 9, 88 (2023), arXiv:gr-qc/0109035 .
  * Carballo-Rubio [2018] R. Carballo-Rubio, Stellar equilibrium in semiclassical gravity, Phys. Rev. Lett. 120, 061102 (2018), arXiv:1706.05379 [gr-qc] .
  * Arrechea  _et al._ [2021] J. Arrechea, C. Barceló, R. Carballo-Rubio, and L. J. Garay, Semiclassical constant-density spheres in a regularized Polyakov approximation, Phys. Rev. D 104, 084071 (2021), arXiv:2105.11261 [gr-qc] .
  * Kaplan and Rajendran [2019] D. E. Kaplan and S. Rajendran, Firewalls in General Relativity, Phys. Rev. D 99, 044033 (2019), arXiv:1812.00536 [hep-th] .
  * Mathur [2005] S. D. Mathur, The Fuzzball proposal for black holes: An Elementary review, Fortsch. Phys. 53, 793 (2005), arXiv:hep-th/0502050 .
  * Mathur [2008] S. D. Mathur, Fuzzballs and the information paradox: A Summary and conjectures, (2008), arXiv:0810.4525 [hep-th] .
  * Ikeda  _et al._ [2021] T. Ikeda, M. Bianchi, D. Consoli, A. Grillo, J. F. Morales, P. Pani, and G. Raposo, Black-hole microstate spectroscopy: Ringdown, quasinormal modes, and echoes, Phys. Rev. D 104, 066021 (2021), arXiv:2103.10960 [gr-qc] .
  * Cardoso and Pani [2019b] V. Cardoso and P. Pani, Testing the nature of dark compact objects: a status report, Living Rev. Rel. 22, 4 (2019b), arXiv:1904.05363 [gr-qc] .
  * Cardoso  _et al._ [2014] V. Cardoso, L. C. B. Crispino, C. F. B. Macedo, H. Okawa, and P. Pani, Light rings as observational evidence for event horizons: long-lived modes, ergoregions and nonlinear instabilities of ultracompact objects, Phys. Rev. D 90, 044069 (2014), arXiv:1406.5510 [gr-qc] .
  * Keir [2016] J. Keir, Slowly decaying waves on spherically symmetric spacetimes and ultracompact neutron stars, Class. Quant. Grav. 33, 135009 (2016), arXiv:1404.7036 [gr-qc] .
  * Cunha  _et al._ [2023] P. V. P. Cunha, C. Herdeiro, E. Radu, and N. Sanchis-Gual, Exotic Compact Objects and the Fate of the Light-Ring Instability, Phys. Rev. Lett. 130, 061401 (2023), arXiv:2207.13713 [gr-qc] .
  * Zhong  _et al._ [2022] Z. Zhong, V. Cardoso, and E. Maggio, On the instability of ultracompact horizonless spacetimes, (2022), arXiv:2211.16526 [gr-qc] .
  * Maggio  _et al._ [2020] E. Maggio, L. Buoninfante, A. Mazumdar, and P. Pani, How does a dark compact object ringdown?, Phys. Rev. D 102, 064053 (2020), arXiv:2006.14628 [gr-qc] .
  * Vellucci  _et al._ [2022] V. Vellucci, E. Franzin, and S. Liberati, Echoes from backreacting exotic compact objects, (2022), arXiv:2205.14170 [gr-qc] .
  * Jackiewicz and Renaut [2002] Z. Jackiewicz and R. Renaut, A note on stability of pseudospectral methods for wave propagation, Journal of Computational and Applied Mathematics 143, 127 (2002).
  * Cardoso and Gualtieri [2016] V. Cardoso and L. Gualtieri, Testing the black hole ‘no-hair’ hypothesis, Class. Quant. Grav. 33, 174001 (2016), arXiv:1607.03133 [gr-qc] .
  * Cardoso and Pani [2017] V. Cardoso and P. Pani, Tests for the existence of black holes through gravitational wave echoes, Nature Astron. 1, 586 (2017), arXiv:1709.01525 [gr-qc] .
  * Dolan [2007] S. R. Dolan, Instability of the massive Klein-Gordon field on the Kerr spacetime, Phys. Rev. D 76, 084001 (2007), arXiv:0705.2880 [gr-qc] .
  * Rosa and Dolan [2012] J. G. Rosa and S. R. Dolan, Massive vector fields on the Schwarzschild spacetime: quasi-normal modes and bound states, Phys. Rev. D 85, 044043 (2012), arXiv:1110.4494 [hep-th] .
  * Cardoso  _et al._ [2011] V. Cardoso, S. Chakrabarti, P. Pani, E. Berti, and L. Gualtieri, Floating and sinking: The Imprint of massive scalars around rotating black holes, Phys. Rev. Lett. 107, 241101 (2011), arXiv:1109.6021 [gr-qc] .
  * Pani  _et al._ [2012] P. Pani, V. Cardoso, L. Gualtieri, E. Berti, and A. Ishibashi, Perturbations of slowly rotating black holes: massive vector fields in the Kerr metric, Phys. Rev. D 86, 104017 (2012), arXiv:1209.0773 [gr-qc] .
  * Destounis [2019b] K. Destounis, Superradiant instability of charged scalar fields in higher-dimensional Reissner-Nordström-de Sitter black holes, Phys. Rev. D 100, 044054 (2019b), arXiv:1908.06117 [gr-qc] .
  * Vieira and Kokkotas [2021] H. S. Vieira and K. D. Kokkotas, Quasibound states of Schwarzschild acoustic black holes, Phys. Rev. D 104, 024035 (2021), arXiv:2104.03938 [gr-qc] .
  * Vieira  _et al._ [2022] H. S. Vieira, K. Destounis, and K. D. Kokkotas, Slowly-rotating curved acoustic black holes: Quasinormal modes, Hawking-Unruh radiation, and quasibound states, Phys. Rev. D 105, 045015 (2022), arXiv:2112.08711 [gr-qc] .
  * Vieira  _et al._ [2023] H. S. Vieira, K. Destounis, and K. D. Kokkotas, Analog Schwarzschild black holes of Bose-Einstein condensates in a cavity: Quasinormal modes and quasibound states, Phys. Rev. D 107, 104038 (2023), arXiv:2301.11480 [gr-qc] .
  * Cunningham  _et al._ [1980] C. T. Cunningham, R. H. Price, and V. Moncrief, Radiation from collapsing relativistic stars. III - Second order perturbations of collapse with rotation, ApJ 236, 674 (1980).
  * Gleiser  _et al._ [1996] R. J. Gleiser, C. O. Nicasio, R. H. Price, and J. Pullin, Colliding black holes: How far can the close approximation go?, Phys. Rev. Lett. 77, 4483 (1996), arXiv:gr-qc/9609022 .
  * Miller  _et al._ [2016] J. Miller, B. Wardell, and A. Pound, Second-order perturbation theory: the problem of infinite mode coupling, Phys. Rev. D 94, 104018 (2016), arXiv:1608.06783 [gr-qc] .
  * Speri  _et al._ [2023] L. Speri, A. Antonelli, L. Sberna, S. Babak, E. Barausse, J. R. Gair, and M. L. Katz, Probing Accretion Physics with Gravitational Waves, Phys. Rev. X 13, 021035 (2023), arXiv:2207.10086 [gr-qc] .
  * Speeney  _et al._ [2022] N. Speeney, A. Antonelli, V. Baibhav, and E. Berti, Impact of relativistic corrections on the detectability of dark-matter spikes with gravitational waves, Phys. Rev. D 106, 044027 (2022), arXiv:2204.12508 [gr-qc] .
  * Akiyama  _et al._ [2019] K. Akiyama _et al._ (Event Horizon Telescope), First M87 Event Horizon Telescope Results. I. The Shadow of the Supermassive Black Hole, Astrophys. J. Lett. 875, L1 (2019), arXiv:1906.11238 [astro-ph.GA] .
  * Akiyama  _et al._ [2022] K. Akiyama _et al._ (Event Horizon Telescope), First Sagittarius A* Event Horizon Telescope Results. I. The Shadow of the Supermassive Black Hole in the Center of the Milky Way, Astrophys. J. Lett. 930, L12 (2022).
  * Lynden-Bell [1969] D. Lynden-Bell, Galactic nuclei as collapsed old quasars, Nature 223, 690 (1969).
  * Gurzadian and Ozernoi [1979] V. G. Gurzadian and L. M. Ozernoi, Accretion on massive black holes in galactic nuclei, Nature 280, 214 (1979).
  * Armitage [2004] P. J. Armitage, Theory of disk accretion onto supermassive black holes, Astrophys. Space Sci. Libr. 308, 89 (2004), arXiv:astro-ph/0405170 .
  * Dittmann and Miller [2020] A. J. Dittmann and M. C. Miller, Star formation in accretion discs and SMBH growth, Mon. Not. Roy. Astron. Soc. 493, 3732 (2020), arXiv:1911.08685 [astro-ph.HE] .
  * Sofue  _et al._ [2009] Y. Sofue, M. Honma, and T. Omodaka, Unified Rotation Curve of the Galaxy – Decomposition into de Vaucouleurs Bulge, Disk, Dark Halo, and the 9-kpc Rotation Dip –, Publ. Astron. Soc. Jap. 61, 227 (2009), arXiv:0811.0859 [astro-ph] .
  * Catena and Ullio [2010] R. Catena and P. Ullio, A novel determination of the local dark matter density, JCAP 08, 004, arXiv:0907.0018 [astro-ph.CO] .
  * Weber and de Boer [2010] M. Weber and W. de Boer, Determination of the local dark matter density in our galaxy, Astronomy and Astrophysics 509, A25 (2010).
  * Mark  _et al._ [2017] Z. Mark, A. Zimmerman, S. M. Du, and Y. Chen, A recipe for echoes from exotic compact objects, Phys. Rev. D 96, 084002 (2017), arXiv:1706.06155 [gr-qc] .
  * Krivan  _et al._ [1997] W. Krivan, P. Laguna, P. Papadopoulos, and N. Andersson, Dynamics of perturbations of rotating black holes, Phys. Rev. D 56, 3395 (1997), arXiv:gr-qc/9702048 .
  * Pazos-Ávalos and Lousto [2005] E. Pazos-Ávalos and C. O. Lousto, Numerical integration of the teukolsky equation in the time domain, Physical Review D 72, 10.1103/physrevd.72.084022 (2005).
  * Zenginoglu and Khanna [2011] A. Zenginoglu and G. Khanna, Null infinity waveforms from extreme-mass-ratio inspirals in Kerr spacetime, Phys. Rev. X 1, 021017 (2011), arXiv:1108.1816 [gr-qc] .
  * Panosso Macedo [2020b] R. Panosso Macedo, Hyperboloidal framework for the Kerr spacetime, Class. Quant. Grav. 37, 065019 (2020b), arXiv:1910.13452 [gr-qc] .