11institutetext:  Universität Heidelberg, Zentrum für Astronomie, Institut für
Theoretische Astrophysik, Albert-Ueberle-Straße 2,
D-69120 Heidelberg, Germany 22institutetext: Institut für Astronomie und
Astrophysik, Universität Tübingen, Auf der Morgenstelle 10, 72076 Tübingen,
Germany 33institutetext: Universität Heidelberg, Interdisziplinäres Zentrum
für Wissenschaftliches Rechnen, Im Neuenheimer Feld 205,
D-69120 Heidelberg, Germany

# The rotational disruption of porous dust aggregates from ab-initio kinematic
calculations

Stefan Reissl  The rotational disruption of porous dust aggregates from ab-
initio kinematic calculationsThe rotational disruption of porous dust
aggregates from ab-initio kinematic calculations Philipp Nguyen The rotational
disruption of porous dust aggregates from ab-initio kinematic calculationsThe
rotational disruption of porous dust aggregates from ab-initio kinematic
calculations Lucas M. Jordan The rotational disruption of porous dust
aggregates from ab-initio kinematic calculationsThe rotational disruption of
porous dust aggregates from ab-initio kinematic calculations Ralf S. Klessen
The rotational disruption of porous dust aggregates from ab-initio kinematic
calculationsThe rotational disruption of porous dust aggregates from ab-initio
kinematic calculationsThe rotational disruption of porous dust aggregates from
ab-initio kinematic calculationsThe rotational disruption of porous dust
aggregates from ab-initio kinematic calculations

###### Abstract

Context. The sizes of dust grains in the interstellar medium follows a
distribution where most of the dust mass is in smaller grains. However, the
re-distribution from larger grains towards smaller sizes especially by means
of rotational disruption is poorly understood.

Aims. We aim to study the dynamics of porous grain aggregates under
accelerated ration. Especially, we determine the deformation of the grains and
the maximal angular velocity up to the rotational disruption event by caused
by centrifugal forces.

Methods. We pre-calculate porous grain aggregate my means of ballistic
aggregation analogous to the interstellar dust as input for subsequent
numerical simulations. In detail, we perform three-dimensional N-body
simulations mimicking the radiative torque spin-up process up to the point
where the grain aggregates become rotationally disrupted.

Results. Our simulations results are in agreement with theoretical models
predicting a characteristic angular velocity $\omega_{\mathrm{disr}}$ of the
order of ${10^{8}-10^{9}\ \mathrm{rad\ s^{-1}}}$, where grains become
rotationally disrupted. In contrast to the theoretical predictions, we show
that for large porous grain aggregates ($\gtrapprox 300\ \mathrm{nm}$)
$\omega_{\mathrm{disr}}$ does not strictly decline but reaches a lower
asymptotic value. Hence, such grains can withstand an accelerated ration more
efficiently up to a factor of 10 because the displacement of mass by
centrifugal forces and the subsequent mechanical deformation supports the
build up of new connections within the aggregate. Furthermore, we report that
the rapid rotation of grains deforms an ensemble with initially 50:50 prolate
and oblate shapes, respectively, preferentially into oblate shapes. Finally,
we present a best fit formula to predict the average rotational disruption of
an ensemble of porous dust aggregates dependent on internal grain structure,
total number of monomers, and applied material properties.

###### Key Words.:

…, …, ….

##  1 Introduction

Dust is a key component of the interstellar medium (ISM). It is important for
regulating the properties of astrophysical objects across a wide range of
scales: from the cooling of collapsing molecular clouds and subsequent star-
formation down to the formation of planetary systems (Spitzer & Arny 1978;
Dorschner & Henning 1995). However, the origin of dust, its initial physical
properties and the redistribution of grain sizes is still a field of ongoing
research (O’Donnell & Mathis 1997; Ormel et al. 2009; Birnstiel et al. 2010;
Guillet et al. 2018; Draine & Hensley 2021a).
Dust composition and grain size distribution may be derived from the observed
interstellar extinction curve and starlight polarization (Mathis et al. 1977;
Draine & Lee 1984; Guillet et al. 2018; Draine & Hensley 2021a). Initially,
the ISM is enriched by intermediate-mass stars at the asymptotic giant branch
(AGB) supernova ejecta coming with a certain grain size distribution (Nozawa
et al. 2007; Gail et al. 2009; Barlow et al. 2010; Matsuura 2011; Karovicova
et al. 2013; Zhukovska et al. 2015; Bevan & Barlow 2016). Later, the grains
grow in dense molecular clouds by accretion of abundant elements and
coagulation (Spitzer & Arny 1978; Chokshi et al. 1993; O’Donnell & Mathis
1997). Naturally, this process results in porous dust aggregates rather than
solid bodies (Ossenkopf 1993; Dominik & Tielens 1997; Wada et al. 2007). Dust
destruction processes such as gas-grain sputtering or grain-grain collision
(shattering) may redistribute the grown grains towards smaller sizes (Dwek &
Scalo 1980; Tielens et al. 1994; Hirashita & Yan 2009).
More recently, the pioneering work presented by Hoang et al. (2019) describes
a new dust destruction mechanism. Here, a radiation field causes radiative
torques (RAT) acting on dust grains which lead to an angular acceleration(see
e.g. Lazarian & Hoang 2007a; Hoang et al. 2014). Given a sufficiently luminous
environment the grains would inevitably be disrupted by the emerging
centrifugal force (Silsbee & Draine 2016; Hoang et al. 2019; Hoang 2020). The
disruption process is usually quantified by the maximal tensile strength
$\mathcal{S}_{\mathrm{max}}$ which is a measure how materials respond to
stretching (Silsbee & Draine 2016; Tatsuuma & Kataoka 2021). For porous
materials and dust grain analogs the tensile strengths may be determined by
numerical N-body simulations (Kataoka et al. 2013b; Seizinger et al. 2013b;
Tatsuuma et al. 2019). However, what simulations of
$\mathcal{S}_{\mathrm{max}}$ are missing is that individual building blocks
(so called monomers) do not just feel the stretching force between its
neighbours, but additionally the global centrifugal force acting on the entire
aggregate. Hence, up to this point, it remains unclear if the parameter
$\mathcal{S}_{\mathrm{max}}$ describes accurately the internal processes of
material displacement within the rotating aggregate.
The rotational disruption of fractal grain aggregates was already studied
indirectly in (see Reissl et al. 2022, RMK22 hereafter) by evaluating
$\mathcal{S}_{\mathrm{max}}$ in the context of rapid grain ration caused by a
differential gas-dust velocity. However, in this paper we aim to simulate the
dynamics of rapidly rotating interstellar grains directly by taking the time
evolution of the internal aggregate structure into account. The aim is to
develop a model of the average disruption process of large ensembles of porous
grains. This paper is structured as follows: In Sect. 2 we discuss the most
likely composition of elements and minerals to be present in dust aggregates.
An algorithm to mimic the growth of porous dust aggregates is outlined in
Sect. 3 in detail. Here, we also introduce the methods to quantify the grain
shape and porosity. In Sect. 4 we discuss the processes and forces acting
between the monomers connected within the aggregates. The spin-up process of
grains by means of RATs is outlined in Sect. 5. In Sect. 6 we discus the
numerical implementation of the set of equations that governs the internal
aggregate dynamics under rapid rotation. In Sect. 7 we present and discuss our
N-body simulation results. Finally, in Sect. 8 we summarize our findings.

$\mathrm{Mg/O}$ | $\mathrm{Si/O}$ | $\mathrm{Fe/O}$ | $\mathrm{Mg/Fe}$ | $\mathrm{(Mg+Fe)/Si}$ | References
---|---|---|---|---|---
0.31 | 0.26 | 0.15 | 2.03 | 1.82 | Min et al. (2007)
X | X | X | 1.09 | 2.25 | Voshchinnikov & Henning (2010)
0.23 | 0.19 | 0.19 | 1.25 | 2.25 | Compiègne et al. (2011)
0.19 | 0.15 | 0.17 | 1.07 | 2.34 | Hensley & Draine (2021)
0.23 | 0.27 | 0.21 | 1.08 | 1.60 | this work
Table 1:  Ratio of different element abundances as observed in the ISM in comparison with the composition of the co-S material applied in this work.  | $a_{\mathrm{eff}}=200\ \mathrm{nm}$ | $a_{\mathrm{eff}}=350\ \mathrm{nm}$ | $a_{\mathrm{eff}}=500\ \mathrm{nm}$
---|---|---|---
$\mathrm{BA}$ |  |  |
| $a_{\mathrm{out}}=543.1\ \mathrm{nm}$ | $a_{\mathrm{out}}=882.6\ \mathrm{nm}$ | $a_{\mathrm{out}}=1385.7\ \mathrm{nm}$
|  $N_{\mathrm{mon}}=160$, $N_{\mathrm{con}}=162$ |  $N_{\mathrm{mon}}=521$, $N_{\mathrm{con}}=521$ |  $N_{\mathrm{mon}}=1773$, $N_{\mathrm{con}}=1772$
$\mathrm{BAM1}$ |  |  |
| $a_{\mathrm{out}}=457.4\ \mathrm{nm}$ | $a_{\mathrm{out}}=901.2\ \mathrm{nm}$ | $a_{\mathrm{out}}=1214.7\ \mathrm{nm}$
|  $N_{\mathrm{mon}}=117$, $N_{\mathrm{con}}=235$ |  $N_{\mathrm{mon}}=636$, $N_{\mathrm{con}}=1275$ |  $N_{\mathrm{mon}}=1885$, $N_{\mathrm{con}}=3774$
$\mathrm{BAM2}$ |  |  |
| $a_{\mathrm{out}}=320.1\ \mathrm{nm}$ | $a_{\mathrm{out}}=603.3\ \mathrm{nm}$ | $a_{\mathrm{out}}=845.2\ \mathrm{nm}$
|  $N_{\mathrm{mon}}=157$, $N_{\mathrm{con}}=475$ |  $N_{\mathrm{mon}}=649$, $N_{\mathrm{con}}=1950$ |  $N_{\mathrm{mon}}=1814$, $N_{\mathrm{con}}=5546$
Figure 1: Exemplary selection from the total ensemble of porous BA (top row),
BAM1 (moddle row), and BAM2 (bottom row) aggregates for the effective grain
radii $a_{\mathrm{eff}}=200\ \mathrm{nm}$ (left column),
$a_{\mathrm{eff}}=350\ \mathrm{nm}$ (middle column), and
$a_{\mathrm{eff}}=500\ \mathrm{nm}$ (right column) with the corresponding
numbers of monomers $N_{\mathrm{mon}}$ and neighbourhood connections
$N_{\mathrm{con}}$. Monomers are sampled to guarantee an exact effective
radius $a_{\mathrm{eff}}$. The radius $a_{\mathrm{out}}$ is associated to the
smallest sphere enclosing the entire aggregate. The axes
$\hat{a}_{\mathrm{1}}$, $\hat{a}_{\mathrm{2}}$, and $\hat{a}_{\mathrm{3}}$ are
defined by the grain’s moments of inertia
${I_{\mathrm{a1}}>I_{\mathrm{a2}}>I_{\mathrm{a3}}}$, where
$\hat{a}_{\mathrm{1}}$ is the designated axis of grain rotation.

Figure 2: The distribution of the characteristic quantities of total number of
monomers $N_{\mathrm{mon}}$ (top left), number of connected neighbours
$N_{\mathrm{con}}$ (top right), porosity $\mathcal{P}$ (bottom left), and
fractal dimension $D_{\mathrm{f}}$ (bottom right), for all BA (red), BAM1
(green), and BAM2 (blue) aggregates, respectively, dependent on effective
radius $a_{\mathrm{eff}}$. Dots are the ensemble average over all shapes and
materials while vertical bars represent the minima and maxima. Note that the
data points have a small offset for better visibility.

##  2 Dust grain composition

