# Discrete element method model of soot aggregates

Egor V. Demidov  Department of Chemistry and Environmental Science, New Jersey
Institute of Technology, 161 Warren Street, Newark, NJ 07103, United States of
America  Gennady Y. Gor  Department of Chemical and Materials Engineering, New
Jersey Institute of Technology, 161 Warren Street, Newark, NJ 07103, United
States of America  Alexei F. Khalizov111Corresponding author,
<EMAIL_ADDRESS>Department of Chemistry and Environmental Science, New
Jersey Institute of Technology, 161 Warren Street, Newark, NJ 07103, United
States of America  Department of Chemical and Materials Engineering, New
Jersey Institute of Technology, 161 Warren Street, Newark, NJ 07103, United
States of America

## Abstract

Soot is a component of atmospheric aerosols that affects climate by scattering
and absorbing the sunlight. Soot particles are fractal aggregates composed of
elemental carbon. In the atmosphere, the aggregates acquire coatings by
condensation and coagulation, resulting in significant compaction of the
aggregates that changes the direct climate forcing of soot. Currently, no
models exist to rigorously describe the process of soot restructuring,
reducing prediction accuracy of atmospheric aerosol models. We develop a
discrete element method contact model to simulate restructuring of fractal
soot aggregates, represented as assemblies of spheres joined by cohesion and
by sintered necks. The model is parametrized based on atomic force
spectroscopy data and is used to simulate soot restructuring, showing that the
fraction of necks in aggregates determines the restructuring pathway.
Aggregates with fewer necks undergo local compaction, while aggregates with
nearly-full necking prefer global compaction. Additionally, full compaction
occurs within tens of nanoseconds, orders of magnitude faster than the time
scale of soot aging through condensation. An important implication is that in
atmospheric soot aggregates, the rate of condensation determines how many
necks are fractured simultaneously, affecting the restructuring pathway, e.g.,
producing highly compact, thinly-coated soot as observed in recent studies.

##  1 Introduction

Atmospheric aerosols affect climate directly by scattering and absorbing light
[1] and indirectly by enhancing cloud formation [2]. One of the components of
atmospheric aerosols and a strong light absorber is soot [3]. Soot is formed
upon combustion of carbon-containing fuels in a multi-step process involving
thermolysis of the fuel, PAH formation, nucleation, and oxidation [4, 5]. As a
result, fractal-like structures consisting mostly of elemental carbon are
produced [6]. The micro-structure of carbonaceous phase can be either
amorphous, fullerenic, or a combination of both, depending on conditions,
under which the aggregate was formed [7].

While in the atmosphere, soot aggregates are subject to aging through a number
of pathways, including vapor condensation to form coatings around aggregates
[8]. In a process akin to adsorption-induced deformation [9, 10], capillary
forces induced by the liquid coatings cause initially fractal aggregates to
become more compact [11, 12]. This process of restructuring changes the
radiative forcing properties of the aggregates [13, 14, 15]. In order to
improve the representation of soot in aerosol models, it is necessary to
develop a model that can predict the morphology of soot aggregates based on
their surrounding conditions and time spent in the atmosphere. The
implications of mechanical properties of soot aggregates extend far beyond the
study of soot restructuring for environmental considerations. Carbon black
(CB), a material similar in nature to soot [16], is widely used as an additive
in manufacturing of rubbers, polyurethane coatings, oxide cathodes, inks, etc.
Understanding of mechanical properties of carbon black is important for its
industrial applications [17, 18].

The only attempt to model the mechanics of soot has been undertaken by
Schnitzler et al. [19], who studied the effect of surface tension of the
coating material on the extent of aggregate restructuring. Particle bonding is
not explicitly considered in that study; adhesive and liquid-originating
interactions are modeled with simple, pairwise force fields. Additionally,
Schnitzler’s model does not consider rotation of particles, hence it cannot
differentiate necked and point-touch contacts that are present in a real soot
aggregate (Figure 1). While Schnitzler’s model qualitatively supports the
discoveries made in their experimental study, it lacks an accurate
representation of necks that are usually present in fresh soot and resist the
restructuring process. Another noteworthy study is the discrete element method
(DEM) model for humidified TiO2 nanoparticle aggregate films under mechanical
load by Laube et al. [20], who consider the impact of capillary and solvation
forces that occur in pairs of wetted particles. Capillary and solvation forces
are modeled using force fields parametrized by the authors based on molecular
dynamics simulations. The model also considers friction and rotational degrees
of freedom of individual particles. Despite the rigor of Laube’s model, their
methodology is not directly transferable to modeling soot restructuring
because they assumed rigid aggregates where displacement of primary particles
in an aggregate relative to each other is prohibited, which is a valid
assumption for agglomerates of aggregates in their study, but does not allow
for any restructuring within aggregates by definition, and cannot be used in a
model for soot where all contacts are necked. Kelesidis et al. [21] used DEM
to model formation of soot by agglomeration of fine unit clusters into larger
spherules, which in turn agglomerate into fractal aggregates of primary
spherules. Kelesidis et al. did not consider the mechanical behavior of
aggregates.

Our goal is to develop a multi-scale framework for simulation of soot aging
and restructuring. Previously, our group introduced a technique to account for
the impact of non-spherical geometry of soot on vapor condensation kinetics
[22, 12]. Now, we present a DEM contact model for soot aggregate mechanics.
The contact model allows for creation of elastic bonds (necks) and non-bonded
contacts within a soot aggregate. Our goal was to create a model robust enough
to reproduce an experiment yet computationally efficient to be used to make
predictions by integrating it in larger aerosol models. The developed model is
used in this work to simulate an atomic force microscopy (AFM) spectroscopy
experiment [23], as well as aggregate restructuring by capillary forces (with
a simplistic representation of coating materials). The AFM force spectroscopy
simulation is used to parametrize and validate the model. Several
restructuring simulations are conducted to investigate the time scale of soot
restructuring, the impact of inter-particle bonds on the outcome of
restructuring, and to find an optimal parameter to track the progress of a
restructuring simulation.

##  2 Developed framework and methodology

###  2.1 Numerical model

Geometrically, soot aggregates can be approximated as assemblies of spherical
particles [24]. On the scale of this model, we disregard the micro-structure
of primary particles and assume that they are solid spheres. Hence, discrete
element method (DEM) is an appropriate simulation technique. DEM was
introduced by Cundall and Strack [25]. It is a numerical technique for
simulating mechanical behavior of granular systems. This method is based on
computing the resultant force, $\mathbf{f}$, and torque,
$\boldsymbol{\uptau}$, acting on each particle comprising an aggregate and
integrating Newton’s second law with respect to time to obtain the time-
evolution of the system of particles,

$m\ddot{\mathbf{x}}=\mathbf{f}$ (1)
$I\ddot{\boldsymbol{\uptheta}}=\boldsymbol{\uptau}$ (2)

