# Active doping controls the mode of failure in dense colloidal gels

Tingtao Zhou  John F. Brady  Division of Chemistry and Chemical Engineering
and Division of Engineering and Applied Science, California Institute of
Technology

###### Abstract

Mechanical properties of disordered materials are governed by their underlying
free energy landscape. In contrast to external fields, embedding a small
fraction of active particles within a disordered material generates non-
equilibrium internal fields, which can help to circumvent kinetic barriers and
modulate the free energy landscape. In this work, we investigate through
computer simulations how the activity of active particles alters the
mechanical response of deeply annealed polydisperse colloidal gels. We show
that the ’swim force’ generated by the embedded active particles is
responsible for determining the mode of mechanical failure, i.e., brittle vs.
ductile. We find, and theoretically justify, that at a critical swim force the
mechanical properties of the gel decrease abruptly, signaling a change in the
mode of mechanical failure. The weakening of the elastic modulus above the
critical swim force results from the change in gel porosity and distribution
of attractive forces among gel particles. Below the critical swim force, the
ductility enhancement is caused by an increase of gel structural disorder.
Above the critical swim force, the gel develops a pronounced heterogeneous
structure characterized by multiple pore spaces, and the mechanical response
is controlled by dynamical heterogeneities. We contrast these results with
those of a simulated monodisperse gel that exhibits a non-monotonic trend of
ductility modulation with increasing swim force, revealing a complex interplay
between the gel energy landscape and embedded activity.

## Significance

Adaptive materials that can alter their mechanical properties in response to
an external stimulus are highly desirable. Specifically, yielding of
disordered materials is a nonlinear process governed by their complex free
energy landscape. Using Active Brownian Dynamics simulations, we demonstrate
precise control of colloidal gel structure, elasticity and ductility by tuning
the properties of a remarkably small fraction of embedded active particles,
specifically their propulsive swim force. Our findings show that active doping
is an effective approach to modulating the mechanical properties of the host
material, providing guidance for designing novel reconfigurable materials.

## Introduction

Biological matter has the amazing capability of adapting to its environment.
For instance, cells alter their stiffness and viscosity in response to
different substrate rigidities Lo  _et al._ (2000) by polarizing their actin
cytoskeleton network Doss  _et al._ (2020). Adaptations generally render more
robust response to environmental changes and enhance the chances of survival
for living organisms. Inspired by biology, can one endow passive synthetic
materials with similar controlled mechanical adaptivity? In this work, we
consider a simple system of passive colloidal particles where a small fraction
of them can be activated by external stimuli to self-propel, i.e., become
active Brownian particles (ABPs) Fily and Marchetti (2012). The persistent
motion of ABPs distinguishes them from passive Brownian particles and can be
understood as a source of mechanical force or stress, i.e. the swim stress
Takatori  _et al._ (2014); Omar  _et al._ (2020). In fact, it is possible to
extract mechanical work from self-propelling particles and transmit free
energy to the surrounding or to other components of the system Angelani  _et
al._ (2009); Reichhardt and Reichhardt (2017); Ro  _et al._ (2022). This
mechanical perspective suggests that embedded activity may be harnessed to
modulate the mechanical properties of a host material.

Recent advances in colloidal physics and chemistry provide the ability to
engineer colloidal particle shapes Damasceno  _et al._ (2012) and interaction
potentials Wang  _et al._ (2012), which have led to the self-assembly of
various structured colloidal materials with different symmetries and
hierarchies Vlasov  _et al._ (2001); Liu  _et al._ (2014); Engel  _et al._
(2015); Wang  _et al._ (2017); He  _et al._ (2020). As an alternative to this
passive self-assembly, active doping – namely embedding a small fraction of
active particles into the host material – can provide a non-equilibrium route
for directed self-assembly that may overcome kinetic arrest and frustration in
passive self-assembly. For example, adding a small amount of active particles
can promote the crystallization of colloidal hard spheres Ni  _et al._ (2014)
or attractive particles Mallory  _et al._ (2020), and anneal defects in dense
colloidal monolayers Ramananarivo  _et al._ (2019); Saud and Solomon (2023).
In disordered colloidal gels, Omar et al. Omar  _et al._ (2018) showed that
the dynamical free energy landscape can be engineered by changing the self-
propulsive swim force and the length and time scales of the persistent motion
of active particles, resulting in widely different gel structures.

The mechanics and rheology of dense gels and glasses are governed by the
inherent structures of the energy landscape Sastry (2001); Sastry  _et al._
(1998); Heuer (2008) and by dynamical heterogeneities Berthier (2011);
Berthier  _et al._ (2011). Linear response properties in disordered hard
sphere packings are known to relate to the glass transition Wittmer  _et al._
(2013); Yoshino and Zamponi (2014) and follow universal scalings Goodrich  _et
al._ (2016) at the jamming point. As a nonlinear response to external loading,
yielding and plasticity can be viewed as stress and thermal activation
processes on a complex energy landscape Cao  _et al._ (2019). The ductility of
yield-stress colloidal materials, such as cement Masoero  _et al._ (2012);
Ioannidou  _et al._ (2016), is a property that determines the integrity of
human-made structures, such as buildings and bridges, and thus is relevant to
their life-time and durability. At the same time, structural features, such as
pore network connectivity and topology, have fundamental implications for
multiple degradation mechanisms of heterogeneous materials Zhou  _et al._
(2019, 2020); Monfared  _et al._ (2020). It remains to be understood how the
internal energy sources in active matter interplay with the underlying energy
landscape of the host material and consequently modulate its structural and
mechanical properties.

Recent experiments on dilute gels Szakasits  _et al._ (2019); Saud  _et al._
(2021) observed a dramatic threefold change in the yield stress due to a tiny
fraction (0.1 %) of active particles disrupting the gel structure. The
activity-induced structural changes were long-lasting even after activity was
switched off Wei  _et al._ (2023), demonstrating a strong impact of activity
on the aging behavior of dilute gels. A wide range of practical colloidal
materials is dense, glassy and much stiffer, such as cement paste. How does
activity alter the mechanical behaviors of such dense gels?

In this work, we consider a simple system of passive colloidal particles with
a small fraction ($\sim 10\%$) of embedded ABPs. The ABPs self-propel at a
fixed swim force. We develop an active Brownian dynamics (BD) model to
numerically demonstrate precise control of the structure and mechanical
properties of dense disordered colloidal gels by altering the activity of the
embedded active particles. Specifically, we show a fine control of the gel’s
brittle-ductile transition by varying the swim force, accompanied by a
nonlinear change of the gel pore morphology and elastic properties. Our
results demonstrate a protocol for modifying colloidal gels in response to
external loading so that they deform in a ductile manner instead of
catastrophic brittle failure. These results pave the way to adaptive materials
design by active doping.

## Model of deeply annealed gels doped with active particles

We model dense active-doped gels as random polydisperse packings of spheres,
where only $\phi_{A}=$10% of the particles are Active Brownian Particles
(ABPs), as shown in Fig. 1(A). ABPs’ motion is overdamped and a particle moves
at a velocity $\bm{U}$ that satisfies the force balance

$-\zeta\bm{U}+\zeta\bm{F}^{S}+\bm{F}^{B}+\bm{F}^{P}=0,$ (1)

where $\zeta$ is the fluid drag, and $\bm{F}^{S}=\zeta U_{0}\bm{q}$ is the
self-propulsive swim force, with ${U}_{0}$ the swim speed in the direction
$\bm{q}$. $\bm{F}^{P}$ is the interparticle force and $\bm{F}^{B}$ is the
random Brownian force. An ABP is also subject to a random torque that changes
the orientation $\bm{q}$ over the reorientation timescale $\tau_{R}$ (see
Materials and Methods). The passive colloidal particles follow the same
dynamics except that they do not self-propel ($\bm{F}^{S}=0$.).

Activity of the ABPs will be switched on and then off, similar to recent
experiments for quasi-2D gels Wei  _et al._ (2023). The passive gel particles
strongly attract each other (attraction potential well depth
$\epsilon=10~{}k_{B}T$; see Materials and Methods), but the ABPs are not
attracted to the gel particles (or each other) and behave as soft spheres.
Before switching on the activity of the ABPs (their swim force is
$\bm{F}^{S}=0$ at this moment), the active-doped gel samples are annealed to
low effective temperature via the swap Monte Carlo algorithm Grigera and
Parisi (2001); Ninarello  _et al._ (2017). As a result, we obtain stable dense
gels that are brittle in nature Ozawa  _et al._ (2018). After annealing to a
low temperature ($\frac{\epsilon}{k_{B}T}\gg 1$), the interaction energy
dominates over mixing entropy and favors micro-phase separation where the ABPs
aggregate to form pockets naturally, as shown in Fig. 1.

Figure 1: Dense active-doped gel model. (A) A typical initially annealed
configuration. Green for polydisperse passive gel particles, darker particles
have larger coordination numbers. Orange for ABPs. Yellow arrows denote the
randomized initial orientation of the ABPs. The fraction of ABPs in each
sample is about $10\%$. (B) $g(r)$ for the gel-gel particles and ABP-ABP
particles show a micro-phase separation of ABPs.

## Protocol for activated gels: nonlinear dynamics of constant pressure
activation

The activation process proceeds by turning on the active swim force of the
ABPs to a finite value $F^{S}$ instantaneously, with randomized initial
orientations $\bm{q}$ and reorientation time $\tau_{R}=2~{}\tau_{B}$, where
$\tau_{B}$ is the Brownian timescale (see Methods). Activation causes the gel
to expand, and, as in an experiment, the gel expands at constant pressure. The
pressure, or the active stress of the gel is the sum of the kinetic pressure
$nk_{B}T$ and the virial part
$\Pi^{a}=\frac{1}{3V}\sum_{\alpha}\bm{x}_{\alpha}\cdot\bm{F}_{\alpha}^{t}$.
Here, $\bm{F}_{\alpha}^{t}$ is the total force exerted on particle $\alpha$
from the other passive and active particles, including the swim force, for
which we use the impulse formula:
$\bm{x}_{\alpha}\cdot\bm{F}_{\alpha}^{S}=\tau_{R}\bm{U}_{\alpha}\cdot\bm{F}^{S}_{\alpha}$.

Once the ABPs are activated, the unit cell expands isotropically according to
the non-equilibrium constant pressure barostat Wang and Brady (2015)

$\frac{\dot{L}}{L}=\frac{1}{\kappa_{T}}\left(nk_{B}T+\frac{1}{3}\left<\bm{x}\cdot\bm{F}^{t}\right>-\Pi_{0}\right),$
(2)

where $\kappa_{T}$ is an isothermal bulk viscosity of the gel and $\Pi_{0}$ is
the ambient pressure set to balance the thermal pressure before activating the
ABPs. Immediately after the activity is turned on, the active stress is large
and positive (Fig. 2A inset), leading to an instantaneous volumetric expansion
(Fig. 2A) that is far from equilibrium in nature. The excess stress quickly
decays due to the volume expansion, and the entire system then expands quasi-
statically over timescales longer than the Brownian time, as shown in Fig.
2(B).

To make a rough estimate of the effect of activation, a homogeneous gel will
rupture at a pressure
$P_{C}\sim\frac{N}{V}Z_{\text{eff}}^{\prime}F_{\text{mean}}a$, where the
typical attractive force between gel particles $F_{\text{mean}}\sim
10k_{B}T/a$, originating from their interaction potential. Here,
$Z_{\text{eff}}^{\prime}\sim O(1)$ is the effective per particle nearest
neighbor pairs across a plane, and $a$ is the average particle size. The
homogeneous collisional pressure from the ABPs is
$\Pi_{c}\sim\phi_{A}F^{S}a\frac{N}{V}$, which gives the critical value to
rupture the gel through isotropic expansion: $F^{S}\sim
Z_{\text{eff}}^{\prime}F_{\text{mean}}/\phi\sim 10^{2}k_{B}T/a$. Hence, within
a reasonably mild activity range, the gel will not be torn apart by the ABPs.

Figure 2: Nonlinear dynamics during the activation processes. (A) Typical
volume expansions at the early stage immediately after the ABPs’ activity is
turned on. Inset shows the instantaneous active stress is large at $t=0$ but
quickly decays. (B) Quasi-static gel volume expansion during activation;
$F^{S}$ increases from 10 to 30 $k_{B}T/a$ as line colors become darker. (C)
Both the steady state and final gel volumes increase as $F^{S}$ increases. (D)
The self-diffusivity of the gel particles and ABPs during activation.

On the other hand, at a local scale the ABPs can break apart passive particle
pairs and penetrate into the gel. The critical swim force $F^{S}_{c}$ to break
a pair of attractive passive particles can be estimated by considering the
energy of quasi-statically breaking a bond:

$F^{S}_{c}\sim\frac{\epsilon}{a}Z_{\text{eff}},$ (3)

where $Z_{\text{eff}}$ is the effective coordination number of passive gel
particles. Simulation results give $F^{S}_{c}\sim 15k_{B}T/a$, or
$Z_{\text{eff}}\approx 1.5$, which is reasonable. This critical swim force
separates qualitatively different regimes of the gel structure and mechanics,
as discussed below.

After steady state, we turn off the activity of ABPs ($F^{S}=0$), and let the
system relax to the initial ambient pressure $\Pi_{0}$. The total volume
shrinks as the active collisional pressure is eliminated. The final volume is
always larger than the initial value and grows as swim force increases, as
shown in Fig. 2(C), indicating irreversible structural changes. We then
quantify the structural changes compared to the initial annealed state.