The exact composition of dust remains an open question, with a large variety
of observed materials and sizes within the galaxy being similarly possible.
Interstellar dust is usually modelled by grains with a spectrum of different
sizes with silicate and carbonaceous components (Mathis et al. 1977;
Weingartner & Draine 2001; Zhukovska et al. 2008; Voshchinnikov & Henning
2010; Guillet et al. 2018). Small spherical monomers condensate by depletion
of the most abundant elements C, Mg, Si, Fe, and O from their immediate
surrounding (see e.g. Kim et al. 2021). Elements such Ti, Al, S Ca, Ni,
respectively, are less abundant and may only contribute a few percent of the
total dust mass (Hensley & Draine 2021). The abundance of molecules within the
grains may be determined by observing the spectral dust absorption features.
Silicate minerals from the olivine and pyroxene group are the most likely
candidates to be present in order to account for the observed characteristic
features (Mathis et al. 1977; Wada et al. 1999; Li & Draine 2001; Hensley &
Draine 2021). However, it remains yet inconclusive if these minerals form dust
in a crystalline or amorphous structure and to what extend iron is present in
its pure form (Draine & Li 2007; Rogantini et al. 2019; Do-Duy et al. 2020).
Carbonaceous grains may consist of a regular graphite lattice or amorphous
structures and be partly hydrogenated (Wada et al. 1999; Goto et al. 2003;
Mennella 2006). Possibly, carbonaceous and silicate grains are not even
separate dust populations but baked into a single composite material (Draine &
Hensley 2021a).
For mimicking the optical properties of dust grains models based on the
refractive indices of various materials are developed and the data is publicly
available111For the interested reader we refer to the website e.g. of Bruce
Draine https://www.astro.princeton.edu/ draine/dust/dust.diel.html and the
THEMIS model https://www.ias.u-psud.fr/themis/ (see e.g. Draine & Li 2007;
Jones 2012; Draine & Hensley 2021a). However, complementary well constrained
models of the mechanical properties particularly for the composition of
interstellar dust grain materials are still missing. Commonly, grain analogs
consisting exclusively of icy or pure quartz ($\mathrm{SiO_{2}}$) materials
are utilized to simulate grain growth processes (e.g. Dominik & Tielens 1997;
Wada et al. 2007; Seizinger et al. 2012). This is simply due to the fact that
quartz is an easily available material on the market with well constrained
properties by laboratory experiments (Kendall et al. 1987; Heim et al. 1999;
Israelachvili 2011; Krijt et al. 2013).
In this paper, however, we explore the rotational disruption of three distinct
grain materials labeled a-C, q-S, and co-S, respectively. Carbonaceous grains
are represented by the mechanical parameters of amorphous carbon (a-C) and for
compression silicate grains are considered built of pure quartz (q-S). In
addition, we aim to approximate the mechanical properties of composite
silicate (co-S) grains more precisely. Here, we assume the minerals forsterite
($\mathrm{Mg_{2}SiO_{4}}$) and fayalite ($\mathrm{Fe_{2}SiO_{4}}$) of the
olivine series as well as enstatite ($\mathrm{MgSiO_{3}}$) and ferrosilite
($\mathrm{FeSiO_{3}}$) of the pyroxene series to be among the most likely
ingredients of silicate grains (Petrovic 2001; Zolensky et al. 2006; Zhukovska
et al. 2008; Gail et al. 2009; Takigawa & Tachibana 2012; Min et al. 2007;
Kimura et al. 2015; Fogerty et al. 2016; Hoang et al. 2019; Escatllar et al.
2019; Kimura et al. 2020; Hensley & Draine 2021; Draine & Hensley 2021a). In
order to get an approximation of the material mixture, we match the abundance
of individual elements within the considered minerals with the abundance of
elements typical for the ISM (Min et al. 2007; Voshchinnikov & Henning 2010;
Compiègne et al. 2011; Hensley & Draine 2021).
In Table 1 we present the relative abundances of elements within the ISM in
comparison with the composition of our co-S model. For our best fit co-S model
we get that each individual silicate monomer consist of a mixture of $31\ \%$
forsterite, $29\ \%$ fayalite, $20\ \%$ enstatite, $20\ \%$ ferrosilite,
respectively, and average the mechanical material properties accordingly. The
abundance of elements within the co-S agrees well with the overall
observations within the Milky Way but our co-S model shows a slight
overabundance of Si. Naturally, the composition of co-S may locally be vastly
different e.g. in the vicinity of oxygen or carbon rich AGB stars where dust
grains are newly formed Zhukovska & Henning (2013). Thus, we remain agnostic
concerning the actual composition of interstellar dust but consider our co-S
model to be an improvement compared to simulations using pure quartz monomers.

##  3 Dust grain growth and aggregation

Dust grain aggregates may grow by ballistic hit-and-stick processes of
monomers onto a grain’s surface. This process is usually called ballistic
particle-cluster aggregation (BCPA) (Kozasa et al. 1992; Bertini et al. 2007).
Newly formed grains may then grow to even larger aggregates by ballistic
cluster-cluster aggregation (BCCA) Ossenkopf (1993). Such grain-grain
collisions result significantly compressed aggregates and, subsequently, gas
pressure may compress an aggregate even more (see e.g Dominik & Tielens 1997;
Kataoka et al. 2013a; Michoulier & Gonzalez 2022, and references therein).
In our study, we model such compression effects by the ballistic aggregation
with migration (BAM) model introduced in Shen et al. (2008). Here, grains
simply grow by means of ballistic aggregation (BA) 222BA and BCPA are used
synonymous in literature. of monomers hitting the aggregate from random
directions. In the BAM model grain model monomers may migrate along the
surface once (BAM1) or twice (BAM2). Consequently, for BAM1 each monomer has
at least two connections with the aggregate and for BAM2 each monomer has at
least three Shen et al. (2008); Seizinger et al. (2013a).
Commonly, dust aggregation models utilize only a constant monomer size (see
e.g. Kozasa et al. 1992; Shen et al. 2008; Wada et al. 2007; Bertini et al.
2007; Seizinger et al. 2012). However, it seems unlikely that a condensation
process of elements in nature would lead to exactly one monomer size. In fact,
numerous laboratory experiments clearly indicate that grown aggregates consist
of monomers with a variable radius (Karasev et al. 2004; Chakrabarty et al.
2007; Slobodrian et al. 2011; Kandilian et al. 2015; Salameh et al. 2017; Paul
et al. 2017; Baric et al. 2018; Kelesidis et al. 2018; Bauer et al. 2019; Wu
et al. 2020; Zhang et al. 2020; Kim et al. 2021). In this study we focus on a
polydisperse system of monomers, where the exact distribution of the monomer
radii within an aggregate may be approximated by a log-normal distribution
(Köylü & Faeth 1994; Lehre et al. 2003; Slobodrian et al. 2011; Bescond et al.
2014; Kandilian et al. 2015; Liu et al. 2015; Bauer et al. 2019; Wu et al.
2020; Zhang et al. 2020).
We realize the BAM model with a Monte-Carlo approach in order to create an
ensemble of pre-calculated grain analogs resembling the observed parameters of
dust in the circumstellar and interstellar medium. The radius of the i-th
monomer $a_{\mathrm{mon,i}}$ is sampled from a range of
$a_{\mathrm{mon,i}}\in[10\ \mathrm{nm},100\ \mathrm{nm}]$. The log-normal
distribution has a typical average of ${\left<a_{\mathrm{mon}}\right>=20\
\mathrm{nm}}$ and a standard deviation of $1.65\ \mathrm{nm}$. Successively,
each newly sampled monomer is shot on a random trajectory into the simulation
domain until it hits the aggregate. For simplicity we assume that each monomer
sticks onto the surface when colliding. For BAM1 and BAM2, respectively, the
i-th monomer migrates along the the surface of the initially hit monomer in a
random direction to establish additional connections. In order to create an
aggregate in equilibrium the monomer position $\@vec{X}_{\mathrm{i}}$ is
corrected in such a way that each overlap between connected monomers agrees
with the material dependent equilibrium compression length
$\delta_{\mathrm{0}}$ (see Sect. 4 for details). The aggregation process is
repeated until the dust aggregate reaches a certain volume of

$V_{\mathrm{agg}}=\frac{4\pi}{3}\sum_{\mathrm{i}=1}^{N_{\mathrm{mon}}}a_{\mathrm{mon,i}}^{3}\,.$
(1)

Ballistic aggregates are usually quantified by the total number of monomers
$N_{\mathrm{mon}}$ (see e.g. Wada et al. 2007; Shen et al. 2008; Seizinger et
al. 2012). However, dust observations are tightly connected to the effective
size of the grains (Mathis et al. 1977; Weingartner & Draine 2001). Hence, in
our study we rather control for an exact effective radius of

$a_{\mathrm{eff}}=\left(\frac{3V_{\mathrm{agg}}}{4\pi}\right)^{\frac{1}{3}}$
(2)

by an biased sampling of the last three monomer radii instead of getting the
grain size indirectly from ${a_{\mathrm{eff}}\approx
N_{\mathrm{mon}}^{1/3}\left<a_{\mathrm{mon}}\right>}$. Finally, we calculate
the inertia tensor for each aggregate in order to determine the characteristic
moments of inertia ${I_{\mathrm{a1}}>I_{\mathrm{a2}}>I_{\mathrm{a3}}}$, along
the unit vectors $\hat{a}_{\mathrm{1}}$, $\hat{a}_{\mathrm{2}}$, and
$\hat{a}_{\mathrm{3}}$ (we refer t to RMK22 for the exact procedure).
In order to connect the rotational disruption of the grain ensemble to
distinct quantities associated with the internal structure of fluffy dust
grains we introduce the porosity $\mathcal{P}$, volume filling factor $\phi$,
and fractal dimension $D_{\mathrm{f}}$ as well as the semi major axes $a<b<c$
unique for each individual aggregate.
The porosity quantifies the empty space within an aggregate where
$\mathcal{P}=0$ is for a solid object and $\mathcal{P}=1$ for vacuum. Its
value for an individual object depends on the definition of the surface that
envelopes the aggregate. In our study we follow the procedure of determining
$\mathcal{P}$ by utilizing the moments of inertia of an aggregate as outlined
in Shen et al. (2008). Here, the quantity

$\alpha_{\mathrm{i}}=\frac{5}{4\pi}\frac{I_{\mathrm{ai}}}{\rho_{\mathrm{mat}}a_{\mathrm{eff}}^{5}}$
(3)

is the ratio of the moment of inertia $I_{\mathrm{ai}}$ to that of an sphere
with equivalent volume where the index ${\mathrm{i}\in\\{1,2,3\\}}$ denotes
the three spatial directions. The corresponding semi major axes are

$a=a_{\mathrm{eff}}\sqrt{\alpha_{\mathrm{2}}+\alpha_{\mathrm{3}}-\alpha_{\mathrm{1}}}\,,$
(4)
$b=a_{\mathrm{eff}}\sqrt{\alpha_{\mathrm{3}}+\alpha_{\mathrm{1}}-\alpha_{\mathrm{2}}}\,,$
(5)

and

$c=a_{\mathrm{eff}}\sqrt{\alpha_{\mathrm{1}}+\alpha_{\mathrm{2}}-\alpha_{\mathrm{3}}}\,,$
(6)

respectively. Finally, the porosity of an aggregate may then be written as

$\mathcal{P}=1-\frac{a_{\mathrm{eff}}^{3}}{abc}\,,$ (7)

whereas the complementary quantity $\phi=1-\mathcal{P}$ is the volume filling
factor Shen et al. (2008).
The fractal dimension is a measure of the shape of an aggregate where
$D_{\mathrm{f}}=1$ represents a one dimensional line and $D_{\mathrm{f}}=3$ a
compact sphere. We determine the fractal dimension $D_{\mathrm{f}}$ of each
dust aggregate by the correlation function

$C\chi^{D_{\mathrm{f}}-3}=\frac{n(\chi)}{4\pi\chi^{2}lN_{\mathrm{mon}}}$ (8)

as outlined in Skorupski et al. (2014). Here, $C$ is a scaling factor, $\chi$
is the distance from the center of mass, $l$ is a length with $l\ll
a_{\mathrm{eff}}$, and $n(\chi)$ is the number density of connected monomers
within the shell ${[\chi-l/2;\chi+l/2]}$.
We create dust grains with effective radii in the range ${a_{\mathrm{eff}}=50\
\mathrm{nm}-550\ \mathrm{nm}}$ in steps of $50\ \mathrm{nm}$. For each
$a_{\mathrm{eff}}$ we repeat the MC dust growth process with 30 random seeds
for the BA, BAM1, and BAM2 configurations and the a-C, q-S, co-S materials. In
total we pre-calculate an ensemble of 2970 individual grains as input for our
three-dimensional N-body simulations.
We note that the grain growth by BA may heavily be impacted by grain charge
Matthews et al. (2012), high impact velocities (Dominik & Tielens 1997; Ormel
et al. 2009), a preferential impact direction by means of grain alignment with
the magnetic field (Lazarian & Hoang 2007a; Hoang 2022) or a gas-dust drift
(RMK22). Hence, the resulting shapes in our grain ensemble may not be
representative when compared for grain growth processes e.g. in the vicinity
of AGB stars (Zhukovska & Henning 2013) or in protostellar envelopes (Galametz
et al. 2019). However, our grain model cover a sufficiently broad variety of
grain shapes to allow conclusions about their rotational stability even though
each individual shape is not equally likely to be realize in nature.
An exemplary selection of the grains is shown in Fig. 1. In Fig. 2 we present
the characteristic quantities of the entire grain ensemble as introduced
above. The number of monomers $N_{\mathrm{mon}}$ as well as the number of
connections $N_{\mathrm{con}}$ within each aggregate increase with effective
radius $a_{\mathrm{eff}}$. By design BA grains have generally less connections
compared to BAM grains. Compared to the BAM grains with a fixed monomer size
(Shen et al. 2008) the porosity $\mathcal{P}$ is not strictly increasing but
stagnates for higher $a_{\mathrm{eff}}$ because smaller monomers may easily
migrate towards the center making our aggregates overall more compact. The
fractal dimension $D_{\mathrm{f}}$ shown in Fig. 2 increases slightly towards
larger grains but the resulting grain BA, BAM1, and BAM2 grains of different
sizes shapes in this work are not well correlated with the fractal dimension
$D_{\mathrm{f}}$. This is in contrast to the dust models of RMK22 where the
grains are explicitly constructed to get an exact pre-determined
$D_{\mathrm{f}}$ and are not to be compared with the BAM grains as depicted in
Fig. 1.

##  4 Inter-monomer contact effects