where $m$ is mass, $I$ is moment of inertia, $\bf x$ is position, and
$\boldsymbol{\uptheta}$ is orientation of the particle. The key assumption
here that makes DEM different from other discrete techniques, such as
molecular dynamics, is that individual particles are not point masses and can
undergo meaningful rotation. Moment of inertia is $I=2/5mr^{2}$ for solid
spherical primary particles and $I=2/3mr^{2}$ for hollow spherical primary
particles (descriptive of internally oxidized soot [26]). Hence, torques and
rotational degrees of freedom need to be considered. The source of torques are
tractions applied to inter-particle contacts. The tractions are either a
result of deformation of an elastic inter-particle bond (in the case of a
particle pair joined with a neck, see Figure 1) or friction and cohesion (in
the case of a particle pair engaged in a point-touch contact, see Figure 1).
By using different contact force models, we can simulate aggregates with both
necked and point-touch contacts and vary the fraction of necks to see the
impact of necks on the outcome of restructuring. The model for force applied
to particles in contact determines whether the contacting pair behaves as a
bonded union or as non-bonded grains.

Figure 1: The two possible inter-particle contact types shown on a TEM image
of a soot aggregate. Particles that are bonded contain a solid neck between
the particles. Particles that are not bonded are touching at one point.

A non-bonded contact is a contact that does not have any tensile strength.
Such a contact only results in a force when two particles collide with each
other. A bonded contact is a contact that does have tensile strength. A pair
of particles in which a bonded contact has been established behaves
approximately as a rigid body. A small motion of one particle in the pair
relative to the other is resisted by the bonding forces. The resultant forces
seek to restore the un-deformed configuration of the pair. In addition to
contact forces, there are also body forces that can occur between particles
that are not in direct contact with each other, i.e. gravity or van der Waals
attraction. A body force is applied to the center of mass of a particle and
does not directly cause rotation. Gravity is not included in this model
because it is irrelevant at the scale of these simulations and cannot have an
impact on the restructuring process, because particles in an aggregate
accelerated by gravity retain constant velocity relative to each other. In the
next section, we develop the general concepts of the contact model that are
used in both the rigid bond model and the frictional contact model.

####  2.1.1 Development of the contact model

In a system of $N$ particles we can have up to $\frac{1}{2}N(N-1)$ binary
interactions. Let us pick two particles $i$ and $j$ from system and consider
their common reference frame. There are four types of relative motion that
particles $i$ and $j$ can undergo, corresponding to the four degrees of
freedom (DOFs) that we would like to control (Figure 2). These DOFs are (a)
normal translation, (b) torsion, (c) tangential translation, and (d) rolling.

Figure 2: Degrees of freedom of two particles in a common reference frame and
illustrations of springs that are used to constrain them

In order to define these degrees of freedom, consider a unit normal vector
$\mathbf{n}$ that is pointing from the center of particle $i$,
$\mathbf{x}_{i}$, towards the center of particle $j$, $\mathbf{x}_{j}$:

$\mathbf{n}=\frac{\mathbf{x}_{j}-\mathbf{x}_{i}}{\lVert\mathbf{x}_{j}-\mathbf{x}_{i}\rVert}$
(3)

Then the relative velocity at the contact point, $\mathbf{v}_{ij}$, is:

$\mathbf{v}_{ij}=\mathbf{v}_{j}-\mathbf{v}_{i}+\boldsymbol{\upomega}_{j}\times
a\mathbf{n}+\boldsymbol{\upomega}_{i}\times a\mathbf{n}$ (4)

where $\mathbf{v}_{i}$ and $\mathbf{v}_{j}$ are translational velocities of
the particles, $\boldsymbol{\upomega}_{i}$ and $\boldsymbol{\upomega}_{j}$ are
angular velocities of the particles, $r$ is particle radius (all particles are
assumed to have equal radii in this study), and $a$ is particle radius
corrected for overlap/separation $\delta$ between particles $i$ and $j$:

$a=r+\frac{1}{2}\delta$ (5)
$\delta=\lVert\mathbf{x}_{j}-\mathbf{x}_{i}\rVert-2r$ (6)

Then $\mathbf{v}_{ij}$ can be decomposed into normal and residual (tangential)
components:

$\mathbf{v}_{ij,\rm n}=\left(\mathbf{v}_{ij}\cdot\mathbf{n}\right)\mathbf{n}$
(7) $\mathbf{v}_{ij,\rm
t}=\mathbf{v}_{ij}-\left(\mathbf{v}_{ij}\cdot\mathbf{n}\right)\mathbf{n}$ (8)

Similarly, the relative angular velocity of the particles in the pair,

$\boldsymbol{\upomega}_{ij}=\boldsymbol{\upomega}_{j}-\boldsymbol{\upomega}_{i}$
(9)

can be decomposed into normal (torsional) and residual (rolling) components:

$\boldsymbol{\upomega}_{ij,\rm
n}=\left(\boldsymbol{\upomega}_{ij}\cdot\mathbf{n}\right)\mathbf{n}$ (10)
$\boldsymbol{\upomega}_{ij,\rm
r}=\boldsymbol{\upomega}_{ij}-\left(\boldsymbol{\upomega}_{ij}\cdot\mathbf{n}\right)\mathbf{n}$
(11)

Different subscripts are used for $\mathbf{v}_{ij,\rm t}$ and
$\boldsymbol{\upomega}_{ij,\rm r}$ to emphasize the fact that their directions
are not necessarily the same.

In order to control each DOF listed in Figure 2, we insert four springs
between particles $i$ and $j$. The rates of stretching/contraction of these
four (normal, torsional, tangential, and rolling) springs are given by
$\mathbf{v}_{ij,\rm n}$, $r\boldsymbol{\upomega}_{ij,\rm n}$,
$\mathbf{v}_{ij,\rm t}$, and $\boldsymbol{\upomega}_{ij,\rm r}\times
a\mathbf{n}$ respectively. Length of the normal spring can be computed
directly at every time step from positions of particles $i$ and $j$. As to the
remaining three springs, their lengths start at zero when the contact is
established and need to be incremented throughout the duration of the contact.
Let us define a generic term for a vector representing one of the three
springs (torsional, tangential, or rolling), $\boldsymbol{\upxi}$, and its
rate of change, $\dot{\boldsymbol{\upxi}}$. Then the algorithm at every time
step is as follows:

  1. 1.

Load the spring $\boldsymbol{\upxi}^{\prime}$ stored in memory since the
previous time step

  2. 2.

Compute $\boldsymbol{\upxi}$ by rotating $\boldsymbol{\upxi}^{\prime}$ to be
orthogonal to $\mathbf{n}$ and re-scaling the rotated spring to match the
magnitude of $\boldsymbol{\upxi}^{\prime}$

  3. 3.

Increment $\boldsymbol{\upxi}$ by $\dot{\boldsymbol{\upxi}}\Delta t$

  4. 4.

Use the value of $\boldsymbol{\upxi}$ to compute the forces and torques acting
on particles $i$ and $j$

  5. 5.

Store $\boldsymbol{\upxi}$ in the variable $\boldsymbol{\upxi}^{\prime}$ to be
used at the next iteration

####  2.1.2 Elastic bond model