## Activity-modulated elasticity

To determine the effect of active doping on the mechanical properties of the
gel, we perform athermal quasi-static (AQS) simple shear simulations on the
final gel configurations. During shear, the ABPs remain deactivated
($F^{S}=0$). Below, we denote gel samples by the $F^{S}$ (in units of
$k_{B}T/a$) values of the ABPs during the activation processes.

First, we examine the linear response of gels activated at different $F^{S}$.
We quantify the shear modulus using the initial slope of the stress-strain
curve at shear strain $\varepsilon<1\%$. The initial gels are strong compared
to the previous experimental dilute gels Szakasits  _et al._ (2019); Saud  _et
al._ (2021) reflected by their relatively high shear modulus $G$. As shown by
the blue line in Fig. 3(A), $G(F^{S}=0)\approx 0.25~{}G_{\text{FCC}}$, where
$G_{\text{FCC}}$ is the modulus of a FCC crystal calculated from the same AQS
protocol. For these inactivated gels, their elastic regime is quite ideal,
evidenced by the moving-averaged frequency of stress drops as shear strain
increases, as shown by the orange line in Fig. 3(D). Very few stress drop
appear at the early stage of small $\varepsilon<1\%$.

After activation, the shear modulus $G$ is weakened as $F^{S}$ increases.
Noticeably, the decrease in $G$ is nonlinear in $F^{S}$, with a drastic
decline around $F^{S}\sim F^{S}_{c}$. Although $G$ decreases rapidly around
$F^{S}_{c}\sim 15~{}k_{B}T/a$, it remains finite and the same order of
magnitude even at $F^{S}=30~{}k_{B}T/a$ (from initial $G_{0}\sim
0.25~{}G_{\text{FCC}}$ to $G_{0}\sim 0.1~{}G_{\text{FCC}}$). Hence,
$F^{S}_{c}$ is below the force necessary to “fluidize” the gel, where the
shear modulus decreases to zero. According to the estimate in the last
section, ABPs with swim forces above $F^{S}_{c}$ can break local bonds, and
cages formed during gel annealing can be broken in this regime. As a result,
the elastic regime of highly activated gels deviates from the ideal scenario,
and becomes populated with microscopic plastic deformations even at small
strain $\varepsilon<1\%$. These localized irreversible deformations are
recorded in the higher initial frequency of stress drops, as shown by the blue
and green lines in Fig. 3(D); the system loses rigidity as $F^{S}$ increases.

The structural origin of the rigidity loss can be understood from the
perspective of the Cauchy-Born theory of amorphous solids Alexander (1998);
Zaccone  _et al._ (2009). To clarify this, we define the mean number of
contacts per particle $N_{c}$, as shown in Fig. 3(B) (see Materials and
Methods). For gels below $F^{S}_{c}$, $N_{c}$ stays large, which is consistent
with the rigid gel response to loading. For gels activated at
$F^{S}>F^{S}_{c}$, $N_{c}$ decreases rapidly, indicating the onset of loss of
isostaticity Maxwell (1864). As one expects from the Cauchy-Born theory, the
shape of $G(F^{S})$ resembles that of $N_{c}(F^{S})$, comparing Fig. 3(A) and
(B). The correlation between $N_{c}$ and $G$ observed here is consistent with
previous experiments on passive colloidal gels Whitaker  _et al._ (2019);
Hsiao  _et al._ (2012). Overall, our findings here show that tuning $F^{S}$
can control the gel network topology and elastic modulus.

Figure 3: Mechanical properties of the gels after activation-deactivation
processes at different swim forces. (A) Both elastic modulus ($G$, blue line,
normalized by $G_{\text{FCC}}$) and maximum stress drop size max($\sigma$)
decrease with swim force. (B) Mean contact number per particle $N_{c}$ as a
function of $F^{S}$. (C) Stress-strain curves for the initial gel and
activated gels during athermal-quasi-static shear tests, each curve is
averaged over 5 samples. (D) Typical stress drop frequencies (moving averaged)
during shear.

## Activity-induced brittle-to-ductile transition

Preluded by the increasing abundance of microscopic plastic deformations at
small strain, for large deformations that yield the gels, the stress-strain
responses undergo a continuous change from brittle to ductile behavior as
$F^{S}$ increases, as shown by Fig. 3(A) red line and (C). The yield point is
taken as when the macroscopic shear stress saturates or starts to decrease as
shear strain increases. The yield strain is always around $\varepsilon\sim
8\%$, regardless of different $F^{S}$ values. The yield stress, however,
decreases continuously as $F^{S}$ increases. To quantify the gel ductility, we
measure the stress drops $\Delta\sigma$ (see Materials and Methods) during the
simple shear processes. As expected, the maximum stress drop
$\Delta\sigma_{\text{max}}$ coincides with the yielding transition and the
magnitude of $\Delta\sigma_{\text{max}}$ is used as a scalar measure of
brittleness.

The inactivated gels are quite brittle, with a significant stress overshoot in
their stress-strain curves, as shown by the orange curve of Fig. 3(C).
Correspondingly, $\Delta\sigma_{\text{max}}$ is largest at $F^{S}=0$. For gels
activated at different $F^{S}$, $\Delta\sigma$ follow similar truncated power-
law distributions (SI Fig. 1), but $\Delta\sigma_{\text{max}}$ decreases as
$F^{S}$ increases, as seen in Fig. 3(A). In contrast to the drastic decrease
of $G$ around $F^{S}_{c}$, $\Delta\sigma_{\text{max}}$ decreases mildly with
increasing $F^{S}$. Meanwhile, the total number of stress drops, or plastic
events, increases when $F^{S}>F^{S}_{c}$, as shown by the total areas under
different curves in Fig. 3(D). Overall, these statistics of plastic events
indicate that the gels become more ductile as $F^{S}$ increases.

We notice that the $\Delta\sigma_{\text{max}}$ at $F^{S}=30~{}k_{B}T/a$ is
still finite ($\sim O(1)$), raising the question of whether the yielding
behavior can be classified as truly ductile or rather a series of small yet
brittle avalanches. Current continuum models give different predictions on the
nature (discontinuous or continuous) of the brittle-ductile transition Ozawa
_et al._ (2018); Singh  _et al._ (2020); Pollard and Fielding (2022); Rossi
_et al._ (2022). To further clarify our quantification of gel ductility, we
examine the spatial distributions of shear transformations during shear. For
all gels activated at $F^{S}>F^{S}_{c}$, there is no shear banding that
percolates the entire sample; in contrast, such percolating shear bands always
show up in the deeply annealed gels without activation (Fig. 4A-D). The
percolation behavior of shear transformation zones is consistent with the
observed trends of $\Delta\sigma_{\text{max}}$ and the avalanche frequency
$\nu$ (see Materials and Methods) with respect to $F^{S}$. The absence of
shear bands suggests that the gels activated above $F^{S}_{c}$ undergo truly
ductile yielding. Below $F^{S}_{c}$, the characteristic behaviors of
$\Delta\sigma_{\text{max}}$ and $\nu$ indicate that the gels continuously
enhance ductility as $F^{S}$ increases.

Figure 4: Typical shear transformations for $F^{S}=20$ during stress drop
events at shear strains of (A) $\varepsilon=5.47\%$, (B) $\varepsilon=9.90\%$,
(C) $\varepsilon=10.16\%$ and (D) shear banding for $F^{S}=0$ at
$\varepsilon=10.15\%$. No shear banding is observed for $F^{S}=20$. Yellow
particles have large non-affine displacements (¿1 particle diameter) during
the stress drop event. (E) Nonaffine displacements Falk and Langer (1998);
Shimizu  _et al._ (2007) accumulated during the shear strain range
$\varepsilon\in[0,10\%]$ (before yielding) for groups of particles with
different excess entropy.

## Activity-modulated structural entropy change and its connection to brittle-
to-ductile transition

To understand the underlying structural changes responsible for the enhanced
ductility, we examine the degree of disorder of the passive gel particle
configurations by calculating the local excess entropy
$S_{\text{ex}}=S-S_{\text{ig}}$ Rosenfeld (1977); Dyre (2018) (see Materials
and Methods) around each particle and average over the sample, where
$S_{\text{ig}}$ is the ideal gas entropy. Below $F^{S}_{c}$, the gel structure
becomes more disorder as $F^{S}$ increases, as seen by the increase of
$S_{\text{ex}}$ with $F^{S}$ shown in Fig. 5(A). The enhancement in ductility
accompanied by the increase of disorder is consistent with recent observations
in other disordered materials, such as silica or metallic glasses Tang  _et
al._ (2021). In that study, structural disorder was assumed to originate from
irradiation. In fact, ABPs share similarities with irradiation bombardment in
terms of their ballistic motions at short times. $S_{\text{ex}}$ saturates
when $F^{S}$ is slightly above $F^{S}_{c}$ and then decreases with increasing
$F^{S}$, indicating that the local arrangement of gel particles is more
ordered, even though the gel is more ductile.

Figure 5: Local degree of structural disorder measured by the excess entropy
$S_{\text{ex}}$. (A) the excess entropy (left y-axis) averaged over all
particles ($S_{all}$ blue line) and gel particles ($S_{gel}$ black line), and
the final total volume $V$. The initial total volume is $V_{0}=1.03\times
10^{4}~{}a^{3}$. (B) Probability distribution of $S_{\text{ex}}$. (C)(D) show
snapshots of the final gels activated at $F^{S}=10,30$, respectively, with gel
particles colored by their local excess entropy and the ABPs colored in
orange.  Figure 6: Structural changes of the passive gel induced by the ABPs.
(A)(B)(C) compare the final surface morphology of the passive gels activated
at $F^{S}=10$, $F^{S}=20$, and $F^{S}=30$. ABPs are orange spheres with yellow
arrows indicating their orientation, gel particles not shown. (D) shows the
total number of pores in the sample $N_{\text{pore}}$ (blue line, left y-axis)
and the total volume of pore space $V_{\text{pore}}$ (black line, right
y-axis), as a function of the swim force of ABPs.

## Activity-modified pore networks enhance dynamical heterogeneity

To rationalize the reversed trend of $S_{\text{ex}}$ at large $F^{S}$, we
notice that above $F^{S}_{c}$, the gel pore volume increases significantly
with $F^{S}$ (Fig. 6D). The irreversible volume expansion induced by activity
increases the internal pore space volume $V_{\text{pore}}$, as shown in Fig.
6(D). Below the critical swim force, both $V_{\text{pore}}$ and number of
pores $N_{\text{pore}}$ increase mildly with $F^{S}$. In this regime, the ABPs
are not able to penetrate into the gel and disrupt its structure, but rather
they locally strain the gel and create a multitude of small pores. Here,
$S_{\text{ex}}$ remains a good indicator for local degree of order in the
particle configuration and gradually increases.

However, above $F^{S}_{c}$ the ABPs penetrate through the gel backbone,
enhancing the pore connectivity. As a result, $N_{\text{pore}}$ decreases
quickly, while $V_{\text{pore}}$ increases rapidly, as shown in Fig. 6(C).
Typical final gel pore structures are shown in Fig. 6 (A-C) at swim forces
below, near and above the critical force $F^{S}_{c}$. Hence, the gel develops
a larger pore network and exhibits stronger structural heterogeneities, as
supported by the broadened distribution of $S_{\text{ex}}$ in Fig. 5(B). In
this regime, the decrease of $S_{\text{ex}}$ emphasizes only small nucleations
of crystalline ordered clusters, but does not reflect the larger scale pore
topology. A closer examination of the most orderly arranged particles with
$S_{\text{ex}}<{\bar{S}_{\text{ex}}}-\sigma_{S}$, where $\bar{S}_{\text{ex}}$
is the mean of $S_{\text{ex}}$ and $\sigma_{S}$ is the standard deviation,
shows that these particles do not form a percolating backbone to sustain the
globally applied external load; hence, they do not dominate the gel’s response
to shear. Instead, the enhanced overall porosity and heterogeneity now take
over and increase gel ductility. Correspondingly, the system behaves
dynamically more heterogeneous, as shown by the multitude of localized shear
transformations with small stress drops before yielding (Fig. 4).

Figure 7: Active-doped monodisperse gels. (A) Local excess entropy
$S_{\text{ex}}$ averaged over all particles (blue line, left y-axis) and
passive particles only (black line, left y-axis), and the final system volume
$V$ (red line, right y-axis). (B) Averaged stress-strain curves during
athermal quasi-static shear. (C) Elastic moduli $G$ (blue line, left y-axis,
normalized by $G_{\text{FCC}}$) and the maximum stress drop size
$\Delta\sigma_{\text{max}}$ (red line, right y-axis). Typical snapshot for
final gel configurations at (D) $F^{S}=10$ and (E) $F^{S}=30$, where ABPs are
colored orange, passive particles with local FCC order colored blue, HPC order
colored yellow. Other passive particles, i.e. disordered ones, are not
displayed.

## Activity drives the energy landscape

