myprojects-hg-reborn

Fork

%\documentclass{elsart}
\documentclass[preprint,12pt]{elsarticle}





\usepackage{graphicx}% Include figure files
% \usepackage{dcolumn}% Align table columns on decimal point
% \usepackage{bm}% bold math

\usepackage{url}
\usepackage{natbib}
 \usepackage{amsmath}
 \usepackage{amssymb}

\newcommand{\Teleia}{\quad .}
\newcommand{\Comma}{\quad ,}

%\newcommand{\ol}{\overline}
\renewcommand{\i}{\int}
\newcommand{\n}{\nabla}
\newcommand{\x}{\vec x\; }
\renewcommand{\d}{\dag}
\newcommand{\h}{\hat}
\newcommand{\p}{\partial}
\renewcommand{\v}{\vert}
\renewcommand{\l}{\langle}
\renewcommand{\r}{\rangle}
\newcommand{\f}{\frac}
\newcommand{\s}{\sum}
\newcommand{\lm}[1]{\lim_{#1\to\infty}}
\renewcommand{\in}{\infty}
\newcommand{\rro}{\right)}
\newcommand{\lro}{\left( }
\newcommand{\lsq}{\left[}
\newcommand{\rsq}{\right]}
\newcommand{\rcu}{\right\}}
\newcommand{\lcu}{\left\{}
\newcommand{\be}{\begin{equation}}
\newcommand{\ee}{\end{equation}}
\newcommand{\bi}{\begin{itemize}}
\newcommand{\ei}{\end{itemize}}
\newcommand{\ben}{\begin{enumerate}}
\newcommand{\een}{\end{enumerate}}
\newcommand{\esp}{ESPResSo }
\newcommand{\rmin}{r_{\textrm{min}}}
\newcommand{\rcut}{r_{\textrm{cut}}}
\newcommand{\umin}{u_{\textrm{min}}}
\newcommand{\usigma}{u_{\sigma}}
\newcommand{\umod}{u_{\textrm{mod}}}

\newcommand{\pra}{{\it Physical Review A}}
\newcommand{\prb}{{\it Physical Review B}}
\newcommand{\prl}{{\it Physical Review Letters}}

\newcommand{\jc}{{\it Journal of Colloid and Interface Science}}
\newcommand{\jas}{{\it Journal of Aerosol Science}}
%\newcommand{\pra}{{\it Physical Review A}}
%\newcommand{\prb}{{\it Physical Review B}}
%\newcommand{\pre}{{\it Physical Review E}}
%\newcommand{\prl}{{\it Physical Review Letters}}

%Fine preambolo

\newcommand{\unit}{\hat{\bf n}}
%  \newcommand{\pol}{\hat{\bf e}}
\newcommand{\rv}{{\bf r}}
\newcommand{\Ev}{{\bf E}}
\newcommand{\Bv}{{\bf B}}
\newcommand{\Ec}{{\cal E}}
\newcommand{\Rc}{{\cal R}}
\newcommand{\Pc}{{\cal P}}
\newcommand{\Pcv}{\bbox {\cal P}}
\newcommand{\dv}{{\bf d}}
\newcommand{\Dc}{{\cal D}}
\newcommand{\Dcv}{\bbox {\cal D}}
\newcommand{\Hc}{{\cal H}}
\newcommand{\kappav}{\bbox \kappa}
\newcommand{\Dkappav}{\Delta {\bbox\kappa}}
\newcommand{\qv}{{\bf q}}
\newcommand{\kv}{{\bf k}}
\newcommand{\eo}{\epsilon_0}
\newcommand{\ej}{\epsilon_j}
\newcommand{\beq}{\begin{equation}}
\newcommand{\eeq}{\end{equation}}
\newcommand{\bea}{\begin{eqnarray}}
\newcommand{\eea}{\end{eqnarray}}
\newcommand{\up}{\uparrow}
\newcommand{\down}{\downarrow}
\newcommand{\<}{\langle}
\renewcommand{\>}{\rangle}
\renewcommand{\(}{\left(}
\renewcommand{\)}{\right)}
\renewcommand{\[}{\left[}
\renewcommand{\]}{\right]}
\newcommand{\dagg}{d_{\rm{agg}}}
\newcommand{\vagg}{V_{\rm{agg}}}
\newcommand{\nagg}{n_{\rm{agg}}}
\newcommand{\df}{d_{f}}
\newcommand{\ragg}{\rho_{\rm{agg}}}
\newcommand{\reff}{\rho_{\rm{eff}}}
\newcommand{\re}{{\rm{Re}}}
\newcommand{\pr}{{\rm{Pr}}}
\newcommand{\sh}{{\rm{Sh}}}
\newcommand{\Kn}{{\rm{Kn}}}
\newcommand{\ra}{{\rm{Ra}}}
\renewcommand{\sc}{{\rm{Sc}}}
\newcommand{\nusselt}{{\rm{Nu}}}
\newcommand{\magg}{m_{\rm{agg}}}
\newcommand{\tres}{\tau_{\rm{res}}}
\newcommand{\gdif}{{\gamma_{\rm{dif}}}}
\newcommand{\vdep}{{v_{\rm{deb}}}}
\newcommand{\gth}{{\gamma_{\rm{th}}}}
\newcommand{\vth}{{v_{\rm{th}}}}
% \newcommand{\tres}{{\tau_{\rm{res}}}}
\newcommand{\kt}{{K_{\rm{T}}}}
\newcommand{\kair}{{k_{\rm{air}}}}
\newcommand{\vdif}{{v_{\rm{dif}}}}
\newcommand{\kp}{{k_{\rm{p}}}}
\newcommand{\commentout}[1]{{}}
%\newcommand{\half}{\hbox}
% \newcommand{\half}{\hbox{$1\over2$}}
\newcommand{\nv}{{\vec\nabla}}
\renewcommand{\c}{{\cdot}}
\newcommand{\hv}{\harvarditem}
   

\journal{Journal of Colloid and Interface Science}

         
\begin{document}    

\begin{frontmatter}     

\title{Nanoparticle collisional dynamics by Langevin simulations}

\author[jrc]{Lorenzo Isella}\fnref{LorenzoISI}
\ead{lorenzo.isella@gmail.com}
\fntext[LorenzoISI]{Present address: ISI Foundation, Turin 10133,
Italy}

\author[jrc,unewcastle]{Yannis Drossinos\corref{dross}}
\ead{ioannis.drossinos@jrc.it}
\address[jrc]{European Commission Joint Research Centre}
\address[unewcastle]{School of Mechanical \& Systems Engineering,
Newcastle University, Newcastle upon Tyne NE1 7RU, United Kingdom}
\cortext[dross]{Corresponding author}

%\author{L. Isella\corref{isella}\fnref{isellaref}}
%\author{L. Isella\fnref{isellaref}}
%\ead{lorenzo.isella@gmail.com}
%\address[isellaref]{ISI Foundation, Turin 10133, Italy}

%\author{Y. Drossinos\corref{dross}\fnref{drossref}}
%\ead{ioannis.drossinos@jrc.it}
%\address[drossref]{European Commission Joint Research Centre\fnref{label3}}
%\cortext[dross]{Corresponding author}

\date{\today;
\texttt{agglo\_ydv5a.tex}}


\begin{abstract}
Nanoparticle agglomeration as a function of monomer-monomer
interaction potential is simulated by solving the Langevin
equations of motion for a set of interacting monomers in three dimensions.
The numerical simulations are analyzed to investigate the
static structure of the generated clusters, their dynamic transport
properties, and the collision frequency kernel between small clusters.
Cluster restructuring is also observed and discussed.
We identify a time-dependent average fractal dimension whose time
evolution is linked
to the kinetics of two cluster populations. The absence of shielding
in Langevin equations is discussed and its effect on
cluster translational and rotational properties is quantified.
\end{abstract}

\begin{keyword}
% keywords here, in the form: keyword \sep keyword
Langevin, fractal, aggregate, agglomeration.
% PACS codes here, in the form: \PACS code \sep code
\PACS 
\end{keyword}
\end{frontmatter}



% \keywords{Langevin equation, fractal dimension, aggregate,
% agglomeration.}
% PACS codes here, in the form: \PACS code \sep code
%\pacs{}

%\maketitle

% main text
\section{Introduction}
\label{introduction}

Nanoparticle aggregates are of significant importance in
technological and industrial processes such as, for example,
combustion, filtration, and gas-phase particle synthesis.
The fractal nature of these aggregates
has profound implications on their
transport  \citep{filippov_drag} and thermal
\citep{filippov_thermal} properties.

Fractal aggregates arise from the agglomeration of smaller units,
herein taken to be spherical and hereafter referred to as monomers, that
do not coalesce, but rather retain their identity in the resulting aggregate.
In a quiescent fluid the main agglomeration driving mechanism
is diffusion: accordingly we model individual monomers
as interacting Brownian particles whose dynamics is described by a Langevin
equation \citep{risken_book}. Langevin simulations have been employed in aerosol
science to investigate aggregate agglomeration \citep{mountain}, aggregate
collisional properties \citep{pratsinis_kernel}, the limits of validity of Smoluchowski
equation \citep{pratsinis_ld}, and aggregate films
\citep{friedlander_deposition}.

In this study we investigate nanoparticle agglomeration and the dynamics
of the resulting aggregates by considering their diffusive growth by relying solely on the
Langevin equations of motion for a set of interacting monomers in three
dimensions. The monomer-monomer interaction potential is taken to be
either a model potential, composed of short range repulsive and attractive parts that
generate a deep minimum, or a Lennard-Jones intermolecular potential integrated
over the particle volumes. Unlike the works mentioned above,
no assumptions are made about the structure or mobility of the
aggregates generated by the dynamics of the individual monomers.
Our approach is reminiscent of other applications of Langevin simulations in
dilute colloidal suspensions to study diffusion-induced agglomeration and the
structures it gives rise to  \citep{videcoq}.

 An important
limitation of this approach is the neglect of the decrease of the
 surface area of monomers  in an aggregate accessible to surrounding fluid
molecules. Hence, for a $k$-monomer cluster ($k$-mer) generated  in this study, the
accessible area equals  $k$ times the accessible area of an
isolated monomer (sphere). This equality never holds for a real
cluster, hence we  call ideal the clusters generated in this study. 
The decrease of the accessible area for a monomer in a cluster
is  due to the shielding mainly by its first neighbours.
 In fact, the fluid-molecule concentration boundary layers of neighbouring
monomers overlap, thus hindering the access of fluid molecules
to the monomer surface.
%Shielding plays a role even in open-structured aggregates, like linear chains.   
%  The decrease of the accessible area, or equivalently of
% the area exposed to fluid molecules, arises either from shielding of
% inner monomers in the aggregate by neighboring monomers or from
% limited access of fluid molecules to the monomers, as would occur for
% example in a linear chain of monomers where the fluid-molecule concentration boundary
% of neighboring monomers layers overlap.
The decrease of the accessible area
modifies the friction coefficient a monomer feels, decreasing the friction
coefficient of the aggregates, and it modifies the effect of the thermal noise.
This is an inherent limitation of these simulations since Langevin
equation by itself does not provide any means to estimate the
accessible area of a cluster  (see for instance
Ref. \cite{filippov_thermal} for the calculation of the shielding
factor of a generic cluster).   Nevertheless, we  quantify the effect of
this approximations on the dynamics of the aggregate, in particular
the diffusive properties, while we argue
that it has a limited effect on the static structure of the generated
aggregates.
% In the following, we will call the cluster generated in our numerical
% experiments ideal clusters to mean 

Since we solve numerically the Langevin equations for interacting monomers
rather than the aggregate equations of motion,
the simulations output does not contain direct information
on aggregate formation. We infer this datum, together with the
detailed structure of each aggregate, the record of the collisions
it underwent and its eventual restructuring,
by using techniques borrowed from graph theory \citep{book_algorithms}.

The paper is organized as follows: Section \ref{model} provides the
theoretical framework for the Langevin simulations of nanoparticle
agglomeration. Emphasis is placed on the description and justification of the
monomer-monomer interaction potentials used in the numerical experiments.
Section \ref{simulations} offers an overview of the numerical method, and
it introduces the quantities monitored to investigate the statics and
dynamics of the generated aggregates. Section \ref{sec:Statics}
presents the static properties of the aggregates, namely their
fractal dimension, their morphology (specifically, the monomer coordination number),
and some indications of cluster restructuring. Section \ref{sec:Dynamics} summarizes
the cluster dynamical properties, in particular, their translational
and rotational diffusion coefficients, their number evolution in time, and
the determination of a limited number of
kernel elements. The final remarks in Section \ref{conclusions}
summarize the main results and conclude the paper.

\section{Model description}
\label{model}

\subsection{Monomer Langevin equation of motion}
\label{sec:lang-equat-mesosc}

%\label{sec:comp-technique}
% The dynamics of Diesel exhaust nanoparticles is assumed to be
% described by a generalized Langevin equation, where particles are
% allowed to interact via a short-ranged potential consisting of a
% highly-repulsive part

We investigate the non-equilibrium dynamics of interacting clusters in
the continuum regime  by
solving the Langevin  equations of motion of a system of $N$ interacting monomers
in three dimensions. We use the word cluster with the same meaning as the term
aggregate in Ref. \cite{konstandopoulos}, i.e. a set of physically bound
spherules (monomers). The $i$-th monomer obeys the Langevin equation
\begin{equation}
\label{eq:Langevin}
m_1\ddot{\bf r}_i={\bf F}_i-\beta_1 m_1\dot{\bf r}_i+
{\bf W}_i(t) \Comma
\end{equation}
where $m_1$ in the monomer mass, ${\bf r}_i$ its position in
three-dimensional space, ${ \bf F}_i$ an external force acting on
the monomer, $\beta_1$ the (per-unit-mass) friction coefficient,
%is the inverse of the monomer relaxation time $\tau_1$
and ${\bf W}$ is a noise term that
%${\bf W}_i=(W_{i}^x,W_{i}^y,W_{i}^z)$
models the effect of collisions between the monomer
and the molecules of the surrounding quiescent fluid. In the
following, we will always use boldface for vectorial or matricial quantities unless
otherwise stated.
%along the $x$, $y$ and $z$ direction.
We consider that each monomer feels a Stokes drag: the friction
coefficient becomes $\beta_{1} = 3 \pi \mu_f \sigma$ where $\mu_f$
is the fluid viscosity, and $\sigma$ the monomer diameter. Equivalently,
the friction coefficient is the inverse monomer  relaxation
time $\tau_1 = \rho_p \sigma^2/(18 \mu_f)$ where
$\rho_p$ is the monomer material density. The use of Stokes drag implies that
the monomer surface area accessible to fluid molecules equals the
accessible   area of an isolated monomer irrespective of the
state of aggregation of the monomer. We also neglect hydrodynamic interactions
that are important for the dynamics of dense colloidal suspensions,
where the scalar friction coefficient is replaced by a microscopic
friction tensor dependent on the positions of all other monomers.
The noise is assumed to be Gaussian and white (delta-correlated in
time)
\begin{equation}
\label{eq:gaussian_noise}
\langle  W^{k}_i(t) \rangle=0 \quad \quad {\rm{and}} \quad \quad
\langle W^{k}_i(t) W_{l}^{j}(t')\rangle = \Gamma \delta_{ij}
\delta_{k l}\delta(t-t') \Comma
\end{equation}
where subscripts denote monomer number, $i, j = 1, \ldots N$, and
superscripts denote cartesian coordinates $k,l = x,y,z$. The noise strength
$\Gamma$ is determined from the Fluctuation Dissipation Theorem
applied to each individual monomer to be $\Gamma=2\beta_1 m_1k_BT$
with $k_B$ Boltzmann's constant and $T$ the system temperature.
In the case of nanoparticles (mesoscopic system) the damping and noise
terms in Eq.~(\ref{eq:Langevin}) express
the effects of the collisions of small fluid molecules with
much larger solid nanoparticles.
Irrespective of the choice of the friction coefficient of each monomer,
as determined from the approximation of the monomer
accessible surface area, the assumption that the
fluctuation dissipation theorem applied to each monomer gives
the noise strength of that monomer is one of the
main assumptions of Langevin simulations of Brownian agglomeration
via Eq.~(\ref{eq:Langevin}). Alternatively, the Fluctuation Dissipation
Theorem could be applied on the aggregate and not on its
monomer constituents, rendering however difficult the identification
of the noise felt by each monomer.
In future sections we will discuss
in depth the implicit assumptions of using a Stokes drag and a
corresponding noise strength on all monomers at all times, i.e.,
irrespective of the formation of aggregates.

The force acting on a monomer is considered to be conservative
\begin{equation}
\label{eq:potential_pairwise}
{\bf F}_i=- \mbox{\boldmath$\nabla$}_{{\bf r}_i} U_i \Comma
\end{equation}
arising from its interaction with all the other monomers in the system.
We consider that the total intermonomer interaction potential $U_i$ monomer
$i$ feels derives from a pair-wise additive radial potential $u_{ij}(r_{ij})$
between monomer $i$ and $j$
\begin{equation}
\label{eq:pairwise_pot_explicit}
U_i=\f{1}2{}\s_{j\neq i}u_{ij}(r_{ij}) \Comma
\end{equation}
where the radial distance between the two monomers
is $r_{ij}=|{\bf r}_i-{\bf r}_j|$.

Equation~(\ref{eq:Langevin}) may be cast in dimensionless form
that is more convenient for its numerical solution. By fixing
length, $l^*$, time, $\tau^*$, and  mass, $m^*$, characteristic
scales we introduce the dimensionless distance, time and mass by
\begin{equation}
\label{eq:dimensionless-scales}
r\equiv l^*\tilde r,\;\;\;\; t\equiv \tau^*\tilde t, \;\;\;\;
m_1\equiv m^*\tilde m_1 \Teleia
\end{equation}
% The quantities in Eq. \ref{eq:dimensionless-scales} also fix aa
% characteristic temperature in the system
% $k_BT^*=m^*(l^*)^2/(\tau^*)^2$ which can be used to scale the system
% temperature $T=\tilde TT^*$.
% This way
Equations~(\ref{eq:Langevin}) and~(\ref{eq:gaussian_noise}) can be re-written
component-wise as
\begin{equation}
\label{eq:Langevin-dimensionless-componentwise}
\f{d^2\tilde r_i^l}{d\tilde t^2} =
-\f{1}{\tilde m_1}\f{\p\tilde U}{\p \tilde r_i^l}
-\tilde\beta_1\f{d\tilde r_i^l}{d\tilde t}
+\f{1}{\tilde m_1}\tilde W_i^l(\tilde t) \ , \quad l=1,2,3 \ , \quad i = 1, \ldots N
\end{equation}
\begin{equation}
\langle \tilde W_i^l(\tilde t)\rangle=0 \quad \quad {\rm{and}} \quad \quad
\langle \tilde W_i^l(\tilde t) \tilde W_j^k(\tilde t')
\rangle=\delta_{ij}\delta_{kl}\delta(\tilde t-\tilde t')2\tilde\beta\tilde m_1\tilde T
\end{equation}
where we introduced
\begin{equation}
\label{eq:derived-dimensionless-quantities}
\tilde\beta_1\equiv \beta_1\tau^* \, , \
k_BT^*\equiv\f{m^*(l^*)^2}{(\tau^*)^2} \, , \
\tilde T=\f{T}{T^*} \, , \
\tilde U\equiv \f{U}{k_BT^*} \, , \
\tilde W_i^l\equiv\f{(\tau^*)^2}{m^*l^*}W_i^l \, .
\end{equation}
The characteristic system temperature $T^*$ is
fixed by the chosen characteristic time, length, and mass scales.
Equation~(\ref{eq:Langevin-dimensionless-componentwise}) shows that
there are three independent dimensionless variables in the system, $\tilde m_1$,
$\tilde \beta_1$ and $\tilde T$, that render dimensionless the $i$-monomer
Langevin equation of motion.
A natural, but by no means unique, choice of the units in
Eq.~(\ref{eq:dimensionless-scales}) is: $m^*=m_1$, $\tau^*=\tau_1$, and
$l^*=\sigma$, where $\sigma$ is a characteristic length scale of the intermonomer
interaction potential (to be taken to be the hard-core monomer diameter).  This
amounts to setting $\tilde\beta_1=\tilde m_1=1$ in
Eq.~(\ref{eq:Langevin-dimensionless-componentwise}).
The characteristic temperature $T^*$ then can be evaluated: for example
consider a soot monomer with $\sigma=20$nm,
material density $\rho_p=1.3$g/cm$^3$ and suspended in air at $T=300K$
with dynamic viscosity $\mu_f=1.85\cdot 10^{-4}$g/(cm$\cdot$sec). The Stokes
particle relaxation time is $\tau_1=\sigma^2\rho_p/(18\mu_f) \sim 1.6 \times 10^{-9}$sec
leading to
\begin{equation}
\label{eq:T_star}
T^*=\f{18^2\pi\mu_f^2\sigma}{6k_B\rho_p}\simeq 650{\rm K}
\end{equation}
Therefore, $T=300$ corresponds to a dimensionless temperature
$\tilde T\simeq 0.5$.



\subsection{Monomer-monomer interaction potential}
\label{explain_potential}

The monomer-monomer interaction potential is chosen to mimic the
interaction of hard-core monomers sticking upon collision. The
condition of perfect sticking of two monomers when their relative
distance falls below the monomer diameter was used in the
early studies by \cite{mountain} and \cite{meakin_cluster_models}.
In these works the authors did not use a
monomer-monomer interaction potential, but they examined the system
frequently during its time evolution to identify agglomeration events
(monomer-cluster or cluster-cluster) by calculating the relative distances of the
monomers in the system. After, e.g., two clusters had touched, the
relative distances of all the monomers in the resulting cluster were
``frozen'', and the cluster was allowed to diffuse with a diffusion
coefficient that had to be prescribed \textit{a priori}.

Herein, we investigate aggregate formation as it arises from the collisions
of monomers that interact via a predefined interaction potential. We used two
intermonomer interaction potentials to describe the interaction of two
macroscopic bodies: a model potential $u_{mm}^{\textrm{mod}}$
(a short-ranged potential well) a potential $u_{mm}$ that arises from the integration
of the intermolecular interaction potential over
the volume of the macroscopic bodies The latter was taken to be
an integrated intermolecular Lennard-Jones potential
to model the attractive Van der Waals interactions between two
non-charged soot nanoparticles.

The model potential $u_{mm}^{\textrm{mod}}$ exhibits a deep and narrow
attractive well to represent monomer sticking upon
collision. It tends smoothly to zero at
$r= r_{\rm cut}$ where $r_{\rm cut}$ is a cut-off distance
such that $r_{\rm cut}-\sigma\ll\sigma$. This avoids the introduction
of the so-called impulsive forces in the system \citep{md_book}.
It is attractive on a length scale much smaller than
the monomer diameter $\sigma$, whereas $u_{mm}$ is a more long-ranged, though
quickly decaying, interaction.
We chose an analytical expression for the
monomer-monomer model interaction potential  to be
used in the numerical experiments as required by the numerical solver.
The functional form is
\beq
\label{eq:potential_well}
\umod(r) = \begin{cases}
\frac{\pi}{2} \, \Big ( \frac{\usigma - \umin}{\rmin - \sigma} \Big )
\left ( r - \sigma \right ) + \usigma & \text{if} \quad 0<r \le \sigma \Comma \\
\usigma - \left ( \usigma - \umin \right ) \, \cos \Big [ \frac{\pi}{2} \,
\left ( \frac{\rmin -r}{\rmin - \sigma} \right ) \Big ] & \text{if} \quad \sigma<r \le r_{\rm min} \Comma \\
\umin \cos^2 \Big [ \frac{\pi}{2} \, \left ( \frac{r -\rmin}{\rcut - \rmin} \right ) \Big ]
& \text{if} \quad r_{\rm min}<r \le r_{\rm cut} \Comma \\
0 & \rm{elsewhere} \Comma \\
\end{cases}
\eeq
where $\usigma = \umod (\sigma)$, $\rmin$ is the location of the potential minimum,
and $\umin = \umod (\rmin)$. The model potential depends on five parameters
($\rmin$, $\sigma$, $r_{\textrm{cut}}$, $\usigma$, and $\umin$)
chosen as follows in our numerical simulations. The potential minimum
is located at $r_{\rm min}=1.05\sigma$, where it
evaluates to $u_{\rm min}=-100k_BT$, ensuring a monomer sticking
probability of unity. At monomer separation $\sigma$ the potential evaluates to
$u (\sigma) =60k_BT$ with a steep gradient. For monomer separations in
the range $(0,\sigma]$ the potential is extrapolated linearly
with the slope it has at $r=\sigma$.
Hence, for separations smaller than $\sigma$, the monomers feel a
strong constant repulsive force matching the one at $r=\sigma$.
The cut-off distance for the potential is $r_{\textrm{cut}}=1.1\sigma$

The interaction potential between two equal-sized
nanoparticles is obtained from the integration of the
intermolecular interactions over the nanoparticle volumes
as in Ref.~\cite{yannis_potential}. We assume pairwise additivity
of the intermolecular interactions, continuous medium, and  constant material
properties. Elastic flattening of the monomer is neglected.
As usual the intermolecular interaction potential is
taken to be the Lennard-Jones potential
\begin{equation}
\label{eq:LJ}
u_{LJ}(r)=4\epsilon\lsq\lro\f{\sigma_{LJ}}{r}\rro^{12}-\lro\f{\sigma_{LJ}}{r}\rro^6  \rsq,
\end{equation}
where $\epsilon$ is the depth of the attractive potential,
the maximum attractive energy between two molecules, and $\sigma_{LJ}$
the distance at which the potential evaluates to zero,
the distance of closest approach of two molecules which collide
with zero initial relative kinetic energy. The first
term expresses the repulsive part (Born repulsion),
the second the attractive (van der Waals attraction).
Integration of Eq. \ref{eq:LJ} over two
equal-sized spheres yields, see, for example
Refs.~\cite{yannis_potential} and  \cite{ruckenstein},
\beq
\label{eq:Umm}
\begin{split}
 u_{mm} (r) = & -A\lcu \sigma_{LJ}^{6}\lsq\f{\sigma}{7560r}\lro \f{1}{(r-\sigma)^6}-
\f{1}{(r+\sigma)^6}\rro \right. \right. \\
&  +  \f{1}{37800r}\lro\f{2}{r^5}-\f{1}{(r-\sigma)^5}-\f{1}{(r+\sigma)^5}\rro
  \\
& \left. -\f{\sigma^2}{2520r}\lro\f{1}{2(r-\sigma)^7}
  +\f{1}{2(r+\sigma)^7}+\f{1}{r^7}\rro\rsq \\
& + \f{1}{6}
\lsq\log\lro\f{r^2-\sigma^2}{r^2}\rro
  \left. \left. +\f{\sigma^2}{2(r^2-\sigma^2)}+\f{\sigma^2}{2r^2}\rro\rsq
      \rcu \, ,
\end{split}
\eeq
where $A=4\pi\epsilon\sigma_{LJ}^6n^2$ is the Hamaker constant and
$n$ is the molecular number density in the solid. Note that in the
limit $r \gg \sigma$ the intermonomer interaction potential decays as
$u_{mm} \sim r^{-6}$, i.e. the long-range van der Waals interaction between
two macroscopic bodies decays as the Lennard-Jones intermolecular potential. The Hamaker constant
for a typical example of a soot nanoparticle, for example for a n-hexane soot nanoparticle,
was estimated to be $A=2.38\cdot 10^{-19}$J \citep{kasper_hamaker}
with a corresponding $\sigma_{LJ}=0.5949$nm \citep{substance_book}. We used these
values in our numerical simulations, along with $\sigma = 20$nm.

The repulsive part of $u_{mm}$ diverges for $r\to\sigma$ thus
potentially giving rise to numerical difficulties in the simulations.
We therefore modified the potential in Eq.~(\ref{eq:Umm}) at
distances $r_{\rm match}=1.015\sigma$  (hence less than the position
of the potential minimum at $r_{\min}\simeq 1.017\sigma$) by
extrapolating it linearly with the same slope it exhibits at $r=r_{\rm
sep}$ until $r=0.995\sigma$ (where it evaluates to about $60 k_B T^*$).
(A similar matching condition was used to determine the coefficient
of the model potential linear term.) We set the monomer-monomer interaction
potential $u_{mm}$ to zero at smaller monomer-monomer separations. This modification of the repulsive part
is not expected to affect the dynamics of the system, since it
involves only monomer-monomer separations are energetically unlikely.
Hence, intermonomer interaction potential we use is
\beq
\label{eq:van_der_waals}
u(r)_{mm}^{\textrm{sim}} =
\begin{cases}
0 & \text{if} \quad 0 \le r \le 0.995 \sigma \Comma \\
\f{\p u_{mm}}{\p r}\bigg{|}_{r_{\rm match}} (r_{\rm match}-r)+u_{mm}(r_{\rm match})
& \text{if} \quad 0.995 \sigma <r \le r_{\rm match} \Comma \\
u_{mm}(r) & \text{if}  \quad r_{\rm match} < r \Teleia \\
\end{cases}
\end{equation}
We truncate monomer-monomer interaction potential at distances
$r \geq 7 \sigma$ where it is negligible with respect to the thermal energy,
$u(7\sigma)/(k_BT^*)\sim 10^{-4}$.

In both cases, for distances smaller than $\sigma$, the monomers
feel a very strong repulsive force much larger than the other energy
scale in the system, namely the thermal energy $k_BT$, where  $T$ is
the system temperature. Although neither potential  diverges at
separations $r<\sigma$ and one should call $\sigma$ the soft-core
monomer diameter, monomer separations  below $\sigma$ are energetically
unfavorable and extremely unlikely occur in the system dynamics. This
justifies the identification of $\sigma$ with the monomer hard-core diameter.

In the following, we will be using only dimensionless quantities, but
we will drop the tilde to simplify the notation. Furthermore,
unless specified otherwise, the results presented will be those
obtained with $u_{mm}^{\textrm{sim}}$, Eq.~(\ref{eq:van_der_waals}).
Figure~\ref{plot_potential} presents and compares the dimensionless
radial potential used in the
numerical simulations, namely the model potential, Eq.~(\ref{eq:potential_well}), and the
integrated Lennard-Jones potential, Eq.~(\ref{eq:van_der_waals}).

\section{Numerical method}
\label{simulations}

\subsection{Numerical solver}

The numerical simulations were performed with the software
package ESPResSo \citep{espresso}, a versatile package for generic Molecular Dynamics
simulations in condensed-matter physics. We used the Molecular Dynamics
program with a Langevin thermostat to simulated the Langevin dynamics of
Brownian particles in a quiescent fluid. The Langevin thermostat
was construed as formal method to perform Molecular Dynamics simulations
in a constant temperature canonical (NVT) ensemble. The Langevin
thermostat introduces a fictitious viscous force to model the coupling
the system to a thermal bath according to Eq.~(\ref{eq:Langevin}):
the fictitious, random force slows down hot molecules and speeds up
cold molecules so that the system temperature fluctuates
about the temperature of the thermal bath. The molecular equations of motion
with a Langevin thermostat are formally identical to the Langevin equations
of interacting Brownian particles: the physical scales and their interpretation
differ. The molecular dynamics simulations provide a
microscopic description of the system, whereby all the degrees of freedom are
explicitly resolved (except for the bathe degrees of freedom), whereas the Langevin simulations provide a
mesoscopic description
of the system, whereby the degrees of freedom associated to the fluid molecules
are model as a stochastic term in the equations of motion.
Moreover, the Langevin thermostat does not conserve momentum in
molecular-dynamics simulations of Brownian particles in a bath of solvent molecules
i.e. in simulations where the equations of motion of both the solvent and
the Brownian particles are explicitly solved and both are coupled to
a bath (to ensure a constant temperature simulation). Artificial
non conservation of global momentum leads to an artificial
shielding of hydrodynamic interactions, and hence coupling of the
motion of the Brownian particles: however, this limitation of the Langevin
thermostat is not relevant for the Langevin simulation of interacting
Brownian particles. this is not a problem for the
mesoscopic description we adopt.

As a check of the numerical solver we simulated the motion of a
single Brownian particle in a quiescent fluid. The simulations for
the mean-squared displacement $\langle r_{1}^2(t) \rangle$,
the mean-squared velocity $\langle v_1^2 \rangle$, and
the velocity autocorrelation function $\langle {\bf v}_1(t){\bf v}_1(0) \rangle$
agreed with analytical expressions~\citep{risken_book}.

The initial state is created by placing randomly in a cubic box of size $L$
(measured, like all distances, in units of monomer diameter
$\sigma$),  $N_{\infty}(0)V_{\rm box}=5000$ monomers, where $V_{\rm
  box}=L^3$ is the box volume and $N_{\infty}(0)$ is the initial
monomer concentration.  The box size $L$ is chosen to obtain a given initial monomer
density according to $L=\lro{5000}/{N_{\infty}}\rro ^{1/3}$.
% \begin{equation}
% \label{eq:box_size}
% L=\lro\f{5000}{N_{\infty}}\rro ^{1/3},
% \end{equation}
We chose $N_{\infty}(0)=0.01$ that leads to $L\simeq 80$.
The initial random placement of monomers in the box could give rise
to numerical problems if two or more monomers happen to be placed in
almost overlapping positions, thus experiencing immediately a very
strong repulsive force. We use a well-known technique of
MD simulations to avoid these numerical instabilities by ramping up the repulsive force
monomers feel when $r<1$ to its constant final value  during a time
$\tau_{\rm ramp}=1$. The simulation time step was chosen to be $\delta t_{\textrm{sim}} = 1.25 \times 10^{-3}$.
After the initialization the system is evolved until
a final dimensionless time $3000$, when the cluster concentration is almost two orders of
magnitude smaller than the  initial.
% We notice that the minimum of the
% potential and the value it assumes when the monomer separation is
% equal to the diameter $\sigma$ are both much larger than $k_bT$ which
% is the kinetic energy scale for the monomers.
The results we show in this study were obtained 
using the 
output of $10$ simulations with different initial conditions, i.e. random placement of
the monomers within the computational box, 
  to ensure that enough data for the statistical analysis
are available.
%  It is worth pointing out that the results in Ref. \cite{mountain} were
% based on a system containing one tenth of the monomers in this study.

\subsection{Cluster identification}
\label{sec:cluster-calculation}

The identification of clusters is one of the most time-consuming
tasks of post-processing the simulation results. The \esp simulations
return individual monomer positions and velocity. Unlike the previously
mentioned works \citep{meakin_cluster_models, mountain},
agglomeration events are not identified during the time evolution of
the system, but cluster formation is determined \textit{a posteriori}. Sampling
of the  simulation output was performed at time intervals $\delta
t_{\rm sam}=2$ during the whole simulation.

We resort to an approach based on graph theory.
A perfectly rigid cluster is a set of connected monomers such that all monomer-monomer distances
are constant in time. Since both intermonomer potentials Eq.~(\ref{eq:potential_well})
and Eq.~(\ref{eq:van_der_waals}) are not infinitely narrow,
monomers perform tiny oscillations about their equilibrium position.
Despite of that, the bottom of the potential well is so deep that when
two monomers collide they remain bound (no agglomerate breakup was noticed).
As a consequence, one may consider two monomers bound when their
relative distance is less than a threshold distance $d_{\rm thr}$.
The choice of $d_{\rm thr}$ depends on the position of the minimum of
monomer interaction potential $r_{\rm min}$. For the simulations with
the integrated Lennard-Jones potential we found $d_{\rm thr}=1.04$ to provide accurate
results: for the model potential, $d_{\rm thr}=1.06$ is more
appropriate. Small variations of $d_{\rm thr}$ about
these values do not affect the identification of clusters.
Our definition of a cluster is reminiscent of the liquid cluster definition
proposed by Stillinger and used in gas-liquid
nucleation studies (see, for example, Ref.~\cite{clusterReguera}).
The objections and criticism raised in the Stillinger cluster definition,
e.g., that molecules close to the cluster but not belong to it, are not relevant
for our simulations as the monomer density is low and the potential well
much deeper that the intermolecular well.

Computationally, the first step is the calculation of the distance
matrix ${\bf D}$ among all the monomers in the system.
For a periodic system, the distance between the $i$-th and $j$-th
monomer is the distance between the $i$-th monomer and the nearest
image of the $j$-th monomer \citep{md_book}. For instance, the
(ordered) distance between two monomers along the $l$-th  co-ordinate is
\begin{equation}
\label{eq:distance_calc}
D_{ij}^{l}= r_i^{l}-r_j^{l}-L\cdot{\rm nint}\lro\f{r_i^{l}-r_j^{l}}{L} \rro
\quad , \quad l=1,2,3 \Comma
\end{equation}
where``nint" is the nearest integer function. The periodicity of the
box imposes a cutoff on the maximum distance between two monomers along each axis,
namely $L/2$. Furthermore, $D_{ij}^{l}=-D_{ji}^{l}$, so the distance along $x$ in
Eq. \ref{eq:distance_calc} is actually the $x$ component of the vector
pointing from the $j$-th to the $i$-th monomer. The three-dimensional distance becomes
\begin{equation}
\label{eq:distance_3D}
D_{ij} = \Big \{ \sum_{l=1}^3 \big [ D_{ij}^{l} \big ]^2 \Big \}^{1/2} \Comma
%D_{ij}=\sqrt{{D_{ij}^{(x)}}^2+{D_{ij}^{(y)}}^2+{D_{ij}^{(z)}}^2}.
\end{equation}
always a  non-negative quantity. The distance matrix $\bf{D}$
is subsequently used to calculate the adjacency
matrix $\bf{A}$, defined as
\begin{equation}
\label{eq:adjacency_matrix}
A_{ij} =
\begin{cases}
1 & \quad \text{if } \quad D_{ij} \le d_{\rm thr}\\
0 & \quad \text{otherwise }
\end{cases}
\end{equation}

The adjacency matrix is usually introduced in graph theory \citep{book_algorithms}
as a convenient way to represent uniquely a graph. Monomers in a cluster can be formally regarded as
graph vertices (nodes) and bonds due to the interaction potential become
graph edges (links). The problem of cluster identification, given the distance matrix
${\bf D}$, is then re-formulated as the identification of
the connected components in a non-directed graph expressed by the
(symmetric) adjacency matrix ${\bf A}$. They can be determined using a
standard breadth-first search (BFS) algorithm
\citep{book_algorithms}. Due to its speed and scalability, we resort
to the BFS algorithm  implementation in the igraph library \citep{igraph}.

\subsection{Cluster radius of gyration}

The radius of gyration is a fundamental property of the aggregate; indeed
the dependence of the number of monomers in a cluster on its radius of
gyration can be used to determine its fractal dimension. The radius of gyration is
defined as the root-mean-square distance mass displacement with respect to
the aggregate center of mass
\begin{equation}
\label{eq:r_gyr_definition}
R_g^2 = \frac{1}{k} \, \s_i^k{ ({\bf r}_i -{\bf R}_{CM})^2} + a_1^2,
\end{equation}
where $k$ is the number of monomers in the cluster and $a_1$
the radius of gyration of a monomer. The monomer radius of gyration,
calculated as the radius of gyration of a single three-dimensional sphere,
is included in the definition of the cluster radius of gyration
so that, if an average fractal exponent is used, the large $k$ behaviour
could persist even for small clusters. It evaluates to
$a_1=\sqrt{3/5}(\sigma/2)$, as in Ref.~\cite{sorensen_prefactor}.
We stress that the cluster radius of gyration is a static property
since it does not depend on diffusive properties of the
cluster, but only on its mass distribution.

% However, we do not calculate the radius of gyration according to
% Eq.~(\ref{eq:r_gyr_definition}).

The calculation of $R_g$ according to Eq.~(\ref{eq:r_gyr_definition})
must take into consideration the periodicity of the box.
We adopt a reference frame centered on the position of the first
monomer in the cluster, i.e. we work in a reference frame where
${\bf  r}_1=(0,0,0)$. Since all monomer-monomer distances along each axis are
known from Eq.(\ref{eq:distance_calc}), the position of all the other monomers in a cluster
with respect to the first monomer may be easily calculated;
for instance, the $j$-th monomer position along the $l$-th co-ordinate is
\begin{equation}
\label{eq:reconstruct_cluster}
 r_j^{l} =r_1^{l} +D_{1j}^{l} = D_{1j}^{l} \quad ,
 \quad l=1,2,3 \Teleia
\end{equation}
Once all the new monomer coordinates have been obtained
from Eq~(\ref{eq:reconstruct_cluster}), the coordinates of the
aggregate center of mass ${\bf R}_{CM}=\s_i{\bf r}_i/k$  
may be calculated to evaluate
straightforwardly Eq.~(\ref{eq:r_gyr_definition}).
We note that this procedure is independent on the specific
monomer whose coordinates are taken as the origin of the new reference frame
the radius of gyration  is independent of the cluster spatial position) nor
does it depend on the sign convention chosen for $D_{1j}^{l}$.

\subsection{Collision kernel}

The package \esp allows addressing each monomer individually
at all times during the simulation, i.e. a permanent label, or color, may
be associated with each monomer. Each cluster becomes an unordered
collection of monomers where no monomer label is repeated. This amount
to the mathematical definition of a set.
Viewing clusters as sets of monomers provides a computational method
to investigate the process of cluster collisions even for sampling times considerably
longer than the simulation time step if aggregates do not break up. A monomer label
allows unequivocal identification of the clyster it belong to.
Since the two intermonomer potentials used are much deeper than $k_BT$ the
a sticking probability upon collision is unity.
To be specific, consider two clusters at time $t$.
If during the (sampling) time interval $(t, t+\delta t]$ they collide, they will be
part of the same newly-formed cluster at time $t+\delta t$.
An equivalent description of the collision is that the two monomer sets that
identify the pre-collision clusters at time $t$ will become proper subsets of
the monomer set that identifies the cluster formed by their collision, and detectable at time
$t+\delta t$. The number of collisions that took occurred during $(t, t+\delta t]$
and the clusters involved in the collisions may be recorded simply by comparing different sets.
The collision kernel between $i$ and $j$-mers during the sampling time interval
$\delta t$ is, then, estimated from
 \begin{equation}
\label{eq:kernel_calc}
\f{N_{ij}(t)}{\delta t}=(2-\delta_{ij}) K_{ij} \f{N_i(t)}{V_{\rm
 box}}\f{N_j(t)-\delta_{ij}}{V_{\rm box}}\cong (2-\delta_{ij}) K_{ij}n_i(t)n_j(t),
\end{equation}
where $N_{ij}(t)$ is the number of collisions between $i$- and $j$-mers
that took place during the interval $(t, t+\delta t]$, $V_{\rm box}$ is the
volume of the computational box,  $N_i(t)$ is the number of
$i$-mers at time $t$, $n_i(t)=N_i(t)/V_{\rm box}$ the cluster concentration,
and $\delta_{ij}$ the Kroenecker delta. The collision kernel $K_{ij}$ is estimated
from Eq.~(\ref{eq:kernel_calc}) by treating it as the fitting parameter
that minimizes (in a least-square sense) the distance between the
number of detected collisions
$N_{ij}(t)/\delta_t$ and $(2-\delta_{ij}) K_{ij}n_i(t)n_j(t)$.

\section{Static cluster properties}
\label{sec:Statics}

\subsection{Cluster fractal dimension}
\label{calculate_df}

Examples of 50-monomer clusters at the end of the simulation
are shown in Fig.~\ref{snapshot_end_simulation}. Note their relatively
compact structure.

The fractal dimension of the generated clusters is determined by
fitting the power-law dependence of the clusters size $k$ on its radius of
gyration $R_g$ according to
\begin{equation}
\label{eq:determination_df}
\overline{R}_g = a k^{1/d_f},
\end{equation}
where $d_f$ is the average, time-independent fractal dimension of the aggregates
and $a$ the prefactor. The brackets denote an ensemble average, a time-independent average
over generated cluster configurations. In particular, the average radius of
gyration of a k-monomer aggregate is taken over all k-monomer aggregates
that had been recorded at least $200$ times. This requirement aims at
eliminating outliers when calculating the average radius of gyration. It
is not particularly stringent since for each simulation $1500$ system configurations (snapshots)
were stored, and $10$ different simulations, corresponding to different initial
conditions, were performed. Our results are not sensitive to reasonable choices
of the threshold of minimum number of cluster occurrences: raising the occurrence threshold
$400$ does not change the calculated average, time-independent fractal dimension.
We point out that, unlike e.g. \cite{pratsinis_ld}, where Langevin
simulations are used to investigate the limits of validity of
the continuum kernel, here the aggregate fractal dimension is a simulation
output and not a tunable input parameter.

Figure~\ref{two_df} presents the calculated average $\overline{R}_g $ as a function of the
aggregate size $k$. The double logarithmic plot shows that the power-law scaling breaks down at
cluster sizes $k<5$. Hence, we set $k=5$ as the minimum cluster size
in the fits presented in Fig.~\ref{two_df}. All larger than 5 cluster configurations
were fitted to a single line leading to an average fractal exponent $d_f=1.62$.
Inspection of the figure shows that clusters with $k\le 15$ deviate
significantly from the single-line fit. Instead, simulation results
are fitted better by considering two cluster population (and hence
two different slopes): we identify small ( $k\le 15$)  and
large ($k>15$) clusters. Refitting the data gives
a fractal dimension $d_f^{mc}=2.25$ for small clusters
and $d_f^{cc}=1.56$ for large clusters.

The existence of two cluster populations obeying different scaling laws
is attributed to different agglomeration mechanisms. Small clusters
are generated mainly by monomer-cluster agglomeration (also called
diffusion-limited agglomeration [DLA]); they are more compact and spherical than
large clusters that are mainly generated by cluster-cluster agglomeration.
Therefore, it is expected, and numerically confirmed, that
$d_f^{mc}>d_f^{cc}$: in fact, the calculated fractal dimension for
small and large clusters is indeed  close enough to the
values $2.4-2.5$ and $1.7-1.8$ usually reported in literature
for DLA and cluster-cluster agglomeration,
respectively~\citep{mountain, smirnov_review}.

The existence of two different cluster populations that arise from
the predominance of different agglomeration mechanism suggests
the calculation of a time dependent average fractal dimension
$d_f (t)$. We calculate it by considering the instantaneous radius of gyration
of $k$-monomer cluster for each of the ten initial configuration:
only cluster with $k>5$ without imposing
an occurrence threshold. The fractal dimension is obtained
through the power-law scaling of Eq.~(\ref{eq:determination_df}) for each recorded
configuration at time $t$. As the number of cluster is considerably smaller we
do not calculate first the instantaneous average radius of gyration and
then fit the data (as done for the calculation presented in Fig.~\ref{two_df}):
instead we fit all the date simultaneously (no averaging). The time-evolution of the time-dependent
fractal dimension is shown in Fig.~\ref{evolution_df} (top).
The evolution of the fractal dimension may be linked to the kinetics of
the small and large cluster populations. At early times, the system is almost entirely made
up of small clusters ($5 < k \le 15$), hence $d_f (t) \simeq d_f^{mc}$, whereas at late times
the contribution of large clusters to the overall cluster population
is dominant, hence $d_f (t) \to d_f^{cc}$.

These findings are in the spirit of the study by Kostoglou and
Konstandopoulos~\cite{konstandopoulos} who used a distribution of fractal
dimensions to characterize the morphology of the aggregates: they showed that
the mean fractal dimension relaxes to an asymptotic value that had been
fixed \textit{a priori} by the agglomeration mechanism.
We remark again that in this study no assumption is made
about the asymptotic fractal dimension: in fact, we show explicitly the importance of
the agglomeration kinetics on the evolution of $d_f (t)$.
The properties of the diffusion coefficient are not expected to have a major effect on
the calculated fractal dimension(s). For instance, Refs.~\cite{smirnov_review, mountain}
show that the aggregate fractal dimension depends mainly on the aggregation mechanism only.
Indeed, we remark that the calculated average fractal dimensions do not
depend sensitively on the choice of the intermonomer interaction potential (for the two potential
previously described). This seems in contrast with
the findings in Ref.~\cite{videcoq}, where the fractal dimension is found
to be a function of the interaction range. Here however, both potentials,
and specifically even Eq.~(\ref{eq:van_der_waals}),
are quickly decaying functions of $r$: thus, the interaction range does
not enhance significantly the cluster re-organization, which leads to
the increase of the calculated fractal dimension in the systems studied
in \cite{videcoq}.

\subsection{Cluster Morphology and Coordination number}
\label{sec:coordination-number}

An important quantity to characterize cluster morphology is the cluster coordination
number defined as the mean number of first neighbours of a monomer in the cluster.
In the following, we consider the mean coordination number,
namely the time-dependent, ensemble-averaged coordination number, the
coordination number of each cluster averaged over all clusters at a given time.
This quantity provides valuable information on the openness of the
aggregates, their compactness, and the presence of cavities in their structure, their
porosity. The mean coordination number provides additional
information on the aggregate structure in addition to the fractal dimension.
The coordination number is calculated in the post-processing of the
simulations again resorting  the igraph library
\citep{igraph}.

The time-dependent mean coordination number
is plotted in Fig.~\ref{evolution_coord_number}. The
mean coordination number goes well above five, suggesting that
the aggregates are relatively compact. This relatively high value provides
new information about the cluster morphology: the fractal dimension, as determined
from the scaling law in Eq.~(\ref{eq:r_gyr_definition}), is not sufficient
to characterize fully the aggregate structure.
Indeed, the large clusters generated at late times have a
low average fractal dimension $d_f^{cc}<2$ due to their
tubular, elongated shape (examples are the second
and fourth cluster, clockwise in Fig. \ref{snapshot_end_simulation}): their low fractal
dimension does not arise from the presence of holes or cavities in their structure. These
clusters are compact and not porous. This conclusion (valid in a statistical sense) cannot
be reached on the basis of the fractal dimension alone: hence the use of the fractal dimension
and the coordination number are required for a better characterization of aggregate
morphology. The generation of clusters with a high coordination number (as most
likely those in Ref.~\cite{videcoq}) is an intrinsic property of
Langevin simulations  due to their neglect of interparticle hydrodynamic
interactions, as demonstrated in Ref.~\cite{prl_hydro}.

\subsection{Cluster Restructuring}
\label{sec:clust-restructuring}

Cluster restructuring has been investigated in a 3d 
lattice model in Ref.~\cite{meakin_reorganization_3D}.
The novelty of our approach is that cluster
restructuring does not have to be implemented as an additional feature
of cluster dynamics. The monomers interact via a deep interaction \emph{radial}
potential. Its effect is to lock the distances between neighbouring monomers,
but the bonds between two monomers, once formed, may move on the surface of
the monomers. The mechanism responsible for this motion is collisions
of the molecules of the surrounding fluid with the
monomers, as modelled by the noise term in Eq.~(\ref{eq:Langevin}).
% As in Ref. \cite{meakin_reorganization_3D}, we observe cluster
% restructuring to take place on a smaller time-scale than the

We use two statistical quantities as indicators of the modification of
the structure of a cluster: the radius of gyration and the mean
monomer-monomer distance. For a $k$-mer the mean monomer-monomer distance
$D_{mm}$ is the average monomer-monomer distance over all the monomers,
\beq
D_{mm} = \frac{1}{k^2} \, \sum_{i,j}^k D_{ij} \Comma
\eeq
where the monomer-monomer distance was defined in Eq.~(\ref{eq:distance_3D}).
Specifically, we select a monomer at the beginning of the simulation to
identify a specific cluster (which will of course undergo collisions
and increase its size as time progresses).
In Fig.~\ref{single_cluster_size_evolution} (left) we show the
evolution of the mean monomer-monomer distance $D_{mm}$ (the mean was
taken over all the monomer-monomer distances in the cluster) during a
time interval when the chosen cluster increases its size from $k=14$ to $k=43$
due to a collision at time  $t=568$. We notice that most of the time
before and after this collisions, the cluster $D_{mm}$ is practically constant,
a clear sign of its rigidity. When the collision between the two aggregates occurs,
$D_{mm}$ increases, but on a time scale of a few monomer relaxation
times $\tau_{1}$ it 
settles to a lower value, which remains relatively constant.
We consider this behaviour a strong indication of cluster reorganization
resulting from the collision of two clusters.
A perfectly analogue behavior is observed when the radius
of gyration of the cluster on the same time interval (see the right
plot in Fig.~\ref{single_cluster_size_evolution}) is considered.

\section{Dynamic cluster properties}
\label{sec:Dynamics}

The motion of a cluster in three-dimensional space is a coupled
translational and rotational motion. In the following section,
Sec.~\ref{sec:clust-diff-coeff}, we investigate the
translational motion of the cluster centre-of-mass, while
its rotational properties are summarized
in Sec.~\ref{sec:clust-diff-rotation}.

\subsection{Cluster translational diffusion coefficient and cluster thermalization}
\label{sec:clust-diff-coeff}

Equation~(\ref{eq:Langevin}) determines the dynamics of each monomer and indirectly
the cluster diffusive properties. The numerical determination of the diffusion coefficient
of a $k$-mer (i.e. a cluster consisting of $k$ monomers) in the
continuum regime and its scaling with $R_g$ and the number
of monomers $k$ provide valuable insights into the dynamics of the cluster.
The diffusion coefficient of a $k$-mer $D_k$ may
be expressed as a generalization of the
Stokes-Einstein diffusion coefficient of a
single, spherical, monomer. Specifically, the latter is \citep{risken_book}
\begin{equation}
\label{eq:monomer-diffusion-coefficient}
D_1=\f{k_BT}{m_1\beta_1}.
\end{equation}
% where the Cunningham slip factor~\cite{Hinds} $C_s(Kn)$ was included
% to account for non-continuum effects. It is a function
% of the \textit{gas} Knudsen number defined as
% $\textrm{Kn} =2\lambda_f/d_{agg}$, where $\lambda_f$ is the fluid mean-free path and
% $d_{agg}$ is the aggregate diameter (here taken to be $2R_g$). The Knudsen
% number is used to estimate whether the system is in
% the continuum regime ($\textrm{Kn} \ll 1$), in the free-molecular
% regime $\textrm{Kn} \gg 1$), or in the transition regime.
The Langevin equations~(\ref{eq:Langevin}) characterize the fluid only via
its viscosity, modelled by the damping and noise terms, providing no
information on the fluid mean-free path.
Equation~(\ref{eq:monomer-diffusion-coefficient}) can be extended for a
$k$-mer by introducing the cluster mass $M_k = km_{1}$
and the per-unit-mass friction coefficient for a $k$-mer 
$\beta_k$
\begin{equation}
\label{eq:cluster-diffusion-coefficient}
 D_k=\f{k_BT}{km_1\beta_k} \Teleia
\end{equation}

The cluster diffusion coefficient is determined from numerical simulations by
considering the late-time dependence of the variance of the cluster center-of-mass
position at a given time
displacement~\citep{risken_book}
\begin{equation}
\label{eq:diffusion_simu_cluster}
\langle \delta{ R}^2_{CM}(t) \rangle =
\langle \big [ {\bf R}_{CM} (t) - \langle {\bf R}_{CM} (t) \rangle \big ]^2 \rangle \,
\xrightarrow{t\to\infty} 6D_{k}t \Comma
\end{equation}
where %  the position of the cluster centre of mass is
% \beq 
% \label{centre-of-mass-diffusion}
% {\bf R}_{CM} (t) = \frac{1}{k} \, \sum_i^k {\bf r}_i \Comma
% \eeq
% and
${\bf r}_i$ is the position of the $i$-th monomer in the $k$-mer
and $\langle \rangle$ denotes the average
over many cluster trajectories of the cluster mean square displacement.
Accordingly, aggregate diffusive transport is investigated by placing a given
cluster at rest in  the simulation box with centre-of-mass
position ${\bf R}_{CM}=(0,0,0)$. The simulation box size was chosen to have $L=10,000$,
a size much larger than the box size used for the agglomeration simulations
because in the investigation of aggregate diffusion no displacement of the aggregate from its
initial position larger than $L/2$ is admissible due to the
periodicity of the system. Hence the necessity of choosing a large box to track
the aggregate center of mass for a (potentially) long time. 
For the diffusion simulations described above, one does not need to
resort to 
 Eq.~(\ref{eq:reconstruct_cluster}) but can use directly  the  monomer
 positions in each cluster 
as returned by \esp to evaluate ${\bf R}_{CM}$. Indeed this time no
cluster ever crosses the boundary of the box, therefore box periodicity does
not play any role. Furthermore no cluster detection
based on Eq.~(\ref{eq:adjacency_matrix}) is needed here, since no
cluster undergoes any collisions.
 %  hence the cluster centre-of-mass
% position $R_{CM}$ is no longer determined in an arbitrary coordinate
% frame chosen for calculating the cluster radius of gyration.

We followed the motion of the cluster
centre-of-mass along $800$ trajectories up to time $t=400$ for clusters with $k=4,10,50,98$.
Hence, the different cluster dynamics does not arise from different initial conditions,
as was the case in the calculation of the cluster radius of gyration, but from the
noise. Each trajectory started with identical cluster.

Fig. \ref{diffusion_plot} shows our simulation results for
$\langle \delta{R}^2_{CM}(t) \rangle$ for a $50$-monomer cluster,
chosen to be the cluster at the bottom right of Fig. \ref{snapshot_end_simulation}.
The time dependence of the ensemble-averaged mean-square displacement quickly becomes linear,
in agreement with Eq.~(\ref{eq:diffusion_simu_cluster}).
The inset of Fig. \ref{diffusion_plot} magnifies (on a double
logarithmic scale) the early-time behavior of
$\langle \delta{R}^2_{CM}(t) \rangle$. For $t\le 1$ cluster dynamics
is ballistic, and the mean-square displacement
exhibits a power-law dependence on time,
$\langle \delta{R}^2_{CM}(t) \rangle\propto t^\gamma$ with $\gamma \simeq 3$.
This is the expected behaviour of a single Brownian
monomer with an initial zero velocity \citep{risken_book}.
Similar behaviour was found for the diffusive properties of a cluster
with $k=4$. %, thus ruling out the possibility that only large clusters are in the continuum regime.

The results of numerical simulations for the four clusters
are summarized in Table~\ref{table_diffusion},
where we also report the cluster radius of gyration $R_g$. The numerically
evaluated diffusion coefficients behave as $D_k\propto 1/k$
for $k$ spanning almost two orders of magnitude. % Consequently, the ratio of
% the Cunningham slip factor to the cluster friction coefficient (per cluster mass)
% is independent of the radius of gyration. If we take $C_s(\textrm{Kn} = 0) =1$, i.e.,
% if we assume that the system is in the continuum regime,
The estimated value of $\beta_k$
equals the monomer friction coefficient  $\beta_1$ within a few percents.
The small deviations is attributable to the relatively limited number of
simulated stochastic trajectories. Therefore, our Langevin simulations suggest
that the cluster diffusion coefficient is inversely proportional to the cluster
mass and to the monomer friction coefficient. This behaviour is a direct
consequence of neglecting shielding of inner monomer by outer monomers, or
equivalently by assuming that the cluster accessible area in the sum of the
monomer accessible areas. This is an inherent assumptions of all
simulations of monomer agglomeration that start from the unshielded 
monomer Langevin equations Eqs.~(\ref{eq:Langevin}).
%  Say something more
% and explain why this is so. They do not differentiate between exposed or
% accessible surface area that is screened (the Langevin equation does not imply
% that the accessible area is the geometric area).

In the continuum regime, the diffusion coefficient is frequently expressed
in terms of a mobility radius $R_m$ defined through
\begin{equation}
\label{eq:mobility_radius}
D_{k}=\f{k_BT}{6\pi\mu_f R_m} \Comma
\end{equation}
where the fluid viscosity $\mu_f$ expressed in terms of monomer properties (using the
Stokes drag law  for a monomer of hard-sphere diameter $d$,
$m_1 \beta_1 U_{\infty} = 3 \pi \mu_f \sigma U_{\infty}$, where $U_{\infty}$ is
the fluid velocity far from the monomer) is
\begin{equation}
\label{eq:fluid_viscosity}
\mu_f=\f{m_1\beta_1}{3\pi\sigma} \Teleia
\end{equation}
Thus, given a numerically calculated diffusion coefficient the mobility radius $R_m$
is obtained by inverting Eq.~(\ref{eq:mobility_radius}). The cluster mobility radii
are also reported in Table~\ref{table_diffusion}.

The cluster diffusion coefficient is frequently estimated by
taking the aggregate mobility radius $R_m$ approximately
equal to the cluster radius of gyration $R_g$.
Our simulations (cf. Table~\ref{table_diffusion}) show that the radius of
gyration tends to be of the same order of magnitude for cluster sizes in the range
 $k=4$ to $k=98$, whereas the calculated mobility radius is considerably higher
(the radius of gyration severely underestimates the  mobility radius).
As mentioned earlier, this is an inherent problem of unshielded
Langevin simulations since cluster mobility is determined by the
dynamics of its  monomers, which are always at
least partially shielded by the other monomers from the surrounding fluid.
% The recent study by \cite{moskal} suggests that the approximation may
% be extremely poor, but this approximation is often deemed reasonable
% for clusters having $d_f\le 2$ (see for instance \cite{friedlander_deposition}).

Finally, we also investigated cluster centre-of-mass velocity fluctuations
to gain insight into cluster thermalization. For a single Brownian monomer, the
late-time mean-square velocity fluctuations are \citep{risken_book}
\begin{equation}
\label{eq:delta_v_squared_monomer}
\langle \delta v_{\infty}^2 \rangle =
\lim_{t \rightarrow \infty} \langle \big [ {\bf v}(t) - \langle {\bf v}(t) \rangle \big ]^2 \rangle
 =\f{k_B T}{m_1} \Comma
\end{equation}
where ${\bf v}(t)$ is the monomer instantaneous velocity. For a $k$-mer we generalize the
previous expression, as we did for the diffusion coefficient, to
\begin{equation}
\label{eq:delta_v_squared}
\langle \delta  V_{CM, \infty}^2 \rangle =
\lim_{t \rightarrow \infty} \langle \big [ {\bf V}_{CM} (t) - \langle
{\bf V}_{CM}(t) \rangle \big ]^2 \rangle
 =\f{k_B T}{k m_1} \Comma
\end{equation}
where the cluster centre-of-mass velocity is ${\bf V}_{CM} (t) = \sum_i {\bf v}_i/k$,
% \beq \label{v-cm}
% {\bf V}_{CM} (t) = \frac{1}{k} \, \sum_1^k {\bf v}_i \Comma
% \eeq
  with ${\bf v}_i$ the $i$-th monomer velocity.
In Figure \ref{diffusion_plot} (right) we plot
$\langle \delta V_{CM}^2 (t) \rangle$ as a function of time again for the same 50-monomer
cluster considered in the diagram on the left.
The calculation of $ \delta V_{CM}^{2}$ in
Eq.~(\ref{eq:delta_v_squared}) is straightforward
  since the velocity of each
monomer is a simulation output exactly like the monomer position.
We notice that the late-time behaviour of
$\langle \delta V_{CM, \infty}^2 (t) \rangle$
fluctuates about $0.03$, the value predicted by
Eq~(\ref{eq:delta_v_squared}) for $k_B T=0.5$, $k=50$ and $m_1=1$. As in the
discussion of the cluster diffusion coefficient, the trivial cluster thermalization
arises from neglect of shielding.

The main result of this section is that from Eq.~(\ref{eq:Langevin}),
i.e. neglecting monomer shielding, a $k$-mer  is
found to have the same relaxation time of a monomer
and the same diffusive and thermal properties of a sphere  with a
mass $k$ times the one of a single monomer.
This  can be seen as a consequence of the ideality of the clusters:
the absence of shielding makes the cluster relaxation time (inverse of
the per-unit-mass cluster friction coefficient) identical to the one
of an isolated monomer and the cluster diffusional properties thus deviate from
those  of a monomer only because of the cluster larger mass.

% Hence, we can define our clusters as ideal cluster where the monomers
% are not screened by other monomers, they are transparent.
% This part of the manuscript has to be rewritten keeping this
% cluster ideality in mind.

\subsection{Cluster rotation}
\label{sec:clust-diff-rotation}

As argued, the generated cluster may be treated as a rigid bodies,
since the monomer-monomer interaction potential is very deep and cluster
re-structuring is very limited. Rigid body rotation in three dimensions may be
described by the three Euler angles $\theta$, $\phi$ and $\psi$ \citep{goldstein}.
We evaluate the cluster Euler angles during post-processing
by eliminating the cluster translational motion via rigidly
translating the aggregate centre-of-mass to the origin of the
computational box coordinate system at all times.
The Euler angles may, then be calculating by recording
the time-dependent positions
${\bf r_A}(t)$, ${\bf r_B}(t)$ and ${\bf r_C}(t)$  of three non-coplanar monomers in each cluster.
We define a 3 by 3 matrix ${\bf X}(0) =[{\bf r_A}(0), {\bf r_B}(0), {\bf r_C}(0)]$
containing the initial positions of the three reference monomers  and a
matrix ${\bf X'}(t) =[{\bf r_A}(t), {\bf r_B}(t), {\bf r_C}(t)]$ with
their positions at time $t$. The rotation matrix ${\bf A}$ such that
$\bf X'= AX$ can then be calculated as
\begin{equation}
\label{eq:rotation_matrix}
\bf{A} (t) =\bf{X}'(t) \bf{X}^T (0) \, \Big [ \bf{X} (0) \bf{X}^T(0)  \Big ]^{-1}
= \bf{X}'(t) \bf{X}^{-1} (0) \Comma
\end{equation}
where the superscript $T$ denotes matrix transpose.
Once the rotation matrix $\bf A$ is known,  the Euler angles may be determined.
As they are not uniquely defined, their values depend on the order of the three required
rotations, we chose to calculate them via the algorithm presented in \cite{rotation_alghoritm}.
The Euler angles so determined vary as follows: $\phi$ and $\psi$ both
range in the interval $[-\pi,\pi]$, whereas $\theta$ in the interval $[-\pi/2,\pi/2]$.
We investigate whether there is a preferential orientation of the
clusters by examining the distributions of the recorded Euler angles.
For uniform random rotation matrices \citep{random_euler}, both $\psi$
and $\phi$ are random variables with a uniform probability distribution
in the interval $[-\pi,\pi]$:
hence $\langle\phi\rangle=\langle\psi\rangle=0$ and
$\langle \delta \phi^2(t) \rangle = \langle [ \phi(t) - \langle \phi (t) \rangle ]^2 \rangle =
\langle \delta \psi^2 (t) \rangle={\pi}/{\sqrt3}\simeq 1.81$.
On the other hand, the Euler angle $\theta$ is distributed according to \citep{random_euler}
\begin{equation}
\label{eq:theta_distr}
\theta=\arccos(1-2Z(0,1))-\f{\pi}{2},
\end{equation}
where $Z(0,1)$ is a random variable with a uniform probability
distribution in the interval $[0,1]$. Hence, given the distribution
Eq.~(\ref{eq:theta_distr}) it is easy to show numerically that $\langle \theta (t)\rangle=0$ and
$\langle\ \delta \theta^2 (t) \rangle \simeq 0.68$.

Figure \ref{cluster_rotation} (left) shows the results for the mean
values of the Euler angles calculated from an ensemble of $800$ identical
clusters containing $10$ monomers each. As argued earlier the
numerical simulation confirm that the
expectation value of each Euler angle fluctuates around the  zero
value. On the right diagram of Fig. \ref{cluster_rotation} we
calculate the mean-square fluctuation of Euler angles for the same ensemble of
clusters. We note that they rapidly tend to the expected values
for random rotation matrices. Another calculation of
$\l\delta\theta\r$ for an ensemble of $800$ identical cluster with
$k=50$ shows saturation to the same value but on a longer time-scale.
This is not surprising since we already know from the results of the
previous section that increasing the size of a cluster reduces its mobility.
Once again, the perfect isotropy of cluster rotational properties are
due its ideality: due to the absence of shielding, each monomer in the
cluster feels the same friction and noise and this prevents any
preferential orientation/rotation angle for the cluster.
%Hence the clusters are ideal \ldots.

\subsection{Evolution of the number of clusters}

The total cluster concentration in the system,
$N_{\infty} (t) = \sum_k^{\infty} N_k(t)$, where $N_k(t)$ is
the number of $k$-mer at time $t$ is an important quantity
since it shows the progress of agglomeration and it can be used to compare the results of
the simulations with the numerical solution of the agglomeration
equation \citep{friedlander book}
\begin{equation}
\label{eq:agglomeration_equation}
\f{dn_k}{dt}=\f{1}{2}\s_{i+j=k} K_{ij}n_in_j-n_k\sum_i K_{ik}n_i,
\end{equation}
where $n_k = N_k(t)/V$ is the concentration of clusters consisting of $k$
monomers and $K_{ij}$ is the collisional kernel between
 $i$ and $j$-mers, respectively.

The standard treatment of the collision kernel for fractal aggregates
in the continuum limit (hereafter called continuum kernel), is given by \citep{friedlander book}
\begin{equation}
\label{eq:cont_kernel_general}
K_{ij}=4\pi (D_i+D_j)(R_{i}+R_{j}) \Comma
\end{equation}
where $R_i=R_{g,i}$ is the radius of gyration of an $i$-mer
(we drop the subscript $g$ to simplify the notation) and $D_i$ the fractal
aggregate diffusion coefficient. Equation~(\ref{eq:cont_kernel_general}) is
often expressed in terms of the aggregate volume of solids
$\nu_i\propto(R_{i})^{d_f}$ \cite{LorenzoJAS} ($d_f$ being a single, time-independent fractal dimension
determined in Sec.~\ref{calculate_df}),
while the diffusion coefficient is given by Eq.~(\ref{eq:mobility_radius}) with the radius
of gyration taken to be the mobility radius. The continuum kernel then reads
\begin{equation}
\label{eq:cont_kernel_standard}
K_{ij}^{Sm}=\f{2k_BT}{3\mu_f}
\lro \f{1}{\nu_i^{1/d_f}}+\f{1}{\nu_j^{1/d_f}}  \rro
\lro \nu_i^{1/d_f}+\nu_j^{1/d_f} \rro \Comma
\end{equation}
where the superscript $Sm$ refers to the so-called Smoluchowski kernel.
The Smoluchowski kernel is not expected to model our numerical
results accurately since it is based on the diffusion coefficient
Eq.~(\ref{eq:mobility_radius}) with $R_m=R_g$, an equality not
respected by our simulations, cf. Table~\ref{table_diffusion}).
Instead, the numerical simulations may be used to derive a
kernel that models more accurately the aggregate
collisional dynamics. Such a kernel is obtained by
using the numerically validated expression for the cluster
diffusion coefficient Eq.~(\ref{eq:cluster-diffusion-coefficient})
with $\beta_k=\beta_1$ in Eq.~(\ref{eq:cont_kernel_general}).
For the radius of gyration, we use the fitted value for
$k\ge 5$ and the average values directly calculated from the
simulations, i.e.
\begin{equation}
\label{eq:rules_for_r_gyr}
R_i = \begin{cases}
\textrm{simulation data}, & i< 5 \Comma \\
\alpha_{mc} \, i^{1/d_f^{mc}},  & 5\le i\le 15 \Comma \\
\alpha_{cc} \, i^{1/d_f^{cc}}, &  15 < i \Teleia \\
\end{cases}
\end{equation}

Furthermore, non-continuum effects (specify, related to the \textit{particle}
Knudsen number) can play a role in cluster collisions and
they are usually accounted for by introducing Fuchs
correction  $\beta_F$ in the kernel as in Ref. \citep{LorenzoJAS}.
Hence the kernel derived for Langevin simulations (hereafter called
Langevin-Dynamics kernel) reads
\begin{equation}
\label{eq:kernel_complete}
K_{ij}^{LD}=\f{4\pi k_BT}{m_1\beta_1} \lro\f{1}{i}+\f{1}{j}\rro \lro
R_{i}+R_{j}\rro \beta_F.
\end{equation}

The late-time dependence of the total number of cluster
depends on kernel homogeneity exponent. A kernel $K_{ij}$ is said to be homogeneous of
order $\lambda$ if $K_{\gamma i,\gamma j}=\gamma^{\lambda} K_{ij}$.
The asymptotic time-decay of the total number concentration
is $N_\infty (t) \sim t^{-1/(1-\lambda)}$.
For the Smoluchowski kernel $\lambda=0$ and, therefore, $N_\infty\sim t^{-1}$,
whereas for the Langevin-Dynamics kernel $\lambda=(1/d_f^{\rm large})-1$ and
$N_{\infty}\sim t^{-0.72}$.

The results of the simulations (with both potentials) and the numerical
solution of the agglomeration equation Eq.~(\ref{eq:agglomeration_equation})
with $K^{Sm}$ and $K^{LD}$ are shown in
Fig.~\ref{Smoluchowsky  comparison}. The Langevin-Dynamics kernel,
Eq.~(\ref{eq:kernel_complete}), leads to improved agreement at
short times between the simulations and the numerical solution of the
agglomeration equation Eq.~(\ref{eq:agglomeration_equation}).
We, also, fitted the time decay of the total cluster number to a power law,
as expected from the kernel homogeneity, $N_\infty\sim t^{-\xi}$, at
times $ 2500 \le t\le 3000$. From Eq.~(\ref{eq:agglomeration_equation}) with
the Langevin-Dynamics kernel we find $\xi=0.77$,  a value close to
the values for the discrete simulations ($\xi=0.78$ and $\xi=0.79$ for
model potential and Van der Waals interaction, respectively), whereas
the continuum kernel leads to the rather different result $\xi=1$.
The slight difference between calculated and theoretical exponents
is attributed to the slow approach of the asymptotic limit of
the agglomeration equation: the time asymptotic limit is reached
at times about one order of magnitude longer than the
duration of the simulations, primarily due to the Fuchs factor.
Finally, we note that, as expected, the Van der Waals interaction
enhances the agglomeration rate with respect to the model potential
since the former is a long-ranged interaction, though quickly decaying.

\subsection{Collision kernel elements}

Collision kernel elements are obtained from the simulation
by averaging Eq.~(\ref{eq:kernel_calc}) over the 10 different
initial conditions
\begin{equation}
\label{eq:kernel_nonlin_aver}
\f{\langle N_{ij}(t)\rangle}{\delta t}=
(2-\delta_{ij}) K_{ij}\langle n_i(t)n_j(t) \rangle,
\end{equation}
where $\langle\dots\rangle$ denotes the ensemble average.
We found that the average of the nonlinear term decouples,
namely  $\langle n_in_j\rangle = \langle n_i\rangle \langle n_j\rangle$
to the accuracy of our simulations. We recorded the collisions between
aggregates in the computational box rather frequently, every
$\delta t_{sam}$ to avoid counting ternary
collisions occurring during sampling time interval. We have enough data only
for the evaluation of a few kernel elements
$K_{ij}$ for low $i$ and $j$ indexes, but they suffice to draw some conclusions.
Figure~\ref{calculated_beta_ij}  shows an example of the fitting
procedure to determine $K_{13}$ from $N_{13}/\delta t$ and
$n_1 n_3$ that were evaluated directly from the simulations.
We determined a few kernel elements $K_{ij}$, reported in
Table~\ref{table_kernel},
 that are compared to the
analytical values of $K_{ij}^{LD}$: kernel elements expressed are scaled
to the diagonal of the continuum kernel $K^{Sm}_{ii}=8k_BT/3\mu_f$.
We note an excellent agreement between numerical and analytical
estimates for the $\{(1,1), (1,2), (2,2)\}$
kernel elements. These kernel elements are important in the early-time
evolution of $N_\infty(t)$, and therefore are responsible for
the agreement of the numerical simulations and the numerical
solution of the agglomeration equation at early times, as shown
in Fig.~\ref{Smoluchowsky comparison}.

\section{Concluding remarks}
\label{conclusions}

In this study we employed MD dynamics techniques (and in particular a
MD research code) to investigate aggregate collisional dynamics as a
function of the monomer-monomer interaction potential.
The problem of cluster identification is tackled in the
post-processing of the simulation results   by considering each cluster
as the connected components of a graph.
The time evolution of the system fractal dimension is linked to the
kinetics of two cluster populations, namely  small and large
clusters.



The diffusion coefficient of a $k$-mer is found to scale as
$D_k\propto k^{-1}$. In other words,
aggregates diffuse like massive monomers as a consequence of the
absence of any shielding effect in the Langevin equation for
interacting monomers. This is an inevitable drawback of the chosen
method, unless a shielding coefficient is explicitly introduced in
Langevin equation \ref{eq:Langevin}.

We also successfully calculate numerically and compare to the
analytical prediction the kernel elements $\beta_{ij}$ for low
indexes. An extensive numerical investigation of the collisional  is
beyond the purpose of the present study, but it can be performed using
the methodology described here.

% The Appendices part is started with the command \appendix;
% appendix sections are then done as normal sections
% \appendix

% \section{}
% \label{}

% \clearpage

\begin{thebibliography}{00}

% \bibitem{label}
% Text of bibliographic item

% notes:
% \bibitem{label} \note

% subbibitems:
% \begin{subbibitems}{label}
% \bibitem{label1}
% \bibitem{label2}
% If there is a note, it should come last:
% \bibitem{label3} \note
% \end{subbibitems}

%\bibitem{jegfrje}

% \hv{Burtscher}{2004}{burtscher review} Burtscher, H. (2004).
%Physical characterization of particulate emissions from diesel engines: a review. {\it Journal of Aerosol Science, 36}, 896-932.


\hv{Filippov}{2000}{filippov_drag} Filippov, A.V. (2000). Drag and
Torque on Clusters of N Arbitrary Spheres at Low Reynolds Number.
  {\jc} {\it
  229}, 184-195.

\hv{Filippov {\it et al.}}{2000}{filippov_thermal} Filippov,
A.V., Zurita, \& M., Rosner, D.E. (2000). Morphology and physical
properties of fractal-like aggregates.  {\jc} {\it 229}, 261-273.




\hv{Risken}{1989}{risken_book} Risken, H. (1989). {\it The
  Fokker-Planck Equation}, 2nd edition, Springer Series in
Synergetics, Berlin.


\hv{Mountain {\it et al.}}{1986}{mountain} Mountain, R.D., Mulholland, G.W., \&
Baum, H. (1986). Simulation of Aerosol Agglomeration in the Free
Molecular and Continuum Flow Regimes. {\it Journal of Colloid and
  Interface Science, 114}, 67-81.


\hv{Gutsch {\it et al.}}{1995}{pratsinis_kernel} Gutsch, A.,
Pratsinis, S.E., \& Loeffler, F. (1995).  Agglomerate structure and
growth rate by trajectory calculations of monomer-cluster collisions. {\jas} {\it 26}, 187-199.


\hv{Heine \& Pratsinis}{2007}{pratsinis_ld} Heine, M.C., \& Pratsinis,
S.E., (2007). Brownian Coagulation at High Concentration, {\it Langmuir, 23}, 9882-9890.

\hv{Maedler {\it et al.}}{2006}{friedlander_deposition}
Maedler, L., Lall, A.A., \& Friedlander S.K. (2006). One-step aerosol
synthesis of nanoparticle agglomerate films: simulation of film
porosity and thickness.  {\it
  Nanotechnology 17}, 4783-4795.

\hv{Videcoq {\it al.}}{2007}{videcoq} Videcoq, A., Han, M., Ab\'elard,
P., Pagnoux, C., Rossignol, F. Ferrando, R. (2007). Influence of the
potential range on the aggregation of colloidal particles. {\it
  Physica A, 374}, 507-516.











\hv{Cormen {\it {et al.}}}{2001}{book_algorithms} Cormen, T.H.,
Leiserson, C.E., Rivest, R.L., Stein, C. (2001). {\it Introduction to
algorithms}, 2nd edition, MIT Press, Cambdridge (Massachussets).


\hv{Kostoglou \& Konstandopoulos}{2001}{konstandopoulos} Kostoglou, M., \&
Konstandopoulos, A.G. (2001). Evolution of aggregate size and fractal
dimension during Brownian coagulation. {\it Journal of Aerosol
  Science, 32}, 1399-1420.





\hv{Meakin}{1984}{meakin_cluster_models} Meakin,
P. (1984). Diffusion-controlled aggregation on two-dimensional square
lattices: results from a new cluster-cluster aggregation
model. {\prb}, {\it 29}, 2930-2942.





\hv{clusterReguera}{2007}{clusterReguera} Wedekind, J., \& Reguera, D. (2007).
What is the best definition of a liquid cluster at the molecular scale?.
\textit{Journal of Chemical Physics}, \textit{127}, 154516.





\hv{Frenkel \& Smit}{2002}{md_book} Frenkel, D., \& Smit,
B. (2002). {\it Understanding Molecular Simulation}, 2nd edition,
Academic Press, San Diego.






\hv{Lazaridis \& Drossinos}{1998}{yannis_potential} Lazaridis, M., \&
Drossinos, Y. (1998). Multilayer Resuspension of Small Identical
Particles by Turbulent Flow.
 {\it Aerosol Science and Technology, 28}, 548-560.





\hv{Narsimhan \& Ruckenstein}{1985}{ruckenstein} Narsimhan, G., \&
Ruckenstein, E., (1985). {\jc}, {\it 104}, 344-369.
  

\hv{Rothenbacher {\it et al.}}{2008}{kasper_hamaker}
Rothenbacher, S., Messerer, A., \& Kasper, G., (2008). Fragmentation and bond strength of airborne diesel soot agglomerates. {\it Particle
  and Fiber Toxicology, 5}, 1-7.

  

\hv{Poling {\it et al.}}{2000}{substance_book} Poling,
B.E., Prausnitz, J.M., \& O'Connell, J.P. (2000). {\it The properties
  of gases and liquids}, McGraw-Hill, New York.



\hv{Limbach {\it {et al.}}}{2006}{espresso} Limbach, H.J., Arnold A.,
Mann, B.A., \& Holm, C. (2006). ESPResSo--an extensible simulation
package for research on soft matter systems. {\it Computer Physics
  Communications, 174}, 704-727.



\hv{Csardi \& Nepusz}{2006}{igraph}  Csardi, G., \&  Nepusz,
T. (2006). The igraph software package for complex network
research. {\it Inter Journal Complex Systems}, 1695.





\hv{Sorensen \& Roberts}{1997}{sorensen_prefactor} Sorensen, C.M. \&
Roberts, G.C., (1997). {\jc}, {\it 186}, 447-452.


\hv{Smirnov}{1990}{smirnov_review} Smirnov, B.M. (1990). The
properties of fractal clusters. {\it Physics Reports, 188}, 1-78.

\hv{Tanaka \& Araki}{2000}{prl_hydro} Tanaka, H., \& Araki, T.,
(2000). Simulation Method of Colloidal Suspensions with Hydrodynamic
Interactions: Fluid Particle Dynamics.  \prl, {\it 85}, 1338-1341.




\hv{Jullien \& Meakin}{1989}{meakin_reorganization_3D} Jullien, R., \&
Meakin, P. (1989). Simple models for
the restructuring of three-dimensional ballistic aggreagates. \jc,  {\it
127}, 265-272.







\hv{Goldstein}{1950}{goldstein} Goldstein, H. (1950). {\it Classical
  mechanics}, Cambridge, Massachussets.




\hv{Shoemake}{1994}{rotation_alghoritm} Shoemake, K.  (1994)
\emph{Euler Angle Conversion}, from "Graphics Gems IV",  222-229. Morgan Kaufmann.

\hv{Kuffner}{2004}{random_euler} Kuffner, J.,J., (2004). {\it Effective Sampling and
  Distance Metrics for 3D Rigid Body Path Planning}, Proceeding of the
2004 IEEE International Conference on Robotics and Automation (ICRA 2004).

\harvarditem{Friedlander}{2000}{friedlander book} Friedlander S.K. (2000).
{\it Smoke, Dust and Haze},  2nd Edition,      Oxford University Press, New York.





\hv{LorenzoJAS}{2008}{LorenzoJAS} Isella, L., Giechaskiel, B., \& Drossinos, Y. (2008).
Diesel-exhaust aerosol dynamics from the tailpipe to the dilution tunnel.
\textit{Journal of Aerosol Science}, {\it 39}, 737-758.


% \hv{POV-Ray}{2004}{povray} Persistence of Vision Pty. Ltd. (2004).
%    \url{http://www.povray.org/}.


\hv{mayavi2}{2004}{mayavi2} Ramachandran, P., \& Varoquaux, G. (2008). Mayavi: making 3D data 
visualization reusable, \textit{SciPy08: Proceedings of the 7th Python in 
Science Conference}, Caltech, Pasadena, CA, 19--24 August.








%%%%%%%%%%%%%%%%%%%%%%%%

% \hv{Hinds}{1999}{Hinds} Hinds, C. H. (1999). {\it Aerosol Technology},
% 2nd edition, John Wiley \& Sons, New York.

% \hv{Meakin}{1984}{meakin_trajectories} Meakin, P. (1984). Computer
% simulations of cluster-cluster aggregation using linear trajectories:
% results for three-dimensional simulations and a comparison with
% aggregates formed using Brownian trajectories. {\it Journal of Colloid
% and Interface Science, 102}, 505-512.

% \hv{Meakin}{1984}{meakin_trajectories2} Meakin,
% P. (1984). Effects of cluster trajectories on cluster-cluster
% aggregation: a comparison of linear and Brownian trajectories in two-
% and three-dimensional simulations. \pra, {\it 29}, 997-999.

% \hv{Mulholland {\it et al.}}{1988}{mulholland_free_molecular}
% Mulholland, G.W., Samson, R.J., Mountain, R.D., \& Ernst,
% M.H. (1988). Cluster size distribution for free molecular
% agglomeration. {\it Energy \& Fuels, 2}, 481-486.



% \hv{Meakin}{1986}{meakin_reorganization_2D} Meakin, P. (1986). The
% effects of reorganization processes on two-dimensional cluster-cluster aggregation.
%  {\it Journal of Colloid and Interface Science, 112}, 187-194.



% \hv{Vemury \& Pratsinis}{1995}{pratsinis_self_preserving} Vemury,
% S. \& Pratsinis, S.E. (1995). Self-preserving size distributions of
% agglomerates. {\jas}, {\it 26}, 175-185.


% \hv{Wu \& Friedlander}{1993}{cont_kern_fractal} Wu, M.K., \&
% Friedlander, S.K. (1993). Enhanced power law agglomerate growth in the
% free molecule regime.
%  {\jas}, {\it 24}, 273-282.

% \hv{Moskal \& Payatakes}{2006}{moskal} Moskal, A., \& Payatakes,
% A.C. (2006). Estimation of the diffusion coefficient of aerosol particle aggregates
%        using Brownian simulation in the continuum regime.
%  {\jas}, {\it 37}, 1081-1101.

% \hv{Brasil {\it {et al.}}}{2001}{brasil_numerical} Brasil, A.M.,
% Farias, T.L., Carvalho, M.G., Koylu, U.O. (2001). Numerical characterization of the morphology
%            of aggregated particles.
%  {\jas}, {\it 32}, 489-508.


% \hv{Fuchs}{1989}{fuchs book} Fuchs, N.A. (1989). {\it The mechanics of
% aerosol}, Dover edition, Dover Publications, New York.



% \hv{Lattuada {\it et al.}}{2003}{lattuada_hydro} Lattuada, M.,
% Wu, H., \& Morbidelli, M. (2003). Hydrodynamic radius of fractal clusters. {\jc}, {\it 268}, 96-105.




% \hv{Garcia-Ybarra {\it et al.}}{2006}{ybarra_drag}
% Garcia-Ybarra, P.L, Castillo J.L., \& Rosner, D.E. (2006). {\jas} {\it
% 37}, 413-428.


% \hv{Biswas, \& Kulkarni}{2004}{deposition_interaction} Kulkarni, P.,
% \& Biswas, P. (2004). A Brownian Dynamics Simulation to Predict
% Morphology of Nanoparticle Deposits in the Presence of Interparticle
% Interactions.  {\it Aerosol Science and Technology, 38}, 541-554.

% \hv{Hutter}{2000}{hutter_langevin} Hutter, M. (2000). Local Structure
% Evolution in Particle Network Formation Studied by Brownian Dynamics
% Simulation. {\jc} {\it 231}, 337-350.





% \hv{Filippov {\it et al.}}{2000}{filippov-review} Filippov, A.V., Zurita, M.,  \& Rosner,
% D.E. (2000).  Fractal-like Aggregates: Relation between Morphology and Physical Properties.  
%  {\it Journal of  Colloid and Interface Science}, {\it 229} 261-273.        



\end{thebibliography}

 \clearpage

\begin{table}
\caption{Overview of the cluster diffusive properties}
\bigskip

\begin{center}
\begin{tabular}{|c|c|c|c|c|c|c|} \hline
& $k=4$ & $k=10$ & $k=50$ & $k=98$  \\ \hline \hline
$D_k$ &$1.30\cdot 10^{-1}$ & $4.92\cdot 10^{-2}$ & $9.48\cdot 10^{-3}$
& $4.84\cdot 10^{-3}$  \\ \hline
$R_m$ &$1.92$ & $5.08$ & $26.35$ & $51.56$ \\ \hline
$R_g$ &$7.5\cdot 10^{-1}$ & $1.14$ & $2.12$ & $4.70$  \\ \hline
%$R_H$ &$1.05$ & $1.34$ & $2.22$ & $3.65$  \\ \hline
$\beta_k$ &$9.6\cdot 10^{-1}$ & $1.01$ & $1.05$ & $1.05$  \\ \hline
\hline
\end{tabular}
\label{table_diffusion}
\end{center}
\end{table}


\clearpage

\begin{table}
\caption{Comparison of kernel elements. All kernel elements are given
  in units of $8k_{B}T/3\mu_{f}$.}
\bigskip

\begin{center}
\begin{tabular}{|c|c|c|c|c|c|c|} \hline
$i,j$ & {$K_{ij}^{LD}\; \rm  [Eq.\; (\ref{eq:kernel_complete})]$} & {$K_{ij}\; \rm  
  [Eq. \; (\ref{eq:kernel_nonlin_aver})]$}   \\ \hline \hline
${1,1}$ &$0.213$ & $0.214$   \\ \hline
${1,2}$ &$0.294$ & $0.298$  \\ \hline
${1,3}$ &$0.309$ & $0.397$   \\ \hline
%$R_H$ &$1.05$ & $1.34$ & $2.22$ & $3.65$  \\ \hline
${1,4}$ &$0.309$ & $0.449$  \\ \hline
${2,2}$ &$0.317$ & $0.317$  \\ \hline
${2,3}$ &$0.303$ & $0.349$  \\ \hline
\hline
\end{tabular}
\label{table_kernel}
\end{center}
\end{table}


%\end{document}

% Figure 1
\clearpage   
\begin{figure}
\includegraphics[width=\columnwidth]{NewFigs/presentation_potential_inset.pdf}
\caption{Model potential and monomer-monomer interaction potential used in the
simulations (energies are in units of $k_BT^*$ and distances in  units of $\sigma$).}
\label{plot_potential}
\end{figure}

% Figure 2
\clearpage
\begin{figure}
\includegraphics[width=0.5\columnwidth,height=0.5\columnwidth]{NewFigs/final2.png}
\includegraphics[width=0.5\columnwidth,height=0.5\columnwidth]{NewFigs/50_monomers_neighbor1.png}
%\vspace*{0.4cm}
\includegraphics[width=0.5\columnwidth,height=0.5\columnwidth]{NewFigs/50_monomers_neighbor2.png}
\includegraphics[width=0.5\columnwidth,height=0.5\columnwidth]{NewFigs/50_monomers_neighbor3.png}
%\includegraphics[width=0.5\columnwidth, height=6cm]{figure8_b.pdf}
\caption{Clockwise from top: snapshot of the system at the end of a
simulation, followed by three different  aggregates each
consisting of $50$ monomers. Images created with Mayavi2
software for three-dimensional rendering  \cite{mayavi2}.}
\label{snapshot_end_simulation}
\end{figure}

% Figure 3
\clearpage
\begin{figure}
\includegraphics[width=\columnwidth]{NewFigs/figure_two_slopes_extended_complete.pdf}
\caption{Time-independent mean radius of gyration versus cluster
size. Linear fits performed on a double-logarithmic scale
for $k>5$ (long-dashed line), $5\le k\le 15 $ (short-dashed line)
and $k>15$ (solid line).}
\label{two_df}
\end{figure}


% Figure 4
\clearpage
\begin{figure}
\includegraphics[width=\columnwidth]{NewFigs/evolution_df_and_number_clusters.pdf}
% \includegraphics[width=0.5\columnwidth]{evolution_df.pdf}
%\includegraphics[width=0.5\columnwidth, height=6cm]{figure8_b.pdf}
\caption{Top: Time-dependent fractal dimensions as a function of
time. Bottom: evolution of the total number of clusters and
contributions from small and large clusters.}
\label{evolution_df}
\end{figure}

% Figure 5
\clearpage
\begin{figure}
\includegraphics[width=\columnwidth]{NewFigs/evolution_coordination_number.pdf}
% \includegraphics[width=0.5\columnwidth]{evolution_df.pdf}
%\includegraphics[width=0.5\columnwidth, height=6cm]{figure8_b.pdf}
\caption{Evolution of the mean (ensemble-averaged) coordination number.}
\label{evolution_coord_number}
\end{figure}

% Figure 6
\clearpage
\begin{figure}
\includegraphics[width=0.45\columnwidth]{NewFigs/mean_sep_single_cluster_select.pdf}
\includegraphics[width=0.45\columnwidth]{NewFigs/r_gyr_single_cluster_select.pdf}
\caption{Left: evolution of the mean monomer-monomer distance for a
particular cluster in the simulations. Right: evolution of the
 radius of gyration for the same cluster.}
\label{single_cluster_size_evolution}
\end{figure}


% Figure 7, a and b
\clearpage
\begin{figure}
\includegraphics[width=0.45\columnwidth]{NewFigs/inset_plot.pdf}
\includegraphics[width=0.45\columnwidth]{NewFigs/delta_v_squared.pdf}
\caption{Left: time-dependence of the ensemble-averaged mean-square
displacement (crosses) and linear fit (solid line). Inset:
early-time behaviour of $\langle \delta R^2_{CM}  \rangle$ (crosses),
linear fit (solid line)
and power-law fit (dashed line). Right: velocity fluctuations $\langle \delta V^2_{CM}  \rangle$  as a
function of time (crosses) and analytical expression in
Eq.~(\ref{eq:delta_v_squared}) (solid line). }
\label{diffusion_plot}
\end{figure}

%Figure 8,a  and b
\clearpage
\begin{figure}
\includegraphics[width=0.45\columnwidth]{NewFigs/mean_angle.pdf}
\includegraphics[width=0.45\columnwidth]{NewFigs/delta_angle_alpha.pdf}
%\includegraphics[width=0.5\columnwidth, height=6cm]{figure8_b.pdf}
\caption{Left: mean values of the Euler angles $\l\theta_{10}\r$ (solid line),
$\l\phi_{10}\r$ (dashed line) and
$\l\psi_{10}\r$ (dot-dashed line) calculated for an ensemble of
clusters with $k=10$. Right: standard deviation of Euler
angles $\l\delta\theta_{10}\r$ (solid line), $\l\delta\phi_{10}\r$
(long-dashed line), $\l\delta\psi_{10}\r$ (short-dashed line) for the
same clusters as in the left diagram and $\l\delta\theta_{50}\r$
(dot-dashed line) for a
simulation of clusters with $k=50$. The horizontal lines are the
theoretical values
 of $\l\delta\psi(\delta\phi)\r$
and $\l\delta\theta\r$ for random rotation matrices.}
\label{cluster_rotation}
\end{figure}

% Figure 9
\clearpage
\begin{figure}
\includegraphics[width=\columnwidth]{NewFigs/smolu_comparison_final.pdf}
%\includegraphics[width=0.5\columnwidth]{coagulation_coefficients.pdf}
%\includegraphics[width=0.5\columnwidth, height=6cm]{figure8_b.pdf}
\caption{Evolution of the total number of clusters calculated with 
the monomer-monomer interaction potential (diamonds), model potential (filled circles),
solution of the agglomeration equation with Smoluchowski kernel
(long-dashed line) and with Langevin kernel (short-dashed line).}
\label{Smoluchowsky comparison}
\end{figure}

% Figure 10
\clearpage
\begin{figure}
\includegraphics[width=\columnwidth]{NewFigs/plot_kernel_elements_13.pdf}
%\includegraphics[width=0.45\columnwidth]{figs/kernel_matrix_normalized.pdf}
\caption{Determination of kernel element  $K_{13}$. Diamonds: $N_{13}/\delta
t$. Circles: $2K_{13}n_in_3$  with the fitted
value of  $K_{13}$.}
\label{calculated_beta_ij}
\end{figure}


\end{document}