For a pair of particles that is connected with an elastic bond (i.e., two
monomers connected by a neck), we would like to approximate a rigid–body
motion. In other words, the common reference frame of particles $i$ and $j$
can rotate and translate, but any translation or rotation of particle $i$
relative to particle $j$ should be restricted. The distance between all points
in the pair of particles should be approximately preserved over time. It can
be shown that when using the springs defined in the previous subsection to
restrict the motion of particles $i$ and $j$, then the union of particles $i$
and $j$ will undergo rigid body motion as the stiffness of inserted springs
approaches infinity. In the simulation we need to use a finite stiffness, but
as long as the amplitude of oscillations is much smaller than the length scale
of particles in the simulation, the motion will, approximately, be rigid.

To stabilize the system over time and dissipate any vibrational kinetic energy
in the bonds, each spring is supplemented by a dashpot element. The force,
$\mathbf{f}$, exerted on particle $i$ by each spring in the $ij$ bond is given
by:

$\mathbf{f}=k\boldsymbol{\upxi}+\eta\dot{\boldsymbol{\upxi}}$ (12)

where $k$ is stiffness and $\eta$ is the damping coefficient of the respective
spring. Forces exerted on particle $j$ are equal in magnitude and opposite in
direction. Forces arising from the normal and tangential springs are applied
to the particles. The force associated with the tangential spring will also
give rise to torques because the tangential force is not collinear with
$\mathbf{n}$. Forces computed from the torsion and rolling resistance springs
are quasi-forces that are not applied to particles $i$ and $j$, but are only
used to compute torques that will be applied to the particles.

Soot restructuring is limited by the strength of elastic bonds between primary
particles [27]. It is known that most-to-all neighboring primary particles in
a soot aggregate are initially bonded with rigid necks [28], but some of the
bonds break under the pressure applied by the AFM tip or the stress imposed by
the surface tension of a condensed liquid. In lieu of a bond strength
parametrization, we prescribe a certain fraction of broken bonds at the start
of every simulation. Then, aggregate fragments can fold, slide, and roll
relative to each other only at the non-bonded contacts. If the fraction of
broken bonds is zero, then all adjacent primary particles in the aggregate are
bonded with necks and no folding, sliding, or rolling is possible.

####  2.1.3 Frictional contact model

The model described by Luding [29] is used to simulate frictional contacts
between non-bonded particles. A brief description is provided here and the
reader is referred to Luding [29] for more details. Luding’s model uses the
same four springs described in Section 2.1.1 to compute normal and tangential
forces, rolling and torsion resistance torques. Instead of directly setting
force proportional to spring elongation, a certain degree of slip is allowed
between particles in contact. That is done by computing a test force,
$\mathbf{f}_{0}$, from current spring vector, $\boldsymbol{\upxi}$, and rate
of spring stretching, $\dot{\boldsymbol{\upxi}}$, using Equation 12. Then,
based on the magnitude of the test force, a decision is made whether static or
dynamic friction should be used,

$\text{if}\ \lVert\mathbf{f}_{0}\rVert\leq f_{\rm C,s}\ \text{use static
friction}$ (13) $\text{if}\ \lVert\mathbf{f}_{0}\rVert>f_{\rm C,s}\ \text{use
dynamic friction}$ (14)

where $f_{\rm C,s}$ is static friction force. Static and dynamic friction
forces are computed using Coulomb’s law of friction,

$f_{\rm C,s}=\mu_{\rm s}f_{\rm n}$ (15) $f_{\rm C,d}=\mu_{\rm d}f_{\rm n}$
(16)

where $\mu_{\rm s}$ is the static friction coefficient, $\mu_{\rm d}$ is the
dynamic friction coefficient, and $f_{\rm n}$ is the magnitude of the normal
force between the two particles.

In case the contact is determined to be in the state of static friction, the
tangential spring $\boldsymbol{\upxi}$ is incremented as described in
Subsection 2.1.1 and the test force computed using Equation 12 is used to
compute torques and, in the case of the tangential force, is also applied to
the contacting particles. In case the contact is determined to be in the state
of dynamic friction, the spring is not allowed to stretch anymore. Instead, a
certain extent of slipping is allowed. The spring to be used at the next
iteration, $\boldsymbol{\upxi}^{\prime}$, is set to

$\boldsymbol{\upxi}^{\prime}=-\frac{1}{k}\left(f_{\rm
C,d}\frac{\mathbf{f}_{0}}{\lVert\mathbf{f}_{0}\rVert}+\eta\dot{\boldsymbol{\upxi}}\right)$
(17)

and the magnitude of the Coulomb’s force, $f_{\rm C,d}$, is used to compute
torques and, if applicable, accelerate the particles in contact.

The model is only enabled when the normal force is repulsive. Once particles
are not overlapping, the frictional force is set to zero and accumulated
springs are reset.

####  2.1.4 Non-contact forces

Primary particle diameter of soot from different sources ranges between 10 to
50 nm [30, 31]. Primary particles that constitute soot aggregates in this
study are assumed to be 28 nm in diameter, which corresponds to the size of
primary particles particles in soot used by Enekwizu et al. [32] in an
experimental study of soot restructuring and by Chen et al. [12] and Ivanova
et al. [22] in the study of vapor condensation on soot. At that scale, van der
Waals attraction is significant between two particles that are in proximity.
Hamaker [33] derived the potential energy, $U$, due to van der Waals
attraction between two particles of radii $r_{1}$ and $r_{2}$ whose surfaces
are separated by distance $\delta$ from each other to be:

$U=-\frac{A}{6}\left[\frac{2r_{1}r_{2}}{(2r_{1}+2r_{2}+\delta)\delta}+\frac{2r_{1}r_{2}}{(2r_{1}+\delta)(2r_{2}+\delta)}+\ln\frac{(2r_{1}+2r_{2}+\delta)\delta}{(2r_{1}+\delta)(2r_{2}+\delta)}\right]$
(18)

Equation 18 simplifies when all particles are of equal radius $r$:

$U=-\frac{A}{6}\left[\frac{2r^{2}}{(4r+\delta)\delta}+\frac{2r^{2}}{(2r+\delta)^{2}}+\ln\frac{(4r+\delta)\delta}{(2r+\delta)^{2}}\right]$
(19)

where $A$ is the Hamaker constant – a material property. The magnitude of
force acting on the particles can be derived by differentiating potential
energy, $U$, with respect to separation distance, $\delta$, and the direction
of the force will coincide with the normal unit vector $\bf n$ defined
earlier:

$\mathbf{f}=-\frac{A}{6}\left[\frac{(4r+2\delta)}{(4r+\delta)\delta}-\frac{2}{(2r+\delta)}-\frac{4r^{2}}{(2r+\delta)^{3}}-\frac{2r^{2}(4r+2\delta)}{(4r+\delta)^{2}\delta^{2}}\right]\mathbf{n}$
(20)

Since in the limit as $\delta$ approaches $0$ the magnitude of force (and
potential energy) becomes infinite, a saturation distance $\delta_{0}$ is
introduced. The magnitude of $\delta_{0}$ varies between 0.4 and 1 nm [34].
Equation 20 is used to simulate van der Waals attraction between two primary
particles in the aggregate.