From the perspective of inherent structures in glassy states, the changes in
mechanical properties discussed above should be reflected in the different
local environments on the energy landscape. Indeed, during activation, the gel
particles significantly enhance their mobility as $F^{S}$ increases, as shown
by the self-diffusivities in Fig. 2(D). Their steady-state mean-squared-
displacements (MSD) exceed a particle diameter when $F^{S}>F^{S}_{c}$ (SI
Fig.1), indicating cage-breaking events during activation, consistent with the
bond-breaking argument and the observed change in inter-particle contacts.
Hence, we argue that the embedded activity drives the system to effectively
tunnel through energy barriers and explore the landscape.

## Effects of active doping in monodisperse gels

To contrast the above results obtained from strong glass formers made of
polydisperse particles, we perform the same numerical protocols to 3D gels
made of monodisperse particles, again with 10% ABPs. After thermalization at
$1~{}k_{B}T$, they form disordered gels similar to the polydisperse case.

A similar critical value of swim force $F^{S}_{c}\approx 15$ is observed,
where the system volume increases mildly below $F^{S}_{c}$ and drastically
above $F^{S}_{c}$ (Fig. 7A). However, the effects of activity on the structure
and mechanics show qualitative differences from the polydisperse gels. The
averaged local excess entropy $S_{\text{ex}}$ stays at low values below
$F^{S}_{c}$ and increases quickly for $F^{S}>F^{S}_{c}$. The increase of
$S_{\text{ex}}$ merely reflects the decrease in packing density — locally
particles form clearly visible crystallites, as shown in Fig. 7(D,E), and the
fraction of crystalline particles increases with $F^{S}$. Our result is
consistent with previous studies showing the activity-induced enhancement of
crystalline order in monodisperse packings Mallory and Cacciuto (2019);
Ramananarivo  _et al._ (2019), distinguishing the role of embedded intrinsic
activity from thermal activation.

Intriguingly, the change in mechanics of the monodisperse gels shows a non-
monotonic trend: Below $F^{S}_{c}$, the activated gels become tougher, with
higher yield stress, as shown in Fig. 7(B), and enhanced elastic modulus $G$
(Fig. 7(C) blue line, normalized by $G_{\text{FCC}}$). The ductility of the
gel is marginally reduced (Fig. 7(C) red line). As $F^{S}$ increases further
above $F^{S}_{c}$, the gel enters a weakening regime, where both the elastic
modulus and yield stress decrease. Meanwhile, it becomes more ductile with
smaller yield stress overshoot and $\Delta\sigma_{\text{max}}$. Similar to the
polydisperse gels, in this regime, the monodisperse gel’s mechanical response
to external loading develops multiscale features and enhanced dynamical
heterogeneity since the gel expands and significantly changes its pore
structure. Local clusters of crystallites are no longer good predictors for
the global stress relaxation mechanism, while the globally averaged degree of
disorder $S_{\text{ex}}$ is consistent with the enhanced ductility.

We attribute the different trends in mechanical behaviors between activated
monodisperse and polydisperse gels to their underlying energy landscapes and
different glass-formation propensity. The polydisperse gels are prepared at a
very stable, deeply annealed state. Under the effect of weak embedded
activity, statistically the system always finds pathways to less stable glassy
states near the current local minimum. In the monodisperse gels, nucleation of
small crystalline clusters serves as low energy barrier pathways towards the
global minimum of a crystal. Hence, even when embedded ABPs perturb the system
weakly, it is both thermodynamically and kinetically favorable to first
enhance local crystalline order, before a stronger swim force disrupts the gel
structure and expands the total volume significantly. As a result, below
$F^{S}_{c}$, embedded ABPs reduce the polydisperse gel toughness by increasing
disorder, but enhance the monodisperse gel toughness by slightly increasing
local crystalline order. Above $F^{S}_{c}$, the gel structure is disrupted
with many bonds broken, and the effect of increasing pore volume dominates to
continuously weaken both polydisperse and monodisperse gels, regardless of the
formation of locally ordered clusters.

## Conclusion and Discussions

In this work we demonstrate the possibility to engineer the mechanical
response and mode of failure of colloidal gels by doping them with active
swimmers of different activity strength. With a prototypical computer model of
polydisperese colloidal gels doped with Active Brownian Particles, we show
that activity dramatically alters the structure and mechanical response of
doped colloidal materials. By controlling the particle activity protocol,
i.e., on/off timescales, along with the activity strength, $F^{S}$, we can
enhance the degree of irreversible change and the dynamical heterogeneities
present in the mechanical response and dynamics of the passive colloidal gels.
Through theoretical arguments, we predict that there is a critical swim force
that an active particle needs to break the bonds between the constituents of
the gel. Above this critical swim force, active swimmers enhance the
mesoscopic degree of disorder by developing larger pore networks, and thus
lead the structure to undergo a brittle-to-ductile transition. This scenario
is revealed by our simulations, which show the existence of such a critical
force and that the activity-enhanced porous network softens the colloidal
structure. Such behavior is reflected by the abrupt decrease of the observed
elastic shear modulus and change in degree of disorder, quantified by excess
local entropy and pore connectivity measures.

Recent numerical investigations Tang  _et al._ (2021) revealed connections
between the structural disorder of amorphous materials and their energy
landscape via the statistics of local energy barriers. In this sense, the
embedded activity allows passive particles to “tunnel” through the energy
barriers. Above the critical swim force $F^{S}_{c}$, gels become more porous,
with enhanced connectivity. Porosity and pore connectivity strongly influence
the transport properties of porous media Zhou  _et al._ (2019); Pinson  _et
al._ (2018), and our results suggest that active doping can be an effective
technique to shape the pore networks in disordered materials.

This work is motivated, in part, by the use of dense colloidal materials like
cement Masoero  _et al._ (2012); Ioannidou  _et al._ (2016); Zhou  _et al._
(2019), where recent experiments have explored the possibility of embedding
microbes that can endure the salinity and basic environment in cement paste
Simkiss and Wilbur (2012); Nodehi  _et al._ (2022). Typical manufacturing
methods for cementitious materials (“curing”) involve complex hydration and
aging (“hardening”) processes Bullard  _et al._ (2011) over the timescale of
days. We speculate that our results can be applicable to such strong
materials, yet the timing and duration to stimulate the embedded active
constituents needs further investigation. The effects are different for
polydisperse and monodisperse disordered gels, reflecting their different
underlying energy landscapes. Our investigation here also provides guidelines
for designing reconfigurable materials. For example, we can identify optimal
$F^{S}$, near the critical force $F^{S}_{c}$ to maximally enhance the gel
ductility in real-time with a minimal weakening of the elastic modulus. In
experiments, such tunability can be realized via changing the light intensity,
or concentration of fuel chemicals for the active particles. A possible
variation of the protocol is to maintain activity during the shearing process,
which may be a useful tool to control material response in situ. We hope these
results motivate further work, both in continuum theory on the ductility of
heterogeneous materials and in experiments of active doped gels in the dense
regime.

## Appendix

### Simulation details

The dynamics of the ABPs and gel particles are integrated according to the
over-damped Langevin equation

$\begin{split}-\zeta\bm{U}+\zeta\bm{q}U_{0}+\bm{F}^{B}+\bm{F}^{P}&=0,\\\
\left<\bm{F}^{B}\right>&=0,\\\
\left<\bm{F}^{B}(0)\bm{F}^{B}(t)\right>&=2k_{B}T\zeta\delta(t)\bm{I},\\\
\frac{d\bm{q}}{dt}&=\bm{\Omega}^{B}\times\bm{q},\\\
\left<\bm{\Omega}^{B}\right>&=0,\\\
\left<\bm{\Omega}^{B}(0)\bm{\Omega}^{B}(t)\right>&=2\delta(t)\bm{I}/\tau_{R},\\\
\end{split}$ (4)

where $\zeta$ is the translational drag coefficient, and $\bm{q}$ is a unit
vector denoting the orientation of the ABP. $U_{0}=F^{S}/\zeta$ is the free
space swim speed, which is 0 for passive gel particles. $\bm{U}$ is the
particle velocity. $\bm{F}^{B}$ is the Brownian force with zero mean and
distribution satisfying the fluctuation-dissipation theorem. $\bm{F}^{P}$ is
the force arising from particle-particle interactions. $\bm{\Omega}^{B}$ is
the random angular velocity resulting from random reorientations and we set
$\tau_{R}=2~{}\tau_{B}$ as the persistent runtime. The simulations are in
reduced units: energy is normalized by $k_{B}T$, length by average particle
diameter $a$, and timescale by the Brownian time $\tau_{B}=\zeta
a^{2}/k_{B}T$.

The interaction potential between particle $\alpha$ and $\beta$ is

$V_{\alpha\beta}(r)=C\epsilon\left[\left(\frac{a_{\alpha\beta}}{r}\right)^{n}-\left(\frac{a_{\alpha\beta}}{r}\right)^{m}\right],$
(5)

where the potential depth is $\epsilon=10~{}k_{B}T$, and $n=20$, $m=10$,
$C=\frac{n}{n-m}(\frac{n}{m})^{m/(n-m)}=4$. The effective size is the average
of the two interacting particle diameters
$a_{\alpha\beta}=(a_{\alpha}+a_{\beta})/2$. For passive-passive pairs, the
interaction is cut off at $r_{p}=2.5\sigma_{\alpha\beta}$. While for either
$\alpha$ or $\beta$ being an ABP, the interaction is cut off at
$r_{a}=2^{1/10}a_{\alpha\beta}$ and shifted to only retain the repulsive part.
The gel particle diameters are drawn from a uniform distribution $U(0.8,1.2)$.
All ABPs have diameter $a_{A}=1.0$. Each simulation sample has on average
11000 particles.

### Structural characterizations

The pore surface is determined by the $\alpha$-shape algorithm Stukowski
(2014), with a probe sphere radius 0.9, close to the nearest neighbor distance
between gel particles. The local excess entropy for each particle $\alpha$ is
estimated by the two-body contribution Dyre (2018)

$S^{\alpha}_{ex}\approx
S^{\alpha}_{2}=-2\pi\rho_{\alpha}\int_{0}^{R}\left[g_{\alpha}(r)\ln
g_{\alpha}(r)-g_{\alpha}(r)+1\right]r^{2}dr,$ (6)

where $\rho_{\alpha}$ is the local particle number density averaged over a
sphere of radius $R=5$ around particle $\alpha$. The integration is cut off at
$R=5$, where the $g_{\alpha}(r)$ converges to 1. Including ABPs for computing
$S_{\text{ex}}$ shows similar trends, and we focus on the gel $S_{\text{ex}}$
because in all the gels (activated or not), the passive particles’
rearrangements dominate the stress drops as they constitute the load-bearing
backbone of the gels. The local crystalline order for the monodisperse gels is
identified by the Ackland-Jones bond angle method Ackland and Jones (2006).

The mean number of contacts per particle $N_{c}$ is defined by including all
repulsive pairs of particles, as well as attractive pairs whose separation is
$r_{12}<a_{12}+\delta r$. We show the result for $\delta r=0.01$ in Fig. 3B.
The numeric value of $N_{c}$ depends on $\delta r$, but the trend with
increasing $F^{S}$ is similar, insensitive to the choice of $\delta r$.

### Characterization of stress avalanches

A stress drop event is defined as a continuous interval of shear strain where
the shear stress monotonically decreases. The difference of stress values at
the beginning and the end of this stress drop event, $\Delta\sigma$, is
denoted as the stress drop magnitude of this event. We then define the
frequency of such stress drop events as the number of events per unit amount
of incremental strain. For visual clarity, the $\nu(\varepsilon)$ curves in
Fig. 3(D) are denoised by moving average with window size
$d\varepsilon=0.5\%$.

##  I Acknowledgment

We acknowledge financial support from the Department of Energy under Grant DE-
SC0022966. The simulations in this work are performed on the CFN computing
facility at the Brookhaven National Laboratory (proposal ID: GUP-310204). T.Z.
thanks S. Chen, H. Row, S. Mallory, D. Frag, X. Wan, and X. Cheng for their
helpful discussions. All study data are included in the article.