Figure 3: Schematic representation of the forces and torques acting in between the i-th and the j-th monomer at the positions $\@vec{X}_{\mathrm{i}}$ and $\@vec{X}_{\mathrm{j}}$, respectively. The common contact surface with radius $r_{\mathrm{ij}}$ is shaded in red. (a) Each individual monomer experiences external centrifugal $\@vec{F}_{\mathrm{cent}}$, Coriolis $\@vec{F}_{\mathrm{cor}}$, and Euler forces $\@vec{F}_{\mathrm{eul}}$ as a result of the aggregate’s accelerated rotation with angular velocity $\@vec{\omega}_{\mathrm{agg}}$. (b) The normal force $\@vec{F}_{\mathrm{N,ij}}$ acts normal to the contact surface of two monomers because of surface attraction and mechanical deformation with the contact pointers remaining anti-parellel ${\@vec{n}_{\mathrm{i}}=-\@vec{n}_{\mathrm{j}}}$. (c) The sliding force $\@vec{F}_{\mathrm{S,ij}}$ is parallel to the sliding displacement $\zeta$. Both $\@vec{F}_{\mathrm{S,ij}}$ and the corresponding sliding torque $\@vec{\Gamma}_{\mathrm{S,ij}}$ work tangential to the contact surface. (d) The same for the rolling of monomers withing the aggregate: The rolling displacement $\@vec{\zeta}$, the force $\@vec{F}_{\mathrm{R,ij}}$, and the $\@vec{\Gamma}_{\mathrm{R,ij}}$ are tangential to the contact surface. (e) The twisting of monomers in contact is not associated with a net force. The resulting torque $\@vec{\Gamma}_{\mathrm{T,ij}}$ points is parallel to the normal vector $\@vec{\phi}$ of the twisting. We note that the depicted quantities are not to scale since $r_{\mathrm{ij}}\ll a_{\mathrm{mon}}$. | carbon: | silicate: | silicate: | forsterite | fayalite | enstatite | ferrosilite
---|---|---|---|---|---|---|---
| amorphous (a-C) | quartz (q-S) | composite (co-S) | $31\ \%$ | $29\ \%$ | $20\ \%$ | $20\ \%$
$\gamma\ [\mathrm{mJ\ m}^{-2}]$ | $50_{\ (11)}$ | $20_{\ (1)}$ | 70 | $70_{\ (4)}$ | $70_{\ (4)}$ | $70_{\ (4)}$ | $70_{\ (4)}$
$E\ [\mathrm{GPa}]$ | $168_{\ (9,10)}$ | $54_{\ (1,3)}$ | 169 | $204_{\ (6,8)}$ | $140_{\ (6)}$ | $180_{\ (6)}$ | $142_{\ (7)}$
$\nu$ | $0.16_{\ (9)}$ | $0.17_{\ (1,3)}$ | 0.27 | $0.24_{\ (2,6,8)}$ | $0.32_{\ (2,6)}$ | $0.21_{\ (2,6)}$ | $0.31_{\ (7)}$
$\rho_{\mathrm{mat}}\ [\mathrm{kg\ m}^{-3}]$ | $2368_{\ (9,11)}$ | $2650_{\ (1,3)}$ | 3707 | $3213_{\ (2,5,14)}$ | $4393_{\ (2,5,14)}$ | $3209_{\ (2,5,14)}$ | $4014_{\ (5,7,14)}$
Table 2: Material parameters of the surface energy $\gamma$, Young’s modulus
$E$, Poisson number $\nu$, and material density $\rho_{\mathrm{mat}}$ for the
monomer materials of amorphous carbon (a-C), pure quartz (q-S), and composite
silicate (co-S) considered in our N-body simulations. The co-S material is
assumed to consist of a mixture of the different minerals of forsterite,
fayalite, enstatite, and ferrosilite, respectively. For all the materials we
assume a critical rolling displacement of $\xi_{\mathrm{crit}}=0.2\
\mathrm{nm}_{\ (1,3)}$ and a viscous dumping time of $T_{\mathrm{vis}}=5\
\mathrm{ps}_{\ (1,12,13)}$.
References: (1) Seizinger et al. (2012), (2) Christensen (1996), (3) Dominik &
Tielens (1997), (4) Bogdan et al. (2020), (5) Williams & Jadwick (1980), (6)
Petrovic (2001), (7) Mitchell (2004), (8) Gouriet et al. (2019), (9) Jensen et
al. (2015), (10) Remediakis et al. (2007), (11) Zebda et al. (2008), (12)
Krijt et al. (2013), (13) Seizinger et al. (2013a), (14) Cardarelli (2008)

In this section we outline the forces and torques acting between the monomers
of an aggregate in detail. In Fig. 3 we provide a schematic illustration of
all considered monomer interactions. Two individual monomers in physical
contact establish a common contact surface area and experience an attraction
because of the van der Waals force. For some materials stronger attractions
such as Coulomb forces between charged monomers, dipole-dipole interaction
within ices, or metallic binding between iron pallets may become of relevance.
The attraction is quantified by the material dependent energy per surface area
$\gamma$. Assuming monomers act like elastic spheres with an radius of
$a_{\mathrm{mon,i}}$ and $a_{\mathrm{mon,j}}$, respectively, their elastic
deformation causes an repulsive force. An analytical description of these
forces was first presented in Johnson et al. (1971) with the so called JKR
model (see also Hertz 1896) where the equilibrium radius of the contact
surface is

$r_{\mathrm{0}}=\left(\frac{9\pi\gamma R^{2}}{E^{*}}\right)^{1/3}$ (9)

when no external force are acting between monomers i.e. the attractive and
repulsive forces are in balance. Here, the quantity
${R=a_{\mathrm{mon,i}}a_{\mathrm{mon,j}}/(a_{\mathrm{mon,i}}+a_{\mathrm{mon,j}})}$
is the reduced monomer radius whereas the elastic parameter
${E^{*}=E/(2-2\nu^{2})}$ is determined by the material specific constants of
the Poisson number $\nu$ and Young’s modulus $E$, respectively (we refer to
Johnson 1987, for further details). Utilizing that a monomer contact breaks
for the characteristic pulling force of

$F_{\mathrm{C}}=3\pi\gamma R\,,$ (10)

as outlined in (Johnson et al. 1971), the normal force along the unit vector

$\@vec{n}_{\mathrm{c}}=\frac{\@vec{X}_{\mathrm{i}}-\@vec{X}_{\mathrm{j}}}{\left|\@vec{X}_{\mathrm{i}}-\@vec{X}_{\mathrm{j}}\right|}\,.$
(11)

may be written as

$\@vec{F}_{\mathrm{N,ij}}=4F_{\mathrm{C}}\left[\left(\frac{r}{r_{\mathrm{0}}}\right)^{3}-\left(\frac{r}{r_{\mathrm{0}}}\right)^{3/2}\right]\@vec{n}_{\mathrm{c}}$
(12)

Consequently, for $F_{\mathrm{N,ij}}>0$ the i-th monomer is exerting a pushing
force onto its j-th neighbour. Otherwise, for $F_{\mathrm{N,ij}}<0$, the j-th
monomer is pulling on the i-th monomer.
Later, the JKR model was extended by Dominik & Tielens (1997) considering the
mechanics of rolling (Dominik & Tielens 1995, 1997), sliding (Dominik &
Tielens 1996, 1997) , and twisting motions (Dominik & Tielens 1997) in between
connected monomers. In order to track the relative motion of monomers over we
use the formulation of the contact pointers $\@vec{n}_{\mathrm{i}}$ and
$\@vec{n}_{\mathrm{j}}$ as outlined in Dominik & Nübold (2002). These vectors
point initially towards the centers of neighboring monomers when a new contact
is established (see Fig. 3).
The force acting on the i-th monomer because of the sliding motion of the j-th
monomer may the be written as

$\@vec{F}_{\mathrm{S,ij}}=-8r_{\mathrm{0}}G^{*}\@vec{\zeta}\frac{\left(a_{\mathrm{mon,j}}\@vec{n}_{\mathrm{j}}+a_{\mathrm{mon,i}}\@vec{n}_{\mathrm{i}}\right)\@vec{n}_{\mathrm{c}}}{\left|\@vec{X}_{\mathrm{i}}-\@vec{X}_{\mathrm{j}}\right|}$
(13)

where

$\@vec{\zeta}=a_{\mathrm{mon,i}}\@vec{n}_{\mathrm{i}}+a_{\mathrm{mon,j}}\@vec{n}_{\mathrm{j}}-\left(a_{\mathrm{mon,i}}\@vec{n}_{\mathrm{i}}\@vec{n}_{\mathrm{c}}-a_{\mathrm{mon,j}}\@vec{n}_{\mathrm{c}}\@vec{n}_{\mathrm{i}}\right)\@vec{n}_{\mathrm{c}}$
(14)

is the sliding displacement and

$\@vec{\Gamma}_{\mathrm{S,ij}}=-8r_{\mathrm{0}}G^{*}a_{\mathrm{mon,i}}\@vec{n}_{\mathrm{i}}\times\@vec{\zeta}$
(15)

is the associated sliding torque (Johnson 1987; Dominik & Tielens 1996, 1997).
Here, the material constant ${G^{*}=G/(2-2\nu_{\mathrm{i}}^{2})}$ depends on
the shear modulus ${G=E/(2+2\nu)}$ (see e.g. Cardarelli 2008).
In contrast to sliding a rolling motion does not result in a force (see e.g.
Dominik & Tielens 1997; Wada et al. 2007) but only in a torque of

$\@vec{\Gamma}_{\mathrm{R,ij}}=-4F_{\mathrm{C}}\@vec{n}_{\mathrm{i}}\times\@vec{\xi}\,,$
(16)

where the rolling displacement depends on the contact pointers via
${\@vec{\xi}=R(\@vec{n}_{\mathrm{i}}+\@vec{n}_{\mathrm{j}})}$.
The same for the twisting between monomers with a torque of

$\@vec{\Gamma}_{\mathrm{T,ij}}=\frac{16}{3}G^{*}r_{0}^{3}\@vec{\phi}\,.$ (17)

whereas the motion of the twisting follows the direction of the vector

$\@vec{\phi}=\@vec{n}_{\mathrm{c}}(t)\int_{0}^{t}{\left(\@vec{\omega}_{\mathrm{i}}(t^{\prime})-\@vec{\omega}_{\mathrm{j}}(t^{\prime})\right)\@vec{n}_{\mathrm{c}}(t^{\prime})\mathrm{d}t^{\prime}}\,.$
(18)

Here, the angular velocities $\@vec{\omega}_{\mathrm{i}}(t)$ and
$\@vec{\omega}_{\mathrm{j}}(t)$, respectively, describe the relative rotation
of two monomers in contact (Dominik & Tielens 1997).
In this study we extend the monomer contact physics pioneered by Johnson et
al. (1971), Dominik & Tielens (1995, 1996, 1997), and Wada et al. (2007) by
introducing additional forces emerging from the accelerated rotation of the
grain aggregate. In the notation of our grain model the centrifugal force
acting on each monomer may be write as

$\@vec{F}_{\mathrm{cent,i}}=-m_{\mathrm{mon,i}}\
\@vec{\omega}_{\mathrm{agg}}\times\left(\@vec{\omega}_{\mathrm{agg}}\times\@vec{X}_{\mathrm{i}}\right)$
(19)

where ${m_{\mathrm{mon,i}}=4\pi/3\ \rho_{\mathrm{mat}}a_{\mathrm{mon,i}}^{3}}$
is the mass of the i-th monomer.
We emphasize that our numerical setup operates in a co-rotating coordinate
system with its origin coinciding with the center of mass of the most massive
fragment. Hence, the Coriolis effect acts as an additional fictive force on
each individual monomers via

$\@vec{F}_{\mathrm{cor,i}}=-2m_{\mathrm{mon,i}}\
\@vec{\omega}_{\mathrm{agg}}\times\frac{\mathrm{d}\@vec{X}_{\mathrm{i}}}{\mathrm{d}t}\,.$
(20)

Furthermore, our aggregates are initially at rest and are gradually spun-up.
Consequently, the acceleration of each monomer within the aggregate leads to
an Euler force of

$\@vec{F}_{\mathrm{eul,i}}=-m_{\mathrm{mon,i}}\frac{\mathrm{d}\@vec{\omega}_{\mathrm{agg}}}{\mathrm{d}t}\times\@vec{X}_{\mathrm{i}}$
(21)

acting on each monomer where
${\mathrm{d}\@vec{\omega}_{\mathrm{agg}}/\mathrm{d}t}$ is the angular
acceleration.
When monomers establish a new contact they start to oscillate around the
equilibrium position driven by the surface attraction and monomer deformation.
In nature it is expected that the oscillation becomes dampened because the
deformation dissipates energy. A weak damping force may be introduced based on
the relative velocity of the monomers (Kataoka et al. 2013b; Seizinger et al.
2012) to artificially reduce oscillations. However, as shown in Seizinger et
al. (2012), such a force may come with numerical instabilities. However, the
damping model presented in Krijt et al. (2013) shows that the elastic
dampening can significantly increase the dissipation of energy. In this study
we apply the viscoelastic damping force

