myprojects-hg-reborn

Fork

% Template article for preprint document class `elsart'
% SP 2006/04/26

%  \documentclass{elsart}

% Use the option doublespacing or reviewcopy to obtain double line spacing
  \documentclass[reviewcopy]{elsart}

% if you use PostScript figures in your article
% use the graphics package for simple commands
% \usepackage{graphics}
% or use the graphicx package for more complicated commands
%  \usepackage{graphicx}
% or use the epsfig package if you prefer to use the old commands
% \usepackage{epsfig}


 %\usepackage[pdftex]{thumbpdf}
%\usepackage[pdftex]{hyperref}



\usepackage{graphicx}
\usepackage{url}
\usepackage{natbib}
%\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{\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{dep}}}}
\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}

\usepackage{amsmath}

% The amssymb package provides various useful mathematical symbols
\usepackage{amssymb}

% The lineno packages adds line numbers. Start line numbering with
% \begin{linenumbers}, end it with \end{linenumbers}. Or switch it on
% for the whole article with \linenumbers.
% \usepackage{lineno}

% \linenumbers


\begin{document}

\begin{frontmatter}

% Title, authors and addresses

% use the thanksref command within \title, \author or \address for footnotes;
% use the corauthref command within \author for corresponding author footnotes;
% use the ead command for the email address,
% and the form \ead[url] for the home page:
% \title{Title\thanksref{label1}}
% \thanks[label1]{}
% \author{Name\corauthref{cor1}\thanksref{label2}}
% \ead{email address}
% \ead[url]{home page}
% \thanks[label2]{}
% \corauth[cor1]{}
% \address{Address\thanksref{label3}}
% \thanks[label3]{}

\title{Nanoparticle collisional dynamics by Langevin simulations}

% use optional labels to link authors explicitly to addresses:
% \author[label1,label2]{}

% \address[label1]{}
% \address[label2]{}

\author{L. Isella and Y. Drossinos}

\address{European Commission Joint Research Centre}

\begin{abstract}
% Text of abstract
The process of nanoparticle agglomeration as a function of the
monomer-monomer interaction potential
 is simulated by solving via Langevin
equations for a set of interacting monomers in three dimensions.   The simulation output  is used to investigate the
structure of the generated clusters and the collision
frequency between small clusters. Cluster restructuring is also
observed and discussed.
We identify a time-dependent fractal dimension whose evolution is linked
to the kinetics of two cluster populations.
The absence of screening  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}

% main text
\section{Introduction}
\label{introduction}
Nanoparticle aggregates are of paramount importance in 
technological and industrial processes such as combustion, filtration,
 gas-phase-particle synthesis and many more.
 The fractal nature of these aggregates 
 has profound implications on their
transport  \citep{filippov_drag, ybarra_drag, moskal} and thermal
\citep{filippov_thermal} properties.


Fractal aggregates arise from the agglomeration of smaller spherules,
hereafter called monomers, which do not coalesce, but rather retain
their identity in the resulting aggregate.
Individual monomers in a quiescent fluid are Brownian particles whose
dynamics is described by 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, deposition_interaction}.

In this study we investigate nanoparticle dynamics relying solely on
Langevin equations for a set of interacting monomers in three
dimensions.
Unlike the  works mentioned above, no assumptions are made about the
structure or mobility of the aggregates generated by the dynamics.
This is reminiscent of other applications of Langevin simulations in
the field of dilute colloidal suspensions to study agglomeration and the
structures it gives rise to  \citep{videcoq, hutter_langevin}. 
The main limitation of this approach is the lack of screening of the
inner monomers in an aggregate. 
We  quantify the effect of this
approximation  on the aggregate diffusional properties, while we argue
that it has a limited effect on the structure of the generated aggregates.

Since we solve numerically
Langevin equations for interacting monomers rather than aggregates,
the raw output of the simulations 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  Langevin nanoparticle simulations. 
Emphasis is given to the description and justification of the
monomer-monomer interaction potentials employed in the numerical
experiments.
Section \ref{simulations} offers an overview of the numerical work and  introduces
 the quantities of the interest monitored to investigate
the system dynamics.
Section \ref{results} contains the results and discussion  of
Langevin simulations, whereas the final remarks in Section \ref{conclusions} conclude
the paper.

 
% For instance, the calculation of the drag force felt by an aggregate
% in a Stokes flow typically requires the solution of Stokes equation
% for the flow surrounding the aggregate to be matched with the solution
% of Brinkman equation for the creeping flow inside the aggregate
% \citep{filippov_drag, ybarra_drag}, which can be a numerically
% demanding task for which analytical  treatments are available only in
% the case of spherically-symmetric fractal aggregates.





\section{Model Formulation}\label{model}

\subsection{Langevin Equation for Mesoscopic systems}
\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 Brownian particle dynamics by
solving the  equations for a system of interacting clusters
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),
\end{equation}
where $m_1$ in the  is the monomer mass, ${\bf r}_i$ its position in
three-dimensional space, ${ \bf F}_i$ the force acting on
the monomer, $\beta_1$ is the inverse of the monomer relaxation time $\tau_1$ and
${\bf W}_i=(W_{i}^x,W_{i}^y,W_{i}^z)$ is a noise term expressing the collisions between the monomer
and the molecules of the surrounding fluid along the $x$, $y$ and $z$
direction.
The noise is assumed to be Gaussian and hence it has the properties:
\begin{equation}
  \label{eq:gaussian_noise}
\langle  W^{j}_i(t) \rangle=0\;\;\;\; {\rm{and}} \;\;\;\; \langle
 W^{j}_i(t) W_{i'}^{j'}(t')\rangle=\Gamma\delta_{ii'}\delta_{j
j'}\delta(t-t'),
\end{equation}
where $j=x,y,z$ and the noise strength $\Gamma$ is fixed by the
fluctuation-dissipation theorem to be $\Gamma=2\beta_1 m_1k_BT$ 
%\begin{equation}
%  \label{eq:noise_strength}
%  \Gamma=2\beta mk_BT,
%\end{equation}
 with $k_B$ being the Boltzmann constant and $T$ the system
 temperature.

The force term in Eq. (\ref{eq:Langevin}) is assumed to derive from a
pair-wise additive radial potential:
\begin{equation}
  \label{eq:potential_pairwise}
  {\bf F}_i=-\nabla_{{\bf r}_i} U_i,
\end{equation}
where
\begin{equation}
  \label{eq:pairwise_pot_explicit}
  U_i=\f{1}2{}\s_{j\neq i}u(r_{ij})
\end{equation}
and $r_{ij}=|{\bf r}_i-{\bf r}_j|$.



% Langevin equation is often used both in the context of mesoscopic and
% microscopic systems, but with different physical meanings.
Langevin simulation can arise in different contexts, such as in
molecular dynamics (MD) simulations and in the study of mesoscopic systems.
 In
MD simulations, a microscopic description of the system is provided and    the so-called Langevin
thermostat introduces  a fictitious viscous force according to Eq. \ref{eq:Langevin}  
  to model the  coupling  the system with
  a thermal bath.
This is used to slow down hot molecules and speed up
cold molecules in order to get the temperature of the system to
fluctuate around the temperature of the thermal bath.
%   In the Langevin thermostat it is not necessary to apply
% Eq. \ref{eq:Langevin} at each time step while evolving the system, an
% alternative being 
% to apply a stronger damping less often during the time evolution.

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.
%  Another important difference is that
% in molecular dynamics simulations the effect of Eq. \ref{eq:Langevin}
% is only to slow down hot molecules and speed up cold molecules,
% therefore there is no need to apply it at each time-step.    and introduce
% a fictitious friction in the system to model its coupling to a thermal
% bath.
% Equation \ref{eq:Langevin} is then
% a coarse-grained model for the system dynamics and  has to be used
% at all times.  
\subsection{Dimensionless Formalism}
Equation \ref{eq:Langevin} has to be cast in a convenient
dimensionless form before being solved numerically. By fixing 
length, $l^*$,
time, $\tau^*$, and   mass, $m^*$, units 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.
\end{equation}
% The quantities in Eq. \ref{eq:dimensionless-scales} also fix a
% 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{align}
  \label{eq:Langevin-dimensionless-componentwise}
\nonumber
  \f{d^2\tilde r_i}{d\tilde
    t^2}&=-\f{1}{\tilde
    m_1}\f{\p\tilde U}{\p \tilde r_i}-\tilde\beta_1\f{d\tilde r_i}{d\tilde
    t}+\f{1}{\tilde m_1}\tilde W_i^l(\tilde t) \\
  \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, \;\;\;\;\; \langle
\tilde W_i^l(\tilde t)\rangle=0,
\end{align}
% where we defined $\tilde\beta=\beta\tau^*$,  a characteristic
% temperature of the system $k_BT^*=m^*(l^*)^2/(\tau^*)^2$, 
% the dimensionless monomer-monomer potential $\tilde U=U/(k_BT^*)$ and 
% the dimensionless noise $\tilde W_i^l=W_i^l(\tau^*)^2/(m^*l^*)$ obeying the correlations
where we introduced:
\begin{align}
  \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 
\end{align}



We notice that the characteristic temperature of the system $T^*$ is
fixed by the time, length and mass units.
Equation \ref{eq:Langevin-dimensionless-componentwise} shows that
there are
three independent dimensionless variables in the system, namely $\tilde m_1$,
$\tilde \beta_1$ and $\tilde T$,  which completely characterise Langevin equation.
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 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^*$ can be evaluated for 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
relaxation time can be estimated as $\tau_1=\sigma^2\rho_p/(18\mu_f)$
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$.
%  In the following, we will drop the tilde on the
% dimensionless quantities in order not to burden the notation, unless
% otherwise stated.


\subsection{Interaction potential}
\label{explain_potential}
The monomer-monomer interaction potential is chosen in order to mimick the
interaction of hard-core monomers sticking upon collision. The
condition of perfect sticking of any two monomers when their relative
distance falls below the monomer diameter was also used in the early studies by \cite{mountain} and
\cite{meakin_cluster_models}. In their works the authors did not use a
monomer pairwise interaction potential, but examined the system
frequently during the evolution and looked for
agglomeration events (monomer-cluster or cluster-cluster) by calculating the relative distances among the
monomers in the system. After e.g. two clusters had touched, the
relative distances of the monomers in the resulting cluster were
``frozen'' and the cluster was allowed to diffuse with a diffusion
coefficient which had to be prescribed a priori.
% In this study, we investigate the agglomeration process as a function
% of the monomer-monomer interaction potential.
% Indeed, we use two different potentials: a model potential which is a
% short-ranged potential well and a long-ranged (though quickly
% decaying) interaction modelling the attractive part of Van der Waals
% potential between nanoparticles.
%  occurs by means of the Brownian motion getting monomers and
% clusters in touch  and to a monomer-monomer interaction potential
% which mimicks the condition of sticking upon collision.

In this study investigate the monomer collision dynamics using a model
potential (a short-ranged potential well) and an integrated
Lennard-Jones potential modeling the attractive part of Van der Waals
interaction (from now on simply called Van der Waals potential) between two soot nanoparticles.
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.

 
The model potential also exhibits a deep and  narrow
  attractive part, responsible for the sticking of monomers upon
  collision, while smoothly  going 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}.
 The model potential is attractive on a length scale much smaller than
 the monomer diameter, whereas   
 Van der Waals potential is a long-ranged, though
quickly decaying, interaction.
% A simple potential satisfying the above requirement is a very narrow
% and deep potential well. However, the crude implementation of a potential well would lead
% to a number of numerical problems, mainly due to the discontinuity of
% the potential and of its first derivative.
%  In order to overcome these problems, we introduce a smooth
%  monomer-monomer interaction potential in the following 
 
In the following, we give the analytical expression of the
monomer-monomer model interaction potential $u_M(r)$ used in
the numerical experiments. 

%  and has a constant gradient for
% $r\le \sigma$ to model a strong (constant) repulsive force when the physical
% distance between two monomers falls below $\sigma$.
\begin{equation}
\label{eq:potential_well}
u_M(r) = \left\{
\begin{array}{rl}
 & -ar+b,  \text{if }  0<r< \sigma\\
& u_{\sigma}-2\Delta \lcu \cos^2\lsq \f{\pi}{2\omega_1}\lro r_{\rm
  min}-r \rro-\f{1}{2}  \rsq \rcu,   \text{if } \sigma<r<r_{\rm min}\\
 & u_{\rm min}\cos^2\lsq \f{\pi}{2\omega_2}\lro r-r_{\rm min} \rro \rsq,
  \text{if } r_{\rm min}<r<r_{\rm cut} \\
& 0,  \rm{elsewhere,}
\end{array} \right.
\end{equation}
where  
\begin{align}
  \label{eq:omegas}
  & \omega_1=2(r_{\rm min}-\sigma), \;\;\;\;\; \omega_2=r_{\rm
    cut}-r_{\rm min}, \;\;\;\; b=h_{\rm max}-a\sigma \\
& a=-\f{2\Delta\pi}{\omega_1}\cos\lsq\f{\pi (r_{\rm
    min}-\sigma)}{2\omega_1}\rsq
\sin\lsq \f{\pi(r_{\rm min}-\sigma)}{2\omega_1}  \rsq , \;\;\;\;
\Delta=h_{\rm max} - h_{\rm min} . 
\end{align}

The model potential has a minimum located at $r_{\rm min}=1.05\sigma$, where it
evaluates to $u_{\rm min}=-100k_BT$, thus ensuring a sticking
probability which can be safely taken to be $1$. 
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 prolonged linearly
with the same 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_{cut}=1.1\sigma$ 

We introduce Lennard-Jones (LJ) potential between two molecules  as in \cite{yannis_potential}
\begin{equation}
  \label{eq:LJ_ruckenstein}
  V_{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 LJ potential and $\sigma_{LJ}$ is the distance
at which LJ potential evaluates to zero.
We express the potential obtained from the integration on two
equal-sized spheres of Eq. \ref{eq:LJ_ruckenstein} in an equivalent,
but more compact expression than in Refs.
\cite{ruckenstein} and \cite{yannis_potential} 
\begin{align}
  \label{eq:ruckenstein potential}
\nonumber
%   u_{\rm NR}(R) &=-\f{A}{720\cdot
%     21}\lsq\f{2R/d}{(d-2R)^6}-\f{2R/d}{(d+2R)^6}+\f{1}{720\cdot
%     105}\f{1}{d}\lro\f{2}{d^5}-\f{1}{(d-2R)^5}-\f{1}{(d+2R)^5}\rro \right. \\ \nonumber
%  &-\f{1}{90\cdot
%    28}\lro\f{R^2/d}{(d-2R)^7}+\f{R^2/d}{(d+2R)^7}+\f{2(R^2/d)}{d^7}\rro+\f{1}{6}
% \lro\log\lro\f{(d+2R)(d-2R)}{d^2}\rro\right. \\ 
% & \left. \left. \left. +\f{2R^2}{(d+2R)(d-2R)}+\f{2R^2}{d^2}\rro\rro
%       \rsq , 
  u_{\rm NR}(r) =&-A\lcu
  \sigma_{LJ}^{6}\lsq\f{\sigma}{7560r}\lro \f{1}{(r-\sigma)^6}-
\f{1}{(r+\sigma)^6}\rro \right. \right. \\ \nonumber
& +  \f{1}{37800r}\lro\f{2}{r^5}-\f{1}{(r-\sigma)^5}-\f{1}{(r+\sigma)^5}\rro
  \\ \nonumber
 &
 \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{align}
where  $A=4\pi\epsilon\sigma_{LJ}^6n^2$ is Hamaker constant and
$n$ is the molecular number density in the solid. For a typical $\sigma=20$nm n-hexane soot nanoparticle, one can
estimate $A=2.38\cdot 10^{-19}$J \citep{kasper_hamaker} and a
corresponding $\sigma_{LJ}=0.5949$nm \citep{substance_book}.

The repulsive part of  $u_{\rm NR}$ diverges for $r\to\sigma$ thus
potentially giving rise to numerical difficulties in the simulations. 
 We therefore modify the potential in Eq. \ref{eq:ruckenstein potential} at
distances $r_{\rm sep}=1.15\sigma$  (hence below the  minimum position   of
 $r_{\min}\simeq 1.17\sigma$) by
prolongating it linearly with the same slopes it exhibits at $r=r_{\rm
sep}$ until $r=0.995\sigma$ (where it evaluates to about $60k_bT^*$).
We set Van der Waals potential to zero at lower monomer-monomer
separations.
The modification of the repulsive part of the integrated LJ potential
is not expected to affect the dynamics of the system, since it  
involves only monomer-monomer separations are energetically unfeasible.
  Hence the potential we use to simulate Van der
Waals forces  is given by 
\begin{equation}
\label{eq:van_der_waals}
u_{VdW}(r) = \left\{
\begin{array}{rl}
 & u_{\rm NR}(r),  \text{if }  r>r_{\rm sep}\\
& \lro\f{\p u_{\rm NR}}{\p r}\bigg{|}_{r_{\rm sep}}\rro(r_{\rm
  sep}-r)+u_{\rm NR}(r_{\rm sep}),   \text{if }
0.995\sigma<r\le r_{\rm sep} \\
& 0, {\rm elsewhere}.
  %  & u_{\rm min}\cos^2\lsq \f{\pi}{2\omega_2}\lro r-r_{\rm min} \rro \rsq,
%   \text{if } r_{\rm min}<r<r_{\rm cut} \\
% & 0,  \rm{elsewhere,}
\end{array} \right. 
\end{equation}
We truncate Van der Waals potential at distances $r=7\sigma$ where it
is absolutely negligible with respect to the thermal energy scale
since $u_{VdW}(7\sigma)/(k_BT^*)\sim 10^{-4}$.
In the following, we will be using only dimensionless quantities, but
we will drop the tilde in order not to burden the notation. Furthermore,
unless specified otherwise, the results presented will be those
obtained with $u_{\rm VdW}$.
A plot of the dimensionless radial potentials  given by
Eqs. \ref{eq:potential_well}  is \ref{eq:van_der_waals} is
shown in Fig. \ref{plot_potential}.  

\section{Simulations and Post-Processing}\label{simulations}


\subsection{Preparation of the initial state}
% \label{sec:prep-init-state}

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 ensure a given initial monomer
density $N_{\infty}(0)=0.01$ according to:
\begin{equation}
  \label{eq:box_size}
  L=\lro\f{5000}{N_{\infty}}\rro ^{1/3},
\end{equation}
which leads to $L\simeq 80$.
The initial random displacement of the 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.
In order to avoid these problems, we use a well-known technique of
MD simulations and we ramp up the repulsive force
monomers feel when $r<1$  to its constant final value  during a time $\tau_{\rm
  ramp}=1$.
The system is then 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 by averaging the
output of $10$ simulations,  in
order to ensure that we have enough data for the statistical analysis. 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 Detection}
\label{sec:cluster-calculation}
The determination of the clusters is one of the most time-consuming
tasks of the post-processing of the simulation results.
The simulations we ran with \esp return the individual monomer
positions and velocity obtained by solving Eq. \ref{eq:Langevin}.
  Unlike the previously mentioned works
\citep{meakin_cluster_models, mountain},
we do not have to  look for agglomeration events while
evolving the system.  However, we are left with the problem of
determining a posteriori which clusters formed during the system evolution. 


In this study 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. 
For the specific case of our simulations, since both the model
potential Eq. \ref{eq:potential_well} and Van der Waals potential Eq. \ref{eq:van_der_waals} are
not infinitely narrow, monomers are allowed tiny oscillations around
the equilibrium point.
Despite of that, the bottom of the potential well is so deep that when
two monomers collide, they remain bound together.
As a consequence, one can consider two monomers as bound when their
relative distance falls below
 a threshold distance  $d_{\rm thr}$.
Of course, the choice of $d_{\rm thr}$ depends on the position of the minimum of
monomer interaction potential $r_{\rm min}$. For the parameters of the
Van der Waals
simulations, 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}$ 
around this value do not affect the determination of the clusters in the system.
 

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
distance between two monomers along the $x$ axis  is given by
\begin{equation}
  \label{eq:distance_calc}
  D_{ij}^{(x)}= x_i-x_j-L\cdot{\rm nint}\lro\f{x_i-x_j}{L} \rro,
\end{equation}
and the same holds of course for the distance along the $y$ and $z$ axis.
One notices that the periodicity of the box imposes a cutoff on the
maximum distance between two monomers along each axis, namely $L/2$.
Furthermore, $D_{ij}^{(x)}=-D_{ji}^{(x)}$, 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 calculation needs to take into  account
% the periodicity of the system. \esp return by default particle
% cordinates outside the box length $L$.  Then the interparticle
% distance along e.g. the $x$ axis used for the calculation of the
% interparticle interaction is the distance between the $i$-th particle
% and the nearest image of the $j$-th particle and  is given by:
 
The three-dimensional distance is of course 
\begin{equation}
  \label{eq:distance_3D}
  D_{ij}=\sqrt{{D_{ij}^{(x)}}^2+{D_{ij}^{(y)}}^2+{D_{ij}^{(z)}}^2}.
\end{equation}
and this is always a non-negative quantity.
The distance matrix $\bf{D}$ then is used to calculate the adjacency
matrix $\bf{A}$ which is defined as:
\begin{equation}
\label{eq:adjacency_matrix}
A_{ij} = \left\{
\begin{array}{rl}
 & 1,  \text{if }  D_{ij} \le d_{\rm thr}\\
& 0,   \text{otherwise }
\end{array} \right.
\end{equation}
 

The adjacency matrix is usually introduced in the context of graph
theory \citep{book_algorithms}, as a convenient way of uniquely representing a
graph.
Monomers in a cluster can be formally regarded as the graph  vertices   and the
physical bounds due to the interaction potential are, in this analogy,
the  graph  edges.  
The problem of determining the clusters, given the distance matrix
${\bf D}$, is then re-formulated equivalently as the determination 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{Calculation of the Radius of Gyration and the Hydrodynamic
Radius}
The radius of gyration is a fundamental quantity of interest; indeed
the dependence of the number of monomers in a cluster on its radius of
gyration can be used to determine the fractal dimension of the system.
The radius of gyration is defined as the root-mean-square distance
mass displacement around the aggregate   centre-of-mass 


\begin{equation}
  \label{eq:r_gyr_definition}
R_g^2=\f{\s_i{({\bf r}_i}-{\bf r}_{CM})^2}{k}+a_1^2,   
\end{equation}
where $k$ is the number of monomers in the cluster and $a_1$
represents the radius of gyration of a monomer, which we take as
$a_1=\sqrt{3/5}(\sigma/2)$ as in \cite{sorensen_prefactor}.
 We stress that the cluster radius of gyration is a static property
 since it does not depend on the diffusional properties of the
 cluster, but only on its mass distribution.


In order to evaluate $R_g$ according to Eq. \ref{eq:r_gyr_definition}
we need again to take into account the periodicity of the box.
We adopt a reference frame centred on the position of the first
monomer in the cluster, i.e. we work in a frame where ${\bf
  r}_1=(0,0,0)$. 
 Since all the
monomer-monomer distances along each axis are known from
Eq. \ref{eq:distance_calc}, one can calculate the position of
all the other monomers in a cluster with respect to the first monomer;
for instance the $j$-th monomer position along $x$ will be given by:
\begin{equation}
  \label{eq:reconstruct_cluster}
  x_j=x_1+D_{1j}^{(x)}=D_{1j}^{(x)}
\end{equation}
and similarly along the $y$ and $z$ axis.
Once all the new monomer coordinates have been obtained from Eq \ref{eq:reconstruct_cluster}, one can
safely calculate the coordinates of the 
aggregate centre-of-mass  (e.g. $x_{CM}=\s_ix_i/k$) and apply straightforwardly
Eq. \ref{eq:r_gyr_definition}.
We point out 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 on
the cluster spatial position) nor does it depend on the
sign convention chosen for $D_{1j}^{(x)}$.

Equation \ref{eq:distance_3D} allows one  to also calculate the
hydrodynamic radius of a cluster. In fact the cluster hydrodynamic radius, $R_H$, is defined as \citep{lattuada_hydro} 
\begin{equation}
  \label{eq:hydrodynamic-radius}
\f{1}{R_H}=\f{1}{k(k-1)}\s_{i\neq j}\f{1}{D_{ij}}. 
\end{equation}
The hydrodynamic radius again is another static property of the
cluster and his relevance is given by the fact that it is often used
to characterise the transport properties of a cluster \citep{lattuada_hydro}.
\subsection{Collisional Kernel}
\esp allows one to address
individually each monomer at all times during the simulation, i.e. one
can  ``stick''  a permanent label onto each monomer.
One can then regard each cluster as a not ordered collection of
monomers where no monomer label is repeated. This amounts to the
mathematical definition of a set.
Viewing clusters as sets of monomers provides a computational method
to investigate the process of cluster collisions.
First of all, we should remind ourselves that both the model and Van der Waals
potential used in the simulations are much deeper than $k_BT$ hence we
can safely assume a sticking probability upon collision
equal to one.
To fix the ideas, let us consider two individual clusters  at time
$t$.
If in the time interval $(t, t+\delta t]$ they collide, they will end up in
the same newly-formed cluster observable at time $t+\delta t$.
An equivalent description of the collision process is that  both
sets of monomers at time $t$ (the two initial clusters) end
up being proper subsets of the same set of monomers (the cluster
formed by their collision) detectable  at time
$t+\delta t$. One can then record the number of collisions taking
place in $(t, t+\delta t]$ and the clusters involved in the collisions
simply by comparing different sets.
The collisional kernel between $i$ and $j$-mers is estimated from the
relation
  \begin{equation}  
    \label{eq:kernel_calc}
    \f{N_{ij}(t)}{\delta t}=(2-\delta_{ij})\beta_{ij}\f{N_i(t)}{V_{\rm
        box}}\f{N_j(t)-\delta_{ij}}{V_{\rm box}}\cong (2-\delta_{ij})\beta_{ij}n_i(t)n_j(t),
  \end{equation}
where $N_{ij}(t)$ is the number of collisions between $i$ and $j$-mers
taking place in 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}$ and $\delta_{ij}$
is Kroenecker delta.

From Eq. \ref{eq:kernel_calc}, we choose the optimal value of  $\beta_{ij}$ by
treating it as a fit parameter minimising (in a least-square sense)
the distance between the list of detected collisions $N_{ij}(t)/\delta_t$ and
$(2-\delta_{ij})\beta_{ij}n_i(t)n_j(t)$.

%  A collision is detected whenever two independent
% clusters at time $t$ become a proper subset of the same cluster at a
% later time $t+\delta t$, hence the detection of cluster collisions is
% reduced to comparing different sets of monomers at time $t$ and
% $t+\delta t$ along the evolution of the system. 
 



\section{Results and discussion}\label{results}


\subsection{Cluster diffusion coefficient and cluster thermalization}
\label{sec:clust-diff-coeff}
The complete motion of the cluster in 3D is given by a coupled
rotation and translation. In this section we investigate the
translational motion of the cluster centre-of-mass while leaving the
investigation of its rotational properties to the next section.

Equation \ref{eq:Langevin} determines the dynamics of each monomer in
the system and indirectly also the cluster diffusional properties.
The numerical determination of the diffusion coefficient of a $k$-mer
(i.e. a cluster consisting of
 $k$ monomers)  and its scaling with $R_g$ and the number of monomers $k$ can provide valuable
insight into the dynamics of the system. Furthermore, the knowledge of
a $k$-mer diffusion coefficient $D_k$ 
is useful to investigate the regime where the system finds itself.
In fact, one typically determines whether the system is in the
continuum or free-molecular regime (or in the transition between the
two) by estimating the \emph{gas} Knudsen number defined as
$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$).   
Langevin equation \ref{eq:Langevin} characterises the fluid only via
its viscosity which is modelled by the damping and noise terms and no
information about the fluid molecule mean-free path is available.


We can indirectly extract  information about the Knudsen number
from the numerical evaluation of the diffusion coefficient.
From Einstein relation, the diffusion coefficient of a single (spherical)
monomer is \citep{risken_book}
\begin{equation}
  \label{eq:monomer-diffusion-coefficient}
  D_1=\f{k_BT}{m_1\beta_1}.
\end{equation}
Equation \ref{eq:monomer-diffusion-coefficient} can be extended for a
$k$-mer  by introducing Cunningham
slip factor $C_s(Kn)$ and the inverse of the $k$-mer relaxation time
$\beta_k$
\begin{equation}
  \label{eq:cluster-diffusion-coefficient}
  D_k=\f{k_BT}{km_1\beta_k}C_s(Kn).
\end{equation}
For a system in the continuum regime ($Kn\to 0$) without the screening
effect of  external  monomers on  internal monomers in a cluster
(which is neglected in 
Eq. \ref{eq:Langevin}), $D_k$ is given by
Eq. \ref{eq:cluster-diffusion-coefficient} with $C_s(Kn)=1$ and $\beta_k=\beta_1$.
This can be investigated numerically by evaluating the diffusion
coefficient of a cluster from the time dependence of its
centre-of-mass mean-square displacement \citep{risken_book} 
%  We extend to a cluster the same relationship
% between the mean square displacement and the diffusion coefficient of
% a Brownian monomer 
\begin{equation}
  \label{eq:diffusion_simu_cluster}
  \langle \delta{ r}^2_{CM}(t) \rangle\xrightarrow{t\to\infty} 6D_{k}t,
\end{equation}
where $\langle
\delta{ r}^2_{CM}(t) \rangle$ is the ensemble average of the
cluster
 mean square displacement (variance of the  cluster center-of-mass
 positions at a given time).


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$.
Aggregate transport is investigated by placing the cluster  in the
middle of a box having size $L=10000$. The box size has to be much
larger than for the agglomeration simulations since we now investigate
the aggregate diffusion and 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 a large box to track
the aggregate center of mass for a (potentially) long time.


The results of the numerical simulations are summarised in Table
\ref{table_diffusion}  where we report the diffusion
coefficient calculated according to
Eq. \ref{eq:diffusion_simu_cluster} and the radius of gyration $R_g$ of each cluster used in the
simulations.
We notice that $D_k\propto 1/k$ for $k$ spanning almost two orders of
magnitude, thus confirming that we are in the continuum regime.
Furthermore, the estimated value of $\beta_k$ is equal to the inverse
of the monomer relaxation time  $\beta_1$ within a few percents, the
non-perfect agreement being   
due to the relatively limited number of simulated stochastic trajectories.






Fig. \ref{diffusion_plot} shows the simulation results for $\langle
\delta{ r}^2_{CM}(t) \rangle$ for  a $50$-monomer cluster,
namely the one on the bottom right  in
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$ the cluster finds
itself in the ballistic regime 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 highly reminiscent of the behaviour of a single Brownian
monomer with an initial zero velocity \citep{risken_book}.
A perfectly similar behaviour is found when investigating the
diffusional properties of a cluster with $k=4$, thus ruling out the
possibility that only large clusters are in the continuum regime.
In the continuum regime, one can express the diffusion coefficient
also in terms of a mobility radius $R_m$ defined through the equation
\begin{equation}
  \label{eq:mobility_radius}
  D_{k}=\f{k_BT}{6\pi\mu_f R_m},
\end{equation}


The fluid viscosity can in its turn be expressed in terms of the monomer properties
\begin{equation}
  \label{eq:fluid_viscosity}
  \mu_f=\f{m_1\beta_1}{3\pi\sigma}.
\end{equation}
thus  allowing us to calculate the mobility radius $R_m$
by inverting Eq. \ref{eq:mobility_radius}.



% The diffusion coefficient of a cluster in the continuum limit can be
% expressed through Stokes relations 
% % Langevin equation \ref{eq:Langevin} determines   contains explicitely the monomer
% % relaxation time $\tau_{{\rm mon}}=1/\beta$, which is related to the
% % monomer diffusion coefficient $D_{\rm mon}$ by the relation
% \begin{equation}
%   \label{eq:monomer_diffusion}
%   D_{\rm clu}=\f{k_BT}{6\pi\mu_f R_m},
% \end{equation}
% where $R_m$ is the mobility cluster radius.
% In general, the mobility radius is expected to be of the same order of
% magnitude or larger than the radius of gyration or the hydrodynamic
% radius, the
% latter defined as:
% \begin{equation}
%   \label{eq:hydrodynamic_radius}
%  \f{1}{ R_H}=\f{1}{N_c(N_c-1)}\s_{i\neq j}\f{ 1}{r_{ij}}
% \end{equation}


The cluster radius of gyration $R_g$ or its hydrodynamic radius $R_H$
are often taken to be roughly equal to the aggregate mobility radius
$R_m$ and used to calculate its diffusion coefficient.
Our simulations (cfr. Table \ref{table_diffusion}) show that the hydrodynamic radius and the radius of
gyration tend to be of the same order of magnitude for cluster sizes in the range
 $k=4$ to $k=98$, but that they both severely underestimate the
 mobility radius. 
The major approximation we are making when studying cluster mobility
is that we use  Eq. \ref{eq:Langevin} to describe the dynamics of both an
isolated monomer and  a monomer in a cluster, which is very often 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 the velocity fluctuations in order
to gain some insight  into the cluster
thermalization process.
For a single Brownian monomer,  the velocity fluctuations  $\langle \delta v^2 \rangle$ are given by \citep{risken_book}
\begin{equation}
  \label{eq:delta_v_squared}
  \langle \delta v^2 \rangle =\f{k_BT}{km_1}.
\end{equation}
with $k=1$.
Having determined that $\beta_k=\beta_1$ and that $C_s(kn)= 1$,
we have already tentatively  trivially generalised the well-known
equation for a monomer (see \cite{risken_book}) for
a  $k$-kmer.
In Figure \ref{diffusion_plot} (right) we plot  $\langle \delta v^2
\rangle$ as a function of time again for the same cluster considered
in the diagram on the left. We notice that $\langle \delta v^2 \rangle$   
clearly fluctuates about $0.03$, the value predicted by
Eq. \ref{eq:delta_v_squared} for $k_bT=0.5$, $k=50$ and $m_1=1$.

The main result of this section is that  from Eq.
\ref{eq:Langevin}, i.e. neglecting monomer screening, a $k$-mer  is
found to have the same relaxation time of a monomer
and the same diffusional and thermal properties of a sphere  with a
mass $k$ times the one of a single monomer.



% Other authors like  \cite{mountain}  start by 
% writing Langevin equation directly for a cluster containing
% $N$ monomers whereas \cite{meakin_cluster_models} discusses the
% effects of different assumptions on cluster mobility on their dynamical
% behaviour.
% In the approach developed in this study, the dynamics of a cluster
% arises naturally from monomer-monomer interaction. 
% The main implicit
% approximation in this treatment is that the noise term and particle
% relaxation time in Eq. \ref{eq:Langevin} are clearly different for a
% single monomer and one shielded by other monomers in cluster.



% The interaction potential we use aim at mimicking the condition of
% sticking upon collision between two monomers. It therefore also
% defines the monomer diameter and has to consist of a hard-core
% repulsion and a very deep binding part.


% The analytical expression of the pair-wise radial potential $u(r)$ is
% given in the appendix. Here suffices to say it is a deep smooth  
% potential well. The interaction potential defines the monomer diameter
% $\sigma$ which is the separation where two monomers start experiencing
% a strong constant repulsive force. The potential has a deep minimum
% located at $r_{\rm min}>\sigma$. The potential well is tight ( $r_{\rm
% min} -\sigma\ll \sigma $) and much deeper than $k_BT$, thus leading to
% a capture probability upon collision which can be safely taken as one.

% The numerical simulation are performed using the code \esp
% \citep{espresso} using the Langevin thermostat and the
% user-specified potential described in detail in Appendix \ref{explain_potential}.


%\section{Model results}
% The size of the computational domain (the box containing the monomers)
% is fixed by the density initial particle density $\rho$ according to:


\subsection{Cluster rotation}
Cluster rotational motion is investigated using the results of the
simulations described in the previous section. 
Each cluster can be treated as a rigid body
 and any rotation of a rigid body in 3D can be
described resorting to three angles $\theta$, $\phi$ and $\psi$, usually
called Euler angles \citep{goldstein}.
In the post-processing, we eliminate the contribution of the translational motion by rigidly
  translating the aggregate centre-of-mass to the origin of the
  computational box coordinate system at all times.   
 At this point one can calculate the Euler angles by keeping track of
 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}=[{\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'}=[{\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=X'X}^T\lro {\bf XX}^T \rro^{-1},
\end{equation}
where ${\bf X}^T$ is the transpose of $\bf X$. 
Once the rotation matrix $\bf A$ is known,  Euler angles $\theta$,
$\phi$ and $\psi$ can be determined by the algorithm in \cite{rotation_alghoritm}. 
According to the numerical algorithm employed, $\phi$ and $\psi$ both
range in the interval $[-\pi,\pi]$, whereas $\theta$ ranges in
$[-\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 draws from a uniform probability distribution
in the interval $[-\pi,\pi]$, hence $\langle\phi\rangle=\langle\psi\rangle=0$ and 
$\langle\delta\phi\rangle=\langle\delta\psi\rangle={\pi}/{\sqrt
  3}\simeq 1.81$. 
On the other hand,  the Euler angle $\theta$ distribution
can be expressed as \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 draw from a uniform probability
distribution in the interval $[0,1]$.
We numerically investigated the distribution of $\theta$ angles given
by Eq. \ref{eq:theta_distr} finding that $\langle\theta\rangle=0$ and
$\langle\delta\theta\rangle\simeq 0.68$. 

% For a perfectly random orientation of the clusters, i.e. treating each
% angle as a uniform random variable in its interval of variation, one
% calculates $\langle\phi\rangle=\langle\psi\rangle=0$ and 
% $\langle\delta\phi\rangle=\langle\delta\psi\rangle={\pi}/{\sqrt
%   3}\simeq 1.81$.
%For $\theta$.  random_euler
% These expected results are confirmed by the numerical investigation whose results
% are shown in Fig. \emph{add figure!!!}, where $\langle\phi\rangle$,
% $\langle\psi\rangle$ and $\langle\theta\rangle$ oscillates around zero
% and, at late times, their standard deviation saturates around the
% values calculated above. This clearly shows that, as a consequence of
% the absence of screening in Eq. \ref{eq:Langevin}, each cluster has no
% preferential orientation nor preferential axis of rotation.

Figure \ref{cluster_rotation} (left) shows the results for the mean
values of Euler angles calculated on an ensemble of $800$ identical
clusters consisting of $10$ monomers each. We notice 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 fluctuation of Euler angles for the same ensemble of
clusters and we see that they quickly saturate 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.



% Figure \emph{reference figure} (left) shows the time evolution of the
% ensemble-averaged Euler angles $\langle\theta \rangle$ , $\langle\phi \rangle$,
% $\langle\psi \rangle$ calculated with respect to the initial cluster orientation.
% One notices that their average value at late times fluctuates around
% zero thus indicating that on average there is no preferential orientation of the clusters. 
%  We can actually conclude that there is no preferential rotation axis 
%  for the clusters by looking at the time-dependent values of $\langle\delta\theta \rangle$ , $\langle\delta\phi \rangle$,
% $\langle\delta\psi \rangle$.
% Bearing in mind that both $\phi$ and $\psi$ vary in the interval
% $[0,2\pi]$, then for a uniform random angle distribution of $\phi$ and
% $\psi$ one obtains
% \begin{equation}
%   \label{eq:delta_psi}
%   \langle
%   \delta\psi\rangle=\f{1}{2\pi}\sqrt{\lro\int_0^{2\pi}\lro\psi-\pi\rro^2
%     \rro}=\f{\pi}{\sqrt 3}\simeq 1.81,
% \end{equation}
% and the same result applies of course to
% $\langle\delta\phi\rangle$. For the calculation of
% $\langle\delta\theta\rangle$ one has to consider that, beside ranging
% from $0$ to $\pi$, $\theta$ is invariant for reflection of an
% aggregate with respect to the $xz$ plane. Hence, its variance is
% reduced by a factor two and has to be calculated as

% \begin{equation}
%   \label{eq:delta_theta}
%   \langle
%   \delta\theta\rangle=\f{1}{2\pi}\sqrt{\f{1}{2}\lro\int_0^{\pi}\lro\theta-\f{\pi}{2}\rro^2
%     \rro}=\f{\pi}{2\sqrt 6}\simeq 0.64.
% \end{equation}
% In Fig. \emph{reference figure} we notice that
% $\langle\delta\phi(\psi)\rangle$ and
% $\langle\delta\theta\rangle$ saturate to the asymptotic values
% predicted by Eqs. \ref{eq:delta_psi} and \ref{eq:delta_theta},
% respectively, thus confirming the idea that 

 


\subsection{Determination of the fractal dimension}\label{calculate_df}


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 the power-law:
\begin{equation}
  \label{eq:determination_df}
R_g= \alpha k^{1/d_f},  
\end{equation}
 where $d_f$ is the fractal dimension of the system and $\alpha$ is a
 constant.
% The radius of gyration is a geometric property of the cluster which
% does not depend on the way it diffuses through the medium.


 The average time-independent fractal dimension can be determined by fitting the
dependency of all the recorded cluster sizes on their radia of
gyration, as shown in   Fig. \ref{two_df} 
where we plot the calculated mean $R_g$ as a function of the
aggregate size $k$. 
There we require an aggregate size $k$ to have been recorded
at least $200$ times before the corresponding average $R_g$ is
calculated.
This is not a very stringent requirement if we bear in mind that we
store  $1500$ system configurations for each simulations and that the
results are the average on $10$ different simulations, but it aims at
getting rid of outliers when calculating the mean radius of gyration.
We also point out that our results are not sensitive on the threshold
on the minimum number of occurrences of a cluster size before the
cluster can be included in the statistics. Raising the occurrence threshold from $200$ to
$400$ does not change the calculated time-averaged fractal dimensions.
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 fractal dimension of the system is an
output of the simulations and not a tunable input parameter.

 The double
logarithmic plot shows that the power-law scaling breaks down at
cluster sizes $k<5$, thus setting $k=5$ as the minimum cluster size
included in the fits in Fig. \ref{two_df}.
%  The minimum cluster size
% included in the power-law scaling in Eq. \ref{eq:determination_df} is
% $k=5$, though we report all the calculated radia of gyration down to
% $k=1$.
A fit of the collected
data with a single line, leading to $d_f=1.62$, clearly does not
capture the  clusters  with $k\le 15$.
 Hence the results are better
fitted with two different slopes: we identify small ( $k\le 15$)  and
large ($k>15$) clusters and we calculate a fractal dimension for small clusters,
$d_f^{\rm small}=2.25$ and one for large clusters $d_f^{\rm
  large}=1.56$.
% A careful statistical analysis of the data leads to the conclusion
% that the dependency of the radius of gyration on the cluster size is
% better described
% two rather than one fractal dimension.
The result  $d_f^{\rm small}>d_f^{\rm large}$ can be understood by
considering that small clusters, generated mainly by monomer-cluster
agglomeration (also called diffusion-limited agglomeration, DLA), 
are more compact and spherical than large clusters, which
are mainly brought about by cluster-cluster agglomeration.
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, meakin_trajectories2}.
% We should point out that we believe that the true fractal dimension of
% the clusters generated in these simulations is actually slightly lower
% than this value 
%  but here the
% radial monomer-monomer interaction could introduce multiple contact
% points between the monomers \emph{MAKE  THIS MORE EXPLICIT}.



It is interesting to evaluate a time-dependent ensemble-averaged  fractal
dimension of the system, $d_f^t$, by using the power-law in
Eq. \ref{eq:determination_df} for each ensemble of recorded
configuration at time $t$ and
without taking any average.  The evolution of the time-dependent
fractal dimension is shown in Fig.
\ref{evolution_df} (top).
One can link the evolution of the fractal dimension to the kinetics of
the small and large cluster populations.
At early times, the system is almost entirely made up of small
clusters, hence $d_f^t\simeq d_f^{\rm small}$, whereas at late times
the contribution of large clusters to the overall cluster population
is dominant, hence $d_f^t\to d_f^{\rm large}$.



These findings are in the spirit of the study by
\cite{konstandopoulos}, where a distribution of fractal
dimensions was employed to characterise the aggregate morphology and the
mean fractal dimension was shown to relax to an asymptotic value,
fixed a priori by the agglomeration mechanism. 
We remark again that in this study no assumption is made a priori
about the asymptotic fractal dimension and we show explicitly the importance of the
kinematic 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, \cite{smirnov_review, mountain} show that the fractal dimension is mainly
affected by the aggregation mechanism only. Indeed, we point out that
the corresponding fractal dimensions calculated for the simulations
with the model potential do not differ substantially from the ones
presented for Van der Waals interaction. This seems in contrast with
the findings in \cite{videcoq}, where the fractal dimension is found
to be a function of the interaction range. Here however, $u_{\rm VdW}$
is a quickly decaying function of $r$ and does not enhance
significantly the cluster capability of re-organisation, which leads to
an increase of the calculated fractal dimension in the systems studied
in \cite{videcoq}.

% \subsection{Statistical error on the fractal dimension and the prefactor}
% To be written. Necessary to agree about the formulas used and whether
% we think we understand the meaning of the statistical uncertainties
% (the standard errors are tiny).


\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 each
monomer in the cluster.
The coordination number is calculated in the post-processing of the
simulations again resorting  the igraph library
\citep{igraph}.
In the following, we call mean coordination number the time-dependent
ensemble-averaged coordination number.
This quantity can provide valuable information about the openness of
the aggregates and the presence of cavities in their structure. 
The knowledge of the mean coordination number complements
the information from the fractal dimension about the aggregate
structure.


As one can see from Fig. \ref{evolution_coord_number}
  the mean coordination number goes well above seven.
This relatively high value provides new information about the cluster
morphology.
In fact we argue that the fractal dimension, as determined
from the scaling law in Eq. \ref{eq:r_gyr_definition} is not able
to fully characterise aggregate morphology by itself.
Indeed, the large clusters generated at late times in the simulation
exhibit low fractal dimension $d_f^{\rm large}<2$ due to their
often tubular, elongated shape (exemplified by the second
and fourth cluster counting clockwise in
Fig. \ref{snapshot_end_simulation}) and not because of holes or
cavities in their structure.
This conclusion (valid in a statistical sense) cannot be reached on the basis of the fractal
dimension alone, hence we advocate the use of the fractal dimension
and the coordination number for a better characterisation of aggregate
morphology.
The generation of clusters with a high coordination number (as most
likely those in Ref. \cite{videcoq}) is an intrinsecal property of
Langevin simulation  due to their
neglecting the interparticle hydrodynamic interaction, as demonstrated
in Ref.  \cite{prl_hydro}.

   

% \subsection{Distribution of cluster morphologies}
% \label{sec:distr-clust-morph}
% Most of the cluster generated in our simulations are non-spherical
% objects, as one notices from the system snapshot and the indivual
% clusters shown in Fig. \ref{snapshot_end_simulation}. 
% This result is expected and is a consequence of the short-ranged
% potential (which aims at reproducing an infinitely  narrow potential
% well).



\subsection{Cluster Restructuring}
\label{sec:clust-restructuring}
Cluster restructuring was investigated in the framework of 2 and 3D
lattice models in Refs.
\cite{meakin_reorganization_2D,meakin_reorganization_3D}.
The novelty of the approach that we propose here is that cluster
restructuring does not have to be implemented as an additional feature
of the cluster dynamics. The monomers interact via a deep interaction  \emph{radial}
potential. The effect of the interaction potential is
mainly to lock the relative 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 given by
the collision of the molecules of the surrounding fluid with the
monomers, which is 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  indicator of the modification of
the structure of a cluster, namely the radius of gyration and the mean
monomer-monomer distance.
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 $\l D_{ij}\r$ during a
time interval when the cluster increases its size from $k=14$ to $k=43$
due to a collision taking place at time  $t=568$.
We notice that most of the time before and after  this collisions, the cluster
$\l D_{ij}\r$ is practically constant, a clear sign of its 
rigidity. When the collision between the two aggregates occurs, $\l
D_{ij}\r$ increases, but  on a time scale of a few
$\tau_1$ it settles to another value we can consider again as a
constant.
We see this as a strong indication of reorganization of the
cluster resulting from the collision of two clusters.
A perfectly analogue behavior is observed when looking at the radius
of gyration of the cluster on the same time interval (see the right
plot in Fig. \ref{single_cluster_size_evolution}).

   
% However, shortly  after the collisions at times $t\simeq 1800$ and
% $t\simeq 2700$,  $R_g$ first increases and then quickly decreases
% before reaching a constant value remaining. This happens when the cluster we track collides with
% another one and increases its size by about $30$ monomers in each
% case. 
% The same conclusions can be drawn from studying the mean monomer
% separation (called wrongly mean particle separation) between the
% monomers in the  cluster we track. In
% Fig. \ref{single_cluster_restructuring_mean_distance} we see evidence
% of cluster restructuring this time by the reduction of the mean
% monomer distance after the same two collision events discussed above.









\subsection{Evolution of the number of clusters}

The total  cluster concentration in the system, $N_{\infty}$, is one of
the most important quantities of interest since it shows the
progress of agglomeration and 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}\beta_{ij}n_in_j-n_k\sum_i\beta_{ik}n_i,
\end{equation}
where $n_k$ is the concentration of clusters consisting of $k$
monomers and $\beta_{ij}$ is the collisional kernel between
 $i$ and $j$-mers, respectively.

% Fig. \ref{cluster_decay}  shows the decay of the number of clusters as
% time progresses on a double logarithmic scale. One notices that, after an
% initial transient, $N_\infty(\beta t)$ there are two distinguidhed
% physical regimes. At intermediate times, $t$ in  $[100, 300]$, we find
% $N_{\infty}\sim t^{-1}$, but the late behavior, $t$ in $[1000, 3000]$
% yields  $N_\infty\sim t^{-0.85}$, which is rather
% different  from the $t^{-1}$ behaviour expected for the continuum regime
% (\citep{mountain, smirnov_review, pratsinis_self_preserving}).
% A similar behaviour was found in Ref. \cite{mountain}, where the
% extent of coagulation was studied by g
% The numerical solution of Smoluchowski equation with the continuum
% kernel and a fractal dimension $d_f=1.7$
% (cfr. Eq. \ref{eq:cont_kernel}) is also plotted on the same
% figure. There, the long-time behavior of the total concentration is
% $N_{\infty}\sim t^{-1}$, as expected.
% It is obvious that the two calculations differ quite a lot when one
% takes into 
% On the same figure, we also show the evolution of the number of
% clusters containing more and fewer or equal than $15$ monomers, called
% $N_>$ and $N_<$, respectively. As we will see in the following, small
% and large clusters tend to have a different fractal dimension.
% A power-law decay of the total cluster concentration $N_\infty\sim
% t^{-\alpha}$ together with the fractal law \ref{eq:determination_df}
% leads to the expected time-dependence of the mean radius of gyration
% \begin{equation}
%   \label{eq:evolution_R_gyr}
%   R_g\sim t^{\alpha/d_f}
% \end{equation}
% which is plotted in Fig. \ref{mean_R_gyr_evolution}. The estimated
% exponent of the power-law increase of the radius of gyration $R_g\sim
% t^{0.49}$ is rather close to the ratio $\alpha/d_f\simeq 0.53$.




The standard treatment of the collisional kernel in the continuum
limit (hereafter called continuum kernel),   is given by \citep{friedlander book}
\begin{equation}
  \label{eq:cont_kernel_general}
  \beta_{ij}=4\pi (D_i+D_j)(R_{i}+R_{j}),  
\end{equation}
where $R_i=R_{g,i}$ is the radius of gyration of an $i$-mer 
(we drop the pedix $g$  in order not to burden the notation). 
 Equation \ref{eq:cont_kernel_general} is often expressed in terms of
 the aggregate volume of solid  $\nu_i\propto(R_{i})^{d_f}$
 ($d_f$ being a single, time-independent fractal dimension determined
 in Section \ref{calculate_df} ) while the
 diffusion coefficient is given by
 Eq. \ref{eq:mobility_radius} with the radius of gyration taken as the
 mobility radius. The continuum kernel then reads
\begin{equation}
  \label{eq:cont_kernel_standard}
  \beta_{ij}^{Sm}=\f{2k_BT}{3\mu_f} \lro
  \nu_i^{1/d_f}+\nu_j^{1/d_f} \rro \lro \f{1}{\nu_i^{1/d_f}}+\f{1}{\nu_j^{1/d_f}}  \rro.
\end{equation}
% We already know that Eq. \ref{eq:cont_kernel_standard} is not the most
% appropriate kernel to use in the agglomeration equation since, above
% all for large aggregate size, the approximation $R_g=R_m$ is rather
% poor (see Section \ref{sec:clust-diff-coeff}).

% In Figure \ref{small_and_large_clusters_and_n_infinity}, we plot the results of the discrete simulations
% together with different approaches based on the agglomeration Equation
% \ref{eq:agglomeration_equation}. 
% First of all, we notice that the standard treatment relying on the
% kernel given by Eq. \ref{eq:cont_kernel_standard} does not reproduce
% the system dynamics even at short times. This is not surprising if we
% take into account the fact that the approximation $R_m\simeq R_g$ is
% not justified for the clusters generated in our simulations,
% especially when $k$ is large (cfr. Table \ref{table_diffusion}). 



On the other hand, we can use the information obtained from the numerical simulations to
derive a kernel which describes more accurately the aggregate
collisional dynamics starting from Eq. \ref{eq:cont_kernel_general}.
In order to do this, we use the expression for the diffusion
coefficient as calculated from
Eq. \ref{eq:cluster-diffusion-coefficient} with $\beta_k=\beta_1$ and
$C(kn)=1$. 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 = \left\{
\begin{array}{rl}
&  \text{data from simulations } i< 5    \\
 & \alpha_{\rm small}i^{1/d_f^{\rm small}},  \text{if }  5\le i\le 15 \\
& \alpha_{\rm large}i^{1/d_f^{\rm large}},  \text{if }  i> 15 
  
\end{array} \right.
\end{equation}



Furthermore, non-continuum effects can play a role in cluster collisions and
they are usually accounted for  by   introducing Fuchs
correction \citep{fuchs book} $\beta_f$ in the kernel.
% Furthermore, the presence of two different fractal dimensions in the
% system should be mirrored in the kernel by determining the average
% radius of gyration $R_g$ of a cluster
% from Eq. \ref{eq:determination_df} using $\{\alpha_{\rm
%   small}, d_f^{\rm small} \}$ and $\{\alpha_{\rm large}, d_f^{\rm
%   large} \}$ for small and large clusters, respectively. 
 Hence the kernel derived for Langevin simulations (hereafter called
 Langevin kernel) reads
\begin{equation}
  \label{eq:kernel_complete}
  \beta_{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}
% where 
% \begin{equation}
%   \label{eq:rules_for_r_gyr}
%   R_i = \left\{
% \begin{array}{rl}
%  & \alpha_{\rm small}i^{1/d_f^{\rm small}},  \text{if }  i\le 15 \\
% & \alpha_{\rm large}i^{1/d_f^{\rm large}},  \text{if }  i> 15 
  
% \end{array} \right.
% \end{equation}
It is useful at this point to introduce the kernel homogeneity
exponent. A kernel $\beta_{ij}$ is said to be homogeneous of
order $\lambda$ if $\beta_{\gamma i,\gamma
  j}=\gamma^{\lambda}\beta_{ij}$.
The importance of the kernel homogeneity exponent is that it fixes the
asymptotic time-decay of the total number concentration $N_\infty\sim
t^{-1/(1-\lambda)}$.
For Smoluchowski and Langevin kernel one obtains $\lambda=0$,
$N_\infty\sim t^{-1}$ and $\lambda=(1/d_f^{\rm large})-1$,
$N_{\infty}\sim t^{-0.72}$, respectively.



The results of the simulations (with Van der Waals and model potential) and the numerical
solution of the agglomeration equation \ref{eq:agglomeration_equation}
with $\beta^{Sm}$ and $\beta^{LD}$ are shown in Fig. \ref{Smoluchowsky
  comparison}. Langevin kernel in 
Eq. \ref{eq:kernel_complete} leads to an enhanced agreement at
short times between the simulations and the numerical solution of the
agglomeration equation \ref{eq:agglomeration_equation}.
We point out that the asymptotic limit in the agglomeration equation
is reached at times about one order of magnitude longer than the
duration of the simulations, mainly due to Fuchs factor.

 Nevertheless, we fit the  time-decay of the total cluster
 concentration to a power-law, $N_\infty\sim t^{-\xi}$, at late
 times
$ 2500 \le t\le 3000$.
  From
  Eq. \ref{eq:agglomeration_equation} with  Langevin kernel, we find $\xi=0.77$,  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$. 
Finally, we notice that, as expected, 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.


% Furthermore, since the initial state correspond to a monodisperse
% aerosol, one can try to describe the evolution of $N_\infty(t)$ by
% means of a constant kernel approximation \citep{friedlander book}
% \begin{equation}
%   \label{eq:monodisperse_approx}
%   N_\infty(t)=\f{N_\infty(0)}{1+\mathcal{K}_0t/2},
% \end{equation}
% where $\mathcal{K}_0=2k_BT/(3\mu)$ is the value of the collisional
% kernel for two equally-sized clusters.






%  and
% $\mu$ is the viscosity of the gas.
% The gas viscosity is in turn related to the friction coefficient $\beta$ in
% Eq.~\ref{eq:Langevin} through Stokes relation $\beta=3\pi\mu \sigma/m$.
% One can easily compare the decay of the number of clusters as a
% function of time with a simple approximation for monodisperse aerosols
% which holds at early times \cite{friedlander book}:
% \begin{equation}
%   \label{eq:monodisperse_approx}
%   N_\infty(t)=\f{N_\infty(0)}{1+\mathcal{K}_0t/2},
% \end{equation}
% where $\mathcal{K}_0=2k_BT/(3\mu)$ is the value of the collisional
% kernel for two equally-sized particles.
% We notice that Eq. \ref{eq:monodisperse_approx} has no information
% about the fractal dimension, as a direct consequence of the weak
% dependence of the collisional kernel on $D_f$ in the continuum limit.
% Fig. shows the good agreement between the numerical results and the
% monodisperse approx in Eq. \ref{eq:monodisperse_approx} at early
% times, as expected. 




% One can plot the radius of gyration of the clusters versus the cluster
% sizes. In general, clusters constituted of the same number of monomers
% will have a different radius of gyration, so we want to investigate
% the dependence of the mean radius $\overline R_g$ as a
% function of the cluster size, as in Mountain.
% Calculating naively $\overline R_g$ using the data collected at all
% times during our simulations would overweight the outliers
% (i.e. specific cluster sizes which happen to appear  too few times
% to allow for reliable statistics). In order to avoid that, we require
% a each cluster size to appear at least $200$ times in the collected
% data in order to consider it for the calculation of $\overline R_g$.
% This may seem a lot, but we should bear in mind that we are consider
% the data collected from $10$ different simulations, each starting with
% $5000$ monomers and for each simulation we examine the state of system
% (calculating the number of clusters and their properties) $150$ times
% during the time evolution.
% Furthermore, we find substantially the same results increasing this
% minimum recording  threshold from $200$ to $500$.  
% All this is illustrated in Fig. \ref{cluster_decay}.
% A single fit of $\overline R_g\sim N^{1/d_f}$ spanning all cluster
% sizes would lead to the estimate $d_f\sim 1.7$, which is a value often
% reported in literature \citep{mountain, smirnov_review, meakin_trajectories2}) for
% cluster-cluster aggregation in three dimensions.
% In particular, Mountain included all clusters sizes in his fit, from
% dimers to the largest generated clusters.
% If, on the other hand, we identify small clusters with those having a
% size below $15$ monomers and we consider large the remaining clusters,
% then we can distinguish  two different slopes in the curve in
% Fig. \ref{cluster_decay}. 
% This yields a fractal dimension for small clusters $d_f^{\rm
%   small}=2.2$ and a corresponding value for large clusters $d_f^{\rm large}=1.6$.
% The values of the fractal dimensions for small and large clusters
% calculated above should be meant as time-averaged quantities. By
% looking at the number of small and large clusters in
% Fig. \ref{cluster_decay} one notices that at early times we have
% $N_\infty\simeq N_<$ whereas at late times, after coagulation has
% reduced the initial number concentration by almost two orders of
% magnitude, $N_\infty\simeq N_>$.
% Therefore we expect the time-dependent fractal dimension (calculated
% on all the cluster sizes) to be close to $d_f^{\rm
%   small}$ at early times and to evolve towards  $d_f^{\rm
%   large}$ at late as 
% time progresses.
% The idea of a time-dependent \emph{distribution} of fractal
% dimensions, mirroring the simultaneous existence of different
% cluster morphologies in the system, was suggested by
% \cite{konstandopoulos}.


% Fig. \ref{interaction_potential}  shows the time-evolution of the fractal
% dimension obtained by using Eq. \ref{eq:determination_df} on all the
% clusters sizes for different times. 
% The upper an lower straight lines correspond to the fractal dimensions
% one finds for small and large clusters determined by the scaling law
% Eq. \ref{eq:determination_df} applied to $\overline R_g$.






 \subsection{Comparison between the  Kernel Elements}
Before using Eq. \ref{eq:kernel_calc} to calculate some kernel
elements, we spend a few words clarifying how we evaluate the
non-linear term $n_in_j$.
Since the results presented in this study are ensemble averages taken
on 10 different simulations, we should write Eq. \ref{eq:kernel_calc}
as
\begin{equation}
  \label{eq:kernel_nonlin_aver}
\f{\langle N_{ij}(t)\rangle}{\delta t}=
(2-\delta_{ij})\beta_{ij}\langle n_i(t)n_j(t) \rangle,
\end{equation}
where $\langle\dots\rangle$ denotes the ensemble average.
When evaluating linear quantities from the simulations (such as
e.g. $N_\infty(t)$ or $N_{ij}(t)$) we never introduce explicitly the notation for the
ensemble averages since there is no ambiguity about how the ensemble
average of  a linear quantity is computed.
On the other hand, this has to be made explicit for non-linear
quantities since   in general   $\langle n_in_j\rangle\neq\langle
n_i\rangle \langle n_j\rangle$.
We investigated the behaviour of $\langle n_in_j\rangle$ and $\langle
n_i\rangle \langle n_j\rangle$ finding them to be practically
indistinguishable. 
This justifies (at least for the case of the simulations in this
study) the usual assumption $\langle n_in_j\rangle\simeq\langle
n_i\rangle \langle n_j\rangle$ in the   agglomeration equation.     


We recorded the collisions between aggregates in the computational box rather frequently, every
$\delta t=2/\beta_1$ in order to avoid the occurrence of ternary
collisions in $\delta t$. We have enough data only
for the evaluation of a few kernel elements
$\beta_{ij}$ for low $i$ and $j$ indexes, but they suffice to
 draw some conclusions.
Figure \ref{calculated_beta_ij} (top) shows an example of the fitting
procedure  we use to determine $\beta_{13}$ from the quantities $N_{13}/\delta t$ and
$n_in_j$ evaluated directly from the simulations.
We calculate numerically a few $\beta_{ij}$ which  we compare to the
analytical estimates $\beta_{ij}^{LD}$, both expressed in units of the
diagonal of the continuum kernel $\beta^{Sm}_{ii}=8k_BT/3\mu_f$.
We notice an excellent agreement between numerical and analytical
estimates for the $\{(1,1), (1,2), (1,3), (2,2)\}$
kernel matrix elements, which play an important role in the evolution of $N_\infty(t)$  above all at
early times and are responsible for the agreement between numerical
simulations and agglomeration equation at early times 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 screening 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{table}
% \caption{Values of the dimensionless parameters used in the simulations.}
% \bigskip

% \begin{center}
% \begin{tabular}{|c|c|c|c|c|c|c|c|c|c|c|c|} \hline
% $\rho$ & $\beta$ & $t$ & $N_{\rm mon}$ & $L$ & $\sigma$ & $r_{\rm
%   cut}$ & $r_{\rm min}$ & $k_BT$ & $u_{\rm min}$ & $u_{\sigma}$ \\ \hline \hline
% 0.01 & $1$ & $3000$ & $5000$ & $79.4$ & $1$ & $1.1$ & $1.05$ & $0.5$ &
% $-30$ & $60$ \\
% %120 km/h & $\tau_{\rm ag}\simeq 2$s & $\tau_{\rm dif}\simeq  10^3$s & $\tau_{\rm th}\simeq 30$s & $\tau_{\rm res}\simeq 2$s \\
% \hline
% \end{tabular}
% \label{table_dimensionless}
% \end{center}
% \end{table}


 \clearpage

\begin{table}
\caption{Overview of the cluster diffusional 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



%120 km/h & $\tau_{\rm ag}\simeq 2$s & $\tau_{\rm dif}\simeq  10^3$s & $\tau_{\rm th}\simeq 30$s & $\tau_{\rm res}\simeq 2$s \\
\hline
\end{tabular}
\label{table_diffusion}
\end{center}
\end{table}





\clearpage
\begin{figure}
\includegraphics[width=\columnwidth]{presentation_potential.pdf}
\caption{Model potential and Van der Waals potential used in the
  simulations (energies are in units of $k_BT^*$ and distances in
  units of $\sigma$). }
\label{plot_potential}
\end{figure}




\clearpage
\begin{figure}
\includegraphics[width=0.5\columnwidth]{inset_plot.pdf}
\includegraphics[width=0.5\columnwidth]{cluster_thermalization.pdf}
%\includegraphics[width=0.5\columnwidth, height=6cm]{figure8_b.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  \rangle$ (crosses),
  linear fit (solid line)
  and power-law fit (dashed line). Right: velocity fluctuations as a
  function of time (crosses) and analytical expression in
  Eq. \ref{eq:delta_v_squared} (solid line). }
\label{diffusion_plot}
\end{figure}


\clearpage
\begin{figure}
\includegraphics[width=0.5\columnwidth]{mean_angle.pdf}
\includegraphics[width=0.5\columnwidth]{delta_angle.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$ 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}






\clearpage
\begin{figure}
\includegraphics[width=0.5\columnwidth,height=0.45\columnwidth]{my_output.png}
\includegraphics[width=0.5\columnwidth,height=0.45\columnwidth]{50_monomers.png}
\vspace*{0.4cm}
\includegraphics[width=0.5\columnwidth,height=0.45\columnwidth]{50_monomers-2.png}
\includegraphics[width=0.5\columnwidth,height=0.45\columnwidth]{50_monomers-3.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 \cite{povray}
  software for three-dimensional rendering.}
\label{snapshot_end_simulation}
\end{figure}


\clearpage
\begin{figure}
\includegraphics[width=\columnwidth]{figure_two_slopes_extended_complete.pdf}
% \includegraphics[width=0.5\columnwidth]{evolution_df.pdf}
%\includegraphics[width=0.5\columnwidth, height=6cm]{figure8_b.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}



  
\clearpage
\begin{figure}
\includegraphics[width=\columnwidth]{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}



  
\clearpage
\begin{figure}
\includegraphics[width=\columnwidth]{time_evol_coord_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}


\clearpage
\begin{figure}
\includegraphics[width=\columnwidth]{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 Van
der Waals 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}




\clearpage
\begin{figure}
%\includegraphics[width=0.5\columnwidth]{figure_number_clusters.pdf}
\includegraphics[width=\columnwidth]{plot_kernel_elements_and_kernel_matrix.pdf}
%\includegraphics[width=0.5\columnwidth, height=6cm]{figure8_b.pdf}
\caption{Top: determination of kernel element  $\beta_{13}$. Diamonds: $N_{13}/\delta
  t$. Circles: $2\beta_{13}n_in_3$  with the fitted
  value of  $\beta_{13}$.  Bottom: comparison between some kernel
  elements $\beta_{ij}$ calculated from the simulation and
  corresponding elements $\beta_{ij}^{LD}$, both given in units of $\beta_{ii}^{Sm}$.} 
\label{calculated_beta_ij}
\end{figure}






% \clearpage
% \begin{figure}
% \includegraphics[width=\columnwidth]{number_cluster_vs_time_monodisperse_approx_log-log.pdf}
% %\includegraphics[width=0.5\columnwidth, height=6cm]{figure8_b.pdf}
% \caption{Decay of the total number of clusters from simulations and
%   monodisperse approximation in Eq. \ref{eq:monodisperse_approx}}
% \label{cluster_decay_monodisperse}
% \end{figure}








% \clearpage
% \begin{figure}
% \includegraphics[width=0.5\columnwidth]{figure_number_clusters.pdf}
% %\includegraphics[width=0.5\columnwidth]{coagulation_coefficients.pdf}
% \includegraphics[width=0.5\columnwidth]{evolution_mean_radius_gyration.pdf}

% %\includegraphics[width=0.5\columnwidth, height=6cm]{figure8_b.pdf}
% \caption{Decay of the total number of clusters and contributions from
%   small and large clusters.}
% \label{cluster_decay}
% \end{figure}


% \clearpage
% \begin{figure}
% \includegraphics[width=\columnwidth]{figure_two_slopes.pdf}
% %\includegraphics[width=0.5\columnwidth, height=6cm]{figure8_b.pdf}
% \caption{Fitting of the selected mean $R_g$ as a function of the
%   cluster size for large and small clusters.}
% \label{r_gyr_mean_double_fit}
% \end{figure}





% \clearpage
% \begin{figure}
% %\includegraphics[width=\columnwidth]{coagulation_coefficients.pdf}
% \includegraphics[width=0.5\columnwidth]{figure_two_slopes.pdf}
% %\includegraphics[width=0.5\columnwidth, height=6cm]{figure8_b.pdf}
% \caption{Time-dependent calculation of the fractal dimension together
%   with its value calculated for small clusters, large cluster and its
%   time-averaged value on all cluster sizes.}
% \label{coagulation_coeff}
% \end{figure}





% \clearpage
% \begin{figure}
% \includegraphics[width=\columnwidth]{restructuring.pdf}
% %\includegraphics[width=0.5\columnwidth, height=6cm]{figure8_b.pdf}
% \caption{Evolution of the radius of gyration of a single specific
%   cluster. The vertical lines correspond to the time when the cluster
%   increases its size due to the collision with another cluster or a monomer.}
% \label{single_cluster_size_evolution}
% \end{figure}


\clearpage  
\begin{figure}
\includegraphics[width=0.5\columnwidth]{mean_sep_single_cluster_select.pdf}
\includegraphics[width=0.5\columnwidth]{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}





%radius_gyration_single_cluster_vs_time.pdf



% \clearpage
% \begin{figure}
% \includegraphics[width=\columnwidth]{single_cluster_size_vs_time.pdf}
% %\includegraphics[width=0.5\columnwidth, height=6cm]{figure8_b.pdf}
% \caption{Evolution of the radius of gyration of a single specific
%   cluster. The vertical lines correspond to the time when the cluster
%   increases its size due to the collision with another cluster or a monomer.}
% \label{single_cluster_size_evolution}
% \end{figure}


% \clearpage
% \begin{figure}
% \includegraphics[width=\columnwidth]{radius_gyration_single_cluster_vs_time.pdf}
% %\includegraphics[width=0.5\columnwidth, height=6cm]{figure8_b.pdf}
% \caption{Evolution of the radius of gyration of a single specific
%   cluster. The vertical lines correspond to the time when the cluster
%   increases its size due to the collision with another cluster or a monomer.}
% \label{single_cluster_restructuring}
% \end{figure}


% \clearpage
% \begin{figure}
% \includegraphics[width=\columnwidth]{single_cluster_mean_particle_separation_vs_time.pdf}
% %\includegraphics[width=0.5\columnwidth, height=6cm]{figure8_b.pdf}
% \caption{Evolution of the mean distance between the monomers of a single specific
%   cluster. The vertical lines correspond to the time when the cluster
%   increases its size due to the collision with another cluster or a monomer.}
% \label{single_cluster_restructuring_mean_distance}
% \end{figure}



 






% \clearpage
% \begin{figure}
% %\includegraphics[width=0.5\columnwidth, height=6cm]{figure8_b.pdf}
% \caption{Visualization of the clusters at the end of a simulation.}
% \label{snapshot_end_simulation}
% \end{figure}







%\appendix

% \clearpage
% \begin{figure}
% \includegraphics[width=\columnwidth]{potential.pdf}
% %\includegraphics[width=0.5\columnwidth, height=6cm]{figure8_b.pdf}
% \caption{Plot of the interaction potential used in the simulations.}
% \label{interaction_potential}
% \end{figure}




\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{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{Smirnov}{1990}{smirnov_review} Smirnov, B.M. (1990). The
properties of fractal clusters. {\it Physics Reports, 188}, 1-78.

\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_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_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{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{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{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{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{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{Vemury \& Pratsinis}{1995}{pratsinis_self_preserving} Vemury,
S. \& Pratsinis, S.E. (1995). Self-preserving size distributions of
agglomerates. {\jas}, {\it 26}, 175-185.


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


\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{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{Csardi \& Nepusz}{2006}{igraph}  Csardi, G., \&  Nepusz,
T. (2006). The igraph software package for complex network
research. {\it Inter Journal Complex Systems}, 1695.

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

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

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

\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{Sorensen \& Roberts}{1997}{sorensen_prefactor} Sorensen, C.M. \&
Roberts, G.C., (1997). {\jc}, {\it 186}, 447-452.

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


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


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

\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).


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

\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{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{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{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{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{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{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{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{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{Goldstein}{1950}{goldstein} Goldstein, H. (1950). {\it Classical
  mechanics}, Cambridge, Massachussets. 


\end{thebibliography}


\end{document}

% LocalWords:  Langevin