The attractive potential between a sphere and an infinite plane can be found
by taking the limit of Equation 18 as one of the particle radii approaches
infinity:

$U=-\frac{A}{6}\left[\frac{r}{\delta}+\frac{r}{2r+\delta}+\ln\frac{\delta}{2r+\delta}\right]$
(21)

The expression for attractive force acting on a particle in the vicinity of
the plane is then obtained by differentiation of Equation 21:

$\mathbf{f}=\frac{2A}{3}\left[\frac{r^{3}}{\delta^{2}\left(\delta+2r\right)^{2}}\right]\mathbf{n}$
(22)

Equation 22 is used to model attraction between a primary particle in the
aggregate and either the substrate or the AFM tip in force spectroscopy
simulations. We use $A=1.0\times 10^{-19}\ \rm J$ for particle-particle and
particle-substrate attraction, which is within the range of Hamaker constant
values for graphite reported in literature [35, 36]. We use $A=1.0\times
10^{-18}\ \rm J$ for particle-tip attraction to ensure rapid adhesion of the
aggregate to the tip during indentation.

####  2.1.5 Capillary force representation

It is well-known that capillary forces exerted by liquid menisci on the pore
walls inside porous materials cause adsorption-induced deformation [9, 10].
The same capillary force drives soot restructuring, except the menisci form
between primary particles in a fractal-like aggregate and not between the
walls of concave pores. To simulate the capillary force that occurs in pairs
of primary particles, we chose a sigmoid force-distance function:

$\mathbf{f}=\frac{1}{2}\psi\left[1-\tanh\left(\kappa\left(\lVert\mathbf{x}_{j}-\mathbf{x}_{i}\rVert-
d_{\mathrm{cutoff}}\right)\right)\right]$ (23)

where $\psi$ is the maximum force magnitude (related to surface tension of the
coating), $d_{\mathrm{cutoff}}$ is inter-particle distance at which the force
drops to zero (related to the amount of coating present on the aggregate), and
$\kappa$ is the sharpness of the drop that occurs at $d_{\mathrm{cutoff}}$. As
$\kappa$ becomes large, the expression in Equation 23 becomes a step function.
The reason why we chose a sigmoid dependence to simulate the capillary force
is that it is a good mathematical model for a short-range force. For example,
if $d_{\mathrm{cutoff}}=4r$, then only monomers that are separated by one
monomer can interact. We use $\psi=1.0\times 10^{-9}\ \rm N$ in this article,
so that the magnitude of capillary force is of the same order as reported in
literature [37]. In a more rigorous implementation, capillary force can be
calculated in the model from the shape of a meniscus formed between two
particles (i.e. catenoid or globoid [22]).

####  2.1.6 Generation of aggregates

Aggregation naturally produces self-similar shapes that can be conveniently
described through by the fractal scaling law [38]:

$N=k_{\rm f}\left(\frac{R_{\rm g}}{r}\right)^{D_{\rm f}}$ (24)

where $N$ is the number of primary particles in the aggregate, $k_{\rm f}$ is
the pre-exponential factor, $R_{\rm g}$ is the root mean square distance of
primary particles from the center of mass of the aggregate, and $D_{\rm f}$ is
the fractal dimension. The two fractal parameters, $k_{\rm f}$ and $D_{\rm
f}$, are determined experimentally from imaging [39] and are typically
characteristic for all aggregates in a sample, while $N$ and $R_{\rm g}$ vary
between different aggregates in one sample. Using the notion of fractals then
allows to numerically generate aggregates, geometries of which are
representative of real aggregates observed experimentally. Fresh soot
aggregates have a $D_{\rm f}$ of $1.8$. As the aggregate acquires a coating
and restructures, its $D_{\rm f}$ increases and approaches the limit of $\sim
2.2$ for thickly-coated aggregates [40].

A typical algorithm for numerical generation of fractal aggregates involves
sequentially combining primary particles into clusters and subsequently
combining sub-clusters into larger clusters in a way that the fractal scaling
law (Equation 24) holds until the desired aggregate size, $N$, is attained
[41]. By using this procedure to generate aggregates with $D_{\rm f}$ of
$1.8$, we get structures that are representative of real soot aggregates. In
this study, we used the cluster-cluster aggregation (CCA) algorithm described
by Filippov et al. and implemented by Skorupski et al. [42] as part of FLAGE
software.

##  3 Results and discussion

The contact model described in this article was implemented in an in-house
code. Technical details of the implementation are described in Supporting
Information. Each elementary interaction model (non-bonded contact, bonded
contact, and van der Waals attraction) were validated for compliance with
conservation of energy and momentum. Validation was carried out by setting up
simple simulations involving two or three particles with only the interaction
models of interest enabled. Conserved quantities were recorded throughout the
simulations. Descriptions of validation simulations and plots of conserved
quantities are available in Supporting Information.

To put the developed framework to a test, two types of simulations are set up
with the goal of qualitatively reproducing experimental results. The first
simulation involves only the mechanical properties of aggregates and
corresponds to the AFM indentation experiment conducted by Rong et al. [23]
and also in our lab, where soot samples were collected on a silicon substrate
and AFM indentation measurements were performed to determine the mechanical
properties of those aggregates. The second simulation is a simplistic
representation of restructuring of a thinly-coated soot aggregate, such as in
experiments by Chen et al. [43, 12] (Figure 7).

AFM deposition/indentation simulations were performed with an integration time
step of $5.0\times 10^{-14}\ \rm s$ for a total simulation duration of
$5.0\times 10^{-7}\ \rm s$. Restructuring simulations were performed with an
integration time step of $5.0\times 10^{-15}\ \rm s$ for a total simulation
duration of $5.0\times 10^{-8}\ \rm s$. Both AFM deposition/indentation and
restructuring simulations were performed with a neighbor list update period of
$20$ time steps, particle radius of $14\ \rm nm$, and particle material
density of $1700\ \rm kg/m^{3}$. Verlet radius used to build neighbor lists
was $4.2\times 10^{-8}\ \rm m$ for AFM indentation simulations and $7.0\times
10^{-8}\ \rm m$ for restructuring simulations. Parameters of the force models
used in each simulation type are listed in the diagram in Figure 3.

Figure 3: Force model parameters used in each simulation type

###  3.1 AFM indentation simulation

The AFM indentation simulation consists of two stages: deposition of an
aggregate onto a substrate and indentation of a deposited aggregate with a
solid tip.

####  3.1.1 Aggregate deposition

This simulation stage begins with an aggregate, where all neighboring primary
particles are bonded with necks. The substrate is introduced into the
simulation as a static plane. It can exert forces on the granular system, like
van der Waals attraction, normal force, and friction, but cannot accelerate
itself. AFM tip is not present at this simulation stage. The substrate plane
is normal to the z-axis and is located at $z=0$. A numerically generated
fractal aggregate is loaded and positioned such that the lowest monomer (with
respect to the z-axis) is only slightly above the substrate plane (Figure 4a).
Each primary particle is initialized with $1\ \rm m/s$ velocity in the
negative-z direction.