## References

  * Lo  _et al._ (2000) C.-M. Lo, H.-B. Wang, M. Dembo, and Y.-l. Wang, Biophysical journal 79, 144 (2000).
  * Doss  _et al._ (2020) B. L. Doss, M. Pan, M. Gupta, G. Grenci, R.-M. Mège, C. T. Lim, M. P. Sheetz, R. Voituriez, and B. Ladoux, Proceedings of the National Academy of Sciences 117, 12817 (2020).
  * Fily and Marchetti (2012) Y. Fily and M. C. Marchetti, Physical review letters 108, 235702 (2012).
  * Takatori  _et al._ (2014) S. C. Takatori, W. Yan, and J. F. Brady, Physical review letters 113, 028103 (2014).
  * Omar  _et al._ (2020) A. K. Omar, Z.-G. Wang, and J. F. Brady, Physical Review E 101, 012604 (2020).
  * Angelani  _et al._ (2009) L. Angelani, R. Di Leonardo, and G. Ruocco, Physical review letters 102, 048104 (2009).
  * Reichhardt and Reichhardt (2017) C. O. Reichhardt and C. Reichhardt, Annual Review of Condensed Matter Physics 8, 51 (2017).
  * Ro  _et al._ (2022) S. Ro, B. Guo, A. Shih, T. V. Phan, R. H. Austin, D. Levine, P. M. Chaikin, and S. Martiniani, Physical Review Letters 129, 220601 (2022).
  * Damasceno  _et al._ (2012) P. F. Damasceno, M. Engel, and S. C. Glotzer, Science 337, 453 (2012).
  * Wang  _et al._ (2012) Y. Wang, Y. Wang, D. R. Breed, V. N. Manoharan, L. Feng, A. D. Hollingsworth, M. Weck, and D. J. Pine, Nature 491, 51 (2012).
  * Vlasov  _et al._ (2001) Y. A. Vlasov, X.-Z. Bo, J. C. Sturm, and D. J. Norris, Nature 414, 289 (2001).
  * Liu  _et al._ (2014) Q. Liu, M. G. Campbell, J. S. Evans, and I. I. Smalyukh, Advanced Materials 26, 7178 (2014).
  * Engel  _et al._ (2015) M. Engel, P. F. Damasceno, C. L. Phillips, and S. C. Glotzer, Nature materials 14, 109 (2015).
  * Wang  _et al._ (2017) Y. Wang, I. C. Jenkins, J. T. McGinley, T. Sinno, and J. C. Crocker, Nature communications 8, 14173 (2017).
  * He  _et al._ (2020) M. He, J. P. Gales, É. Ducrot, Z. Gong, G.-R. Yi, S. Sacanna, and D. J. Pine, Nature 585, 524 (2020).
  * Ni  _et al._ (2014) R. Ni, M. A. C. Stuart, M. Dijkstra, and P. G. Bolhuis, Soft Matter 10, 6609 (2014).
  * Mallory  _et al._ (2020) S. Mallory, M. Bowers, and A. Cacciuto, The Journal of Chemical Physics 153, 084901 (2020).
  * Ramananarivo  _et al._ (2019) S. Ramananarivo, E. Ducrot, and J. Palacci, Nature communications 10, 3380 (2019).
  * Saud and Solomon (2023) K. T. Saud and M. J. Solomon, Journal of Colloid and Interface Science (2023).
  * Omar  _et al._ (2018) A. K. Omar, Y. Wu, Z.-G. Wang, and J. F. Brady, ACS nano 13, 560 (2018).
  * Sastry (2001) S. Sastry, Nature 409, 164 (2001).
  * Sastry  _et al._ (1998) S. Sastry, P. G. Debenedetti, and F. H. Stillinger, Nature 393, 554 (1998).
  * Heuer (2008) A. Heuer, Journal of Physics: Condensed Matter 20, 373101 (2008).
  * Berthier (2011) L. Berthier, arXiv preprint arXiv:1106.1739 (2011).
  * Berthier  _et al._ (2011) L. Berthier, G. Biroli, J.-P. Bouchaud, and R. L. Jack, Dynamical heterogeneities in glasses, colloids, and granular media 150, 68 (2011).
  * Wittmer  _et al._ (2013) J. Wittmer, H. Xu, P. Polińska, F. Weysser, and J. Baschnagel, The Journal of chemical physics 138, 12A533 (2013).
  * Yoshino and Zamponi (2014) H. Yoshino and F. Zamponi, Physical Review E 90, 022302 (2014).
  * Goodrich  _et al._ (2016) C. P. Goodrich, A. J. Liu, and J. P. Sethna, Proceedings of the National Academy of Sciences 113, 9745 (2016).
  * Cao  _et al._ (2019) P. Cao, M. P. Short, and S. Yip, Proceedings of the National Academy of Sciences 116, 18790 (2019).
  * Masoero  _et al._ (2012) E. Masoero, E. Del Gado, R.-M. Pellenq, F.-J. Ulm, and S. Yip, Physical review letters 109, 155503 (2012).
  * Ioannidou  _et al._ (2016) K. Ioannidou, K. J. Krakowiak, M. Bauchy, C. G. Hoover, E. Masoero, S. Yip, F.-J. Ulm, P. Levitz, R. J.-M. Pellenq, and E. Del Gado, Proceedings of the National Academy of Sciences 113, 2029 (2016).
  * Zhou  _et al._ (2019) T. Zhou, K. Ioannidou, F.-J. Ulm, M. Z. Bazant, and R.-M. Pellenq, Proceedings of the National Academy of Sciences 116, 10652 (2019).
  * Zhou  _et al._ (2020) T. Zhou, M. Mirzadeh, R. J.-M. Pellenq, and M. Z. Bazant, Physical Review Fluids 5, 124201 (2020).
  * Monfared  _et al._ (2020) S. Monfared, T. Zhou, J. E. Andrade, K. Ioannidou, F. Radjaï, F.-J. Ulm, and R. J.-M. Pellenq, Physical Review Letters 125, 255501 (2020).
  * Szakasits  _et al._ (2019) M. E. Szakasits, K. T. Saud, X. Mao, and M. J. Solomon, Soft Matter 15, 8012 (2019).
  * Saud  _et al._ (2021) K. T. Saud, M. Ganesan, and M. J. Solomon, Journal of Rheology 65, 225 (2021).
  * Wei  _et al._ (2023) M. Wei, M. B. Zion, and O. Dauchot, arXiv preprint arXiv:2302.08360 (2023).
  * Grigera and Parisi (2001) T. S. Grigera and G. Parisi, Physical Review E 63, 045102 (2001).
  * Ninarello  _et al._ (2017) A. Ninarello, L. Berthier, and D. Coslovich, Physical Review X 7, 021039 (2017).
  * Ozawa  _et al._ (2018) M. Ozawa, L. Berthier, G. Biroli, A. Rosso, and G. Tarjus, Proceedings of the National Academy of Sciences 115, 6656 (2018).
  * Wang and Brady (2015) M. Wang and J. F. Brady, Physical review letters 115, 158301 (2015).
  * Alexander (1998) S. Alexander, Physics reports 296, 65 (1998).
  * Zaccone  _et al._ (2009) A. Zaccone, H. Wu, and E. Del Gado, Physical review letters 103, 208301 (2009).
  * Maxwell (1864) J. C. Maxwell, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 27, 294 (1864).
  * Whitaker  _et al._ (2019) K. A. Whitaker, Z. Varga, L. C. Hsiao, M. J. Solomon, J. W. Swan, and E. M. Furst, Nature communications 10, 2237 (2019).
  * Hsiao  _et al._ (2012) L. C. Hsiao, R. S. Newman, S. C. Glotzer, and M. J. Solomon, Proceedings of the National Academy of Sciences 109, 16029 (2012).
  * Singh  _et al._ (2020) M. Singh, M. Ozawa, and L. Berthier, Physical Review Materials 4, 025603 (2020).
  * Pollard and Fielding (2022) J. Pollard and S. M. Fielding, Physical Review Research 4, 043037 (2022).
  * Rossi  _et al._ (2022) S. Rossi, G. Biroli, M. Ozawa, G. Tarjus, and F. Zamponi, Physical Review Letters 129, 228002 (2022).
  * Falk and Langer (1998) M. L. Falk and J. S. Langer, Physical Review E 57, 7192 (1998).
  * Shimizu  _et al._ (2007) F. Shimizu, S. Ogata, and J. Li, Materials transactions 48, 2923 (2007).
  * Rosenfeld (1977) Y. Rosenfeld, Physical Review A 15, 2545 (1977).
  * Dyre (2018) J. C. Dyre, The Journal of chemical physics 149, 210901 (2018).
  * Tang  _et al._ (2021) L. Tang, H. Liu, G. Ma, T. Du, N. Mousseau, W. Zhou, and M. Bauchy, Materials Horizons 8, 1242 (2021).
  * Mallory and Cacciuto (2019) S. A. Mallory and A. Cacciuto, Journal of the American Chemical Society 141, 2500 (2019).
  * Pinson  _et al._ (2018) M. B. Pinson, T. Zhou, H. M. Jennings, and M. Z. Bazant, Journal of colloid and interface science 532, 118 (2018).
  * Simkiss and Wilbur (2012) K. Simkiss and K. M. Wilbur,  _Biomineralization_ (Elsevier, 2012).
  * Nodehi  _et al._ (2022) M. Nodehi, T. Ozbakkaloglu, and A. Gholampour, Journal of Building Engineering 49, 104038 (2022).
  * Bullard  _et al._ (2011) J. W. Bullard, H. M. Jennings, R. A. Livingston, A. Nonat, G. W. Scherer, J. S. Schweitzer, K. L. Scrivener, and J. J. Thomas, Cement and concrete research 41, 1208 (2011).
  * Stukowski (2014) A. Stukowski, Jom 66, 399 (2014).
  * Ackland and Jones (2006) G. Ackland and A. Jones, Physical Review B 73, 054104 (2006).


\times\left.\left(\frac{1}{\rho}\frac{\partial\rho}{\partial
t}\right)\left|\frac{1}{\rho}\frac{\partial\rho}{\partial t}\right|\right]\;,$
(A17)

where

$\frac{1}{\rho}\frac{\partial\rho}{\partial
t}=\frac{2\left(\rho_{j-\frac{1}{2}}^{n+1}-\rho_{j-\frac{1}{2}}^{n}\right)}{\Delta
t^{n+\frac{1}{2}}\left(\rho_{j-\frac{1}{2}}^{n+1}+\rho_{j-\frac{1}{2}}^{n}\right)}\;,$
(A18)

and

$\Delta t^{n+\frac{1}{2}}=t^{n+1}-t^{n}\;.$ (A19)

The time step is chosen to ensure numerical stability. The criterion for
linear stability of the difference equations used in the WONDY code (Hicks,
1978) is

$\begin{multlined}\Delta t<\Delta
x\left[B_{2}c+2B_{1}^{2}\left|\frac{\dot{\rho}}{\rho}\right|\Delta x\right.\\\
\left.+\sqrt{\left(B_{2}c++2B_{1}^{2}\left|\frac{\dot{\rho}}{\rho}\right|\Delta
x\right)^{2}+c^{2}}\,\right]^{-1}\;.\end{multlined}\Delta t<\Delta
x\left[B_{2}c+2B_{1}^{2}\left|\frac{\dot{\rho}}{\rho}\right|\Delta x\right.\\\
\left.+\sqrt{\left(B_{2}c++2B_{1}^{2}\left|\frac{\dot{\rho}}{\rho}\right|\Delta
x\right)^{2}+c^{2}}\,\right]^{-1}\;.$ (A20)

At the end of each cycle, the value of the right hand side of Equation A20 is
evaluated for each cell to determine the maximum time step that can be used to
ensure stability. Since the criteria in Equation A20 is based on a linear
stability analysis, and then assumed to apply generally, the critical value
for each cell is scaled by a constant, $K_{\rm t_{1}}$, to account for any
non-linear effects and ensure stability, such that

$\begin{multlined}\Delta t^{n+\frac{3}{2}}_{j-\frac{1}{2}}=K_{\rm t_{1}}\Delta
x^{n+1}_{j-\frac{1}{2}}\left[B_{2}c_{j-\frac{1}{2}}^{n+1}+2B_{1}^{2}\left|\frac{\dot{\rho}}{\rho}\right|\Delta
x^{n+1}_{j-\frac{1}{2}}\right.\\\
\left.+\sqrt{\left(B_{2}c_{j-\frac{1}{2}}^{n+1}+2B_{1}^{2}\left|\frac{\dot{\rho}}{\rho}\right|\Delta
x^{n+1}_{j-\frac{1}{2}}\right)^{2}+\left(c_{j-\frac{1}{2}}^{n+1}\right)^{2}}\,\right]^{-1}\;.\end{multlined}\Delta
t^{n+\frac{3}{2}}_{j-\frac{1}{2}}=K_{\rm t_{1}}\Delta
x^{n+1}_{j-\frac{1}{2}}\left[B_{2}c_{j-\frac{1}{2}}^{n+1}+2B_{1}^{2}\left|\frac{\dot{\rho}}{\rho}\right|\Delta
x^{n+1}_{j-\frac{1}{2}}\right.\\\
\left.+\sqrt{\left(B_{2}c_{j-\frac{1}{2}}^{n+1}+2B_{1}^{2}\left|\frac{\dot{\rho}}{\rho}\right|\Delta
x^{n+1}_{j-\frac{1}{2}}\right)^{2}+\left(c_{j-\frac{1}{2}}^{n+1}\right)^{2}}\,\right]^{-1}\;.$
(A21)

The next time step is then set as

$\Delta t^{n+\frac{3}{2}}=\min{\left(\Delta
t^{n+\frac{3}{2}}_{j-\frac{1}{2}},K_{\rm t_{2}}\Delta
t^{n+\frac{1}{2}}\right)}\;,$ (A22)

to ensure stability and to limit the increase in the time step to a factor
$K_{\rm t_{2}}$.

Although we will discuss the code in dimensional parameters here, these
equations are implemented in the code in a dimensionless form to reduce
numerical errors. All the parameters, with the exception of temperature, are
non-dimensionalised using combinations of: the ground velocity, $u_{\rm G}$;
the radius of the planet, $R_{\rm p}$; and $p_{0}$.

Note that the gravity field in the 1D calculation is assumed to be fixed and
radial. Since the time for loss is small ($\sim 100$ s) compared to the
timescale of the impact, any change to the field over this period is minimal,
but the field could still be distorted by the deformation of the target and
the presence of the impactor. As long as the gravity field does not vary
significantly from a $1/r^{2}$ dependence, a local escape velocity in the
direction of the ground motion should be sufficient to scale our results to
calculate the loss in 3D impact simulations (Kegerreis et al., 2020a, b).

###  A.2 Implementation of Equations of State