$\@vec{F}_{\mathrm{D,ij}}=\frac{2E^{*}}{\nu_{\mathrm{i}}^{2}}\frac{\mathrm{d}\delta}{\mathrm{d}t}r_{\mathrm{ij}}T_{\mathrm{vis}}\@vec{n}_{\mathrm{c}}\,,$
(22)

as introduced by Seizinger et al. (2013a). Here, the dampening depends not on
the relative velocity, but instead on the time evolution of the compression
lengths

$\delta=a_{\mathrm{i}}+a_{\mathrm{j}}-\left|\@vec{X}_{\textit{i}}-\@vec{X}_{\textit{j}}\right|\,,$
(23)

i.e. the overlap between two spherical monomers with a distance of
$|\@vec{X}_{\textit{i}}-\@vec{X}_{\textit{j}}|$ apart from each other. For an
aggregate in equilibrium the compression length is
${\delta_{\mathrm{0}}=r_{\mathrm{0}}/(3R)}$. The characteristic viscoelastic
timescale $T_{\mathrm{vis}}$ in the order of ${1\ \mathrm{ps}-10\
\mathrm{ps}}$ (Krijt et al. 2013) but poorly constrained for specific grain
materials. For smaller values of $T_{\mathrm{vis}}$ the material effectively
behaves elastically and no damping is expected. For simplicity, we adapt a
fixed value of $T_{\mathrm{vis}}=5\ \mathrm{ps}$ for all considered grains
independent of material.
The full set of material parameters applied in our simulations is listed in
Table 2. We emphasize that all parameters are best estimates for ISM
conditions i.e. assuming a cold and dry environment. In general, these
parameters (especially the surface energy $\gamma$), have a wide range because
they are highly sensitive to temperature (Bogdan et al. 2020), monomer size
(Bauer et al. 2019), and humidity (Heim et al. 1999; Fuji et al. 1999; Kimura
et al. 2015; Steinpilz et al. 2019; Bogdan et al. 2020). Subsequent studies on
this matter must complement the mechanical parameters utilized in this paper.

##  5 Radiative torque disruption (RATD)

In this section we briefly outline the radiative torques (RAT) that lead to
rapid grain rotation and potentially to the disruption by centrifugal forces.
It is well established that grains with an irregular shape acquire a certain
amount of angular velocity when exposed to directed radiation (Dolginov &
Silantev 1976; Draine & Weingartner 1996, 1997; Weingartner & Draine 2003;
Lazarian & Hoang 2007a). Here, the gain of angular momentum over time follows

$\omega_{\mathrm{agg}}(t)=\omega_{\mathrm{RAT}}\left[1-\exp\left(-\frac{t}{\tau_{\mathrm{drag}}}\right)\right]\,,$
(24)

where $\tau_{\mathrm{drag}}$ is the characteristic timescale of the rotational
drag by means of gas collisions (Draine 1996) and photon emission (Draine &
Lazarian 1998) and

$\omega_{\mathrm{RAT}}=\frac{\Gamma_{\mathrm{RAT}}\tau_{\mathrm{drag}}}{I_{\mathrm{a1}}}$
(25)

is the terminal angular velocity. The torque $\Gamma_{\mathrm{RAT}}$ is
characteristic for individual grains and depends on the spectrum of the
radiation field as well as the shape of a particular grain and its material
composition (see e.g. Hoang et al. 2014). A solution of
$\Gamma_{\mathrm{RAT}}$ may be calculated by numerical approximations (Draine
1996; Draine & Flatau 2013; Herranen et al. 2019) or analytical toy models
(Lazarian & Hoang 2007a). However, we emphasize that we do not evaluate
$\omega_{\mathrm{RAT}}$ explicitly within the scope of this paper but assume a
maximal grain rotation to guarantee the disruption of all aggregates. For the
corresponding angular acceleration follows then

$\frac{\mathrm{d}\omega_{\mathrm{agg}}(t)}{\mathrm{d}t}=\frac{\omega_{\mathrm{RAT}}}{\tau_{\mathrm{drag}}}\exp\left(-\frac{t}{\tau_{\mathrm{drag}}}\right)\,.$
(26)

In principle, grains would also spin-up when exposed to an gaseous flow
(Lazarian & Hoang 2007b; Das & Weingartner 2016; Hoang et al. 2018, , RMK22)
leading to an mechanical torque (MET) $\Gamma_{\mathrm{MET}}$. However, the
time evolution and acceleration of the angular velocity would still follow the
same curves governed by Eq. 24 and Eq. 26, respectively, but evaluated with
$\Gamma_{\mathrm{MET}}$ instead of $\Gamma_{\mathrm{RAT}}$.
A dust aggregate cannot simply be modelled as a rigid body. Internal
relaxation processes such as Barnett relaxation (Purcell 1979; Lazarian &
Roberge 1997), nuclear relaxation (Lazarian & Draine 1999) or inelastic
relaxation (Purcell 1979; Lazarian & Efroimsky 1999) would dissipate
rotational energy. Subsequently, the dissipation would re-orient the direction
of the angular velocity $\omega_{\mathrm{agg}}$. For typical interstellar
conditions the total relaxation time is much smaller than the drag time
$\tau_{\mathrm{drag}}$ (Weingartner & Draine 2003) and the grain axis
$\hat{a}_{\mathrm{1}}$ becomes the most likely axis of grain rotation (compare
Fig. 1). For simplicity we assume in this study that the angular velocity
points always in the direction of $\hat{a}_{\mathrm{1}}$ i.e.
${\@vec{\omega}_{\mathrm{agg}}(t)=\hat{a}_{\mathrm{1}}\omega(t)}$ for all of
our three-dimensional N-body simulations.
Potentially, the aggregate becomes ripped apart by centrifugal forces long
before reaching the terminal angular velocity $\omega_{\mathrm{RAT}}$ (Silsbee
& Draine 2016; Hoang et al. 2019). A mathematical framework for the radiative
torque disruption (RATD) of grains was presented in Hoang et al. (2019). Here,
the critical angular velocity for grain suggested for RATD is

$\omega_{\mathcal{S}}=\frac{2}{a_{\mathrm{eff}}}\left(\frac{\mathcal{S}_{\mathrm{max}}}{\rho_{\mathrm{mat}}}\right)^{1/2}\,.$
(27)

Any dust aggregate exceeding the rotational limit
$\omega_{\mathrm{agg}}>\omega_{\mathcal{S}}$ becomes inevitably destroyed.
This theoretical upper limit is derived from the material density
$\rho_{\mathrm{mat}}$, the effective radius $a_{\mathrm{eff}}$, and the
maximal tensile strength $\mathcal{S}_{\mathrm{max}}$ of the aggregate (see
Hoang 2020, for further details). However, the tensile strength of porous dust
aggregates is not well constrained and may be lower by several orders of
magnitude compared to solid bodies because individual substructures are not
fully connected. In Greenberg et al. (1995) a means to estimate the tensile
strength of aggregates was provided by evaluating

$\mathcal{S}_{\mathrm{max}}=\frac{3}{2}\left<N_{\mathrm{con}}\right>\frac{\phi
E}{a_{\mathrm{mon}}h}$ (28)

(see also Li & Greenberg 1997). Here, $E$ is the binding energy, $h$ is
related to the overlap of monomers, $\left<N_{\mathrm{con}}\right>$ is the
average number of connections between all monomers, and $a_{\mathrm{mon}}$ is
the monomer size of a monodisperse aggregate.
Complementary N-body simulations by Seizinger et al. (2013b) suggest that the
tensile strength is not directly related to the initial volume filling factor
$\phi$ as assumed e.g. by Greenberg et al. (1995) or Blum et al. (2006).
Instead, the tensile strength is
${\mathcal{S}_{\mathrm{max}}\propto\phi^{1.6}}$ or
${\mathcal{S}_{\mathrm{max}}\propto\phi^{1.9}}$, respectively, for different
quartz BA, BAM, and hexagonal aggregates. Comparable results are presented by
Tatsuuma et al. (2019) where the relation
${\mathcal{S}_{\mathrm{max}}\propto\phi^{1.8}}$ was suggested in particular
for icy and quartz BCCA aggregates and
${\mathcal{S}_{\mathrm{max}}\propto\phi^{2/(3-D_{\mathrm{f}})}}$ for arbitrary
grain shapes with fractal dimension $D_{\mathrm{f}}$. We emphasize that in
these simulations the aggregates got merely stretched or compressed,
respectively, upon breaking, but do not experience any effects associated with
rotation at all.

##  6 Numerical setup

We implement the physical effects as outlined above into a C++ code that
combines the physics of Seizinger et al. (2012, 2013a) with the code presented
in RMK22 for the growth of aggregates. The time evolution of the net force
acting on each monomer,

$m_{\mathrm{mon,i}}\frac{\mathrm{d}\@vec{\varv}_{\mathrm{i}}}{\mathrm{d}t}=\@vec{F}_{\mathrm{cent,i}}+\@vec{F}_{\mathrm{cor,i}}+\@vec{F}_{\mathrm{eul,i}}+\sum_{\mathrm{j=1,i\neq
j}}^{N_{\mathrm{mon}}}{\left(\@vec{F}_{\mathrm{N,ij}}+\@vec{F}_{\mathrm{S,ij}}\right)}\,,$
(29)

and that of the corresponding torque

$I_{\mathrm{mon,i}}\frac{\mathrm{d}\@vec{\omega}_{\mathrm{i}}}{\mathrm{d}t}=\sum_{\mathrm{j=1,i\neq
j}}^{N_{\mathrm{mon}}}{\left(\@vec{\Gamma}_{\mathrm{S,ij}}+\@vec{\Gamma}_{\mathrm{R,ij}}+\@vec{\Gamma}_{\mathrm{T,ij}}\right)}\,,$
(30)

is calculated with second order accuracy by an symplectic Leap-Frog
integration scheme. Here, the quantity $\@vec{\varv}_{\mathrm{i}}$ is the
velocity of an individual monomer within the simulation domain, and
$\@vec{\omega}_{\mathrm{i}}$ is the angular velocity caused by sliding,
rolling, and twisting motions of its neighbors, whereas
${I_{\mathrm{mon,i}}=2/5\ m_{\mathrm{mon,i}}a_{\mathrm{mon,i}}^{2}}$ is the
moment of inertia of a particular monomer.
In order to account for the dissipation of energy we assume that sliding,
rolling, and twisting operate in the elastic limit until the corresponding
displacement reaches some characteristic critical limit. The theoretical
limits are ${\zeta_{\mathrm{c}}=r_{0}(1-\nu)/(16\pi)}$ for sliding and
${\phi_{\mathrm{c}}=1/(16\pi)}$ for twisting (compare Eq. 14 and Eq. 18).
However, the critical limit of rolling $\xi_{\mathrm{c}}$ (see Eq. 16) is
still highly debated. For silicate monomers of different sizes this limit is
within the range of ${\xi_{\mathrm{c}}\in[0.2\ \mathrm{nm},3.2\ \mathrm{nm}]}$
(Dominik & Tielens 1997; Heim et al. 1999; Paszun & Dominik 2008). In this
study we apply the conservative value of ${\xi_{\mathrm{c}}=0.2\ \mathrm{nm}}$
for the entire ensemble of aggregates independent of size and material. When
the displacements because of the relative motion of monomers in contact exceed
their critical values, energy is dissipated and our code modifies the contact
pointers $\@vec{n}_{\mathrm{i}}$ and $\@vec{n}_{\mathrm{j}}$ to ensure that
${|\@vec{\zeta}|=\zeta_{\mathrm{c}}}$, ${|\@vec{\phi}|=\phi_{\mathrm{c}}}$,
and ${|\@vec{\xi}|=\xi_{\mathrm{c}}}$, respectively (see Dominik & Nübold
2002; Wada et al. 2007, for details).
Individual monomers operate on a characteristic time scale that may be write
as (Wada et al. 2007)

$t_{\mathrm{dis,i}}=\frac{1}{6^{2/3}}\sqrt{\frac{m_{\mathrm{mon,i}}r_{\mathrm{0}}^{2}}{\pi\
\gamma R^{2}}}\,.$ (31)