In the course of the simulation, the aggregate approaches the substrate plane
and collides with it. Adhesion occurs due to van der Waals attraction between
the substrate plane and the primary particles in the aggregate and viscous
dissipation of kinetic energy at the contact points during the collision of
the aggregate with the substrate. The aggregate may continue its motion
towards the substrate and rotate by inertia after one of the branches has
adhered to the substrate. The motion continues until adhesion has occurred at
multiple contact points (Figure 4b). Then, the system becomes strongly
anchored to the substrate.

Figure 4: Initial and final states of an aggregate deposition simulation. (a)
The aggregate is loaded, assigned an initial downward velocity, and positioned
above the substrate. (b) The aggregate then collides with the substrate and
adheres to it. A video file with a visualization of a deposition simulation is
available in Supporting Information

####  3.1.2 Aggregate indentation

This simulation stage consists of an aggregate that is deposited onto the
substrate. A preset fraction of neighboring primary particles are left with
necked bonds, and necked bonds are replaced with non-bonded VdW contacts
randomly in the remaining pairs. The substrate has the same representation and
position as in the deposition simulations. AFM tip is introduced as a
tetrahedron, consisting of solid triangle facets. Each triangle facet can
exert forces on the granular system, like van der Waals attraction, normal
force, and friction, but the facets themselves cannot accelerate. The triangle
facets move along a pre-determined trajectory.

The aggregate from the final data dump of a deposition simulation is loaded
and the necessary fraction of bonds is established. A primary particle that is
not hindered (along the z-axis) is randomly selected as the indentation point.
The AFM tip is then positioned with its lowest vertex $0.5r$ above the
selected primary particle and begins moving towards the primary particle
(Figure 5a). Soon after the tip has collided with the primary particle, it
reverses its direction of motion and starts retraction. The motion of the AFM
tip along the z-axis in indentation simulations is prescribed by following law
(in $\rm m/s$):

$\dot{z}_{\mathrm{tip}}=1.5\tanh\left[{10^{8}(t-1.08\times 10^{-8})}\right]$
(25)

The velocity profile has a sigmoid shape. The tip moves with a constant
velocity of $1.5\ \rm m/s$ towards the substrate, then at time $t=10^{-8}\ \rm
s$ it smoothly reverses its motion and starts moving away from the substrate
with a constant velocity of $1.5\ \rm m/s$. Due to strong adhesion, a branch
of the aggregate attaches to and follows the tip (Figure 5b). The force
exerted by the AFM tip on the granular system is recorded as a function of tip
height.

In the course of the simulation, a part of the aggregate separates from the
rest of the system along the broken necked bonds due to tension from the AFM
tip. The separated part collides with and engages in multiple frictional
contacts with the remaining aggregate pieces, until the piece adhered to the
AFM tip finally separates from the rest of the system (Figure 5c). These
repeated frictional contacts result in the characteristic “sawtooth” profile,
as reported by Rong et al. [23].

Figure 5: Initial, intermediate, and final states of an AFM indentation
simulation. The AFM tip is positioned above a deposited aggregate, (a) is
driven into the aggregate, and (b) retracted away from the aggregate. (c)
Tension causes a part of the aggregate to separate from the system. A video
file with a visualization of an indentation simulation is available in
Supporting Information

A simulated sawtooth curve is shown in Figure 6a. It was obtained from an
indentation simulation with the transfer function described in Section A.1
applied to smooth the recorded force. Figure 6b shows the simulated force-
displacement curve next to an example of our experimental curve. The force
spectroscopy experiment was performed similar to Rong et al. [23], with a
Bruker Dimension Icon AFM under ambient conditions using a cantilever with the
spring constant of $0.4\ \rm N/m$ and retraction velocity of $1\ \rm\mu m/s$.
Soot aggregates were electrostatically size-selected to $240\ \rm nm$, which
corresponds to about 120 primary particles per aggregate. A manuscript
describing the experimental AFM measurements is in preparation.

Figure 6: (a) A simulated force-displacement curve (with the transfer function
filter applied) and (b) the simulated force-displacement curve with a reduced
number of points plotted next to an example experimental force-displacement
curve

The simulated force-displacement curve has the expected sawtooth profile. The
magnitude of force peaks, their number, and the displacement at which the
aggregate is completely severed are comparable between experimental and
simulated data. The outcome of indentation simulations is mostly affected by
(1) initial morphology of the aggregate and (2) the Hamaker constant used to
model van der Waals attraction between monomers. Aggregate morphology controls
the number of peaks seen on a force-displacement curve, as the probability of
multiple collisions between branches occurring is higher in a larger and more
lacey aggregate. Hamaker constant controls the magnitude of the force peaks on
the force-displacement curve. The good agreement between simulations and
experiments can be explained by us using realistic fractal parameters (Section
2.1.6) and a Hamaker constant value for graphite consistent with literature
(Section 2.1.4).

###  3.2 Restructuring simulation

The simulation consists of an aggregate with necked connections between
primary particles and a capillary force acting on it. A preset fraction of
neighboring primary particles are left in non-bonded contacts with each other
to allow for restructuring. In the course of the simulation, the aggregate
undergoes compaction. The shortest branches collapse first, while longer
branches take time to accelerate. The change in morphology of an aggregate
during a restructuring simulation is presented in Figure 7.

Figure 7: Initial through final states of an aggregate restructuring
simulation. (a-d) Capillary forces cause short, peripheral branches to
collapse first. Backbone branches rearrange at a slower rate and the aggregate
reaches a stable configuration eventually, where all the non-bonded contacts
are locked. A video file with a visualization of a restructuring simulation is
available in Supporting Information

Average coordination number of primary particles in the aggregate was chosen
as the metric to track the progress of restructuring. It is defined as the
average number of neighbors that a primary particle has and is computed as,

$c=\frac{1}{N}\sum_{i=1}^{N}N_{\mathrm{nb},i}$ (26)

where $c$ is the average coordination number of the aggregate and
$N_{\mathrm{nb},i}$ is the number of neighbors that particle $i$ has. For the
purpose of finding $N_{\mathrm{nb},i}$, we consider as neighbors two
particles, which are separated by less that $1\%$ of primary particle
diameter. In restructuring simulations, 50 different aggregates with the same
size (150 primary particles) and fractal parameters ($1.8$ fractal dimension
and $1.4$ pre-exponential factor) were restructured with the same coating
force for the same duration. Average coordination numbers of aggregates with
various necking fractions are presented in Figure 8a,b. Coordination number is
descriptive of local morphological changes in an aggregate (inset in Figure
8b). We also compute convexity of the restructuring aggregates, which in three
dimensions is defined as the ratio of the volume of an aggregate to the volume
of the convex hull constructed around the aggregate [44]. Convexity is
descriptive of global morphological changes in the aggregate (inset in Figure
8d). Convexity data is presented in Figure 8c,d.

Figure 8: Evolution in morphology of soot aggregates undergoing restructuring
due to capillary force. (a) Aggregate-average coordination number of 20
restructuring aggregates as a function of time, (b) average terminal
coordination number as a function of aggregate necking fraction. The dashed
line corresponds to the limiting coordination number of 6.0 for bulk random-
jammed packing of frictionless spheres [45]. (c) Average convexity of 50
restructuring aggregates as a function of time, (d) average terminal convexity
as a function of aggregate necking fraction. Inset in sub-plot b with an arrow
and aggregates illustrates local compaction and inset in sub-plot d
illustrates global compaction. Shaded regions and error bars are 95%
confidence interval