In this work we have used a variety of equations of state, in particular for
water. The EOS dictates the method used to solve the finite difference form of
the energy equation (Equation A12). We will now outline the equations of state
used and their implementation.

####  A.2.1 Ideal Gas Equation of State

The ideal gas EOS is the simplest EOS for gases and its implementation is
straightforward. The EOS can be written in the form

$p=(\gamma-1)\rho\epsilon\;,$ (A23)

where $\gamma$ is the ratio $c_{\rm p}/c_{\rm V}$, and $c_{\rm p}$ and $c_{\rm
V}$ are the specific heat capacities at constant pressure and volume
respectively. This linear relation between pressure and energy allows for an
analytical solution to the energy equation,

$\begin{multlined}\epsilon^{n+1}_{j-\frac{1}{2}}=\left[\epsilon^{n}_{j-\frac{1}{2}}+\right.\\\
\left.2\left(p^{n}_{j-\frac{1}{2}}+2q^{n+\frac{1}{2}}_{j-\frac{1}{2}}\right)\frac{\left(\rho^{n+1}_{j-\frac{1}{2}}-\rho^{n}_{j-\frac{1}{2}}\right)}{\left(\rho^{n+1}_{j-\frac{1}{2}}-\rho^{n}_{j-\frac{1}{2}}\right)^{2}}\right]\\\
\left[1-(\gamma-1)\rho^{n+1}_{j-\frac{1}{2}}\right]^{-1}\;,\end{multlined}\epsilon^{n+1}_{j-\frac{1}{2}}=\left[\epsilon^{n}_{j-\frac{1}{2}}+\right.\\\
\left.2\left(p^{n}_{j-\frac{1}{2}}+2q^{n+\frac{1}{2}}_{j-\frac{1}{2}}\right)\frac{\left(\rho^{n+1}_{j-\frac{1}{2}}-\rho^{n}_{j-\frac{1}{2}}\right)}{\left(\rho^{n+1}_{j-\frac{1}{2}}-\rho^{n}_{j-\frac{1}{2}}\right)^{2}}\right]\\\
\left[1-(\gamma-1)\rho^{n+1}_{j-\frac{1}{2}}\right]^{-1}\;,$ (A24)

which exclusively uses quantities that are known from the previous time step.
The pressure at the next time step is given by

$p^{n+1}_{j-\frac{1}{2}}=(\gamma-1)\rho^{n+1}_{j-\frac{1}{2}}\epsilon^{n+1}_{j-\frac{1}{2}}\;.$
(A25)

The only parameter that needs to be specified for an ideal gas EOS is
$\gamma$, although other parameters are needed to initialise the atmosphere
(see Appendix A.4).

####  A.2.2 Tabulated Equations of State

We have also added to WONDY the ability to use a tabulated EOS in the form of
a $T$-$\rho$ table. In this work we have used the water table form Senft &
Stewart (2008), a high quality water EOS based partially on the IAPWS water
EOS (see Appendix A.2.3) but extended to the high pressures relevant for
planetary impacts.

Use of a tabulated EOS requires a numerical solution to the energy equation.
We solve for the temperature at each time step that solves the energy equation
using a Newton-Raphson method. The solution is found iteratively using

$T_{i+1}=T_{i}-\frac{f(T_{i})}{f^{\prime}(T_{i})}\;,$ (A26)

where the subscripts mark the number of the iteration,

$\begin{multlined}f(T_{i})=\epsilon(\rho^{n+1},T_{i})-\epsilon^{n}\\\
-\frac{\left(\rho^{n+1}-\rho^{n}\right)}{\left(\rho^{n+1}-\rho^{n}\right)^{2}}\left(p(\rho^{n+1},T_{i})-p^{n}-2q^{n+\frac{1}{2}}\right)\;,\end{multlined}f(T_{i})=\epsilon(\rho^{n+1},T_{i})-\epsilon^{n}\\\
-\frac{\left(\rho^{n+1}-\rho^{n}\right)}{\left(\rho^{n+1}-\rho^{n}\right)^{2}}\left(p(\rho^{n+1},T_{i})-p^{n}-2q^{n+\frac{1}{2}}\right)\;,$
(A27)

and

$\begin{multlined}f^{\prime}(T_{i})=c_{V}(\rho^{n+1},T_{i})\\\
-\frac{\left(\rho^{n+1}-\rho^{n}\right)}{\left(\rho^{n+1}-\rho^{n}\right)^{2}}\left[\frac{\partial\,p}{\partial
T}\bigg{\rvert}_{\rho}(\rho^{n+1},T_{i})\right]\;.\end{multlined}f^{\prime}(T_{i})=c_{V}(\rho^{n+1},T_{i})\\\
-\frac{\left(\rho^{n+1}-\rho^{n}\right)}{\left(\rho^{n+1}-\rho^{n}\right)^{2}}\left[\frac{\partial\,p}{\partial
T}\bigg{\rvert}_{\rho}(\rho^{n+1},T_{i})\right]\;.$ (A28)

Note we have dropped the spatial subscripts for clarity as all values required
are for the zone under consideration. The values of; $\epsilon$, $p$, $c_{V}$
and $\partial\,p/\partial T\rvert_{\rho}$ are found by linear interpolation of
the tabulated values. The derivatives are pre-calculated at each point in the
table using a second order modified central difference method for unevenly
spaced points. The initial value for the iteration is determined by assuming
an isentropic volume change from the previous time step. The change in
temperature, $\Delta T$, for a given density change, $\Delta\rho$, is given by

$\Delta T=\frac{\partial\,T}{\partial\rho}\bigg{\rvert}_{S}\Delta\rho\;,$
(A29)

which can be written

$\Delta T=\frac{T}{c_{V}\rho^{2}}\frac{\partial\,p}{\partial
T}\bigg{\rvert}_{\rho}\Delta\rho\;,$ (A30)

which is straightforward to calculate from the tabulated values. The solution
is considered converged when

$\bigg{|}\frac{T_{i+1}-T_{i}}{T_{i+1}}\bigg{|}<10^{-10}\;.$ (A31)

A maximum of 1000 iterations are allowed but the routine typically converges
within $2$-$3$ iterations.

####  A.2.3 IAPWS Water Equation of State

The IAPWS water EOS (Wagner, 2002) is a high quality empirically fitted EOS
widely used in industrial applications. The advantages of this EOS is that it
provides a smooth form for the material properties, however the EOS is limited
in its pressure range and therefore it can only be used for lower values of
$u_{\rm G}$. The solution to the energy equation using this EOS is found in
the same way as for the tabulated EOS (Appendix A.2.2), except with the
thermodynamic parameters and their derivatives calculated using the fortran
subroutines provided by the IAPWS rather than by interpolation of an EOS
table.

####  A.2.4 Tillotson Equation of State

The Tillotson EOS (Tillotson, 1962) has been widely used in shock physics and
the study of planetary collisions. The EOS is based on empirical fits along
the hugoniot forced to converge to the Thomas-Fermi limit at high pressures.
Although the EOS is designed for use primarily in compression it includes
parameters to approximate the behaviour of material upon release. For this
work we use the parameters from O’Keefe & Ahrens (1982) for water.

Following Melosh (1989) the EOS implemented here has four different regimes:
compressed states ($\rho/\rho_{0}\geq 1$, where $\rho_{0}$ is the reference
density), cold expanded states ($\rho/\rho_{0}<1$ and
$\epsilon<\epsilon_{iv}$, where $\epsilon_{iv}$ is the energy of incipient
vaporisation), hot expanded states ($\rho/\rho_{0}\leq 1$ and
$\epsilon>\epsilon_{cv}$, where $\epsilon_{cv}$ is the energy of complete
vaporisation), and a transition region between the hot and cold expanded
states (where $\rho/\rho_{0}<1$ and $\epsilon_{iv}<\epsilon<\epsilon_{cv}$).
Firstly, for compressed states and for cold expanded states the pressure is
given by

$p=\left[a+\frac{b}{\frac{\epsilon}{\epsilon_{0}\eta^{2}}+1}\right]\rho\epsilon+A\mu+B\mu^{2}\;,$
(A32)

where: $a$, $b$, $A$, $B$ and $\epsilon_{0}$ are fitted parameters;
$\eta=\rho/\rho_{0}$; $\mu=\eta-1$; and $\rho_{0}$ is the reference density of
the material. In order to ensure numerical stability, a low density cutoff is
applied in the cold expanded state when $\rho/\rho_{0}<0.8$, where a nominal
pressure of $10$ Pa is given to a zone. In hot expanded states:

$\begin{multlined}p=a\rho\epsilon\\\
+\left[\frac{b\rho\epsilon}{\frac{\epsilon}{\epsilon_{0}\eta^{2}}+1}+A\mu\exp\left\\{-\beta\left(\frac{1}{\eta}-1\right)\right\\}\right]\\\
\exp\left\\{-\alpha\left(\frac{1}{\eta}-1\right)^{2}\right\\}\;,\end{multlined}p=a\rho\epsilon\\\
+\left[\frac{b\rho\epsilon}{\frac{\epsilon}{\epsilon_{0}\eta^{2}}+1}+A\mu\exp\left\\{-\beta\left(\frac{1}{\eta}-1\right)\right\\}\right]\\\
\exp\left\\{-\alpha\left(\frac{1}{\eta}-1\right)^{2}\right\\}\;,$ (A33)

where $\alpha$ and $\beta$ are two more fitted parameters that control the
rate at which the EOS tends to the ideal gas law. For expanded states with
internal energies between that of incipient and complete vaporization a hybrid
formula is used to transition smoothly between the cold and hot expanded
states,

$p=\frac{\left(\epsilon-\epsilon_{iv}\right)p_{E}+\left(\epsilon_{cv}-\epsilon\right)p_{C}}{\epsilon_{cv}-\epsilon_{iv}}\;,$
(A34)

where $p_{C}$ is the pressure calculated using Equation A32 and $p_{E}$ is the
pressure calculated using Equation A33.

Using the Tillotson EOS the energy equation can be solved analytically.
However, identifying the relevant regime to use in our code requires knowledge
of the solution. We therefore find the possible solutions for $\epsilon$ in
all four regimes and use the physically meaningful solution. Subbing Equations
A32, A33 and A34 into the energy equation in finite difference form (Equation
A12) produces three separate quadratic equations. These can then be solved
analytically giving two solutions for each equation. In the case that the
material is in compression we choose the solution of the equation derived from
Equation A32 that: 1) gives the correct sign of the change in $\epsilon$ for
the density change; and 2) is closest to the previous value of $\epsilon$. If
the material is in an expanded state solutions from all three equations are
potentially valid and so we choose the solution that: 1) is in the correct
energy range for its regime; 2) gives the correct sign of the change in
$\epsilon$ for the density change; and 3) is closest to the previous value of
$\epsilon$. The relevant equation can then be used to find the pressure. Note
that care must be taken in averaging properties in the cross over region.

In order to calculate artificial viscosity we also need to calculate the
adiabatic sound speed. The sound speed is defined as

$c^{2}=\frac{\partial p}{\partial\rho}\bigg{\rvert}_{S}\;,$ (A35)

where $S$ is specific entropy. Using standard differential relations this can
be expressed as

$c^{2}=\frac{\partial p}{\partial\rho}\bigg{\rvert}_{\epsilon}+\frac{\partial
p}{\partial\epsilon}\bigg{\rvert}_{\rho}\frac{\partial\epsilon}{\partial\rho}\bigg{\rvert}_{S}\;.$
(A36)

Substituting for the final term from the fundamental thermodynamic relation,

$c^{2}=\frac{\partial
p}{\partial\rho}\bigg{\rvert}_{\epsilon}+\frac{p}{\rho^{2}}\frac{\partial
p}{\partial\epsilon}\bigg{\rvert}_{\rho}\;.$ (A37)

All the terms in this equation are known or can be found by differentiating
the relevant equation for pressure in each regime. Equation A37 is used to
calculate the sound speed in each regime with the exception of the low density
cutoff where a nominal sound speed of $0.5u_{\rm G}$ is used.

Due to the limitations of the Tillotson EOS, WONDY can encounter unphysical
solutions in expanded states. In such cases, the code is unable to continue
and most of our calculations using Tillotson did not complete the full model
run. However, unphysical states are usually encountered late in the simulation
when the loss has plateaued and so the total loss calculated is only slightly
underestimated.

###  A.3 Boundary Conditions

The breakout of the shock at the base of the atmosphere/ocean is simulated by
giving the lower boundary an initial velocity, $u_{\rm G}$, and then allowing
the boundary to evolve ballistically. The lower boundary evolves following the
momentum equation (Equation A5) ignoring forces other than gravity, which in
finite difference form is expressed as:

$\begin{multlined}a^{n}_{0}=-\frac{GM_{p}}{\left(x_{0}^{n}\right)^{2}}\;.\end{multlined}a^{n}_{0}=-\frac{GM_{p}}{\left(x_{0}^{n}\right)^{2}}\;.$
(A38)

In this framework, $u_{\rm G}$ is the particle velocity of the surface upon
release of the shock in the planet to the base of the atmosphere or ocean. As
for the case of the ocean/atmosphere interface, $u_{\rm G}$ can be determined
using an impedance-match calculation. Due to the lower impedance of gases as
opposed to water, release of the shock directly to an atmosphere results in
higher ground velocities than in the presence of an ocean and in both cases
the ground velocity is lower than if the shock released to vacuum. When
calculating loss by convolving 1D results with the ground velocity calculated
in 2D or 3D impact simulations it is imperative that corrections are made to
translate the ground velocity.