Furthermore, it is not allowed for an monomer to move a distance larger than
its own radius limiting the time step further by
${t_{\mathrm{vel,i}}=a_{\mathrm{mon,i}}/\varv_{\mathrm{i}}}$. The final time
step for our Leap-Frog integration scheme is then the minimum of these
characteristic times ${\Delta t=0.02\
\min_{\mathrm{i}}\left(t_{\mathrm{dis,i}},t_{\mathrm{vel,i}}\right)}$ of all
the monomers as well as monomer connections within the simulations domain. The
factor of $\Delta t$ is a little smaller than the one suggested by Wada et al.
(2007) that guarantees conservation of energy below an error of $10^{-3}$.
However, monomers may roll along the surface of the rapidly rotating aggregate
and choosing a smaller time step allows for a more accurate spatial resolution
of such displacement of monomer within the aggregate.
The rotational damping acts usually on short timescales $\tau_{\mathrm{drag}}$
in the order of days up to several hundred years compared to typical
astronomical processes in the ISM (see e.g. Weingartner & Draine 2003; Tazaki
et al. 2017; Hoang 2022). A computational burden arises, because integration
time $\Delta t$ is only a few ns. Consequently, some simulations may take up
to $10^{18}$ time steps to terminate. Hence, a much smaller drag time
$\tau_{\mathrm{drag}}$ in the order of hours instead of years is selected to
reduce the total run-time of the code. However, this does not impact the
result of the rotational disruption simulations as long as the change in
centrifugal forces per time step remains much smaller than the displacement
processes for individual monomers to reach a new equilibrium position within
the rotating aggregate. Even for a $\tau_{\mathrm{drag}}$ of a few hours it is
still guaranteed that centrifugal forces are the dominant cause of grains
disruption because the spin-up remains slow enough that the shear within the
aggregate by Euler forces remains negligible. In detail, a
$\tau_{\mathrm{drag}}$ is selected such that the aggregate reaches an angular
velocity where the destruction is guaranteed within a few hours of simulation
time but the Centrifugal force remains always the dominant force exerted on
each monomer i.e. ${F_{\mathrm{eul,i}}\ll F_{\mathrm{cent,i}}}$. The
simulations setup starts at ${\omega_{\mathrm{agg}}=0\ \mathrm{rad\ s^{-1}}}$
and follows the curve of Eq.24 up to a given terminal angular velocity
${\omega_{\mathrm{RAT}}}$ by default. Alternatively the code terminates for
the condition $N_{\mathrm{con}}=1$ i.e. only two connected monomers remain
within the simulation domain.
For detecting monomer collisions with variable sizes and the tracking of
contact breaking events we use a loose octree data structure (see e.g. Ulrich
2000; Raschdorf & Kolonko 2009, for further details) in order to optimize the
runtime of the code even more. A new contact is established when two moving
monomers touch each other, which means $\delta\leq 0$. When individual
monomers or even entire sub-structures of the initial aggregate become
unconnected we only track the evolution and spin-up process of the most
massive remaining fragment while smaller fragments are allowed to leave the
simulation domain. After each time step the remaining connected fragments are
identified with a recursive flood fill algorithm on an undirected graph that
represents the monomer connection relations. For every thousandth time step,
the affiliation of monomer to a certain fragment, the number of monomers
$N_{\mathrm{mon}}$ remaining within the simulation domain and the monomer
positions $\@vec{X}_{\mathrm{i}}$ as well as the forces and connections in
between monomers are recorded.
In order to quantify the total stress in between all connected monomers within
the largest fragment we introduce the quantity of the total stress

$\sigma_{\Sigma}=\sum_{\mathrm{i=1}}^{N_{\mathrm{mon}}}\sum_{\mathrm{j>i}}^{N_{\mathrm{mon}}}{\frac{\left|F_{\mathrm{N,ij}}\right|}{\pi
r_{\mathrm{ij}}^{2}}}\,,$ (32)

where $F_{\mathrm{N,ij}}=0$ for non-connected monomers.

##  7 Results and discussion

We perform numerical N-body simulations for each individual pre-calculated
grain aggregate with our numerical N-body setup upon rotational disruption.
All grains are disrupted at an angular velocity of $\omega_{\mathrm{disr}}$
before the overall spin-up process reaches the terminal angular velocity of
${\omega_{\mathrm{RAT}}=3\cdot 10^{10}\ \mathrm{rad\ s^{-1}}}$. We emphasize
once again that the exact values of $\omega_{\mathrm{RAT}}$ is of minor
relevance for our N-body simulations as long as the disruption of each
aggregate is guaranteed and the centrifugal force and the Euler force follow
strictly the relation ${F_{\mathrm{eul,i}}\ll F_{\mathrm{cent,i}}}$ for the
entire spin-up process as governed by Eq. 24 and Eq. 26, respectively.

###  7.1 Time evolution of the rotational disruption process

Figure 4: Evolution of the mass $M(\omega_{\mathrm{agg}})$ (top panel), number
of connections $N_{\mathrm{con}}(\omega_{\mathrm{agg}})$ (middle panel), and
total stress $\sigma_{\Sigma}(\omega_{\mathrm{agg}})$ (bottom panel) within
the largest fragment dependent on the increasing angular velocity
$\omega_{\mathrm{agg}}$. The blue lines represent the exemplary ensemble of
q-S BAM2 grains with ${a_{\mathrm{eff}}=350\ \mathrm{nm}}$. An arrow points to
the peak value of the number of connections $N_{\mathrm{con}}$. Red dots
indicate the characteristics angular velocity $\omega_{\mathrm{agg}}$ up to
which an individual aggregate has lost $25\ \%$ of its initial mass, initial
number of connections, or reached its peak stress.

Figure 5: Three exemplary snapshots of a BAM2 grain aggregate with
$a_{\mathrm{eff}}=500\ \mathrm{nm}$ of an particular rotational disruption
simulation with the corresponding number of connection $N_{\mathrm{con}}$ for
the angular velocities of ${\omega/\omega_{\mathrm{disr}}=0.1}$ (top row),
${\omega/\omega_{\mathrm{disr}}=0.8}$ (bottom row),
${\omega/\omega_{\mathrm{disr}}=1.0}$ (bottom row). We emphasize that the
exemplary grain is identical to the one in the lower right corner of Fig. 1.
Left column: Color coded is the affiliation of the distinct connected
fragments. Smaller fragments are arbitrarily colored whereas the most massive
fragment is always depicted in blue. Right column: The monomers are colored
according to the largest magnitude of the normal force $F_{\mathrm{N}}$
exerted from its connected neighbours where $F_{\mathrm{N}}>0$ (green)
represents pushing forces while $F_{\mathrm{N}}<0$ (red) are pulling forces,
respectively.  Figure 6: Evolution of grain shapes for a-C BA (red), BAM1
(green), and BAM2 (blue) grains with an effective radius of
${a_{\mathrm{eff}}=200\ \mathrm{nm}}$ (left panel), ${a_{\mathrm{eff}}=350\
\mathrm{nm}}$ (middle panel), and ${a_{\mathrm{eff}}=500\ \mathrm{nm}}$ (right
panel), respectively. Dots represent the initial grain shape i.e.
${\omega_{\mathrm{agg}}/\omega_{\mathrm{disr}}=0}$ whereas crosses are the
grain shape at disruption for
${\omega_{\mathrm{agg}}/\omega_{\mathrm{disr}}=1}$. We note the clear tendency
of the porous grains to become oblate in shape with increasing
$\omega_{\mathrm{agg}}$.

Figure 7:  The angular velocity of disruption $\omega_{\mathrm{disr}}$
dependent on grain size for the ensembles of co-S (top panel), q-S (biddle
panel), and a-C (bottom panel) grain materials dependent on the effective
radius $a_{\mathrm{eff}}$. Color coded are the BA (red), BAM1 (green), and
BAM2 (blue) grains. Vertical bars are the range between minimal and maximal
values of $\omega_{\mathrm{disr}}$ resulting from our N-body disruption
simulations while solid lines represent the best fit model. Figure 8: The same
as Fig. 7 but only for co-S BAM2 grains (blue) in comparison with the
predictions of $\omega_{\mathrm{agg}}$ based on the tensile strengths models
presented in Greenberg et al. (1995) (G95), Seizinger et al. (2013b) (S13),
and Tatsuuma et al. (2019) (T19), respectively (gray). The latter are re-
scaled to match our simulation results for ${a_{\mathrm{eff}}=100\
\mathrm{nm}}$.

A typical simulation result for a q-S BAM2 grain with an effective radius of
$a_{\mathrm{eff}}=350\ \mathrm{nm}$ is shown in Fig. 4. The evolution of the
mass $M(\omega_{\mathrm{agg}})$ during the spin up process remains constant up
to ${\omega_{\mathrm{agg}}\approx 3\cdot 10^{8}\ \mathrm{rad\ s^{-1}}}$. Once
the grain spins even faster, contacts break and individual monomers as well as
smaller fragments start to break from the aggregates surface. Simultaneously,
monomers wander outwards driven by centrifugal forces. For a small period the
total number connections exceeds even the initial value and drops then
steadily. This features is most pronounced for BAM2 grains. Eventually, the
simulations reach the condition of $N_{\mathrm{con}}=1$ and terminate. The
breaking of contacts and the subsequent mass loss is continuous in nature
while the angular velocity $\omega_{\mathrm{agg}}$ increases by roughly one
order of magnitude. To designate a characteristic angular velocity
$\omega_{\mathrm{disr}}$ of rotational disruption based on mass loss or broken
connections would be arbitrary. In Fig. 4 we show also the time evolution of
the total stress $\sigma_{\Sigma}(\omega_{\mathrm{agg}})$ of each aggregate.
In contrast to the mass loss and the braking of contacts the stress
$\sigma_{\Sigma}(\omega_{\mathrm{agg}})$ rises steadily but drops then
abruptly. Hence, we define this unique feature in the time evolution of each
aggregate to be associated with the angular velocity of
$\omega_{\mathrm{disr}}$ of rotational disruption. The value of
$\omega_{\mathrm{disr}}$ is close to but not identical with a mass loss of
$25\ \%$.

###  7.2 The fragmentation of rotating grains

In Fig. 5 we show snapshots of a typical spin-up process and rotational
disruption event for one representative grain. The grain is identical to the
BAM2 with $a_{\mathrm{eff}}=500\ \mathrm{nm}$ depicted to the bottom right of
Fig. 1. Once the grain has spun up to an angular velocity of
${\omega_{\mathrm{agg}}/\omega_{\mathrm{disr}}=0.1}$ the aggregate starts to
experience a stretching force and first connections break. However, all
monomers are still connected to the same aggregate. Compared to its original
configuration as depicted in Fig. 1 the aggregate changed its shape already,
because monomers are re-arranged within the aggregate by means of rolling. The
aggregate fragments for ${\omega_{\mathrm{agg}}/\omega_{\mathrm{disr}}=0.8}$
and new connections are established as monomers move into a new equilibrium
position driven by pulling forces and pushing forces between connected
neighbours. We report that in this phase small disconnected fragments may
rarely reconnect to the most massive fragment. For an angular velocity of
${\omega_{\mathrm{agg}}/\omega_{\mathrm{disr}}=1.0}$ larger fragments are
separated from the most massive cluster. The remaining cluster goes through a
phase of relaxation, where its total stress $\sigma_{\Sigma}$ declines while
simultaneously the mass loss is still an ongoing process. We note that if the
spin-up process of the grain were to stop at an angular velocity
${\omega_{\mathrm{agg}}/\omega_{\mathrm{disr}}\leq 1.0}$ the aggregate would
still fragment, but the mass loss would stop eventually as the remaining most
massive fragment reaches a new equilibrium configuration.
Dust destruction processes such as shattering efficiently redistribute larger
grains sizes towards smaller ones (Dwek & Scalo 1980; Tielens et al. 1994). In
modelling this process it is usually assumed that the new size distribution of
the fragments follows a power-law (Hellyer 1970; Hirashita & Yan 2009;
Kirchschlager et al. 2019). A similar power-law is assumed in models of the
grain redistribution by means of RATD (e.g. Giang et al. 2020). However, the
fragmentation process of rapidly rotating porous dust aggregates remains
poorly constrained. In fact, our N-body simulation results so far indicate
that the considered BAM aggregates almost completely break down into their
individual building blocks. Hence, the resulting size distribution of the
fragments is expected to be almost identical to the initial size distribution
of the monomers. However, we note that our numerical setup does not track
smaller fragments once they become separated from the most massive fragment.
Furthermore, each individual fragment will experience its own own individual
spin-up and drag (see Sect. 5). Smaller fragments may also eventually become
too small for an RAT effective for a further spin-up upon the disruption point
(Lazarian & Hoang 2007a). Such effects need to be taken into account in
forthcoming studies to answer question about the resulting size distribution
of an ensemble of rotationally disrupted grains conclusively.

###  7.3 Rotational deformation of grain shapes

In Fig. 6 we present the time evolution of the shapes of a-C grains with
different sizes up to the breaking point. Initially, for
${\omega_{\mathrm{agg}}/\omega_{\mathrm{disr}}=0.0}$ the ensemble of grains
with an effective radius of $a_{\mathrm{eff}}=200\ \mathrm{nm}$ are almost
oblate and prolate shapes in equal parts. Larger grains are with
$a_{\mathrm{eff}}=350\ \mathrm{nm}$ and $a_{\mathrm{eff}}=500\ \mathrm{nm}$
are slightly more prolate where most of the axis ratios are $c/b<1.7$. At the
breaking point i.e. ${\omega_{\mathrm{agg}}/\omega_{\mathrm{disr}}=1.0}$ the
a-C grains become deformed where a oblate shapes are the most likely outcome
with an axis ratio up to $b/a\leq 5.0$ for BAM grains. The exception is the
ensemble of small BA grains with $a_{\mathrm{eff}}=200\ \mathrm{nm}$ where a
prolate shape seem to be the more favourable configuration. We note that BA
grains reach a smaller axis ratios because they break more easily with a much
shorter deformation phase compared to BAM1 and BAM2 grains. We report similar
trends for the deformation of grains with q-S and co-S materials,
respectively, where oblate grains are the most likely shape of rotational
disruption.