In all simulations, the coordination number is a good indicator of
restructuring progress. As it reaches its maximum value, restructuring slows
down and stops. Maximum coordination number correlates with the fraction of
necks present in the aggregate. Fewer necks allow for more compaction and a
higher coordination number is attained by the aggregate at the end of
restructuring. A slight increase in coordination number (under $10\%$) is
observed for aggregates with $100\%$ necking due to elastic deformation of
necks. Aggregates with no necks attain the highest coordination number,
$c=4.48$. This is well below the limiting coordination number of $6.0$ for
bulk random-jammed frictionless packing, due to inter-particle friction being
enabled and a large fraction of particles being on the surface of the
aggregate [45]. Full compaction occurs within tens of nanoseconds. The
timescale of restructuring depends mostly on the magnitude of the driving
force, which was selected by us based on literature (Section 2.1.5). Convexity
takes longer to reach steady state and exhibits more variance between
different aggregates. Coordination number reaches steady state sooner than
convexity because it characterizes local compactness of the aggregate and
local compaction occurs faster than global movement of long aggregate
branches.

###  3.3 Limitations and potential future developments

One of the limitations of the capillary force model used in this study is that
it does not take into account hysteresis related to formation / breakage of
liquid bridges [46]. This capillary force model cannot be used to model
thickly coated soot, where the liquid layer covers the entire aggregate.
Particularly, in the latter case, the moving air-liquid interface of an
evaporating coating droplet pushes all particles towards the center of the
droplet instead of a causing particle-particle attraction, as described by the
present model. Future work will involve development of a more elaborate
representation of coating materials and parametrization of the restructuring
model based on experiments. An algorithm for brittle fracture of necks will be
implemented and neck strengths will be parametrized. With those additions
implemented on top of the DEM contact model presented in this study, the
restructuring model in conjunction with the capillary condensation model
developed by our group will be used to study evolution of soot composition,
morphology, and optical properties under atmospheric conditions. The DEM
contact model will also be applied to simulating mechanical properties of
industrial carbon blacks incorporated into various media.

##  4 Conclusion

We have developed a DEM contact model for mechanics of soot aggregates.
Elementary interaction models were validated for consistency. The contact
model was parametrized based on AFM spectroscopy experiments and literature
data on carbon Hamaker constant. Restructuring simulations were conducted
using a simplistic model for capillary force. The magnitude of capillary force
used to drive restructuring in this study was on the order of nano-Newtons,
which is comparable to values reported in literature for similar systems [37].
With this order of capillary force, the time scale of restructuring is tens of
nanoseconds. Hence, the rate of restructuring is much faster than the rate of
condensation in saturator-based experiments (milliseconds [22, 12]), chamber-
based experiments (minutes [47, 48]), and the atmosphere (hours [49, 50, 51]).
Therefore, the rate of soot compaction is controlled by the rate of coating
formation by vapor condensation, not the rate of restructuring. This means
that condensation and restructuring simulations can be decoupled, which
significantly simplifies the computational burden for future condensation-
driven restructuring modeling.

Simulations also show that the outcome of the restructuring process is
controlled by the number of necks present in the aggregate. In an aggregate
with near-full necking, where only a few necks are broken, global
restructuring occurs due to movement of bonded branches as a whole. In an
aggregate with most of the necks broken, local restructuring is prevalent,
because exterior particles quickly collapse onto backbone branches, interlock,
and inhibit further restructuring. The fact that major restructuring can occur
with only one or two necks broken in the aggregate can explain the highly
compact thinly-coated soot, where the amount of condensate is only sufficient
to break a few necks, observed in field measurements [52, 53, 54] and in
laboratory studies [55]. Likewise, minor restructuring of thickly coated soot,
where the amount of condensate is sufficient to break most necks in the
aggregate, observed in some studies [56, 57], can be attributed to fast
condensation followed by interlocking of primary particles.

## Authors’ contributions

Egor V. Demidov: Software, Writing – review & editing, Writing – original
draft, Visualization, Validation, Methodology, Investigation, Formal analysis,
Conceptualization. Gennady Y. Gor: Writing – review & editing, Supervision,
Resources, Project administration, Methodology, Funding acquisition, Formal
analysis, Conceptualization. Alexei F. Khalizov: Writing – review & editing,
Supervision, Resources, Project administration, Methodology, Funding
acquisition, Formal analysis, Conceptualization.

## Data availability

Supporting information provided with this manuscript consists of a PDF
document with a description of the implementation details of the model and
validation steps undertaken; and an archive with videos of select deposition,
indentation, and restructuring simulations; an archive with input files that
can be used to reproduce all data used in this article. Software is available
in a GitHub repository. A tutorial on how to run restructuring simulations
using our software is available on Egor Demidov’s web page.

## Acknowledgment

A.K., G.G., and E.D. acknowledge the U.S. National Science Foundation, grant
AGS-2222104. We thank Ali Hasani for the experimental force-displacement curve
used in Figure 6 and for the TEM image of a soot aggregate used in Figure 1.

## Appendix

###  A.1 AFM transfer function

A filter is needed to reduce the high-frequency noise on simulated force-
displacement curves. The real AFM cantilever is a damped oscillator system
[58] with nonzero mass. Due to the inertia of the cantilever, an infinitesimal
force impulse applied to the cantilever will not cause it to deflect by a
detectable distance and will not be measured by the instrument. Here, we
derive a transfer function that relates the applied force to the force
measured by the AFM instrument.

Let the AFM cantilever be a second-order system:

$m_{\rm c}\ddot{u}+\eta_{\rm c}\dot{u}+k_{\rm c}u=f(t)$ (27)

where $m_{\rm c}$ is the effective mass of the cantilever, $\eta_{\rm c}$ is
the damping coefficient of the cantilever, $k_{\rm c}$ is the effective
stiffness of the cantilever, $u$ is the displacement of the cantilever from
the equilibrium position, and $f(t)$ is the external force applied to the
cantilever. The instrument determines force by measuring $u$ by detecting
deflection of a laser beam that is being reflected from the cantilever. Then
instrument software computes the apparent force using Hooke’s law, $g=k_{\rm
c}u$, and records it. We can establish a relationship between the true force
applied to the AFM probe and the apparent force that is recorded by the
instrument as:

$\frac{m_{\rm c}}{k_{\rm c}}\ddot{g}+\frac{\eta_{\rm c}}{k_{\rm
c}}\dot{g}+g(t)=f(t)$ (28)

Assuming that the cantilever is critically damped,

$\eta_{\rm c}=2\sqrt{m_{\rm c}k_{\rm c}}$ (29)

and by taking a Laplace transform of Equation 28, we can arrive at the
expression for the transfer function of the instrument:

$\frac{G(s)}{F(s)}=\frac{\omega_{0}^{2}}{\left(\omega_{0}+s\right)^{2}}$ (30)