Driving the lower boundary in this way creates some numerical complications.
The sudden acceleration of the boundary generates excessively high artificial
viscosities in the lowermost zones. This causes these zones to have
artificially high temperatures and low densities which can be outside the
range of the EOS. We overcome this by limiting the artificial viscosity in the
first four zones to be less than $q_{max}$. The exact value of $q_{max}$,
within a reasonable range, has little effect on the amount of loss so we set
it to $0.2$ of that calculated for the first time step for the first zone
using Equation A17. We also initialise the first zone with $q^{0}_{1}=q_{max}$
and the second and third zones with
$q^{0}_{j-\frac{1}{2}}=0.01q^{0}_{j-\frac{3}{2}}$. This helps to maintain
reasonable values for variables in the lowermost zones with no affect on the
final results.

In some simulations (see Section 6) the velocity of the ground is forced to
increase linearly over some finite rise time, $t_{\rm rise}$, until the
specified ground velocity is reached. The specification of the artificial
viscosity in cells just above the boundary is the same as in the standard
case.

The lower boundary is assumed to stop when it returns to its initial position.
However it is necessary to slow the boundary gradually to avoid numerical
errors. The boundary is prescribed to fall exponentially after it returns to a
position $x_{0}<R_{t}+0.15\left(x_{0}^{max}-R_{t}\right)$ where $x_{0}^{max}$
is the maximum height reached by the boundary. The rate is chosen so that the
velocity of the boundary is continuous from before being slowed. In finite
difference form

$u_{0}^{n+\frac{1}{2}}=u_{0}^{n-\frac{1}{2}}\exp\left\\{\frac{u_{slow}\Delta
t^{n-\frac{1}{2}}}{0.15\left(x_{0}^{max}-R_{t}\right)}\right\\}\;,$ (A39)

where $u_{slow}$ is the velocity of the boundary at the point the boundary
starts to be slowed. Since the slow down happens late in the simulation the
time step can be quite large by the time the slow down begins. We therefore
artificially reduce the time step when the boundary is
$R_{t}+0.15\left(x_{0}^{max}-R_{t}\right)<x_{0}<R_{t}+0.3\left(x_{0}^{max}-R_{t}\right)$
to $\leq 0.1$ s so that the slow down can be properly calculated. After slow
down begins the time step is allowed to evolve normally.

In some simulations (see Section 6) the ground is stopped at some time,
$t_{\rm stop}$. This is achieved by immediately setting the velocity of the
boundary to zero at the specified time. Once the boundary has stopped, it
begins to fall under gravity and its evolution continues as described above.

The boundary at the top of the atmosphere is treated as a stress free boundary
i.e. a vacuum. This is implemented by having an additional zone at the top of
the atmosphere with $(\rho,p)=0$ at all time steps. The finite difference
equations are then applied as normal.

###  A.4 Atmosphere and Ocean Initialisation

The initial conditions for our simulations are a hydrostatic atmosphere and,
optionally, an ocean. We set the extent of the atmosphere by setting the
pressure and temperature at the base of the atmosphere. The atmosphere is then
assumed to be polytropic,

$\left(\frac{p}{p_{0}}\right)^{\frac{\gamma-1}{\gamma}}=\left(\frac{\rho}{\rho_{0}}\right)^{\gamma-1}=\frac{T}{T_{0}}\;,$
(A40)

where $p_{0}$, $\rho_{0}$ and $T_{0}$ are the pressure, density and
temperature at the base of the atmosphere. We then integrate up from the base
of the atmosphere using a 4th order Runge-Kutta method to find the top of the
atmosphere (in this model the pressure drops to zero at a finite height). The
atmosphere is then divided into the specified number of zones, $N_{\rm atm}$,
each of equal thickness. For the ideal gas EOS the structure of the atmosphere
is determined by $p_{0}$, $T_{0}$, $\gamma$ and $m_{\rm a}$ (the molar mass of
the gas) which sets the density at a given $p$, $T$. The values for $\gamma$
and $m_{\rm a}$ for each of the atmospheres used here are given in Table 2.

Atmosphere | $\gamma$ |  $m_{\rm a}$ [g]
---|---|---
H2 | $1.4$ | $2$
H2O | 1.25 | 18
CO2 | $1.29$ | $44$
N2 | $1.4$ | $28$
Earth-like | 1.4 | 29
Table 2: The ideal gas parameters (ratio of specific heat capacities,
$\gamma$, and molar mass, $m_{\rm a}$)for the atmospheres used in this paper.

The extent of the ocean is set by specifying the depth of the ocean, $H_{\rm
oc}$. The temperature and pressure of the bottom of the atmosphere are assumed
to be continuous with the ocean. The exact structure of the ocean however
depends on the EOS being used. For the tabulated EOS an isothermal ocean is
initialised with temperature $T_{0}$. In the case of the IAPWS EOS an
isentropic ocean is used with the entropy set by $p_{0}$ and $T_{0}$. In order
to compare directly with the work of Genda & Abe (2005) when using the
Tillotson EOS we initialise a constant specific internal energy ocean with
$\epsilon=120$ J kg-1. In each case a 4th order Runge-Kutta method is used to
integrate down from the surface to find the properties for each of the $N_{\rm
oc}$ ocean cells. Due to the limited compressibility of water, the ocean is
close to isothermal in both the isentropic and isoenergetic initializations.
Cells are of equal thickness.

The mass, radii, and escape velocity of the planets used in this study are
given in Table 3. For Earth and Mars -mass planets the radii were taken as the
present-day equatorial radii of those planets. For intermediate mass planets,
the radii were calculated using the HERCULES planetary structure code (Lock &
Stewart, 2017; Lock, 2019) using the present-day Earth’s core mass fraction of
0.323 (Yoder, 1995). In HERCULES, a body is modeled as consisting of a series
of nested concentric layers of constant density and a potential field method
is used to calculate the equilibrium structure of the body with a given
thermal state, composition, mass and angular momentum, using realistic
equations of state. For this work, the mantle and core of each planet were
assumed to be forsterite and pure iron, respectively, and were modeled using
the the ANEOS equation of state (EOS) model (Thompson & Lauson, 1972; Canup,
2012; Melosh, 2007). The EOS are documented as ‘aneos-T70’ (iron) and ‘aneos-
gadget’ (forsterite) respectively in two Zenodo repositories (Stewart, 2020;
Stewart et al., 2019). The mantle was assumed to be isentropic with a specific
entropy of 3.2 kJ K-1 kg-1 corresponding to a mantle potential temperature of
around 1900 K. This thermal state approximates that of the Hadean mantle. The
core was also assumed to be isentropic and have a thermal state similar to the
present day. The core specific entropy was set at 1.5 kJ K-1 kg-1,
corresponding to a temperature of 3800 K at the pressure of the present-day
core-mantle boundary. We used the same HERCULES parameters as in previous
studies, that have been shown to accurately model the structure of Earth-like
planets (Lock & Stewart, 2017, 2019; Lock, 2019; Lock et al., 2020). The
planets were modeled by $N_{\rm lay}^{\rm core}=20$ evenly spaced layers in
the core and $N_{\rm lay}^{\rm mantle}=80$ layers in the mantle, with
$N_{\mu}=1000$ points describing the shape of each layer. The expression for
the gravity field was truncated at order $2k_{\rm max}=12$. The minimum
pressure at the surface of the planet was set to 10 bar. The tolerance for the
convergence of the shape of equipotential layers was $\xi_{\rm
toll}^{\mu}=10^{-10}$ and the tolerance for the convergence of the mass of the
planet was $\xi_{\rm toll}=10^{-8}$. The step used for calculating gradients
in the solution algorithm was $\delta\xi=10^{-2}$. For further details of the
definitions of these parameters the reader is referred to the HERCULES user
manual (Lock, 2019).

Mass [$M_{\rm Earth}$] | Radius [km] |  $v_{\rm esc}$ [km s-1]
---|---|---
0.107 | 3396 | 5.02
0.3 | 4607 | 7.20
0.5 | 5366 | 8.62
0.7 | 5907 | 9.72
0.9 | 6337 | 10.64
1.0 | 6371 | 11.18
Table 3: Properties of planets used for hydrodynamic simulations in this study

##  Appendix B Sensitivity Tests for 1D Model

###  B.1 Sensitivity to Code Parameters

We have tested the sensitivity of our 1D models to the intrinsic code
parameters ($B_{1}$, $B_{2}$, $K_{\rm t_{1}}$ and $K_{\rm t_{2}}$) as well as
the number of zones in the ocean and atmosphere ($N_{\rm atm}$ and $N_{\rm
oc}$). To do this we ran a series of calculations varying each of the
parameters in turn and comparing the calculated loss. The final results showed
little variation with any of the parameters over the ranges we explored (e.g.,
Figure 16), although some runs with lower $B_{1}$ values failed. For some
combinations of parameters ringing was observed around the shock front
indicating insufficient numerical viscosity, but such effects did not affect
the final loss. We are therefore confident that our results are not
significantly affected by our choice of these parameters.

Parameter | Standard Value
---|---
$B_{1}$ | $2$
$B_{2}$ | $0.1$
$K_{\rm t_{1}}$ | $0.75$
$K_{\rm t_{2}}$ | $1.05$
$N_{\rm oc}$ | $500$
$N_{\rm atm}$ | $500$
$t_{\rm final}$ |  $5\times 10^{4}$ s
Table 4: Values for constants used in all 1D model runs for which results are
presented here (details in text). $B_{1}$ and $B_{2}$ are used to determine
the strength of the artificial viscosity. $K_{\rm t_{1}}$ and $K_{\rm t_{2}}$
control the size of each time step. $N_{\rm oc}$ and $N_{\rm atm}$ are the
number of zones used for the ocean and atmosphere respectively. $t_{\rm
final}$ is the total runtime for the simulations. Figure 16: The results of
our 1D calculations are not sensitive to the intrinsic parameters in the code.
Shown is the loss as a function of ground velocity for an example simulation
($H_{\rm oc}=3$ km, $p_{0}=100$ bar and a CO2 atmosphere at $300$ K)
calculated varying run parameters $B_{1}$ (1-10), $B_{2}$ (0.1-5), $K_{\rm
t_{1}}$ (1.05-1.2), $K_{\rm t_{2}}$ (0.6-0.9), $N_{\rm oc}$ (250-750), and
$N_{\rm atm}$ (250-750). Each parameter was varied while holding all others
constant at the values given in Table 4.

##  Appendix C Parameterization of loss as a function of ground velocity

In this section, we present a parameterization of the loss due to a given
ground motion from a body with a given escape velocity and atmosphere to ocean
mass ratio. This parameterization is not intended to give physical insight,
but rather to provide expressions for calculating global loss from 3D impact
simulations (e.g., Kegerreis et al., 2019; Denman et al., 2020, 2022). We will
first describe the parameterization for the limiting case of no ocean and then
for the general scenario including oceans and atmospheres of varying masses.
Python functions and an interactive widget that implement this
parameterization are available through GitHub.

###  C.1 The no-ocean case

For the no-ocean case, we describe the loss using a modified logistics
function,

$f_{\rm NO}\left(\frac{u_{\rm G}}{v_{\rm
esc}}\right)=\alpha_{1}\left[1+\exp{\left\\{\alpha_{2}\left(\frac{u_{\rm
G}}{v_{\rm
esc}}\right)+\alpha_{3}\right\\}}\right]^{\alpha_{4}}+\alpha_{5}\;,$ (C1)

where $u_{\rm G}$ is the ground velocity, $v_{\rm esc}$ is the escape
velocity, and $\alpha_{i}$ are constants. By requiring loss to be zero in the
case of zero ground velocity, i.e., $f_{\rm NO}=0$ when $u_{\rm G}=0$, we find

$\alpha_{5}=-\alpha_{1}[1+\exp{(\alpha_{3})}]^{\alpha_{4}}\;.$ (C2)

Similarly, by requiring loss to be complete when $u_{\rm G}=v_{\rm esc}$ we
find

$\alpha_{1}=\frac{1}{[1+\exp{(\alpha_{2}+\alpha_{3}})]^{\alpha_{4}}-[1+\exp{(\alpha_{3})}]^{\alpha_{4}}}\;.$
(C3)

In addition, we force $f_{\rm NO}=0$ when $u_{\rm G}/v_{\rm esc}<0$ and
$f_{\rm NO}=1$ when $u_{\rm G}/v_{\rm esc}>1$.

In the no-ocean case we set $\alpha_{3}=-1$ and performed a least squares fit
to find $\alpha_{2}=-4.32$ and $\alpha_{4}=-3.77$. For the fit we used the
calculations shown in Figure 5A: 0.1, 0.5, 1, 5, 10, 50, 100, and 500 bar
atmospheres on planets with mass 0.107 ($M_{\rm Mars}$), 0.3, 0.5, 0.7, 0.9
and 1 $M_{\rm Earth}$. The surface temperature was 300 K, and the atmosphere
was CO2 with a mean molecular weight of $m_{\rm a}=44$ and a ratio of specific
heat capacities of $\gamma=1.29$. Loss was calculated at 0.05 $v_{\rm esc}$
increments from 0.05 to 0.95 $v_{\rm esc}$.