|  | a-C |  |  | q-S |  |  | co-S |
---|---|---|---|---|---|---|---|---|---
| BA | BAM1 | BAM2 | BA | BAM1 | BAM2 | BA | BAM1 | BAM2
A | 1.97 | 2.28 | 2.25 | 2.13 | 2.18 | 2.22 | 1.50 | 1.65 | 1.57
$\alpha$ | 0.10 | 0.12 | 0.13 | 0.10 | 0.14 | 0.14 | 0.10 | 0.14 | 0.15
Table 3: Best fit parameters of $\omega_{\mathrm{disr}}$ based on our N-body
disruption simulation results. Figure 9: The same as Fig. 7 but only for a-C
BA (deshed lines) and BAM2 (dotted lines) grains. The left panel shows the
impact of critical rolling displacement $\xi_{\mathrm{crit}}\in[0.2\
\mathrm{nm},6.4\ \mathrm{nm}]$ on $\omega_{\mathrm{agg}}$ while the viscous
dumping time $T_{\mathrm{vis}}$ remains constant. The right panel is for a
constant $\xi_{\mathrm{crit}}$ but with $T_{\mathrm{vis}}\in[1\ \mathrm{ps},9\
\mathrm{ps}]$. Red lines represent the default parameters of our N-body
disruption simulations.

Quantifying the grain shape by the fractal dimension $D_{\mathrm{f}}$ during
the spin-up process reveals no clear trend. As depicted in Fig. 2 the fractal
dimension is not well correlated with the grain shapes of the initial
ensemble. Up to the braking point the variation of $D_{\mathrm{f}}$ becomes
even larger. The same for the porosity $\mathcal{P}$. At the beginning of the
spin up process $\mathcal{P}$ starts to increase slightly for BAM1 and BAM2
grains while the BA ones break without an increase in $\mathcal{P}$. Close to
the breaking point the porosity distribution has a large variation and becomes
virtually identical for BAM1 and BAM2 grains. Eventually, the analysis of the
fractal dimension $D_{\mathrm{f}}$ as well as the porosity $\mathcal{P}$ does
no longer apply because the calculation of these quantities fails as the grain
aggregates break down into its individual monomers (see Sect. 3).
The shape and porosity of interstellar dust is still a matter of debate. For
example BA grain growth processes favor roundish grain aggregates with a
fractal dimension of about $D_{\mathrm{f}}=2.0$ where the principle axis are
${a\approx b\approx c}$. Guillet et al. (2018) developed a grain model based
on spheroids of amorphous silicate and amorphous carbon to reproduce both
starlight polarization and polarized sub-millimeter emission. The best fit
model suggest prolate grains with an axes ratio of ${b/a=1/3}$ and a porosity
of $\mathcal{P}=0.2$. More recently, Draine & Hensley (2021b) suggest prolate
grains with an axis ratio of $b/a=0.6$ or oblate grains with $b/a=1.5$ with a
porosity of about $\mathcal{P}=0.4$. Comparing the initial porosity of our
ensemble grown by BAM (see Fig. 2) reveals that only the BAM2 ensemble would
match the limitation in porosity with values $\mathcal{P}=0.4-0.5$.
Elongated grains may be the result of hit and stick processes of BCPA with
preferential direction, e.g. for magneto-hydrodynamic turbulence (Yan &
Lazarian 2003) or for grain aggregates aligned with the magnetic field
direction (Hoang 2022). However, the latter effect requires rapid grain
rotation. Subsequently, the spin-up process would increase the relative
velocity between the surface of the rotating aggregate and impinging monomers
may potentially destroy the entire aggregate in a catastrophic disruption
event (Benz & Asphaug 1999; Morris & Burchell 2017; Schwartz et al. 2018).
What we find with our N-body simulations is that accelerated rotation deforms
initially more roundish aggregates preferentially into oblate shapes.
For RATs the axis ration of the grains is linked to the radiation field via
the maximal angular velocity $\omega_{\mathrm{RAT}}$. The more extreme cases,
as depicted in Fig. 6 with ${b/a=4-5}$ would require for the grains to be in
close vicinity of a supernova or
active galactic nucleus (Hoang et al. 2019; Giang et al. 2020; Giang & Hoang
2021). We emphasize that the torques $\Gamma_{\mathrm{RAT}}$ and
$\Gamma_{\mathrm{MET}}$ (see Sect. 5) are tightly connected to the grain shape
(Lazarian & Hoang 2007a; Das & Weingartner 2016, RMK22). Consequently, the
terminal angular velocity $\omega_{\mathrm{RAT}}$ would be marginally
increased during the deformation phase accelerating the grain rotation even
more (see Eq. 26). However, this effect would be continuous over a long time
span and not considerably impact the internal grain dynamics as far as the
final value of $\omega_{\mathrm{disr}}$ is concerned. More severe is the
change in $\omega_{\mathrm{RAT}}$ during the fragmentation phase. Here, the
RAT would decrease as the aggregate fragments and subsequently the rotation
would slow down and the grain stabilizes by reaching a new equilibrium
configuration. However, considering a dynamical spin-up process in our N-body
simulation would not only require to calculate the torque
$\Gamma_{\mathrm{RAT}}$ but also the rotational drag timescale
$\tau_{\mathrm{drag}}$ for each time step by means of time consuming numerical
approximate methods (Draine & Flatau 2013, RMK22). Potentially, both
quantities $\Gamma_{\mathrm{RAT}}$ and $\tau_{\mathrm{drag}}$ may be
parameterized based on the shape and material parameters of each individual
BAM grain. For now a dynamical spin-up goes beyond the scope of this paper.

###  7.4 Average rational disruption of the grain ensemble

In Fig. 7 we present the characteristic angular velocity
$\omega_{\mathrm{disr}}$ for the entire set of simulation results. The
magnitude of $\omega_{\mathrm{disr}}$ is roughly in the range of about
${5\cdot 10^{8}-5\cdot 10^{9}\ \mathrm{rad\ s^{-1}}}$ for different materials
and grain sizes. This result agrees well with the predictive model presented
in Hoang et al. (2019). However, the exact value of $\omega_{\mathrm{disr}}$
depends on the exact material properties and the internal structure of the
initial grain. While rapidly rotating solid bodies break into few fragments by
driven by centrifugal forces, a porous material may fragment at lower angular
velocity and into a larger number of pieces. The maximal tensile strength
$\mathcal{S}_{\mathrm{max}}$ is usually utilized to quantify the response of a
material during stretching. However, what we find that utilizing the maximal
tensile strength $\mathcal{S}_{\mathrm{max}}$ as introduced in Eq. 28 to
calculate the critical angular velocity $\omega_{\mathcal{S}}$ with Eq. 27
does not reproduce the average rotational disruption $\omega_{\mathrm{disr}}$
resulting from our N-body simulations. We argue that
$\mathcal{S}_{\mathrm{max}}$ is insufficient to describe the dynamical
behavior of rotating aggregates because in a stretching aggregate each
monomers experience only the local forces in between its connected neighbours.
In a rotating aggregate, however, each individual monomer experiences also the
centrifugal force. The maximal tensile strength $\mathcal{S}_{\mathrm{max}}$
is determined by stretching of granular materials performed along a distinct
axis, whereas the centrifugal force is radial with respect to the rotation
axis. Furthermore, considering only stretching produces acting on an aggregate
does not not lead to a large displacement of monomers compared to the
aggregate’s scale. It is important to note that, our N-body simulations reveal
that monomers are rolling outwards within the grain aggregates. This movement
of $N_{\mathrm{mon}}$ monomers allows potentially to establish up to
${N_{\mathrm{mon}}(N_{\mathrm{mon}}-1)\approx N_{\mathrm{mon}}^{2}}$ new
connections. Naturally, some monomers are already connected and only a small
fraction of all unconnected monomers is close enough within the aggregate to
newly connect. Hence, the relation cannot be quadratic but with an much
smaller exponent $\alpha\ll 2$. In order to match our data we suggest to
extend the definition of Eq. 28 by the number of monomers $N_{\mathrm{mon}}$.
By utilizing Eq. 27 the characteristic angular velocity of rotation disruption
for a polydisperse grain aggregate reads then

$\omega_{\mathrm{disr}}=\frac{A}{a_{\mathrm{eff}}}\sqrt{\frac{\gamma}{\rho_{\mathrm{mat}}\left<a_{\mathrm{mon}}\right>}}\Biggl{(}\
\left<N_{\mathrm{con}}\right>N_{\mathrm{mon}}\Phi\ \Biggr{)}^{\alpha}$ (33)

where $\left<a_{\mathrm{mon}}\right>$ is the average monomer radius and $A$
and $\alpha$, respectively, are fit parameters333Note that
${\left<N_{\mathrm{con}}\right>=N_{\mathrm{con}}/N_{\mathrm{mon}}}$ and thus
Eq. 33 may also be written in an equivalent form with a factor
${(\left<N_{\mathrm{con}}\right>N_{\mathrm{mon}}\Phi\
)^{\alpha}=(N_{\mathrm{con}}\Phi\ )^{\alpha}}$ to match the simulation
results. Assigning a separate exponent to each of the individual quantities
$\left<N_{\mathrm{con}}\right>$, $N_{\mathrm{mon}}$, and $\Phi$, respectively,
does not improve the accuracy of the fit. Best fit results of
$\omega_{\mathrm{disr}}$ are plotted in Fig. 7. The fit matches very well with
the ensemble average of the different considered grain sizes.
In Fig. 8 we show a comparison of our best fit of $\omega_{\mathrm{disr}}$ for
co-S BAM 2 grains with that calculated with different parameterizations of the
tensile strength $\mathcal{S}_{\mathrm{max}}$ given in literature (see also
Sect. 5). These parameterizations of $\mathcal{S}_{\mathrm{max}}$ may not
necessarily be evaluated with the parameters provided by our grain models.
Hence, we scale the resulting to match our results of $\omega_{\mathrm{disr}}$
at for an effective radius of ${a_{\mathrm{eff}}=100\ \mathrm{nm}}$. The
comparison reveals that previous attempts of modelling
$\omega_{\mathrm{disr}}$ by a volume filling factor $\phi$ dependent tensile
strength $\mathcal{S}_{\mathrm{max}}$ cannot reproduce the asymptotic behavior
of our simulation results towards larger grain radii $a_{\mathrm{eff}}$.
Finally, we emphasize that the presented results of rotational disruption are
calculated for aggregates of carbonaceous and silicate monomers loosely
connected by van der Waals force. However, materials other than carbon or
silicates may form much stronger bonds (Dominik & Tielens 1997). For instance,
pure iron in the form of small pallets may be present in the interstellar dust
that would create metallic bonds between monomers (Dominik & Tielens 1997;
Draine & Hensley 2021a). The same for ices covering the surface of the grain
aggregate where dipole-dipole interactions may increase the resistance against
rotational disruption. Furthermore, high impact collisions of monomers or
compressive stress acting on the aggregate may lead to sintering, an effect
where the monomers fuse at the contact surface. Subsequently, sintering leads
to a neck between monomers (Maeno & Ebinuma 1983; Blackford 2007; Sirono &
Ueno 2017) strengthening the connection between monomers as well but at the
same time it makes the aggregate a whole more brittle, because as a neck would
not allow for rolling motions between monomers. Moreover, an effect that may
weaken monomer connections is dust heating. A higher dust temperature
decreases the surface energy $\gamma$ up to one order of magnitude (Bogdan et
al. 2020). Any luminous environment with a radiation field strong enough to
drive grain rotation up to $\approx 10^{9}\ \mathrm{rad\ s^{-1}}$ would
inevitably heat the dust grains up to several hundredth of Kelvin. Thus
$\gamma$ and subsequently $\omega_{\mathrm{disr}}$ would decrease. Altogether,
such effects need to be taken into consideration in forthcoming studies in
order to complement our numerical setup of the rotational disruption of of
porous dust aggregates.

###  7.5 Impact of critical rolling displacement and viscous damping

The critical rolling displacement $\xi_{\mathrm{crit}}$ and the viscous
damping time scale $T_{\mathrm{vis}}$ are the least constrained parameters in
our N-body simulations. Laboratory data indicates a huge variation in these
parameters of about one order of magnitude for silicates (Dominik & Tielens
1995; Heim et al. 1999; Krijt et al. 2013).
In order to quantify the impact $\xi_{\mathrm{crit}}$ and $T_{\mathrm{vis}}$,
respectively, on rational disruption we repeat our N-body simulation for a
subset of the pre-calculated grain aggregates varying $\xi_{\mathrm{crit}}$ in
the range $0.2\ \mathrm{nm}$ to $6.4\ \mathrm{nm}$ and $T_{\mathrm{vis}}$ from
$1\ \mathrm{ps}$ to $9\ \mathrm{ps}$. In Fig. 9 we present the exemplary
results of the parameter test of a-C BA and BAM2 grains. Variations in the
rolling displacement $\xi_{\mathrm{crit}}$ show little effect. The resulting
angular velocity of rotational disruption $\omega_{\mathrm{disr}}$ shows small
fluctuations but no clear trend with an increase in $\xi_{\mathrm{crit}}$.
Hence, we attribute the fluctuation to numerical effects (see also appendix
A). In contrast to $\xi_{\mathrm{crit}}$ an increase in $T_{\mathrm{vis}}$
leads to a clear increase of $\omega_{\mathrm{disr}}$. This effect is most
evident for BAM2 grains with $\omega_{\mathrm{disr}}$ being about $11\ \%$
higher for the largest applied viscous time of $T_{\mathrm{vis}}=9\
\mathrm{ps}$. Similar trends can be reported for the q-S and co-S materials.
Further improvements in the predictive accuracy of our rational disruption
simulations cannot be achieved upon forthcoming laboratory data especially for
the viscoelastic damping of oscillations within carbonaceous aggregates.