where $\omega_{0}$ is the natural angular frequency of the undamped oscillator
given by:

$\omega_{0}=\sqrt{\frac{k_{\rm c}}{m_{\rm c}}}$ (31)

By exploring the limits of the transfer function (Equation 30), we can see
that for a large value of $\omega_{0}$, the input is preserved, and for a
small value of $\omega_{0}$, the input is smoothed.

$\lim_{\omega_{0}\rightarrow\infty}\frac{\omega_{0}^{2}}{(s+\omega_{0})^{2}}=1\quad\Rightarrow\quad\text{input
is preserved for large $\omega_{0}$}$ (32) $\lim_{\omega_{0}\rightarrow
0}\frac{\omega_{0}^{2}}{(s+\omega_{0})^{2}}=0\quad\Rightarrow\quad\text{input
is attenuated for small $\omega_{0}$}$ (33)

## References

  * [1] Petr Chylek and J Wong  “Effect of absorbing aerosols on global radiation budget”  In _Geophysical Research Letters_ 22.8 Wiley Online Library, 1995, pp. 929–931
  * [2] Ulrike Lohmann and Johann Feichter  “Global indirect aerosol effects: a review”  In _Atmospheric Chemistry and Physics_ 5.3 Copernicus GmbH, 2005, pp. 715–737
  * [3] James Haywood and Olivier Boucher  “Estimates of the direct and indirect radiative forcing due to tropospheric aerosols: A review”  In _Reviews of Geophysics_ 38.4 Wiley Online Library, 2000, pp. 513–543
  * [4] Murray Thomson and Tirthankar Mitra  “A radical approach to soot formation”  In _Science_ 361.6406 American Association for the Advancement of Science, 2018, pp. 978–979
  * [5] KO Johansson et al.  “Resonance-stabilized hydrocarbon-radical chain reactions may explain soot inception and growth”  In _Science_ 361.6406 American Association for the Advancement of Science, 2018, pp. 997–1000
  * [6] Benjamin Demirdjian et al.  “Heterogeneities in the microstructure and composition of aircraft engine combustor soot: impact on the water uptake”  In _Journal of Atmospheric Chemistry_ 56 Springer, 2007, pp. 83–103
  * [7] William J Grieco, Jack B Howard, Lenore C Rainey and John B Vander Sande  “Fullerenic carbon in combustion-generated soot”  In _Carbon_ 38.4 Elsevier, 2000, pp. 597–614
  * [8] H Saathoff et al.  “Coating of soot and (NH4)2SO4 particles by ozonolysis products of $\alpha$-pinene”  In _Journal of Aerosol Science_ 34.10 Elsevier, 2003, pp. 1297–1321
  * [9] Gennady Y Gor, Patrick Huber and Noam Bernstein  “Adsorption-induced deformation of nanoporous materials—A review”  In _Applied Physics Reviews_ 4.1 AIP Publishing, 2017
  * [10] Gennady Y. Gor and A.. Kolesnikov  “What Drives Deformation of Smart Nanoporous Materials During Adsorption and Electrosorption?”  In _Langmuir_ ACS Publications, 2024, pp. DOI: 10.1021/acs.langmuir.4c00443
  * [11] S Kütz and A Schmidt-Ott  “Characterization of agglomerates by condensation-induced restructuring”  In _Journal of Aerosol Science_ 23 Elsevier, 1992, pp. 357–360
  * [12] Chao Chen et al.  “Single parameter for predicting the morphology of atmospheric black carbon”  In _Environmental Science & Technology_ 52.24 ACS Publications, 2018, pp. 14169–14179
  * [13] CM Sorensen  “Light scattering by fractal aggregates: a review”  In _Aerosol Science & Technology_ 35.2 Taylor & Francis, 2001, pp. 648–687
  * [14] Egor V Demidov et al.  “Differences and similarities in optical properties of coated fractal soot and its surrogates”  In _Journal of Aerosol Science_ 180 Elsevier, 2024, pp. 106392
  * [15] Alexei F Khalizov et al.  “Enhanced light absorption and scattering by carbon soot aerosol internally mixed with sulfuric acid”  In _The Journal of Physical Chemistry A_ 113.6 ACS Publications, 2009, pp. 1066–1074
  * [16] Christopher M Long, Marc A Nascarella and Peter A Valberg  “Carbon black vs. black carbon and other airborne materials containing elemental carbon: Physical and chemical distinctions”  In _Environmental Pollution_ 181 Elsevier, 2013, pp. 271–286
  * [17] Hugh M Smallwood  “Limiting law of the reinforcement of rubber”  In _Journal of Applied Physics_ 15.11 American Institute of Physics, 1944, pp. 758–766
  * [18] Avrom I Medalia  “Effect of carbon black on ultimate properties of rubber vulcanizates”  In _Rubber Chemistry and Technology_ 60.1, 1987, pp. 45–61
  * [19] Elijah G Schnitzler, Jakub M Gac and Wolfgang Jäger  “Coating surface tension dependence of soot aggregate restructuring”  In _Journal of Aerosol Science_ 106 Elsevier, 2017, pp. 43–55
  * [20] Jens Laube et al.  “A new contact model for the discrete element method simulation of TiO2 nanoparticle films under mechanical load”  In _Granular Matter_ 20.2 Springer, 2018, pp. 28
  * [21] Georgios A Kelesidis, Simon Benz and Sotiris E Pratsinis  “Process design for carbon black size and morphology”  In _Carbon_ 213 Elsevier, 2023, pp. 118255
  * [22] Ella V Ivanova, Alexei F Khalizov and Gennady Y Gor  “Kinetic model for competitive condensation of vapor between concave and convex surfaces in a soot aggregate”  In _Aerosol Science and Technology_ 55.3 Taylor & Francis, 2020, pp. 302–315
  * [23] Weizhi Rong et al.  “Complementary TEM and AFM force spectroscopy to characterize the nanomechanical properties of nanoparticle chain aggregates”  In _Nano Letters_ 4.11 ACS Publications, 2004, pp. 2287–2292
  * [24] Patience Gwaze et al.  “Comparison of three methods of fractal analysis applied to soot aggregates from wood combustion”  In _Journal of Aerosol Science_ 37.7 Elsevier, 2006, pp. 820–838
  * [25] Peter A Cundall and Otto DL Strack  “A discrete numerical model for granular assemblies”  In _Geotechnique_ 29.1 Thomas Telford Ltd, 1979, pp. 47–65
  * [26] Georgios A Kelesidis, Patrizia Crepaldi and Sotiris E Pratsinis  “Oxidation dynamics of soot or carbon black accounting for its core-shell structure and pore network”  In _Carbon_ 219 Elsevier, 2024, pp. 118764
  * [27] Alfred P Weber and Sheldon K Friedlander  “In situ determination of the activation energy for restructuring of nanometer aerosol agglomerates”  In _Journal of Aerosol Science_ 28.2 Elsevier, 1997, pp. 179–192
  * [28] Sonja Rothenbacher, Armin Messerer and Gerhard Kasper  “Fragmentation and bond strength of airborne diesel soot agglomerates”  In _Particle and Fibre Toxicology_ 5 Springer, 2008, pp. 1–7
  * [29] Stefan Luding  “Cohesive, frictional powders: contact models for tension”  In _Granular Matter_ 10.4 Springer, 2008, pp. 235–246
  * [30] Kouji Adachi, Serena H Chung, Heiner Friedrich and Peter R Buseck  “Fractal parameters of individual soot particles determined using electron tomography: Implications for optical properties”  In _Journal of Geophysical Research: Atmospheres_ 112.D14 Wiley Online Library, 2007
  * [31] M Wentzel et al.  “Transmission electron microscopical and aerosol dynamical characterization of soot aerosols”  In _Journal of Aerosol Science_ 34.10 Elsevier, 2003, pp. 1347–1370
  * [32] Ogochukwu Y Enekwizu, Ali Hasani and Alexei F Khalizov  “Vapor condensation and coating evaporation are both responsible for soot aggregate restructuring”  In _Environmental Science & Technology_ 55.13 ACS Publications, 2021, pp. 8622–8630
  * [33] Hugo C Hamaker  “The London—van der Waals attraction between spherical particles”  In _Physica_ 4.10 Elsevier, 1937, pp. 1058–1072
  * [34] MB Ranade  “Adhesion and removal of fine particles on surfaces”  In _Aerosol Science and Technology_ 7.2 Taylor & Francis, 1987, pp. 161–176
  * [35] Brad A Krajina, Lakshmi S Kocherlakota and René M Overney  “Direct determination of the local Hamaker constant of inorganic surfaces based on scanning force microscopy”  In _The Journal of Chemical Physics_ 141.16 AIP Publishing, 2014
  * [36] Panagiotis Tolias  “Lifshitz calculations of Hamaker constants for fusion relevant materials”  In _Fusion Engineering and Design_ 133 Elsevier, 2018, pp. 110–116
  * [37] Mahdi Farshchi-Tabrizi et al.  “On the adhesion between fine particles and nanocontacts: an atomic force microscope study”  In _Langmuir_ 22.5 ACS Publications, 2006, pp. 2171–2184
  * [38] R Jullien  “Aggregation phenomena and fractal aggregates”  In _Contemporary Physics_ 28.5 Taylor & Francis, 1987, pp. 477–493
  * [39] Magin Lapuerta, Rosario Ballesteros and Francisco J Martos  “A method to determine the fractal dimension of diesel soot agglomerates”  In _Journal of Colloid and Interface Science_ 303.1 Elsevier, 2006, pp. 149–158
  * [40] Yuanyuan Wang et al.  “Fractal dimensions and mixing structures of soot particles during atmospheric processing”  In _Environmental Science & Technology Letters_ 4.11 ACS Publications, 2017, pp. 487–493
  * [41] AV Filippov, Mauricio Zurita and Daniel E Rosner  “Fractal-like aggregates: Relation between morphology and physical properties”  In _Journal of Colloid and Interface Science_ 229.1 Elsevier, 2000, pp. 261–273
  * [42] Krzysztof Skorupski, Janusz Mroczka, Thomas Wriedt and Norbert Riefler  “A fast and accurate implementation of tunable algorithms used for generation of fractal-like aggregate models”  In _Physica A: Statistical Mechanics and its Applications_ 404 Elsevier, 2014, pp. 106–117
  * [43] Chao Chen et al.  “An unexpected restructuring of combustion soot aggregates by subnanometer coatings of polycyclic aromatic hydrocarbons”  In _Geophysical Research Letters_ 43.20 Wiley Online Library, 2016, pp. 11–080
  * [44] Susan Hert and Stefan Schirra  “3D Convex Hulls”  In _CGAL User and Reference Manual_ CGAL Editorial Board, 2024  URL: https://doc.cgal.org/5.6.1/Manual/packages.html#PkgConvexHull3
  * [45] Leonardo E Silbert  “Jamming of frictional spheres and random loose packing”  In _Soft Matter_ 6.13 Royal Society of Chemistry, 2010, pp. 2918–2924
  * [46] CD Willett, MJ Adams, SA Johnson and JPK Seville  “Effects of wetting hysteresis on pendular liquid bridges between rigid spheres”  In _Powder Technology_ 130.1-3 Elsevier, 2003, pp. 63–69
  * [47] Cerina Wittbom et al.  “Cloud droplet activity changes of soot aerosol upon smog chamber ageing”  In _Atmospheric Chemistry and Physics_ 14.18 Copernicus GmbH, 2014, pp. 9831–9854
  * [48] Alexei F Khalizov et al.  “Role of OH-initiated oxidation of isoprene in aging of combustion soot”  In _Environmental Science & Technology_ 47.5 ACS Publications, 2013, pp. 2254–2263
  * [49] N Riemer, H Vogel and B Vogel  “Soot aging time scales in polluted regions during day and night”  In _Atmospheric Chemistry and Physics_ 4.7 Copernicus GmbH, 2004, pp. 1885–1893
  * [50] J Zhang, J Liu, S Tao and GA Ban-Weiss  “Long-range transport of black carbon to the Pacific Ocean and its dependence on aging timescale”  In _Atmospheric Chemistry and Physics_ 15.20 Copernicus GmbH, 2015, pp. 11521–11535
  * [51] Nicole Riemer, Matthew West, Rahul Zaveri and Richard Easter  “Estimating black carbon aging time-scales with a particle-resolved aerosol model”  In _Journal of Aerosol Science_ 41.1 Elsevier, 2010, pp. 143–158
  * [52] Janarjan Bhandari et al.  “Extensive soot compaction by cloud processing from laboratory and field observations”  In _Scientific Reports_ 9.1 Nature Publishing Group UK London, 2019, pp. 11824
  * [53] Swarup China et al.  “Morphology and mixing state of aged soot particles at a remote marine free troposphere site: Implications for optical properties”  In _Geophysical Research Letters_ 42.4 Wiley Online Library, 2015, pp. 1243–1250
  * [54] Yuanyuan Wang et al.  “Variability in morphology of soot particles during non-cloud and in-cloud processes”  In _Atmospheric Environment_ 326 Elsevier, 2024, pp. 120489
  * [55] Ray P Bambha, Mark A Dansson, Paul E Schrader and Hope A Michelsen  “Effects of volatile coatings and coating removal mechanisms on the morphology of graphitic soot”  In _Carbon_ 61 Elsevier, 2013, pp. 80–96
  * [56] Eben S Cross et al.  “Soot particle studies—instrument inter-comparison—project overview”  In _Aerosol Science and Technology_ 44.8 Taylor & Francis, 2010, pp. 592–611
  * [57] R Ghazi and JS Olfert  “Coating mass dependence of soot aggregate restructuring due to coatings of oleic acid and dioctyl sebacate”  In _Aerosol Science and Technology_ 47.2 Taylor & Francis, 2013, pp. 192–200
  * [58] Hans-Jürgen Butt, Brunero Cappella and Michael Kappl  “Force measurements with the atomic force microscope: Technique, interpretation and applications”  In _Surface Science Reports_ 59.1-6 Elsevier, 2005, pp. 1–152