This parameterization offers a very accurate fit to our model runs. Figure 5C
shows the misfit between our parameterization and the fitted calculations. The
root mean square (RMS) misfit for all runs is 0.0046 and the maximum misfit is
0.0089.

###  C.2 Loss in the presence of an ocean

We chose to parameterize loss in the $\mathcal{R}$-$u_{\rm G}$-loss space
(Figure 10). Loss is described by a modified logisitics function

$f(u_{\rm G},v_{\rm
esc},\mathcal{R})=\alpha_{1}[1+\exp{(\alpha_{2}\log_{10}{\left(\mathcal{R}\right)}+\alpha_{3}})]^{\alpha_{4}}+\alpha_{5}\;,$
(C4)

where $\alpha_{i}$ are functions of the ground velocity and the escape
velocity, $v_{\rm esc}$. We find that loss due to a given ground velocity
tends to the no-ocean case as $\mathcal{R}\rightarrow\infty$ (Figure 10) and
so we enforce $f\rightarrow f_{\rm NO}$ as $\mathcal{R}\rightarrow\infty$ by
setting

$\alpha_{1}=f_{\rm NO}-\alpha_{5}\;,$ (C5)

and requiring that $\alpha_{2}\leq 0$.

The velocity dependence of $\alpha_{i}$ are given by arbitrary function chosen
as they provided a good fit to our simulation results. $\alpha_{2}$ and
$\alpha_{3}$ are described by 3rd-order polynomials:

$\alpha_{i}=a^{i}_{1}+a^{i}_{2}\frac{u_{\rm G}}{v_{\rm
esc}}+a^{i}_{3}\left(\frac{u_{\rm G}}{v_{\rm
esc}}\right)^{2}+a^{i}_{4}\left(\frac{u_{\rm G}}{v_{\rm esc}}\right)^{3}\;,$
(C6)

and $\alpha_{4}$ is given by

$\displaystyle\begin{aligned}
\alpha_{4}=&a^{4}_{1}\exp\left\\{\sin{\left(\frac{u_{\rm G}}{v_{\rm
esc}}\pi\right)}\right.\\\ &\left.\times\left[a^{4}_{2}\frac{u_{\rm G}}{v_{\rm
esc}}+a^{4}_{3}\left(\frac{u_{\rm G}}{v_{\rm
esc}}\right)^{2}+a^{4}_{4}\left(\frac{u_{\rm G}}{v_{\rm
esc}}\right)^{3}\right]\right\\}\;,\end{aligned}$ (C7)

where $a^{i}_{j}$ are coefficients that can depend on $v_{\rm esc}$. The
functional form for the ground velocity dependence of $\alpha_{5}$ is more
complex. The value of $\alpha_{5}$ is the asymptote for loss as
$\mathcal{R}\rightarrow 0$ and is close to the lowest $\mathcal{R}$ loss curve
we calculate (e.g., the yellow line in Figure 6). To describe a loss curve in
$u_{\rm G}$-loss space, we divide the functional form for $\alpha_{5}$ into
two parts for the atmospheric loss ($g_{1}$) and ocean loss ($g_{2}$) regimes
respectively:

$\alpha_{5}=\begin{dcases*}0&\; $u_{\rm G}\leq 0$\\\ g_{1}&\; $0<u_{\rm
G}<u_{\rm tal}$\\\ 1&\; $u_{\rm tal}\leq u_{\rm G}<v_{\rm esc}$ \; \& \;
$g_{2}\leq 1$\\\ g_{2}&\; $u_{\rm tal}\leq u_{\rm G}<v_{\rm esc}$ \; \& \;
$g_{2}>1$\\\ 2&\; $u_{\rm G}\geq v_{\rm esc}$ \; ,\end{dcases*}$ (C8)

where

$g_{1}=a^{5}_{1}\left(\frac{u_{\rm G}}{v_{\rm
esc}}\right)^{a^{5}_{3}}\exp{\left(a^{5}_{2}\frac{u_{\rm G}}{v_{\rm
esc}}\right)}\;,$ (C9) $\displaystyle\begin{aligned}
g_{2}=2&+a^{5}_{4}\left[\frac{u_{\rm G}-u_{\rm tal}}{v_{\rm esc}-u_{\rm
tal}}-1\right]\\\ &+a^{5}_{5}\left[\left(\frac{u_{\rm G}-u_{\rm tal}}{v_{\rm
esc}-u_{\rm tal}}\right)^{2}-1\right]\;,\end{aligned}$ (C10)

and $u_{\rm tal}$ is the velocity at which total atmospheric loss is achieved
(i.e., when $g_{1}=1$) given by

$u_{\rm tal}=v_{\rm
esc}\left(\frac{a^{5}_{3}}{a^{5}_{2}}\right)\min_{k=-1,0}W_{k}\left\\{\frac{a^{5}_{3}}{a^{5}_{2}}\left(a^{5}_{1}\right)^{-1/a^{5}_{3}}\right\\}\;,$
(C11)

where $W_{k}$ is the Lambert $W$ function of order $k$. The velocity regime
between the atmospheric and ocean loss regimes where $\alpha_{5}=1$ loss
emulates a feature seen in our low $\mathcal{R}$ simulations where significant
ocean loss is delayed for a short range of $u_{\rm G}$ after total atmospheric
loss is achieved (Figure 6).

Dependence on planetary mass is captured by a linear dependence of the
coefficients for $\alpha_{2}$ ($a^{2}_{j}$), $\alpha_{4}$ ($a^{4}_{j}$), and
$\alpha_{5}$ ($a^{5}_{j}$) on $v_{\rm esc}$ such that

$a^{i}_{j}=a^{i}_{j,1}+a^{i}_{j,2}\frac{v_{\rm esc}}{v_{\rm esc}^{\rm
Earth}}\;,$ (C12)

where $v_{\rm esc}^{\rm Earth}=11.2$ km s-1 is the escape velocity of Earth.
Note that it is possible to parameterize the $v_{\rm esc}$ dependence on
$\alpha_{5}$ by using a function for the absolute value of $u_{\rm G}$
convolved with the no-ocean loss function, but we found that, because of the
complications of the varying pressure of release and the highly non-linear no-
ocean loss function, a simple escape velocity scaling of the parameters was
able to give a more accurate fit.

To determine the parameters $a^{i}_{j,n}$ we performed a least squares fit on
the results of the simulations shown in Figure 10: all combinations of 1, 5,
10, 50, 100, 300 and 500 bar atmospheres and oceans of depths 0.1, 0.5, 1, 2,
3, 5, 10, 20 and 30 km, on planets with mass $M_{\rm Mars}$ (0.107 $M_{\rm
Earth}$) and 0.3, 0.5, 0.7, 0.9 and 1 $M_{\rm Earth}$. Additional simulations
with 900 bar atmospheres and oceans of 0.1 km depth were also included. The
surface temperature was assumed to be 300 K, and the atmosphere was CO2
($m_{\rm a}=44$, $\gamma=1.29$). Loss was calculated at 0.05 $v_{\rm esc}$
increments from 0.05 to 0.95 $v_{\rm esc}$. The best-fit parameters are given
in Table 5.

$i$ |  $j$ |  $a_{j,1}^{i}$ |  $a_{j,2}^{i}$
---|---|---|---
$2$ | 1 | -16.2584384 | 14.2828078
2 | 2 | 71.5029708 | -86.4465281
2 | 3 | -117.862430 | 164.871407
2 | 4 | 60.1395645 | -95.6142191
3 | 1 | -7.66878923 | -
3 | 2 | 43.9524753 | -
3 | 3 | -105.616622 | -
3 | 4 | 77.4557563 | -
4 | 1 | -0.767532594 | -0.263620386
4 | 2 | -12.0443695 | -11.5024085
4 | 3 | 44.7091911e | 39.8696345
4 | 4 | -38.8346991 | -28.7008315
5 | 1 | 10660.1668 | 13897.9517
5 | 2 | -12.9490638 | 0.640100268
5 | 3 | 7.35306571 | -0.384763832
5 | 4 | 2.12776033 | -0.324256357
5 | 5 | -1.14723912 | 0.418100341
Table 5: Fitted coefficients for each of the parameters $a_{j,k}^{i}$ (see
Equation C12).

We find that our parameterization offers an accurate representation of the
dependence of loss on $\mathcal{R}$, $u_{\rm G}$ and $M_{\rm p}$. Figure 10
shows the fit in $\mathcal{R}$-loss space for example ground velocities. The
global RMS is 0.041, the RMS at any given ground velocity does not exceed
0.062, and the maximum misfit of any simulation result is 0.25.
Unsurprisingly, the misfit is largest in the transition between the low and
high $\mathcal{R}$ regimes where the gradient of the function in the
$\mathcal{R}$-loss regime is greatest and the slight sensitivity to parameters
such as absolute velocity, atmospheric pressure, ocean height etc. are
greatest.

The parameterization also reproduces realistic loss curves in $u_{\rm G}$-loss
space. Figure 17 shows our simulation results for the same four example loss
curves as in Figure 6A, but for CO2 atmospheres, (solid lines), with
corresponding curves calculated using our parameterization (dashed lines). The
grey band gives the full range of possible loss for an Earth-mass planet for
any value of $\mathcal{R}$ using our parameterization.

Figure 17: Our parameterization can accurately determine loss as a function of
ground velocity. Shown are our simulation results for the same example surface
conditions as in Figure 6A but for CO2 atmospheres (solid lines) and loss
curves calculated using our parameterization (dashed lines). Grey band shows
the full range of possible loss from Earth-mass bodies determined using our
parameterization.

###  C.3 Connecting shock strength to ground velocity

As discussed in Section 5, the relationship between the strength of the shock
in the planet and the ground velocity varies depending on the properties of
the ocean and/or atmosphere in contact with the ground. It is also worth
noting that the relationship could be affected somewhat by changing the shock
properties of the ground. When calculating the loss due to a given shock, for
example when combining 1D and 3D simulation (Section 7.4), it is therefore
necessary to calculate the relationship between the properties of the shock in
the planet (here we use the particle velocity, $u_{p}$, as the reference
variable) and the ground velocity for the specific combination of ground and
ocean and/or atmosphere that is under consideration. As this relationship can
vary substantially, we do not attempt to provide a comprehensive set of
$u_{p}$-$u_{\rm G}$ relations. Instead, we provide a python package and a
widget that workers can use to calculate individual or sets of $u_{p}$-$u_{\rm
G}$ relations (see data availability) for different surface conditions and a
range of EOS.

We do make one exception. The $u_{p}$-$u_{\rm G}$ relationship for breakout of
a shock into the base of the ocean from a forsterite mantle is relatively
constant over the range of conditions we have considered (oceans of depths
0.1-30 km, and atmospheric pressure up to 900 bar). In this case, the
$u_{p}$-$u_{\rm G}$ relationship can be well described by

$u_{\rm G}=\beta_{1}u_{p}^{\beta_{2}}\;,$ (C13)

where $\beta_{1}=2.6441$, and $\beta_{2}=0.9252$.