##  8 Summary

We aimed to study the rotational disruption of interstellar dust aggregate
analogs. An ensemble of porous dust aggregates by means of ballistic
aggregation and migration (BAM) is pre-calculated where for BA, BAM1, and BAM2
grains each monomer has at least one, two, or three connections, respectively,
with its neighbors. Here, we modified the original BAM algorithm presented in
Shen et al. (2008) to work with variable monomer sizes. We estimated the
composition of the grain aggregates based on the abundance of elements in the
ISM to approximate their mechanical properties. Numerical three-dimensional
N-body simulations are performed to determine characteristic angular velocity
of rational disruption $\omega_{\mathrm{disr}}$ for each aggregate
individually. The numerical setup is based on the work of Dominik & Tielens
(1997), Wada et al. (2007), Seizinger et al. (2012), and RMK22, respectively.
We modified their setup by intruding additional forces associated with an
accelerated rotation of an aggregate. A subsequent analysis of the disruption
event allows to describe the average rotational disruption of porous grain
ensembles.
The findings of this study are summarized as follows:

  * •

Compared to the original BAM algorithm considering only single a monomer size
we report that a grain growth by BAM with polydisperse monomers results in
aggregates that are more compact and less porous as smaller monomer tend to
migrate deeper into the aggregate. For the same reason the porosity
$\mathcal{P}$ stagnates with an increasing effective radius $a_{\mathrm{eff}}$
with values up to $\mathcal{P}=0.4-0.85$ for BA whereas BAM1 and BAM2 grains
are less porous with $\mathcal{P}=0.35-0.7$ and $\mathcal{P}=0.35-0.5$,
respectively.

  * •

Our simulation results reveal that rotating porous dust aggregates are not
disrupted by a single abrupt event characteristic for a brittle material but
rather by a continuous process where aggregates experience a continuous mass
loss. Subsequently, the initial aggregate breaks ultimately down into
fragments not larger than a few monomers.

  * •

The initial distribution of the pre-calculated polydisperse BAM grains have
oblate and prolate shapes, respectively, almost in equal parts. However, under
accelerated rotation the grain aggregates enter a phase of deformation and,
subsequently, the grain shapes are finally preferentially redistributed
towards oblate shapes.

  * •

We introduce the quantity of the total stress $\sigma_{\Sigma}$ as a measure
for the internal aggregate dynamics and time evolution. We report that
$\sigma_{\Sigma}$ increases up to a peak value characteristic for individual
aggregates while the rotation accelerates and subsequently starts to loose
mass. This peak roughly coincides with a mass loss of $25\ \%$ with respect to
the initial grain mass. For higher angular velocities $\sigma_{\Sigma}$ drops
sharply while the mass loss continues. We utilize this peak in
$\sigma_{\Sigma}$ to define the breaking point on a physical basis and finally
to determine the characteristic angular velocity $\omega_{\mathrm{disr}}$ of
rotational disruption.

  * •

In the deformation phase of BAM aggregates individual monomers are moving
along the grain surface driven by the centrifugal force and new connections
may be established. Consequently, the additional connections stabilize the BAM
aggregates against the disruptive centrifugal forces acting on the monomers.

  * •

Our N-body simulations reveal that the angular velocity
$\omega_{\mathrm{disr}}$ reaches an asymptotic limit towards larger grain
sizes of $a_{\mathrm{eff}}\gtrapprox 300\ \mathrm{nm}$. This finding is in
contrast to previous attempts to describe the rational disruption of porous
aggregates analytically based on the maximal tensile strength where
$\omega_{\mathrm{disr}}$ would continuously decrease for an increasing
$a_{\mathrm{eff}}$.

##  Appendix A Connection breaking test

In this sections we provide a test scenario for the accuracy of our code. Due
to the complexity of N-body simulation it is not feasible to derive an
analytical solution for an entire aggregate. However, the problem can exactly
be solved for two connected monomer with of equal size i.e.
${a_{\mathrm{mon,i}}=a_{\mathrm{mon,j}}}$ where the reduced radius simply
becomes $R=a_{\mathrm{mon,i}}/2$. In this case the centrifugal force (see Eq.
19) may be written as

$F_{\mathrm{cent,i}}=-m_{\mathrm{i}}\omega_{\mathrm{agg}}^{2}a_{\mathrm{mon,i}}$
(34)

by using the criterion of a broken contact $F_{\mathrm{N,ij}}=-F_{\mathrm{C}}$
the corresponding critical contact radius becomes
$r_{\mathrm{C}}=(1/6)^{2/3}r_{0}$. Putting $r_{\mathrm{C}}$ into Eq. 12 allows
the solve for the exact angular velocity

$\omega_{\mathrm{ref}}=\sqrt{\frac{5\pi}{6}\frac{\gamma}{m_{\mathrm{mon,i}}}}$
(35)

where the two equally sized monomers become separated by centrifugal forces.
We use this quantity for reference to evaluate the accuracy of our code. A
number of 80 benchmark runs are perform with randomly selected monomer raddii
${a_{\mathrm{mon,i}}\in[10\ \mathrm{nm},100\ \mathrm{nm}]}$. In Fig. 10 we
compare the resulting angular velocity $\omega_{\mathrm{disr}}$ of our
numerical setup introduced in Sect. 6 with Eq. 35. The resulting error
${(\omega_{\mathrm{ref}}-\omega_{\mathrm{disr}})/\omega_{\mathrm{ref}}}$ is
below $1.8\ \%$ with a trend of underpredicting $\omega_{\mathrm{ref}}$. An
error of $1.8\ \%$ is comparable to the range we observe for the parameters
test of the rolling displacement $\xi_{\mathrm{crit}}$ as presented in Fig. 9.
Hence, we assume a numerical accuracy of about $1.8\ \%$ for our numerical
N-body simulations. Consequently, the uncertainties in our reported results of
the angular velocity of rotational disruption $\omega_{\mathrm{ref}}$ is most
impacted by the lack of exact laboratory material parameters rather than
numerical errors.

Figure 10: The distribution of the angular velocity $\omega_{\mathrm{disr}}$
resulting from a number of runs $N_{\mathrm{run}}$ of N-body simulations in
compression with the corresponding exact analytical solution
$\omega_{\mathrm{disr}}$.

###### Acknowledgements.

Special thanks goes to Wilhelm Kley for numerous fruitfull discussion about
N-body simulations of aggregates. The authors thank Thiem Hoang, Cornelis P.
Dullemond, and Bruce T. Draine for usefull insights into the topic of dust
composition and dynamics. S.R., P.N., and R.S.K. acknowledge financial support
from the Heidelberg cluster of excellence (EXC 2181 - 390900948) “STRUCTURES:
A unifying approach to emergent phenomena in the physical world, mathematics,
and complex data”, specifically via the exploratory project EP 4.4. S.R. and
R.S.K. also thank for support from Deutsche Forschungsgemeinschaft (DFG) via
the Collaborative Research Center (SFB 881, Project-ID 138713538) ’The Milky
Way System’ (subprojects A01, A06, B01, B02, and B08). And we thanks for
funding form the European Research Council in the ERC synergy grant “ECOGAL –
Understanding our Galactic ecosystem: From the disk of the Milky Way to the
formation sites of stars and planets” (project ID 855130). The project made
use of computing resources provided by The Länd through bwHPC and by DFG
through grant INST 35/1134-1 FUGG. Data are in part stored at SDS@hd supported
by the Ministry of Science, Research and the Arts and by DFG through grant
INST 35/1314-1 FUGG.