## References

  * Abe & Matsui (1988) Abe, Y., & Matsui, T. 1988, Journal of the Atmospheric Sciences, 45, 3081, doi: 10.1175/1520-0469(1988)045<3081:EOAIGH$>$2.0.CO;2
  * Asphaug et al. (2021) Asphaug, E., Emsenhuber, A., Cambioni, S., Gabriel, T. S., & Schwartz, S. R. 2021, Planetary Science Journal, 2, 200, doi: 10.3847/PSJ/ac19b2
  * Bower et al. (2019) Bower, D. J., Kitzmann, D., Wolf, A. S., et al. 2019, Astronomy & Astrophysics, 631, A103, doi: 10.1051/0004-6361/201935710
  * Cameron & Ward (1976) Cameron, A. G. W., & Ward, W. R. 1976, in Lunar and Planetary Science Conference Abstracts, Vol. 7, 120
  * Canup (2004) Canup, R. M. 2004, Annual Review of Astronomy and Astrophysics, 42, 441, doi: 10.1146/annurev.astro.41.082201.113457
  * Canup (2012) —. 2012, Science, 338, 1052, doi: 10.1126/science.1226073
  * Canup & Asphaug (2001) Canup, R. M., & Asphaug, E. 2001, Nature, 412, 708, doi: 10.1038/35089010
  * Caracas & Stewart (2023) Caracas, R., & Stewart, S. T. 2023, Earth and Planetary Science Letters, 608, 118014, doi: 10.1016/j.epsl.2023.118014
  * Carter et al. (2018) Carter, P. J., Leinhardt, Z. M., Elliott, T., Stewart, S. T., & Walter, M. J. 2018, Earth and Planetary Science Letters, 484, 276, doi: 10.1016/j.epsl.2017.12.012
  * Carter et al. (2020) Carter, P. J., Lock, S. J., & Stewart, S. T. 2020, Journal of Geophysical Research: Planets, 125, doi: 10.1029/2019JE006042
  * Chen & Ahrens (1997) Chen, G. Q., & Ahrens, T. J. 1997, Physics of the Earth and Planetary Interiors, 100, 21, doi: 10.1016/S0031-9201(96)03228-1
  * Denman et al. (2022) Denman, T. R., Leinhardt, Z. M., & Carter, P. J. 2022, Monthly Notices of the Royal Astronomical Society, 513, 1680, doi: 10.1093/mnras/stac923
  * Denman et al. (2020) Denman, T. R., Leinhardt, Z. M., Carter, P. J., & Mordasini, C. 2020, Monthly Notices of the Royal Astronomical Society, 496, 1166, doi: 10.1093/mnras/staa1623
  * Elkins-Tanton (2012) Elkins-Tanton, L. T. 2012, Annual Review of Earth and Planetary Sciences, 40, 113, doi: 10.1146/annurev-earth-042711-105503
  * Fegley et al. (2023) Fegley, B., Lodders, K., & Jacobson, N. S. 2023, Geochemistry, 125961, doi: 10.1016/j.chemer.2023.125961
  * Genda & Abe (2003) Genda, H., & Abe, Y. 2003, Icarus, 164, 149, doi: 10.1016/S0019-1035(03)00101-5
  * Genda & Abe (2005) —. 2005, Nature, 433, 842, doi: 10.1038/nature03360
  * Grady et al. (1987) Grady, M., Aylmer, D., Kurat, G., et al. 1987, Memoirs of National Institute of Polar Research. Special Issue, 46, 162
  * Grady et al. (1986) Grady, M. M., Wright, I. P., Carr, L. P., & Pillinger, C. T. 1986, Geochimica et Cosmochimica Acta, 50, 2799, doi: 10.1016/0016-7037(86)90228-0
  * Halliday (2013) Halliday, A. N. 2013, Geochimica et Cosmochimica Acta, 105, 146, doi: 10.1016/j.gca.2012.11.015
  * Hartmann & Davis (1975) Hartmann, W. K., & Davis, D. R. 1975, Icarus, 24, 504, doi: 10.1016/0019-1035(75)90070-6
  * Henderson (1989) Henderson, L. F. 1989, Journal of Fluid Mechanics, 198, 365, doi: 10.1017/S0022112089000170
  * Hicks (1978) Hicks, D. L. 1978, Mathematics of Computation, 32, 1123, doi: 10.1090/S0025-5718-1978-0483944-0
  * Hubbard & MacFarlane (1980) Hubbard, W. B., & MacFarlane, J. J. 1980, Journal of Geophysical Research: Solid Earth, 85, 225, doi: 10.1029/JB085IB01P00225
  * Ikoma & Genda (2006) Ikoma, M., & Genda, H. 2006, The Astrophysical Journal, 648, 696, doi: 10.1086/505780
  * Jontof-Hutter (2019) Jontof-Hutter, D. 2019, Annual Review of Earth and Planetary Sciences, 47, 141, doi: 10.1146/annurev-earth-053018-060352
  * Kasting (1988) Kasting, J. F. 1988, Icarus, 74, 472, doi: 10.1016/0019-1035(88)90116-9
  * Kegerreis et al. (2020a) Kegerreis, J. A., Eke, V. R., Catling, D. C., et al. 2020a, The Astrophysical Journal, 901, L31, doi: 10.3847/1538-4357/ab9810
  * Kegerreis et al. (2019) Kegerreis, J. A., Eke, V. R., Gonnet, P., et al. 2019, Monthly Notices of the Royal Astronomical Society, 487, 5029, doi: 10.1093/mnras/stz1606
  * Kegerreis et al. (2020b) Kegerreis, J. A., Eke, V. R., Massey, R. J., & Teodoro, L. F. A. 2020b, The Astrophysical Journal, 897, 161
  * Kegerreis et al. (2018) Kegerreis, J. A., Teodoro, L. F. A., Eke, V. R., et al. 2018, The Astrophysical Journal, 861, 52, doi: 10.3847/1538-4357/aac725
  * Kerridge (1985) Kerridge, J. F. 1985, Geochimica et Cosmochimica Acta, 49, 1707
  * Kipp & Lawrence (1982) Kipp, M. E., & Lawrence, R. J. 1982, NASA STI/Recon Technical Report N, 83, 10396
  * Lammer et al. (2014) Lammer, H., Stokl, A., Erkaev, N. V., et al. 2014, Monthly Notices of the Royal Astronomical Society, 439, 3225, doi: 10.1093/mnras/stu085
  * Lemmon et al. (2000) Lemmon, E. W., Jacobsen, R. T., Penoncello, S. G., & Friend, D. G. 2000, Journal of Physical and Chemical Reference Data, 29, 331, doi: 10.1063/1.1285884
  * Lock (2023a) Lock, S. 2023a, sjl499/Lock_Stewart_2023_PSJ_atmospheric_loss: Supporting resources for Lock & Stewart (2023), Zenodo, doi: 10.5281/zenodo.8368279
  * Lock (2023b) —. 2023b, sjl499/simples: Shock Impedance Match Package and Loss Event Simulator, Zenodo, doi: 10.5281/zenodo.8367762
  * Lock (2019) Lock, S. J. 2019, Zenodo, doi: 10.5281/zenodo.3509365
  * Lock & Stewart (2017) Lock, S. J., & Stewart, S. T. 2017, Journal of Geophysical Research: Planets, 122, 950, doi: 10.1002/2016JE005239
  * Lock & Stewart (2019) —. 2019, Science Advances, 5, eaav3746, doi: 10.1126/sciadv.aav3746
  * Lock et al. (2018) Lock, S. J., Stewart, S. T., Petaev, M. I., et al. 2018, Journal of Geophysical Research: Planets, 123, 910, doi: 10.1002/2017JE005333
  * Lock et al. (2020) Lock, S. J., Stewart, S. T., & Ćuk, M. 2020, Earth and Planetary Science Letters, 530, 115885, doi: 10.1016/j.epsl.2019.115885
  * Lodders (2003) Lodders, K. 2003, The Astrophysical Journal, 591, 1220, doi: 10.1086/375492
  * Marty (2012) Marty, B. 2012, Earth and Planetary Science Letters, 313-314, 56, doi: 10.1016/j.epsl.2011.10.040
  * Marty et al. (2016) Marty, B., Avice, G., Sano, Y., et al. 2016, Earth and Planetary Science Letters, 441, 91 , doi: http://dx.doi.org/10.1016/j.epsl.2016.02.031
  * Melosh (1989) Melosh, H. J. 1989, Impact cratering: A geologic process (Oxford University Press)
  * Melosh (2003) —. 2003, Journal of Applied Physics, 94, 4320, doi: 10.1063/1.1602566
  * Melosh (2007) —. 2007, Meteoritics & Planetary Science, 42, 2079, doi: 10.1111/j.1945-5100.2007.tb01009.x
  * Moore & Gibson (1969) Moore, C. B., & Gibson, E. K. 1969, Science, 163, 174, doi: 10.1126/science.163.3863.174
  * Moore & Lewis (1966) Moore, C. B., & Lewis, C. F. 1966, Earth and Planetary Science Letters, 1, 376, doi: 10.1016/0012-821X(66)90029-X
  * Mukhopadhyay & Parai (2019) Mukhopadhyay, S., & Parai, R. 2019, Annual Review of Earth and Planetary Sciences, 47, 389, doi: 10.1146/annurev-earth-053018-060238
  * Nakajima & Stevenson (2015) Nakajima, M., & Stevenson, D. J. 2015, Earth and Planetary Science Letters, 427, 286, doi: 10.1016/j.epsl.2015.06.023
  * Odert et al. (2018) Odert, P., Lammer, H., Erkaev, N. V., et al. 2018, Icarus, 307, 327, doi: 10.1016/j.icarus.2017.10.031
  * O’Keefe & Ahrens (1982) O’Keefe, J. D., & Ahrens, T. J. 1982, Journal of Geophysical Research, 87, 6668, doi: 10.1029/JB087iB08p06668
  * Olson & Sharp (2019) Olson, P. L., & Sharp, Z. D. 2019, Physics of the Earth and Planetary Interiors, 294, 106294, doi: 10.1016/j.pepi.2019.106294
  * O’Brien et al. (2014) O’Brien, D. P., Walsh, K. J., Morbidelli, A., Raymond, S. N., & Mandell, A. M. 2014, Icarus, 239, 74, doi: 10.1016/j.icarus.2014.05.009
  * Quintana et al. (2016) Quintana, E. V., Barclay, T., Borucki, W. J., Rowe, J. F., & Chambers, J. E. 2016, The Astrophysical Journal, 821, 126, doi: 10.3847/0004-637X/821/2/126
  * Raymond & Izidoro (2017) Raymond, S. N., & Izidoro, A. 2017, Icarus, 297, 134, doi: 10.1016/j.icarus.2017.06.030
  * Raymond et al. (2007) Raymond, S. N., Quinn, T., & Lunine, J. I. 2007, Astrobiology, 7, 66, doi: 10.1089/ast.2006.06-0126
  * Reufer et al. (2012) Reufer, A., Meier, M. M., Benz, W., & Wieler, R. 2012, Icarus, 221, 296, doi: 10.1016/j.icarus.2012.07.021
  * Robert (2003) Robert, F. 2003, Space Science Reviews, 106, 87, doi: 10.1023/A:1024629402715
  * Rufu et al. (2017) Rufu, R., Aharonson, O., & Perets, H. B. 2017, Nature Geoscience, 10, 89, doi: 10.1038/ngeo2866
  * Schaefer & Fegley (2007) Schaefer, L., & Fegley, B. 2007, Icarus, 186, 462, doi: 10.1016/j.icarus.2006.09.002
  * Schaefer & Fegley (2017) —. 2017, The Astrophysical Journal, 843, 120, doi: 10.3847/1538-4357/aa784f
  * Schaller et al. (2018) Schaller, M., Gonnet, P., Draper, P. W., et al. 2018, Astrophysics Source Code Library, ascl:1805.020
  * Schlichting & Mukhopadhyay (2018) Schlichting, H. E., & Mukhopadhyay, S. 2018, Space Science Reviews, 214, 34, doi: 10.1007/s11214-018-0471-z
  * Schlichting et al. (2015) Schlichting, H. E., Sari, R., & Yalinewich, A. 2015, Icarus, 247, 81, doi: 10.1016/j.icarus.2014.09.053
  * Senft & Stewart (2008) Senft, L. E., & Stewart, S. T. 2008, Meteoritics & Planetary Science, 43, 1993, doi: 10.1111/j.1945-5100.2008.tb00657.x
  * Sossi et al. (2020) Sossi, P. A., Burnham, A. D., Badro, J., et al. 2020, Science advances, 6, eabd1387, doi: 10.1126/sciadv.abd1387
  * Span et al. (2000) Span, R., Lemmon, E. W., Jacobsen, R. T., Wagner, W., & Yokozeki, A. 2000, Journal of Physical and Chemical Reference Data, 29, 1361, doi: 10.1063/1.1349047
  * Span & Wagner (1996) Span, R., & Wagner, W. 1996, Journal of Physical and Chemical Reference Data, 25, 1509, doi: 10.1063/1.555991
  * Stewart et al. (2020) Stewart, S., Davies, E., Duncan, M., et al. 2020, AIP Conference Proceedings, 2272, 080003, doi: 10.1063/12.0000946
  * Stewart (2020) Stewart, S. T. 2020, Zenodo, doi: 10.5281/ZENODO.3866507
  * Stewart et al. (2018) Stewart, S. T., Lock, S. J., & Caracas, R. 2018, in Lunar and Planetary Science Conference Abstracts, Vol. 49, 1708
  * Stewart et al. (2019) Stewart, S. T., Davies, E. J., Duncan, M. S., et al. 2019, Zenodo, doi: 10.5281/ZENODO.3478631
  * Thompson & Lauson (1972) Thompson, S. L., & Lauson, H. S. 1972, Improvements in the CHART D radiation-hydrodynamic code III: Revised analytic equations of state, Tech. rep., Sandia Laboratories
  * Tillotson (1962) Tillotson, J. H. 1962, Report No. GA-3216, General Atomic, San Diego, CA
  * Tucker & Mukhopadhyay (2014) Tucker, J. M., & Mukhopadhyay, S. 2014, Earth and Planetary Science Letters, 393, 254, doi: 10.1016/j.epsl.2014.02.050
  * Wagner (2002) Wagner, W. 2002, Journal of Physical and Chemical Reference Data, 31, 387, doi: 10.1063/1.1461829
  * Wiik (1956) Wiik, H. B. 1956, Geochimica et Cosmochimica Acta, 9, 279, doi: 10.1016/0016-7037(56)90028-X
  * Yalinewich & Schlichting (2018) Yalinewich, A., & Schlichting, H. E. 2018, Monthly Notices of the Royal Astronomical Society, 486, 2780, doi: 10.1093/mnras/stz1049
  * Yoder (1995) Yoder, C. F. 1995, in AGU Reference Shelf 1, ed. T. J. Ahrens (American Geophysical Union), 1–31
  * Young et al. (2019) Young, E. D., Shahar, A., Nimmo, F., et al. 2019, Icarus, 323, 1, doi: 10.1016/j.icarus.2019.01.012
  * Zahnle et al. (2010) Zahnle, K., Schaefer, L., & Fegley, B. 2010, Cold Spring Harbor Perspectives in Biology, 2, a004895, doi: 10.1101/cshperspect.a004895
  * Zel’Dovich & Raizer (2002) Zel’Dovich, Y. B., & Raizer, Y. P. 2002, Physics of shock waves and high-temperature hydrodynamic phenomena (Dover Publications Inc)
  * Ćuk & Stewart (2012) Ćuk, M., & Stewart, S. T. 2012, Science, 338, 1047, doi: 10.1126/science.1225542