## References

  * Baric et al. (2018) Baric, V., Grossmann, H. K., Koch, W., & Mädler, L. 2018, Particle & Particle Systems Characterization, 35, 1800177
  * Barlow et al. (2010) Barlow, M. J., Krause, O., Swinyard, B. M., et al. 2010, A&A, 518, L138
  * Bauer et al. (2019) Bauer, F., Daun, K., Huber, F., & Will, S. 2019, Appl. Phys. B, 125, 109
  * Benz & Asphaug (1999) Benz, W. & Asphaug, E. 1999, Icarus, 142, 5
  * Bertini et al. (2007) Bertini, I., Thomas, N., & Barbieri, C. 2007, A&A, 461, 351
  * Bescond et al. (2014) Bescond, A., Yon, J., Ouf, F. X., et al. 2014, Aerosol Science and Technology, 48, 831
  * Bevan & Barlow (2016) Bevan, A. & Barlow, M. J. 2016, MNRAS, 456, 1269
  * Birnstiel et al. (2010) Birnstiel, T., Dullemond, C. P., & Brauer, F. 2010, A&A, 513, A79
  * Blackford (2007) Blackford, J. R. 2007, Journal of Physics D: Applied Physics, 40, R355
  * Blum et al. (2006) Blum, J., Schräpler, R., Davidsson, B. J. R., & Trigo-Rodríguez, J. M. 2006, ApJ, 652, 1768
  * Bogdan et al. (2020) Bogdan, T., Pillich, C., Landers, J., Wende, H., & Wurm, G. 2020, A&A, 638, A151
  * Cardarelli (2008) Cardarelli, F. 2008, Materials handbook: a concise desktop reference, 2nd edn. (London: Springer)
  * Chakrabarty et al. (2007) Chakrabarty, R. K., Moosmüller, H., Arnott, W. P., et al. 2007, Appl. Opt., 46, 6990
  * Chokshi et al. (1993) Chokshi, A., Tielens, A. G. G. M., & Hollenbach, D. 1993, ApJ, 407, 806
  * Christensen (1996) Christensen, N. I. 1996, Journal of Geophysical Research: Solid Earth, 101, 3139
  * Compiègne et al. (2011) Compiègne, M., Verstraete, L., Jones, A., et al. 2011, A&A, 525, A103
  * Das & Weingartner (2016) Das, I. & Weingartner, J. C. 2016, MNRAS, 457, 1958
  * Do-Duy et al. (2020) Do-Duy, T., Wright, C. M., Fujiyoshi, T., et al. 2020, MNRAS, 493, 4463
  * Dolginov & Silantev (1976) Dolginov, A. Z. & Silantev, N. A. 1976, Ap&SS, 43, 337
  * Dominik & Nübold (2002) Dominik, C. & Nübold, H. 2002, Icarus, 157, 173
  * Dominik & Tielens (1995) Dominik, C. & Tielens, A. G. G. M. 1995, Philosphical Magazine A, 72, 783
  * Dominik & Tielens (1996) Dominik, C. & Tielens, A. G. G. M. 1996, Philosophical Magazine, Part A, 73, 1279
  * Dominik & Tielens (1997) Dominik, C. & Tielens, A. G. G. M. 1997, ApJ, 480, 647
  * Dorschner & Henning (1995) Dorschner, J. & Henning, T. 1995, A&A Rev., 6, 271
  * Draine (1996) Draine, B. T. 1996, in Astronomical Society of the Pacific Conference Series, Vol. 97, Polarimetry of the Interstellar Medium, ed. W. G. Roberge & D. C. B. Whittet, 16
  * Draine & Flatau (2013) Draine, B. T. & Flatau, P. J. 2013, arXiv e-prints, arXiv:1305.6497
  * Draine & Hensley (2021a) Draine, B. T. & Hensley, B. S. 2021a, ApJ, 909, 94
  * Draine & Hensley (2021b) Draine, B. T. & Hensley, B. S. 2021b, ApJ, 919, 65
  * Draine & Lazarian (1998) Draine, B. T. & Lazarian, A. 1998, ApJ, 508, 157
  * Draine & Lee (1984) Draine, B. T. & Lee, H. M. 1984, ApJ, 285, 89
  * Draine & Li (2007) Draine, B. T. & Li, A. 2007, ApJ, 657, 810
  * Draine & Weingartner (1996) Draine, B. T. & Weingartner, J. C. 1996, ApJ, 470, 551
  * Draine & Weingartner (1997) Draine, B. T. & Weingartner, J. C. 1997, ApJ, 480, 633
  * Dwek & Scalo (1980) Dwek, E. & Scalo, J. M. 1980, ApJ, 239, 193
  * Escatllar et al. (2019) Escatllar, A. M., Lazaukas, T., Woodley, S. M., & Bromley, S. T. 2019, ACS Earth and Space Chemistry, 3, 2390, pMID: 32055761
  * Fogerty et al. (2016) Fogerty, S., Forrest, W., Watson, D. M., Sargent, B. A., & Koch, I. 2016, ApJ, 830, 71
  * Fuji et al. (1999) Fuji, M., Machida, K., Takei, T., Watanabe, T., & Chikazawa, M. 1999, Langmuir, 15, 4584
  * Gail et al. (2009) Gail, H. P., Zhukovska, S. V., Hoppe, P., & Trieloff, M. 2009, ApJ, 698, 1136
  * Galametz et al. (2019) Galametz, M., Maury, A. J., Valdivia, V., et al. 2019, A&A, 632, A5
  * Giang & Hoang (2021) Giang, N. C. & Hoang, T. 2021, ApJ, 922, 47
  * Giang et al. (2020) Giang, N. C., Hoang, T., & Tram, L. N. 2020, ApJ, 888, 93
  * Goto et al. (2003) Goto, M., Gaessler, W., Hayano, Y., et al. 2003, ApJ, 589, 419
  * Gouriet et al. (2019) Gouriet, K., Carrez, P., & Cordier, P. 2019, Minerals, 9, 787
  * Greenberg et al. (1995) Greenberg, J. M., Mizutani, H., & Yamamoto, T. 1995, A&A, 295, L35
  * Guillet et al. (2018) Guillet, V., Fanciullo, L., Verstraete, L., et al. 2018, A&A, 610, A16
  * Heim et al. (1999) Heim, L.-O., Blum, J., Preuss, M., & Butt, H.-J. 1999, Phys. Rev. Lett., 83, 3328
  * Hellyer (1970) Hellyer, B. 1970, MNRAS, 148, 383
  * Hensley & Draine (2021) Hensley, B. S. & Draine, B. T. 2021, ApJ, 906, 73
  * Herranen et al. (2019) Herranen, J., Lazarian, A., & Hoang, T. 2019, ApJ, 878, 96
  * Hertz (1896) Hertz, H. 1896, Miscellaneous papers (Macmillan)
  * Hirashita & Yan (2009) Hirashita, H. & Yan, H. 2009, MNRAS, 394, 1061
  * Hoang (2020) Hoang, T. 2020, Galaxies, 8, 52
  * Hoang (2022) Hoang, T. 2022, ApJ, 928, 102
  * Hoang et al. (2018) Hoang, T., Cho, J., & Lazarian, A. 2018, ApJ, 852, 129
  * Hoang et al. (2014) Hoang, T., Lazarian, A., & Martin, P. G. 2014, ApJ, 790, 6
  * Hoang et al. (2019) Hoang, T., Tram, L. N., Lee, H., & Ahn, S.-H. 2019, Nature Astronomy, 3, 766
  * Israelachvili (2011) Israelachvili, J. N. 2011, Intermolecular and surface forces (Academic press)
  * Jensen et al. (2015) Jensen, B. D., Wise, K. E., & Odegard, G. M. 2015, Journal of Physical Chemistry A, 119, 9710
  * Johnson (1987) Johnson, K. L. 1987, Contact Mechanics
  * Johnson et al. (1971) Johnson, K. L., Kendall, K., & Roberts, A. D. 1971, Proceedings of the Royal Society of London Series A, 324, 301
  * Jones (2012) Jones, A. P. 2012, A&A, 542, A98
  * Kandilian et al. (2015) Kandilian, R., Heng, R.-L., & Pilon, L. 2015, Journal of Quantitative Spectroscopy and Radiative Transfer, 151, 310
  * Karasev et al. (2004) Karasev, V., Onischuk, A., Glotov, O., et al. 2004, Combustion and Flame, 138, 40
  * Karovicova et al. (2013) Karovicova, I., Wittkowski, M., Ohnaka, K., et al. 2013, A&A, 560, A75
  * Kataoka et al. (2013a) Kataoka, A., Tanaka, H., Okuzumi, S., & Wada, K. 2013a, A&A, 557, L4
  * Kataoka et al. (2013b) Kataoka, A., Tanaka, H., Okuzumi, S., & Wada, K. 2013b, A&A, 554, A4
  * Kelesidis et al. (2018) Kelesidis, G. A., Furrer, F. M., Wegner, K., & Pratsinis, S. E. 2018, Langmuir, 34, 8532, pMID: 29940739
  * Kendall et al. (1987) Kendall, K., McN. Alford, N., & Birchall, J. D. 1987, Proceedings of the Royal Society of London Series A, 412, 269
  * Kim et al. (2021) Kim, T. H., Takigawa, A., Tsuchiyama, A., et al. 2021, A&A, 656, A42
  * Kimura et al. (2015) Kimura, H., Wada, K., Senshu, H., & Kobayashi, H. 2015, The Astrophysical Journal, 812, 67
  * Kimura et al. (2020) Kimura, H., Wada, K., Yoshida, F., et al. 2020, MNRAS, 496, 1667
  * Kirchschlager et al. (2019) Kirchschlager, F., Schmidt, F. D., Barlow, M. J., et al. 2019, MNRAS, 489, 4465
  * Köylü & Faeth (1994) Köylü, U. O. & Faeth, G. M. 1994, Journal of Heat Transfer, 116, 971
  * Kozasa et al. (1992) Kozasa, T., Blum, J., & Mukai, T. 1992, A&A, 263, 423
  * Krijt et al. (2013) Krijt, S., Güttler, C., Heißelmann, D., Dominik, C., & Tielens, A. G. G. M. 2013, Journal of Physics D Applied Physics, 46, 435303
  * Lazarian & Draine (1999) Lazarian, A. & Draine, B. T. 1999, ApJ, 520, L67
  * Lazarian & Efroimsky (1999) Lazarian, A. & Efroimsky, M. 1999, MNRAS, 303, 673
  * Lazarian & Hoang (2007a) Lazarian, A. & Hoang, T. 2007a, MNRAS, 378, 910
  * Lazarian & Hoang (2007b) Lazarian, A. & Hoang, T. 2007b, ApJ, 669, L77
  * Lazarian & Roberge (1997) Lazarian, A. & Roberge, W. G. 1997, ApJ, 484, 230
  * Lehre et al. (2003) Lehre, T., Jungfleisch, B., Suntz, R., & Bockhorn, H. 2003, Appl. Opt., 42, 2021
  * Li & Draine (2001) Li, A. & Draine, B. T. 2001, ApJ, 554, 778
  * Li & Greenberg (1997) Li, A. & Greenberg, J. M. 1997, A&A, 323, 566
  * Liu et al. (2015) Liu, C., Yin, Y., Hu, F., Jin, H., & Sorensen, C. M. 2015, Aerosol Science and Technology, 49, 928
  * Maeno & Ebinuma (1983) Maeno, N. & Ebinuma, T. 1983, The Journal of Physical Chemistry, 87, 4103
  * Mathis et al. (1977) Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425
  * Matsuura (2011) Matsuura, M. 2011, in Astronomical Society of the Pacific Conference Series, Vol. 445, Why Galaxies Care about AGB Stars II: Shining Examples and Common Inhabitants, ed. F. Kerschbaum, T. Lebzelter, & R. F. Wing, 531
  * Matthews et al. (2012) Matthews, L. S., Land, V., & Hyde, T. W. 2012, ApJ, 744, 8
  * Mennella (2006) Mennella, V. 2006, ApJ, 647, L49
  * Michoulier & Gonzalez (2022) Michoulier, S. & Gonzalez, J.-F. 2022, MNRAS, 517, 3064
  * Min et al. (2007) Min, M., Waters, L. B. F. M., de Koter, A., et al. 2007, A&A, 462, 667
  * Mitchell (2004) Mitchell, R. H. 2004, 8th International Kimberlite Conference: The J. Barry Hawthorne volume, Vol. 2 (Gulf Professional Publishing)
  * Morris & Burchell (2017) Morris, A. & Burchell, M. 2017, Icarus, 296, 91
  * Nozawa et al. (2007) Nozawa, T., Kozasa, T., Habe, A., et al. 2007, ApJ, 666, 955
  * O’Donnell & Mathis (1997) O’Donnell, J. E. & Mathis, J. S. 1997, ApJ, 479, 806
  * Ormel et al. (2009) Ormel, C. W., Paszun, D., Dominik, C., & Tielens, A. G. G. M. 2009, A&A, 502, 845
  * Ossenkopf (1993) Ossenkopf, V. 1993, A&A, 280, 617
  * Paszun & Dominik (2008) Paszun, D. & Dominik, C. 2008, A&A, 484, 859
  * Paul et al. (2017) Paul, K. C., Silverstein, J., & Krekeler, M. P. 2017, Environmental Earth Sciences, 76, 1
  * Petrovic (2001) Petrovic, J. J. 2001, Journal of Materials Science, 36, 1573
  * Purcell (1979) Purcell, E. M. 1979, ApJ, 231, 404
  * Raschdorf & Kolonko (2009) Raschdorf, S. & Kolonko, M. 2009, Clausthal University of Technology
  * Reissl et al. (2022) Reissl, S., Meehan, P., & Klessen, R. S. 2022, arXiv e-prints, arXiv:2201.03694
  * Remediakis et al. (2007) Remediakis, I. N., Fyta, M. G., Mathioudakis, C., Kopidakis, G., & Kelires, P. C. 2007, Diamond and Related Materials, 16, 1835, proceedings of the 6th Specialists Meeting in Amorphous Carbon
  * Rogantini et al. (2019) Rogantini, D., Costantini, E., Zeegers, S. T., et al. 2019, A&A, 630, A143
  * Salameh et al. (2017) Salameh, S., van der Veen, M. A., Kappl, M., & van Ommen, J. R. 2017, Langmuir, 33, 2477, pMID: 28186771
  * Schwartz et al. (2018) Schwartz, S. R., Michel, P., Jutzi, M., et al. 2018, Nature Astronomy, 2, 379
  * Seizinger et al. (2013a) Seizinger, A., Krijt, S., & Kley, W. 2013a, A&A, 560, A45
  * Seizinger et al. (2012) Seizinger, A., Speith, R., & Kley, W. 2012, A&A, 541, A59
  * Seizinger et al. (2013b) Seizinger, A., Speith, R., & Kley, W. 2013b, A&A, 559, A19
  * Shen et al. (2008) Shen, Y., Draine, B. T., & Johnson, E. T. 2008, ApJ, 689, 260
  * Silsbee & Draine (2016) Silsbee, K. & Draine, B. T. 2016, ApJ, 818, 133
  * Sirono & Ueno (2017) Sirono, S.-i. & Ueno, H. 2017, ApJ, 841, 36
  * Skorupski et al. (2014) Skorupski, K., Mroczka, J., Wriedt, T., & Riefler, N. 2014, Physica A Statistical Mechanics and its Applications, 404, 106
  * Slobodrian et al. (2011) Slobodrian, R., Rioux, C., & Piche, M. 2011, Open Journal of Metal, 1, 7
  * Spitzer & Arny (1978) Spitzer, L. & Arny, T. T. 1978, American Journal of Physics, 46, 1201
  * Steinpilz et al. (2019) Steinpilz, T., Teiser, J., & Wurm, G. 2019, ApJ, 874, 60
  * Takigawa & Tachibana (2012) Takigawa, A. & Tachibana, S. 2012, ApJ, 750, 149
  * Tatsuuma & Kataoka (2021) Tatsuuma, M. & Kataoka, A. 2021, ApJ, 913, 132
  * Tatsuuma et al. (2019) Tatsuuma, M., Kataoka, A., & Tanaka, H. 2019, ApJ, 874, 159
  * Tazaki et al. (2017) Tazaki, R., Lazarian, A., & Nomura, H. 2017, ApJ, 839, 56
  * Tielens et al. (1994) Tielens, A. G. G. M., McKee, C. F., Seab, C. G., & Hollenbach, D. J. 1994, ApJ, 431, 321
  * Ulrich (2000) Ulrich, T. 2000, in Game Programming Gems, ed. M. DeLoura (Charles River Media), 444–453
  * Voshchinnikov & Henning (2010) Voshchinnikov, N. V. & Henning, T. 2010, A&A, 517, A45
  * Wada et al. (2007) Wada, K., Tanaka, H., Suyama, T., Kimura, H., & Yamamoto, T. 2007, ApJ, 661, 320
  * Wada et al. (1999) Wada, S., Kaito, C., Kimura, S., Ono, H., & Tokunaga, A. T. 1999, A&A, 345, 259
  * Weingartner & Draine (2001) Weingartner, J. C. & Draine, B. T. 2001, ApJ, 548, 296
  * Weingartner & Draine (2003) Weingartner, J. C. & Draine, B. T. 2003, ApJ, 589, 289
  * Williams & Jadwick (1980) Williams, R. J. & Jadwick, J. J. 1980, Handbook of lunar materials (Houston, TX, United States: NASA Reference Publication)
  * Wu et al. (2020) Wu, J., Faccinetto, A., Grimonprez, S., et al. 2020, Atmospheric Chemistry and Physics, 20, 4209
  * Yan & Lazarian (2003) Yan, H. & Lazarian, A. 2003, ApJ, 592, L33
  * Zebda et al. (2008) Zebda, A., Sabbah, H., Ababou-Girard, S., Solal, F., & Godet, C. 2008, Applied Surface Science, 254, 4980
  * Zhang et al. (2020) Zhang, J.-Y., Qi, H., Shi, J.-W., Gao, B.-H., & Ren, Y.-T. 2020, Opt. Express, 28, 37249
  * Zhukovska et al. (2008) Zhukovska, S., Gail, H. P., & Trieloff, M. 2008, A&A, 479, 453
  * Zhukovska & Henning (2013) Zhukovska, S. & Henning, T. 2013, A&A, 555, A99
  * Zhukovska et al. (2015) Zhukovska, S., Petrov, M., & Henning, T. 2015, ApJ, 810, 128
  * Zolensky et al. (2006) Zolensky, M. E., Zega, T. J., Yano, H., et al. 2006, Science, 314, 1735