# Stability of nanoparticle laden aerosol liquid droplets

A.J. Archer<EMAIL_ADDRESS>Department of Mathematical Sciences and
Interdisciplinary Centre for Mathematical Modelling, Loughborough University,
Loughborough LE11 3TU, UK  B.D. Goddard<EMAIL_ADDRESS>School of
Mathematics and the Maxwell Institute for Mathematical Sciences, University of
Edinburgh, Edinburgh EH9 3FD, UK  R. Roth<EMAIL_ADDRESS>Institute for Theoretical Physics, University of Tübingen, 72076 Tübingen,
Germany

###### Abstract

We develop a model for the thermodynamics and evaporation dynamics of aerosol
droplets of a liquid such as water, surrounded by the gas. When the
temperature and the chemical potential (or equivalently the humidity) are such
that the vapour phase is the thermodynamic equilibrium state, then of course
droplets of the pure liquid evaporate over a relatively short time. However,
if the droplets also contain nanoparticles or any other non-volatile solute,
then the droplets can become thermodynamically stable. We show that the
equilibrium droplet size depends strongly on the amount and solubility of the
nanoparticles within, i.e. on the nature of the particle interactions with the
liquid, and of course also on the vapour temperature and chemical potential.
We develop a simple thermodynamic model for such droplets and compare
predictions with results from a lattice density functional theory that takes
as input the same particle interaction properties, finding very good
agreement. We also use dynamical density functional theory to study the
evaporation/condensation dynamics of liquid from/to droplets as they
equilibrate with the vapour, thereby demonstrating droplet stability.

##  I Introduction

An aerosol droplet is a small liquid drop in the colloidal size range of tens
of nanometers up to a few micrometers in diameter, suspended in a gas like
air. The lifetime of such an aerosol droplet is determined by a competition
between gravity and evaporation [1, 2]. Larger droplets sediment from a height
of the order of 2 meters in less than a few seconds, while smaller droplets
evaporate rapidly, as long as the temperature and pressure conditions are such
that the vapour of the volatile liquid is the thermodynamic equilibrium phase.
In the case of water, such small droplets evaporate completely in the order of
seconds. Wells, who was interested in the spreading of diseases through
aerosols [1] showed that the average sedimentation time scales with
$1/R_{0}^{2}$, where $R_{0}$ is the initial droplet radius, while the
evaporation time scales with $R_{0}^{2}$. Thus, a droplet of pure water always
either sediments or evaporates sufficiently quickly, preventing it from
travelling a significant horizontal distance from the source. In contrast to
droplets described by these simple estimates, aerosols produced by people when
they breathe, talk, cough or sneeze can stay air-borne and thus potentially
dangerous for much longer times, of the order of hours instead of seconds.

This is because, if the droplets contain nanoparticles, e.g. pathogenic germs,
suspended within the liquid or some other non-volatile solute, then a
sufficient quantity of these can stabilize the droplets for hours or even
indefinitely [3, 4]. In other words, once such a droplet is formed, there may
be some fast initial evaporation of the solvent liquid, but in the long time
limit, a finite-size droplet remains that can stay air-borne. The droplet is
stabilised by the sufficiently large concentration of nanoparticles within. In
the equilibrium state, the rate at which liquid evaporates from this droplet
is balanced by condensation of the vapour phase onto the droplet. Clearly, the
temperature and relative humidity of the vapour atmosphere have to be taken
into account [3, 4].

Thus, studying such droplets allows to address questions such as how long does
a droplet of saliva survive in the air before it evaporates? How does this
depend on the amount of solute within? These questions are particularly
relevant if any of that solute is an infectious virus, such as COVID or
influenza. A related additional question of particular interest here is: what
influence do the interactions between the non-volatile material within and the
surrounding liquid have on the droplet size? Saliva droplets contain a
significant amount of material other than water, including various mucus
components, proteins, salts and sometimes virus particles [5, 4, 6]. That this
non-volatile material can stabilize the droplets goes some way to explaining
why it was observed that influenza viruses remain stable and infectious in
aerosols across a wide range of relative humidities [7]. Note that the air
humidity is directly related to the vapour pressure or, equivalently, is
determined by the chemical potential of the vapour. Because these quantities
are so closely connected, we refer to them here almost interchangeably.

Aerosol droplets are not just relevant to the spread of disease: They are
ubiquitous in the global environment, playing a crucial role in atmospheric
science and meteorology. For example, aerosol droplets play a role in
determining the how long clouds persist in the sky, before the water they
contain returns to the ground as rain [8]. They are also often the host
locations for chemical reactions between air-borne species [9] and the
lifetime and stability of aerosol droplets must be considered when addressing
the question of how long harmful chemicals and other pollutants remain in the
atmosphere [10].

In what follows, we refer to all of the non-volatile material that may be
found in an aerosol liquid droplet collectively as ‘nanoparticles’. Droplets
of pure water evaporate in the air, but droplets containing sufficient numbers
of nanoparticles are stabilised by the particles within. The recent study by
Netz [3] addressed this and many of the related issues to do with the
evaporation of saliva droplets in air. When considering the thermodynamic
stability of saliva droplets, Netz used a simple thermodynamics for mixtures,
based on assuming ideal-mixing and ideal-volume additivity. Here, we go beyond
this to develop a model to determine the equilibrium size of such droplets and
how the droplet size depends on the humidity and temperature of the
surrounding gas and also on the nature of the interactions of the solute
particles with the solvent liquid. We show that nanoparticles that have a
higher solubility (i.e. have a greater preference for being dispersed in the
liquid) lead to larger droplets. We develop a simple capillarity-approximation
based thermodynamic model that we validate by comparing with the results from
a lattice density functional theory (DFT), which is a theory for the
microscopic density profile of the liquid and nanoparticles within the droplet
[11, 12]. We also study the evaporation/condensation dynamics of the liquid
from/onto non-equilibrium droplets using dynamical density functional theory
(DDFT) [12, 13, 14, 15], thereby demonstrating the stability of the droplets
and elucidating properties of the formation dynamics. By using DDFT, we assume
that the motion of all particles in the system is diffusive, stemming from our
expectation that the droplet dynamics is dominated by diffusive exchange of
solvent molecules from the liquid to the vapour phase and also that all motion
within the droplet can be treated as diffusive. In our work here, we treat the
surrounding gas as being solely made up of the vapour phase and do not treat
explicitly the inert gas molecules that are present, e.g. in the atmosphere.
However, since these largely play the role of passive spectators, the results
presented here do also apply to the case of aerosol water droplets in air.

The lattice DFT that we use here and variants of it have previously been
applied successfully to the study of various properties of liquids condensing
and adsorbed in pores and porous media [16, 17, 18] and to liquid droplets on
surface [19, 20, 21, 22]. The recent study in Ref. [23] is particularly
noteworthy because comparisons with experimental measurements on water
droplets containing nanoparticles drying on various surfaces was made,
demonstrating good agreement between the theory and the experiments. In view
of this, we are confident that the lattice DFT and DDFT used here also
provides a good description of the thermodynamics and dynamics of aerosol
droplets. Additionally, in a future study, we will support our lattice DFT
results by additional off-lattice DFT calculations.

This paper proceeds a follows: In Sec. II we develop a generalised lattice-gas
model for nanoparticle laden droplets and present the lattice DFT that we use
to determine the density profiles of the liquid and nanoparticles within the
droplets. This theory also yields all the relevant thermodynamic quantities
pertaining to the droplets. Bulk thermodynamic quantities are also briefly
discussed in this section. In Sec. III we present results for the liquid and
nanoparticle density profiles in droplets calculated using our DFT model. In
Sec. IV we then present our capillarity-approximation based thermodynamic
model, comparing with DFT results for the total amount of liquid in the
droplets, as quantities like the number of nanoparticles, humidity and
interaction strength between the nanoparticles and the liquid are varied. In
Sec. V we present our DDFT model and results for the formation dynamics of
droplets. Then, in Sec. VI we compare results for droplet dynamics obtained
from DDFT with the dynamics from Picard iteration, which is a fictitious
dynamics used to solve the equations of DFT. This section will be of interest
to DFT practitioners, but may be skipped by those from a more general
audience. Finally, in Sec. VII we make a few concluding remarks.

##  II Model for nanoparticle laden droplets

Figure 1: On the left is an illustration of the system of interest, namely an
aerosol liquid droplet containing nanoparticles that is surrounded by gas
(i.e. the vapour phase). The red circles represent the suspended
nanoparticles. To treat this, we coarse-grain the system onto a square
lattice, as illustrated on the right. We choose the size of each lattice site
to correspond roughly to the diameter of the nanoparticles. Lattice sites are
described as either occupied with a nanoparticle (red circle) or occupied by
liquid (blue square) or empty (white). We also define effective interaction
potentials between pairs of lattice sites that represent coarse-grained
analogues of the inter-particle interaction potentials in the original system
on the left. Thus, the liquid within the drop has the majority of lattice
sites full of either liquid or nanoparticles, while the vapour outside has
most (but not all) sites being empty.

In Fig. 1 we illustrate a cross-section through an aerosol liquid droplet
containing nanoparticles. To model this, we discretise the system onto a
lattice, as also illustrated in Fig. 1. We choose the lattice spacing to be
the diameter of the nanoparticles, $\sigma$. Thus, lattice sites that contain
a nanoparticle have just one of them within, while lattice sites that are
‘full’ of the liquid, have a large number of the liquid molecules within. As
an example of a concrete physical system, consisting of the case where the
solvent is water and for $\sigma\approx 100$nm (roughly the diameter of a
COVID viron), we would find $\approx 3.3\times 10^{7}$ molecules in each
lattice site ‘full’ of the water. Following Ref. [21], this discretization
onto a lattice allows us to map the system onto a two-species lattice-gas
(generalised Ising model) and thereby to investigate the thermodynamics of the
system. The resulting model can also be thought of as a discretised partial
differential equation for the density distributions of the liquid and
nanoparticles [24].

We denote the location of each lattice site by the index $\mathbf{i}$. In
three dimensions, we have $\mathbf{i}=(i,j,k)$, where $i$, $j$ and $k$ are
integers. To simplify our calculations below, we assume instead that our
system is two-dimensional (2D) and therefore we have $\mathbf{i}=(i,j)$. We
introduce two occupation numbers for each lattice site, $p_{\mathbf{i}}^{n}$
and $p_{\mathbf{i}}^{l}$, for the nanoparticles and for the liquid,
respectively. If a lattice site is empty, then both
$p_{\mathbf{i}}^{n}=p_{\mathbf{i}}^{l}=0$. If lattice site $\mathbf{i}$ is
occupied by a nanoparticle, then $p_{\mathbf{i}}^{n}=1$, while
$p_{\mathbf{i}}^{l}=0$. If instead lattice site $\mathbf{i}$ is full of
liquid, then $p_{\mathbf{i}}^{n}=0$ and $p_{\mathbf{i}}^{l}=1$. We assume that
it is impossible for both occupation numbers to simultaneously equal one. The
potential energy of the system can then be approximated as

$\displaystyle
E=-\sum_{\mathbf{i},\mathbf{j}}\left(\frac{1}{2}\varepsilon_{\mathbf{i}\mathbf{j}}^{ll}p_{\mathbf{i}}^{l}p_{\mathbf{j}}^{l}+\varepsilon_{\mathbf{i}\mathbf{j}}^{nl}p_{\mathbf{i}}^{l}p_{\mathbf{j}}^{n}+\frac{1}{2}\varepsilon_{\mathbf{i}\mathbf{j}}^{nn}p_{\mathbf{i}}^{n}p_{\mathbf{j}}^{n}\right)+\sum_{\mathbf{i}}\Phi_{\mathbf{i}}^{l}\,p_{\mathbf{i}}^{l}+\sum_{\mathbf{i}}\Phi_{\mathbf{i}}^{n}p_{\mathbf{i}}^{n},$
(1)

where the first three terms arise from the interactions between pairs of
occupied lattice sites, while the last two terms are the contributions from
any external potentials $\Phi_{\mathbf{i}}^{l}$ and $\Phi_{\mathbf{i}}^{n}$ on
the liquid and nanoparticles, respectively. Here, we assume
$\Phi_{\mathbf{i}}^{l}=\Phi_{\mathbf{i}}^{n}=0$. Note that
$\sum_{\mathbf{i}}=\sum_{i=1}^{M_{x}}\sum_{j=1}^{M_{y}}$, i.e. this denotes a
sum over all lattice sites in the system, where $M_{x}$ is the number of
lattice sites along the Cartesian $x$-direction indexed by $i$ and $M_{y}$ is
the number along the $y$-direction, indexed by $j$. Similarly,
$\sum_{\mathbf{i},\mathbf{j}}$ denotes a sum over all pairs of lattice sites.
The pair interaction terms involved the discretised pair-potentials
$\varepsilon_{\mathbf{i}\mathbf{j}}^{ll}$,
$\varepsilon_{\mathbf{i}\mathbf{j}}^{nl}$ and
$\varepsilon_{\mathbf{i}\mathbf{j}}^{nn}$, which we assume are given by the
matrices
$\varepsilon_{\mathbf{i}\mathbf{j}}^{ll}=\varepsilon^{ll}c_{\mathbf{i}\mathbf{j}}$,
$\varepsilon_{\mathbf{i}\mathbf{j}}^{nl}=\varepsilon^{nl}c_{\mathbf{i}\mathbf{j}}$
and
$\varepsilon_{\mathbf{i}\mathbf{j}}^{nn}=\varepsilon^{nn}c_{\mathbf{i}\mathbf{j}}$,
where the parameters $\varepsilon^{ll}$, $\varepsilon^{nl}$ and
$\varepsilon^{nn}$ determine the overall strength of the pair-interactions and
the matrix

$c_{\mathbf{i}\mathbf{j}}=\begin{cases}1&\text{if
}\mathbf{j}\in{NN\mathbf{i}},\\\ \frac{1}{2}&\text{if
}\mathbf{j}\in{NNN\mathbf{i}},\\\ 0&\text{otherwise},\end{cases}$ (2)

where ${NN\mathbf{i}}$ denotes the lattice sites that are nearest neighbours
to site $\mathbf{i}$, while ${NNN\mathbf{i}}$ indicates the sites that are the
next nearest neighbours. These three short-ranged (truncated) potentials are
chosen to have the above form for the reasons discussed in [24, 21], namely
that this choice greatly reduces the influence of the approximation of having
discretised the system onto a lattice. If other values for the entries of the
matrix $c_{\mathbf{i}\mathbf{j}}$ were used, we would obtain e.g. non-circular
droplets and other similar influences from the underlying grid. There is a
corresponding choice one can make when applying the model in 3D [25, 21, 22].
Note that with the sign convention used in Eq. (1) for the pair interaction
terms, positive values of $\varepsilon^{ll}$, $\varepsilon^{nl}$ and
$\varepsilon^{nn}$ correspond to attractive interactions between neighbouring
particles, while negative values (not considered here) correspond to repulsive
interactions. The parameter $\varepsilon^{nl}$ determines the overall strength
of attraction between a nanoparticle and any liquid surrounding it. Thus,
varying this parameter most directly determines the free energy to insert a
nanoparticle into the liquid [26, 27, 28], which is the quantity that
determines the solubility of the nanoparticles in the liquid.

###  II.1 Lattice DFT and thermodynamics

Having defined the Hamiltonian (1), one could proceed e.g. by performing
Monte-Carlo computer simulations, as was done in Ref. [25]. However, here we
prefer to follow e.g. Refs. [19, 20, 21] to develop a theory for the ensemble-
averaged densities

$\displaystyle\rho_{\mathbf{i}}^{l}=\langle
p_{\mathbf{i}}^{l}\rangle\,\,\,\,\text{and}\,\,\,\,\rho_{\mathbf{i}}^{n}=\langle
p_{\mathbf{i}}^{n}\rangle.$ (3)

Thus, we apply an extension of DFT [12, 11] to lattice-systems. We use the
following approximation for the Helmholtz free energy of the system [29, 21]:

$\displaystyle F(\\{\rho^{l}_{\mathbf{i}},\rho^{n}_{\mathbf{i}}\\})$
$\displaystyle=k_{B}T\sum_{\mathbf{i}}\left[\rho^{l}_{\mathbf{i}}\ln\rho^{l}_{\mathbf{i}}+(1-\rho^{l}_{\mathbf{i}}-\rho^{n}_{\mathbf{i}})\ln(1-\rho^{l}_{\mathbf{i}}-\rho^{n}_{\mathbf{i}})+\rho^{n}_{\mathbf{i}}\ln\rho^{n}_{\mathbf{i}}\right]$
$\displaystyle\quad-\frac{1}{2}\sum_{\mathbf{i},\mathbf{j}}\varepsilon^{ll}_{\mathbf{i}\mathbf{j}}\rho^{l}_{\mathbf{i}}\rho^{l}_{\mathbf{j}}-\sum_{\mathbf{i},\mathbf{j}}\varepsilon^{ln}_{\mathbf{i}\mathbf{j}}\rho^{l}_{\mathbf{i}}\rho^{n}_{\mathbf{j}}-\frac{1}{2}\sum_{\mathbf{i},\mathbf{j}}\varepsilon^{nn}_{\mathbf{i}\mathbf{j}}\rho^{n}_{\mathbf{i}}\rho^{n}_{\mathbf{j}}+\sum_{\mathbf{i}}\left(\Phi^{l}_{\mathbf{i}}\rho^{l}_{\mathbf{i}}+\Phi^{n}_{\mathbf{i}}\rho^{n}_{\mathbf{i}}\right),$
(4)

where $k_{B}$ is Boltzmann’s constant and $T$ is the temperature. To determine
the equilibrium density profiles $\rho^{l}_{\mathbf{i}}$ and
$\rho^{n}_{\mathbf{i}}$ corresponding to an aerosol droplet, we minimise the
free energy $F$ in the semi-grand ensemble, i.e. with fixed total number of
nanoparticles in the system

$N_{n}=\sum_{\mathbf{i}}\rho^{n}_{\mathbf{i}},$ (5)

but with the chemical potential of the liquid $\mu\equiv\mu^{l}$ (or
equivalently the relative humidity) as an input, so that the total amount of
liquid in the system

$N_{l}=\sum_{\mathbf{i}}\rho^{l}_{\mathbf{i}},$ (6)

is an output of our calculations. We calculate the equilibrium density
profiles using a Picard iteration scheme similar to that described in Ref.
[19]. Since we are in the semi-grand ensemble, the equilibrium density
profiles $\\{\rho_{\mathbf{i}}^{l},\rho_{\mathbf{i}}^{n}\\}$ are those that
minimize the semi-grand free energy

$\Omega=F-\sum_{\mathbf{i}}\mu\rho_{\mathbf{i}}^{l},$ (7)

subject to the constraint that $N_{n}$, given by (5), is equal to the desired
value, $\hat{N}_{n}$. In other words, we solve the coupled set of equations

$\displaystyle\frac{\partial\Omega}{\partial\rho_{\mathbf{i}}^{l}}$
$\displaystyle=0$ (8)
$\displaystyle\frac{\partial\Omega}{\partial\rho_{\mathbf{i}}^{n}}$
$\displaystyle=0$ (9) $\displaystyle N_{n}$ $\displaystyle=\hat{N}_{n}.$ (10)

Differentiating $\Omega$ with respect to $\rho_{\mathbf{i}}^{l}$ and
$\rho_{i}^{n}$ and rearranging gives the conditions

$\displaystyle\rho_{\mathbf{i}}^{l}$
$\displaystyle=(1-\rho_{\mathbf{i}}^{l}-\rho_{\mathbf{i}}^{n})\exp\Big{[}\beta\Big{(}\frac{1}{2}\sum_{\mathbf{j}}(\epsilon_{\mathbf{i}\mathbf{j}}^{ll}+\epsilon_{\mathbf{j}\mathbf{i}}^{ln})\rho_{\mathbf{i}}^{l}+\sum_{\mathbf{j}}\epsilon_{\mathbf{i}\mathbf{j}}^{ln}\rho_{\mathbf{j}}^{n}-\Phi_{\mathbf{i}}^{l}+\mu\Big{)}\Big{]}$
(11) $\displaystyle\rho_{\mathbf{i}}^{n}$
$\displaystyle=(1-\rho_{\mathbf{i}}^{l}-\rho_{\mathbf{i}}^{n})\exp\Big{[}\beta\Big{(}\frac{1}{2}\sum_{\mathbf{j}}(\epsilon_{\mathbf{i}\mathbf{j}}^{nn}+\epsilon_{\mathbf{j}\mathbf{i}}^{ln})\rho_{\mathbf{i}}^{n}+\sum_{\mathbf{j}}\epsilon_{\mathbf{i}\mathbf{j}}^{ln}\rho_{\mathbf{j}}^{l}-\Phi_{\mathbf{i}}^{n}\Big{)}\Big{]}$
(12) $\displaystyle N_{n}$ $\displaystyle=\hat{N}_{n}.$ (13)

We solve this set of equations iteratively, starting from an initial guess for
the profiles. A standard Picard solver would take the current state
$\rho_{\mathbf{i}}^{m,\textrm{old}}$ and replace it with the result of
evaluating the right hand sides of Eqs. (11) and (12) with these densities,
denoted $\rho_{\mathbf{i}}^{m,\textrm{rhs}}$. Here, and below,
$m\in\\{l,n\\}$. Since we are working in the semi-grand ensemble where $N_{n}$
is fixed, we enforce this by renormalising the density profile of the
nanoparticles $\\{\rho_{\mathbf{i}}^{n}\\}$ at each step, so that Eq. (13) is
satisfied. Note also that it is often necessary to mix the results from the
previous step $\rho_{\mathbf{i}}^{m,\textrm{old}}$ with the result from
evaluating the right hand sides, $\rho_{\mathbf{i}}^{m,\textrm{rhs}}$:

$\rho_{\mathbf{i}}^{m,\textrm{new}}=\alpha\rho_{\mathbf{i}}^{m,\textrm{rhs}}+(1-\alpha)\rho_{\mathbf{i}}^{m,\textrm{old}},$
(14)

where $\alpha$ is typically small, e.g., $0.01\leq\alpha\leq 0.1$. This mixing
increases the robustness of the scheme, in particular by preventing
$\rho_{\mathbf{i}}^{m,\textrm{new}}$ from lying outside of the range $(0,1)$.

Before presenting results from this lattice-DFT model, we briefly discussing a
few relevant aspects of the bulk thermodynamics and phase behaviour of the
system. From Eq. (4) we find that for a uniform system with
$\rho_{\mathbf{i}}^{l}=\rho^{l}=N_{l}/V$ and
$\rho_{\mathbf{i}}^{n}=\rho^{n}=N_{n}/V$, constants for all $\mathbf{i}$, the
Helmholtz free energy per unit volume, $f=F/V$, where $V$ is the volume of the
system, is given by:

$\displaystyle f=$ $\displaystyle
k_{B}T[\rho^{l}\ln{\rho^{l}}+(1-\rho^{l}-\rho^{n})\ln(1-\rho^{l}-\rho^{n})+\rho^{n}\ln{\rho^{n}}]$
$\displaystyle-3\varepsilon^{ll}(\rho^{l})^{2}-6\varepsilon^{nl}\rho^{l}\rho^{n}-3\varepsilon^{nn}(\rho^{n})^{2}.$
(15)

From this, we can obtain the following relation between the chemical potential
of the liquid and the particle densities,

$\mu=\frac{\partial f}{\partial\rho^{l}}.$ (16)

The pressure (or equivalently minus the semi-grand potential density) is
obtained from the relation

$p=-f+\mu\rho^{l}.$ (17)

For the case where no nanoparticles are present, i.e. where $\rho^{n}\to 0$,
the system exhibits vapour-liquid phase separation. The critical point occurs
for $\beta\varepsilon^{ll}=2/3$, where $\beta=(k_{B}T)^{-1}$; i.e. the
critical point is at
$(\rho^{l},T)=(\rho^{l}_{c},T_{c})=(\frac{1}{2},\frac{3\varepsilon^{ll}}{2k_{B}})$.
Bulk liquid-vapour phase coexistence occurs for $\mu=-3\varepsilon^{ll}$.
Adding nanoparticles to the system can in general completely change the bulk
phase behaviour – see e.g. Ref. [30]. However, for the interaction parameter
values considered here, generally the effect of the nanoparticles is to just
shift somewhat the coexisting density and chemical potential values, but the
overall qualitative phase behaviour remains the same.

##  III Equilibrium droplet density profiles

Figure 2: Top row: the density profiles (liquid on the left, nanoparticles on
the right), for an equilibrium nanoparticle laden droplet. The system
interaction parameters are $\beta\varepsilon^{ll}=1.2$,
$\beta\varepsilon^{nn}=0.9$, $\beta\varepsilon^{nl}=1.5$ and liquid chemical
potential $\beta\mu=-3.8$ (corresponding to moist air with $H_{r}=80\%$),
while the total number of nanoparticles in the system is constrained to be
$N_{n}=300$. Middle row: the corresponding case for $\beta\mu=-4.5$
(comfortable moisture level, with $H_{r}=38\%$), and all other parameters the
same. Bottom row: the corresponding case for $\beta\mu=-5$ (dry air, with
$H_{r}=23\%$), again with all other parameters the same. Note that the droplet
decreases in size with decreasing $\mu$, i.e. with decreasing air humidity.
Moreover, the liquid density within the droplet is also less.

As mentioned above, to simplify our DFT calculations, we treat the system as
being 2D, but our results could fairly easily be repeated in 3D, at the
expense of more computational cost. We consider the system with
$\beta\varepsilon^{ll}=1.2$, where the densities of the coexisting (pure)
liquid and gas are $\rho^{l}=0.034$ and $\rho^{l}=0.966$. For the interaction
strengths with the nanoparticles, we initially set $\beta\varepsilon^{nn}=0.9$
and $\beta\varepsilon^{nl}=1.5$. This choice has
$\varepsilon^{nn}<\varepsilon^{ll}<\varepsilon^{nl}$, which ensures that the
nanoparticles prefer to stay well-dispersed within the liquid and do not
aggregate. If $\varepsilon^{nn}$ were larger and/or if $\varepsilon^{nl}$ were
smaller, this could potentially lead to the nanoparticles and the liquid
demixing [30]. In Sec. IV below, we present results for different values of
$\varepsilon^{nl}$, which is the parameter that most directly determines the
solubility of the nanoparticles in the liquid.

In Fig. 2 we present both the liquid and nanoparticle density profiles for the
case when the total number of nanoparticles in the system is set to be
$N_{n}=300$, in a square box of size $55\sigma\times 55\sigma$. In the top row
are displayed results for the case when $\beta\mu=-3.8$, in the middle row for
$\beta\mu=-4.5$ and in the bottom row for $\beta\mu=-5$. For the temperature
considered here (i.e. for the value of $\beta\varepsilon^{ll}=1.2$, used
here), the chemical potential at gas-liquid coexistence is
$\beta\mu_{coex}=-3.6$ and so these three chemical potential values correspond
to relative humidity values of $H_{r}=80\%$, $H_{r}=38\%$ and $H_{r}=23\%$,
respectively. If this were for the case of water in the air, then these
$H_{r}$ values correspond roughly speaking a moist environment, a typical
comfortable day-to-day level and to a fairly dry level, where one may need
extra moisturising cream on skin. Note that the relative humidity is defined
as $H_{r}(\mu)=100\times p(\mu)/p(\mu_{coex})$, where the corresponding
pressures are calculated via Eqs. (II.1)–(17). We see that when the chemical
potential is closer to the value at phase-coexistence (more humid air), the
equilibrium droplet size is much larger than in the dry air, which of course
is to be expected.

Figure 3: Density profiles (liquid on the left, nanoparticles on the right),
for the case when $\beta\varepsilon^{ll}=1.2$, $\beta\varepsilon^{nn}=0.9$,
$\beta\varepsilon^{nl}=0.6$ and liquid chemical potential $\beta\mu=-4$
(corresponding to $H_{r}=65\%$), while the total number of nanoparticles in
the system is constrained to be $N_{n}=1000$. The smaller value of
$\varepsilon^{nl}$ used here (compared to that used in Fig. 2), results in a
high nanoparticle density within the drop and interestingly the liquid density
profile is highest on the edge of the droplet.

In Fig. 3 we display results for a case where $\varepsilon^{nl}$ is set to be
much lower than for the cases in Fig. 2. This corresponds to the nanoparticles
having a fairly poor solubility in the liquid. As a consequence, we see from
the density profiles that the nanoparticles gather to form a dense clump with
only a relatively small amount of liquid within. Interestingly, the density of
the liquid is highest around the edge of the droplet. We present this
particular result to illustrate that the density profiles are not always as
simple as those presented in Fig. 2. Before discussing any further DFT results
and in particular how the droplet size varies, depending on the number of
nanoparticles within and the value of the parameter $\varepsilon^{nl}$, we
first present our capillarity-approximation based thermodynamic model in Sec.
IV below. Our aim is to compare results from the microscopic DFT model with
this mesoscopic capillarity-approximation model. We show in the following
section that the two are in excellent agreement.

##  IV Capillarity model

Since our DFT calculations in the previous section are for a 2D system, here
we present our capillarity approximation (CA) model in 2D. However, it is no
more difficult to construct it instead in 3D, so at various points in the
following we additionally give the corresponding 3D equations.

Consider first a system of volume $V$ (strictly, area in 2D) containing just
the bulk vapour phase. The grand potential of the system is then just

$\Omega_{\rm vap}=-p_{\rm vap}V,$ (18)

where $p_{\rm vap}$ is the bulk pressure in the vapour phase. In our model we
assume this is given by Eq. (17) in the limit $\rho^{n}\to 0$. Note that the
vapour density $\rho^{l}$ is determined by the chemical potential $\mu$, which
we assume to be specified. Alternatively, one could assume that the vapour
density (i.e. the humidity) is given and then the chemical potential can be
calculated via Eq. (16), with $\rho^{n}\to 0$.

Now consider the case when the system also contains a circular (in 2D) droplet
of radius $R$, surrounded by the vapour phase; see e.g. Fig. 1. The grand
potential can now be approximated as

$\Omega_{\rm drop}=-p_{\rm drop}\pi R^{2}-p_{\rm vap}(V-\pi R^{2})+\gamma 2\pi
R.$ (19)

The first term is the bulk contribution due to the droplet which has volume
(area in 2D) equal to $\pi R^{2}$; c.f. Eq. (18). The second term is the
corresponding contribution due to the vapour filling the remainder of the
system. The final term is the contribution from the interface between the
liquid and gas phases; $\gamma$ is the interfacial tension. The quantity of
relevance in the calculations that follow is the difference between these,
$\Delta\Omega=\Omega_{\rm drop}-\Omega_{\rm vap}$, which is given by

$\Delta\Omega=-(p_{\rm drop}-p_{\rm vap})\pi R^{2}+\gamma 2\pi R.$ (20)

In 3D, the above should be replaced by the following

$\Delta\Omega_{3D}=-(p_{\rm drop}-p_{\rm vap})\frac{4}{3}\pi R^{3}+\gamma 4\pi
R^{2}.$ (21)

Here we assume that the interfacial tension $\gamma$ is that of the pure
vapour-liquid system at the given temperature. For our model, when
$\beta\varepsilon^{ll}=1.2$, we obtain $\beta\gamma\sigma=0.68587$, which is
calculated in the usual way – see e.g. [11, 19]. The fixed number of
nanoparticles in the droplet means that the nanoparticle density in the
droplet is (in 2D)

$\rho_{n}=\frac{N_{n}}{\pi R^{2}},$ (22)

or (in 3D)

$\rho_{n}=\frac{3N_{n}}{4\pi R^{3}}.$ (23)

Substituting Eq. (22) into Eq. (17) [together with Eq. (II.1)], we obtain the
following expression for the pressure in the drop (in 2D)

$\displaystyle p_{\rm drop}(\rho^{l},R)=$ $\displaystyle
k_{B}T\left[\rho^{l}\ln{\rho^{l}}+\left(1-\rho^{l}-\frac{N_{n}}{\pi
R^{2}}\right)\ln\left(1-\rho^{l}-\frac{N_{n}}{\pi
R^{2}}\right)+\frac{N_{n}}{\pi R^{2}}\ln{\frac{N_{n}}{\pi R^{2}}}\right]$
$\displaystyle-3\varepsilon^{ll}(\rho^{l})^{2}-6\varepsilon^{nl}\rho^{l}\left(\frac{N_{n}}{\pi
R^{2}}\right)-3\varepsilon^{nn}\left(\frac{N_{n}}{\pi R^{2}}\right)^{2}.$ (24)

Substituting this into Eq. (20), we obtain an expression for
$\Delta\Omega=\Delta\Omega(\rho^{l},R)$, that is a function of $(\rho^{l},R)$.
The equilibrium droplet radius $R$ is the value that minimises the free
energy, therefore we require

$\frac{\partial\Delta\Omega(\rho_{l},R)}{\partial R}=0.$ (25)

Similarly, substituting Eq. (22) into Eq. (16) yields the following second
equation

$k_{B}T\ln(\rho^{l})-k_{B}T\ln\left(1-\rho^{l}-\frac{N_{n}}{\pi
R^{2}}\right)-6\varepsilon^{ll}\rho^{l}-6\varepsilon^{nl}\left(\frac{N_{n}}{\pi
R^{2}}\right)-\mu=0,$ (26)

that is also a function of the two unknowns $(\rho^{l},R)$. This equation
corresponds to requiring that the liquid density in the drop equals the higher
of the two possible values determined by the selected value of the chemical
potential $\mu$. In other words, we require the value of the chemical
potential within the droplet to be the same as that in the surrounding vapour.
We then simply solve numerically (using fsolve in Maple) the pair of
simultanous Eqs. (25) and (26) for the two unknowns, $\rho^{l}$ and $R$. From
these, we can then easily obtain the total amount of liquid in the droplet as
(in 2D)

$\Gamma_{l}=\rho^{l}\pi R^{2},$ (27)

or (in 3D)

$\Gamma_{l}=\rho^{l}\frac{4}{3}\pi R^{3}.$ (28)

Figure 4: Plots of the amount of liquid $\Gamma_{l}$ in the droplet as a
function of the number of nanoparticles in the droplet $N_{n}$, for various
values of the chemical potential $\mu$, i.e. for varying humidity. We compare
results from the capillarity approximation (CA), where $\Gamma_{l}$ is
calculated via Eq. (27), with those from DFT, where $\Gamma_{l}$ is calculated
via Eq. (29). The results here are for $\beta\varepsilon^{ll}=1.2$,
$\beta\varepsilon^{nn}=0.9$, and $\beta\varepsilon^{nl}=1.5$. Note that
$\beta\mu=-3.8$ corresponds to to a relative humidity of $H_{r}=80\%$;
$\beta\mu=-4$ corresponds to $H_{r}=65\%$; $\beta\mu=-5$ corresponds to
$H_{r}=23\%$; and $\beta\mu=-6$ corresponds to $H_{r}=8\%$. Note also that the
CA lines for $\beta\mu=-5$ and $\beta\mu=-6$ are barely visible, since the
corresponding DFT result is almost identical. Figure 5: The results here are
the same as those in Fig. 4, except now the strength of the attraction between
the nanoparticles and the liquid is much lower, with
$\beta\varepsilon^{nl}=1$; i.e. these are results for nanoparticles with a
much lower solubility. Note the change in the range of the vertical axis
compared to Fig. 4. Also, note again that the CA lines for $\beta\mu=-5$ and
$\beta\mu=-6$ are barely visible, since the corresponding DFT result is almost
identical.

In Fig. 4 we compare the results from this simple CA with results from the DFT
for the case when $\beta\varepsilon^{ll}=1.2$, $\beta\varepsilon^{nn}=0.9$,
$\beta\varepsilon^{nl}=1.5$ and for various values of the chemical potential
$\mu$ (i.e. for various values of the humidity). To compare, we plot
$\Gamma_{l}$ from the CA in Eq. (27) together with the result from the DFT,
where we determine the amount of liquid in the droplet via the following sum
over lattice sites within the droplet

$\Gamma_{l}=\sum_{\rho_{\mathbf{i}}^{l}>0.99\rho_{\rm
max}^{l}}\rho_{\mathbf{i}}^{l},$ (29)

where $\rho_{\rm max}^{l}=\sup(\\{\rho^{l}_{\mathbf{i}}\\})$ is the maximum
liquid density in the system (i.e. in the centre of the drop). Thus, we do not
include the material in the surrounding vapour, when determining the mass of
liquid in the droplet. Choosing the threshold as
$\rho_{\mathbf{i}}^{l}>0.99\rho_{\rm max}^{l}$ is somewhat arbitrary. However,
changing this to e.g. $\rho_{\mathbf{i}}^{l}>0.9\rho_{\rm max}^{l}$ only
slightly changes our results. From the results in Fig. 4, we see that the
agreement between the CA model and the DFT is rather good. We believe the main
source of error in our CA model is in our decision to use the value for the
interfacial tension $\gamma$ obtained for the pure liquid system. In reality,
the presence of the nanoparticles changes the value of $\gamma$. However,
given the agreement we see in Fig. 4, we conclude that, at least for the
present system, this is a reasonable approximation to make.

In Fig. 5 we show results for the case where $\beta\varepsilon^{nl}=1$ (all
other parameters are the same as in Fig. 4). This corresponds to a much lower
value for the strength of attraction between the nanoparticles and the liquid
and therefore corresponds to nanoparticles which have a much lower solubility.
We see that as a result of the weaker attraction between the nanoparticles and
the liquid, the droplets are therefore smaller, as one should expect. Note
again the good agreement between the CA and the DFT. Thus, we conclude the
simple CA is indeed able to capture the effects of varying attraction
strengths on determining the equilibrium droplet size.

##  V Dynamics

We assume that the non-equilibrium dynamics is governed by dynamical density
functional theory (DDFT), as described in detail in [21]. We could include
hydrodynamics, e.g. following [22], but here restrict ourselves to the usual
overdamped DDFT. DDFT has been derived for Brownian particles suspended in a
liquid [13, 31], such as the nanoparticles studied here, and has also been
developed for molecular liquids [14, 32]. The latter case is particularly
accurate when the fluid is close to equilibrium, which is the case here. As
such, DDFT is applicable to both components of the system studied here.

The (lattice) DDFT for the two component system consists of a coupled pair of
equations at each lattice site:

$\displaystyle\frac{\partial\rho_{\mathbf{i}}^{l}}{\partial t}$
$\displaystyle=\nabla\cdot\Big{[}M_{l}\rho_{\mathbf{i}}^{l}\nabla\frac{\partial
F}{\partial\rho_{\mathbf{i}}^{l}}\Big{]}$ (30)
$\displaystyle\frac{\partial\rho_{\mathbf{i}}^{n}}{\partial t}$
$\displaystyle=\nabla\cdot\Big{[}M_{n}\rho_{\mathbf{i}}^{n}\nabla\frac{\partial
F}{\partial\rho_{\mathbf{i}}^{n}}\Big{]}.$ (31)

Note that in contrast to the DDFT presented in [21] we have constant mobility
coefficients $M_{l}$ and $M_{n}$, for both the liquid and nanoparticles. It is
simple to extend the DDFT to the non-constant case, but we do not find it
necessary here. Note also that we set $M_{l}=M_{n}=1$ for simplicity.

The ensemble-averaged densities at site $\mathbf{i}$, $\rho_{\mathbf{i}}^{l}$
and $\rho_{\mathbf{i}}^{n}$ are now functions of time. The differential
operators in (30) and (31) are finite difference approximations on the
lattice. Care is needed when applying these to prevent numerical
instabilities. We take the approach detailed in [21], in particular
alternating the direction of the spatial finite difference. Here we apply this
method in both spatial directions, in contrast to [21] where this was only
required in the direction parallel to their wall. For the time evolution we
use an Euler scheme but note that the results are almost indistinguishable
from those using higher order schemes such as the fourth order adaptive Runge-
Kutta scheme implemented in Matlab’s ode45 [33].

| $\beta\mu=-3.8$ | $\beta\mu=-5$
---|---|---
$t=0$ |  |
$t=4,000$ |  |
$t=20,000$ |  |
Figure 6: Snapshots of DDFT simulations, starting from the DFT equilibrium
with $\beta\mu=-4.5$ at times $t=0$, $t=4000$, and $t=20000$ (top to bottom).
Pairs of plots on the left (right) show the DDFT dynamics for the liquid and
the nanoparticles when the liquid density on the boundary of the box is set to
the value given by the equilibrium DFT computations with $\beta\mu=-3.8$
($\beta\mu=-5$). As expected from the equilibrium calculations, the droplet on
the left (right) grows (shrinks) over time. Figure 7: The mass evolution of
the liquid and nanoparticles over time for the DDFT simulations shown in Fig.
6. Dashed lines denote the liquid masses in the corresponding DFT
computations. Note, in particular, that increasing (decreasing) $\beta\mu$
from the value of $-5.5$ which provides the initial condition causes the
amount of liquid in the box to decrease (increase).

In Figs. 6 and 7 we show the results for two DDFT simulations. Both start from
the DFT equilibrium droplet for $\beta\varepsilon^{ll}=1.2$,
$\beta\varepsilon^{nn}=0.9$, and $\beta\varepsilon^{nl}=1.5$ with
$\beta\mu=-5.5$ (also displayed in Fig. 2). To induce dynamics, we set the
liquid density on the boundary of the box to a constant value, corresponding
to that of a DFT simulation with a different value of $\mu$. This results in
the DDFT simulation equilibrating to the corresponding DFT result. We show two
cases, $\beta\mu=-3.8$ and $\beta\mu=-5$. The former case, displayed in the
two left hand columns of Fig. 6, corresponds to a droplet having moved into a
more humid environment than that in which it was initially formed. This
results in an increase in the amount of liquid in the box, as additional
liquid diffuses into the box from the boundary and then condenses onto the
droplet, making it grow in size. In contrast, the DDFT results in the two
right hand columns of Fig. 6 correspond to the droplet moving to an
environment that has a lower humidity than that where it was formed. This is
the typical case of aerosol droplets that form in a person’s mouth or
respiratory system at a temperature of about 37∘C at 100% relative humidity
that then subsequently move out of the mouth into the air, which is a less
humid environment. The DDFT results show the size of the droplet decreasing
over time, as it reduces to a smaller equilibrium size. Note, however, that
the droplet never completely evaporates; the nanoparticles within mean that it
remains stable, albeit at a smaller size. This can also be seen in Fig. 7,
where we plot the total mass of liquid in the simulation box as a function of
time for these two cases. The relative errors (in the $\ell_{1}$ norm) between
the final dynamic profiles and the corresponding equlibria are below 1% in all
cases after time $4\times 10^{4}$; this can be further reduced by running the
DDFT simulations for longer times.

We note that a crucial aspect of these DDFT simulations relates to the choice
of boundary conditions. In the examples presented in Figs. 6 and 7, the
nanoparticles are largely concentrated in the middle of the simulation domain
and the interparticle attraction prevents any significant diffusion away from
this area and, as such, periodic boundary conditions for the nanoparticles are
an appropriate choice (any other sensible boundary condition will give almost
indistinguishable results). However, the boundary conditions used for the
liquid are much more important. In particular, the distance from the surface
of the droplet to the edge of the boundary where we keep the liquid density
set to the value corresponding to the desired chemical potential value (i.e.
relative humidity value, specified for each simulation), is all-important. If
the box size is increased, i.e. the distance from the droplet to the boundary
is increased, then the DDFT simulations take correspondingly longer to
equilibrate. What determines the time for a droplet to equilibrate is a
combination of two processes: the first relates to the time it takes for
liquid to move out of the droplet, across the liquid-vapour interface. The
second part of the process is that of the liquid diffusing through the vapour
surrounding the droplet, to reach the boundary and be absorbed. The second
process is well understood: for free diffusion from the centre to the edge of
a circular domain, the total amount in the system $N_{l}(t)$, given by Eq.
(6), follows the well-known result

$N_{l}(t)\approx[N_{l}(0)-N_{l}(t\to\infty)]e^{-\lambda t}+N_{l}(t\to\infty),$
(32)

i.e. the amount of liquid decreases exponentially over time with the rate
constant $\lambda=D(j_{1,0}/a)^{2}$, where $D=M_{l}k_{B}T$ is the diffusion
coefficient, $j_{1,0}$ is the first zero of the Bessel $J_{0}(x)$ function,
and $a$ is the radius of the domain. This is the situation our model reduces
to in the limit where the densities of the liquid and nanoparticles are small
everywhere (where the DDFT equations (30) and (31) reduce to diffusion
equations). However, for the cases of interest here, the additional process of
particles crossing the liquid-vapour interface makes the whole equilibration
process much slower. In our DDFT simulations we still observe $N_{l}(t)$
varying over time with the simple exponential decay form in Eq. (32) (see Fig.
7), but the rate constant $\lambda$ that we observe is much smaller than the
result quoted above for the case for simple diffusion from the centre to the
edge of the domain, due to the additional interface crossing process.
Nonetheless, these considerations demonstrate why the distance from the
droplet to the edge of the simulation box (i.e. the size we assume for the
diffusive boundary layer around the droplets) is important for determining the
overall time scale of the equilibration process. That said, we find that as
long as the edge of droplet is $\approx 5$ or more lattice sites away from the
boundary of the box, then snapshots over time from simulations in a small box
and a larger box are almost indistinguishable.

Figure 8: Snapshots from DDFT simulations, starting from two circular
distributions at times $t=0$, $t=20$, $t=1200$, and $t=2000$ (left, top to
bottom) and $t=0$, $t=10$, $t=100$, and $t=1000$ (right, top to bottom). The
simulations differ only in the initial location of the smaller droplet (see
text for details). However, this small difference in initial condition makes a
very significant difference to the manner in which the two droplets coalesce.

Our DDFT model can be used to predict the dynamics of aerosol droplets in a
great variety of different situations. For example, if we included the
external potentials $\Phi_{\mathbf{i}}^{l}$ and $\Phi_{\mathbf{i}}^{n}$ due to
a surface, we could model the slow impact of droplets with a surface and the
subsequent spreading and drying process. Illustrative results for the later
part of this dynamics can be found e.g. in Refs. [21, 23]. Here, we restrict
ourselves to presenting a pair of illustrative results corresponding to the
coalescence of two different sized droplets. In Fig. 8 we show results for
droplets joining for the case when $\beta\varepsilon^{ll}=1.2$,
$\beta\varepsilon^{nn}=0.9$, $\beta\varepsilon^{nl}=1.5$ and $\beta\mu=-5$.
The initial conditions correspond simply to setting all of the lattice sites
within two circular regions to the density values at the centre a single
equilibrium droplet for this set of parameter values, while the density
outside the circles is set to be that of the corresponding vapour in the
single-droplet DFT calculation. The radii of the two circles (i.e. initial
radii of the two droplets) are 5 and 10 with centres at (20,20) and (30,35)
for the left hand simulations and (21,21) and (30,35) for the right hand
simulations in Fig. 8. Thus, the only difference between the two simulations
is that the smaller droplet is moved sightly closer to the bigger droplet for
the right hand set of results. We see however that this very small change
makes a big difference to the dynamics. For the case on the left, the small
droplet shrinks and joins the larger droplet via diffusion through the vapour,
while in the case on the right the whole droplet moves and joins the larger
one. Why it is that one sees one process at one distance and the other at a
slightly different distance was studied in detail in Ref. [34] in the context
of a different DDFT. The mechanism followed by the case on the left is termed
joining via the Ostwald mode, which was first described in Refs. [35, 36] to
understand the process of Ostwald ripening, while the mechanism followed by
the case on the right is termed the translation mode. One can calculate which
mode will dominate by linearising the DDFT equation around the initial state
and then one obtains two distinct eigenfunctions corresponding to each of
these modes. The mode one observes actually occurring is the one with the
largest corresponding eigenvalue [34]. These results illustrate just one
possibility in the hugely complex dynamics of aerosol droplets. Note also that
our DDFT model assumes an over-damped (diffusive) dynamics. If we extended our
theory to include the effects of inertia, then, for example, the evaporation
and coalescence of droplets in the turbulent airflow following a person
sneezing could be investigated [37]. However, we do not pursue that direction
here. The examples in Fig. 8 illustrate that the interplay of dynamics with an
underlying complex free energy landscape can result in rather complicated
dynamics. In the following section, we illustrate this point further, albeit
with examples that correspond to a somewhat unlikely (in Nature) initial
state.

##  VI Comparing DDFT to Picard

Figure 9: Snapshots from a DDFT simulation (left) and from Picard iteration (a
fictitious dynamics, right), starting from an initial condition where the
liquid density is uniform and the nanoparticles are in a square region at the
centre of the box. The initial average density of the liquid is selected so
that the final states are the same as the equilibrium DFT calculation for
$\beta\mu=-5$. The final equilibria from the two dynamics are the same, but
the intermediate states during the evolution are very different. The DDFT
profiles on the left are for the times $t=0$, 10, 30, 120 and $10^{4}$. The
Picard profiles are from iterations 0, 200, 600, 2400 and 4000.

In Fig. 9 we present results corresponding to an initial state where the
density of the liquid is uniform throughout the system, while the the
nanoparticles are gathered within a square region in the centre of the box. We
then perform both the DFT Picard minimisation and the DDFT simulation (with
periodic boundary conditions for both the fluid and nanoparticles). The
initial and final average densities of the liquid in the box are the same for
both systems. For the DDFT, the average density is a conserved quantity
throughout the dynamics. For the Picard iteration, in this case it is roughly
constant, but more generally this fictitious dynamics does not preserve mass
between iterations. Our Picard iteration results are obtained with mixing
parameter $\alpha=0.01$ [38, 19]. Interestingly, the results from these two
approaches have very different paths to (the same) equilibrium. This
demonstrates that in this case the ‘quasi-dynamics’ generated by the Picard
scheme is not a good approximation of the DDFT dynamics. We see some rather
striking transient states during the (realistic) DDFT dynamics displayed on
the left of Fig. 9, where the initial square block of nanoparticles breaks up
into four smaller droplets that then subsequently re-coalesce into the single
final droplet state. This complex dynamics is driven by a competition between
bulk and interfacial contributions to the free energy, with each dominating at
different stages of the dynamics.

##  VII Concluding remarks

We have presented a simple capillarity-approximation based theory for the size
of nanoparticle laden liquid aerosol droplets. Our theory predicts how the
size of the droplets varies depending on the vapour temperature, humidity,
number of nanoparticles within the droplet and also the nature of the
interactions between the nanoparticles and the liquid. We have validated our
simple theory by comparing it with results from DFT. Our lattice DFT theory
yields the density distribution of the particles within the aerosol droplets
in addition to all the relevant thermodynamic quantities, such as the changes
in the liquid-vapour interfacial tension due to varying concentrations of
nanoparticles within the droplet. We have also developed a DDFT model, able to
describe complex dynamical phenomena, such as droplet coalescence. Concrete
examples of the types of airborne aerosol systems that our model can be
applied to include determining the stability (and therefore the lifetime) of
exhaled droplets that can lead to the spread of COVID and other diseases,
aerosol based therapies such as biomolecule inhalation therapy [39], crop
spraying, and the myriad of different aerosols playing important roles in the
world’s atmosphere [8, 9, 10].

Our simple CA model is useful for quick estimation of droplet sizes as a
function of particle loading. For more precise calculations, our CA model
could easily be improved by replacing the simple lattice-gas free energy
(II.1) used here with a more accurate equation of state. For example, the
Mansoori-Carnahan-Starling-Leland equation of state for hard-sphere mixtures
[40, 38] could easily be used instead, or one of the many other accurate bulk
fluid equations of state that are available in the literature – see e.g. Refs.
[41, 42, 43]. The choice of a particular equation of state would be guided by
obtaining additional information about the precise form of the molecular
interactions in the system.

For those interested in a more accurate description of the density
distribution of the liquid and nanoparticles within droplets, as future work
one could replace the lattice DFT used here with a theory based on the lattice
DFT of Refs. [44, 45], which gives an improved approach for dealing with the
nearest-neighbour attractions between particles on a lattice. Alternatively,
one could improve the model by using an accurate continuum DFT [12], such as a
DFT based on fundamental measure theory [46, 38]. Such an approach would give
a much better description of the liquid structure within droplets. This
approach may be needed especially in cases where particles aggregate at the
droplet liquid-vapour interface. For example, size-selectivity can occur in
the drying of colloidal films containing two sizes of nanoparticles [47, 48,
49], so such effects may occur in the drying of aerosol droplets containing
particle mixtures, potentially resulting in highly structured final states
[50].

## Acknowledgements

This research was funded by the London Mathematical Society, the International
Centre for Mathematical Sciences, and Loughborough University Institute of
Advanced Studies. We are grateful to Emiliano Renzi and David Sibley for
valuable discussions.

## References

  * Wells [1934] W. F. Wells, On air-borne infection. Study II. Droplets and droplet nuclei, Am. J. Hyg. 20, 611–618 (1934).
  * Xie  _et al._ [2007] X. Xie, Y. Li, A. T. Y. Chwang, P. L. Ho, and W. H. Seto, How far droplets can move in indoor environments – revisiting the wells evaporation–falling curve, Indoor Air 17, 211 (2007).
  * Netz [2020] R. R. Netz, Mechanisms of airborne infection via evaporating and sedimenting droplets produced by speaking, J. Phys. Chem. B 124, 7093 (2020).
  * Seyfert  _et al._ [2022] C. Seyfert, J. Rodríguez-Rodríguez, D. Lohse, and A. Marin, Stability of respiratory-like droplets under evaporation, Phys. Rev. Fluids 7, 023603 (2022).
  * Poon  _et al._ [2020] W. C. K. Poon, A. T. Brown, S. O. L. Direito, D. J. Hodgson, L. Le Nagard, A. Lips, C. E. MacPhee, D. Marenduzzo, J. R. Royer, A. F. Silva, _et al._ , Soft matter science and the COVID-19 pandemic, Soft matter 16, 8310 (2020).
  * Vejerano and Marr [2018] E. P. Vejerano and L. C. Marr, Physico-chemical characteristics of evaporating respiratory fluid droplets, J. R. Soc. Interface 15, 20170939 (2018).
  * Kormuth  _et al._ [2018] K. A. Kormuth, K. Lin, A. J. Prussin, E. P. Vejerano, A. J. Tiwari, S. S. Cox, M. M. Myerburg, S. S. Lakdawala, and L. C. Marr, Influenza virus infectivity is retained in aerosols and droplets independent of relative humidity, J. Infect. Dis. 218, 739 (2018).
  * Stevens and Feingold [2009] B. Stevens and G. Feingold, Untangling aerosol effects on clouds and precipitation in a buffered system, Nature 461, 607 (2009).
  * Wei  _et al._ [2018] H. Wei, E. P. Vejerano, W. Leng, Q. Huang, M. R. Willner, L. C. Marr, and P. J. Vikesland, Aerosol microdroplets exhibit a stable pH gradient, Proc. Natl. Acad. Sci. U.S.A. 115, 7272 (2018).
  * Von Schneidemesser  _et al._ [2015] E. Von Schneidemesser, P. S. Monks, J. D. Allan, L. Bruhwiler, P. Forster, D. Fowler, A. Lauer, W. T. Morgan, P. Paasonen, M. Righi, K. Katerina Sindelarova, and M. A. Sutton, Chemistry and the linkages between air quality and climate change, Chem. Rev. 115, 3856 (2015).
  * Evans [1979] R. Evans, The nature of the liquid-vapour interface and other topics in the statistical mechanics of non-uniform, classical fluids, Adv. Phys. 28, 143 (1979).
  * Hansen and McDonald [2013] J.-P. Hansen and I. R. McDonald,  _Theory of simple liquids: with applications to soft matter_ (Academic press, 2013).
  * Marconi and Tarazona [1999] U. M. B. Marconi and P. Tarazona, Dynamic density functional theory of fluids, J. Chem. Phys. 110, 8032 (1999).
  * Archer [2006] A. J. Archer, Dynamical density functional theory for dense atomic liquids, J. Phys.: Condens. Matter 18, 5617 (2006).
  * te Vrugt  _et al._ [2020] M. te Vrugt, H. Löwen, and R. Wittkowski, Classical dynamical density functional theory: From fundamentals to applications, Adv. Phys. 69, 121 (2020).
  * Kierlik  _et al._ [2001] E. Kierlik, P. A. Monson, M. L. Rosinberg, L. Sarkisov, and G. Tarjus, Capillary condensation in disordered porous materials: Hysteresis versus equilibrium behavior, Phys. Rev. Lett. 87, 055701 (2001).
  * Woo  _et al._ [2001] H.-J. Woo, L. Sarkisov, and P. A. Monson, Mean-field theory of fluid adsorption in a porous glass, Langmuir 17, 7472 (2001).
  * Schneider  _et al._ [2014] D. Schneider, R. Valiullin, and P. A. Monson, Filling dynamics of closed end nanocapillaries, Langmuir 30, 1290 (2014).
  * Hughes  _et al._ [2014] A. P. Hughes, U. Thiele, and A. J. Archer, An introduction to inhomogeneous liquids, density functional theory, and the wetting transition, Am. J. Phys. 82, 1119 (2014).
  * Hughes  _et al._ [2015] A. P. Hughes, U. Thiele, and A. J. Archer, Liquid drops on a surface: Using density functional theory to calculate the binding potential and drop profiles and comparing with results from mesoscopic modelling, J. Chem. Phys. 142, 074702 (2015).
  * Chalmers  _et al._ [2017a] C. Chalmers, R. Smith, and A. J. Archer, Dynamical density functional theory for the evaporation of droplets of nanoparticle suspension, Langmuir 33, 14490 (2017a).
  * Areshi  _et al._ [2019] M. Areshi, D. Tseluiko, and A. J. Archer, Kinetic monte carlo and hydrodynamic modeling of droplet dynamics on surfaces, including evaporation and condensation, Phys. Rev. Fluids 4, 104006 (2019).
  * Perez  _et al._ [2021] C. M. Perez, M. Rey, B. D. Goddard, and J. H. J. Thijssen, Changing the flow profile and resulting drying pattern of dispersion droplets via contact angle modification, arXiv preprint arXiv:2111.00464 (2021).
  * Robbins  _et al._ [2011] M. J. Robbins, A. J. Archer, and U. Thiele, Modelling the evaporation of thin films of colloidal suspensions using dynamical density functional theory, J. Phys.: Condens. Matter 23, 415102 (2011).
  * Chalmers  _et al._ [2017b] C. Chalmers, R. Smith, and A. J. Archer, Modelling the evaporation of nanoparticle suspensions from heterogeneous surfaces, J. Phys.: Condens. Matter 29, 295102 (2017b).
  * Chipot and Pohorille [2007] C. Chipot and A. Pohorille,  _Free energy calculations_ , Vol. 86 (Springer, 2007).
  * Roth  _et al._ [2006] R. Roth, Y. Harano, and M. Kinoshita, Morphometric approach to the solvation free energy of complex molecules, Phys. Rev. Lett. 97, 078101 (2006).
  * Coe  _et al._ [2023] M. K. Coe, R. Evans, and N. B. Wilding, Understanding the physics of hydrophobic solvation, J. Chem. Phys. 158, 034508 (2023).
  * Woywod and Schoen [2003] D. Woywod and M. Schoen, Phase behavior of confined symmetric binary mixtures, Phys. Rev. E 67, 026122 (2003).
  * Areshi [2020] M. Areshi,  _Mathematical Modelling of Liquids on Surfaces_ , Ph.D. thesis, Loughborough University (2020).
  * Marconi and Tarazona [2000] U. M. B. Marconi and P. Tarazona, Dynamic density functional theory of fluids, J. Phys.: Condens. Matter 12, A413 (2000).
  * Archer [2009] A. J. Archer, Dynamical density functional theory for molecular and colloidal fluids: A microscopic approach to fluid mechanics, J. Chem. Phys. 130, 014509 (2009).
  * Dormand and Prince [1980] J. R. Dormand and P. J. Prince, A family of embedded Runge-Kutta formulae, J. Comput. Appl. Math. 6, 19 (1980).
  * Pototsky  _et al._ [2014] A. Pototsky, U. Thiele, and A. J. Archer, Coarsening modes of clusters of aggregating particles, Phys. Rev. E 89, 032144 (2014).
  * Lifshitz and Slyozov [1961] I. M. Lifshitz and V. V. Slyozov, The kinetics of precipitation from supersaturated solid solutions, J. Phys. Chem. Solids 19, 35 (1961).
  * Wagner [1961] C. Wagner, Theorie der alterung von niederschlägen durch umlösen (ostwald-reifung), Z. Elektrochem. 65, 581 (1961).
  * Renzi and Clarke [2020] E. Renzi and A. Clarke, Life of a droplet: Buoyant vortex dynamics drives the fate of micro-particle expiratory ejecta, Phys. Fluids 32, 123301 (2020).
  * Roth [2010] R. Roth, Fundamental measure theory for hard-sphere mixtures: a review, J. Phys.: Condens. Matter 22, 063102 (2010).
  * Roudini  _et al._ [2023] M. Roudini, J. M. Rosselló, O. Manor, C.-D. Ohl, and A. Winkler, Acoustic resonance effects and cavitation in saw aerosol generation, Ultrasonics Sonochemistry 98, 106530 (2023).
  * Mansoori  _et al._ [1971] G. A. Mansoori, N. F. Carnahan, K. E. Starling, and T. W. Leland Jr, Equilibrium thermodynamic properties of the mixture of hard spheres, J. Chem. Phys. 54, 1523 (1971).
  * Chapman  _et al._ [1989] W. G. Chapman, K. E. Gubbins, G. Jackson, and M. Radosz, SAFT: Equation-of-state solution model for associating fluids, Fluid Ph. Equilibria 52, 31 (1989).
  * Kontogeorgis  _et al._ [2006] G. M. Kontogeorgis, M. L. Michelsen, G. K. Folas, S. Derawi, N. Von Solms, and E. H. Stenby, Ten years with the CPA (cubic-plus-association) equation of state. Part 2. Cross-associating and multicomponent systems, Ind. Eng. Chem. Res. 45, 4869 (2006).
  * Gross and Vrabec [2006] J. Gross and J. Vrabec, An equation-of-state contribution for polar components: dipolar molecules, AIChE journal 52, 1194 (2006).
  * Maeritz and Oettel [2021a] M. Maeritz and M. Oettel, Density functional for the lattice gas from fundamental measure theory, Phys. Rev. E 104, 024124 (2021a).
  * Maeritz and Oettel [2021b] M. Maeritz and M. Oettel, Droplet condensation in the lattice gas with density functional theory, Phys. Rev. E 104, 034127 (2021b).
  * Rosenfeld [1989] Y. Rosenfeld, Free-energy model for the inhomogeneous hard-sphere fluid mixture and density-functional theory of freezing, Phys. Rev. Lett 63, 980 (1989).
  * Fortini  _et al._ [2016] A. Fortini, I. Martín-Fabiani, J. L. De La Haye, P.-Y. Dugas, M. Lansalot, F. D’Agosto, E. Bourgeat-Lami, J. L. Keddie, and R. P. Sear, Dynamic stratification in drying films of colloidal mixtures, Phys. Rev. Lett. 116, 118301 (2016).
  * He  _et al._ [2021] B. He, I. Martín-Fabiani, R. Roth, G. I. Tóth, and A. J. Archer, Dynamical density functional theory for the drying and stratification of binary colloidal dispersions, Langmuir 37, 1399 (2021).
  * Kundu and Howard [2022] M. Kundu and M. P. Howard, Dynamic density functional theory for drying colloidal suspensions: Comparison of hard-sphere free-energy functionals, J. Chem. Phys. 157, 184904 (2022).
  * Liu  _et al._ [2019] W. Liu, J. Midya, M. Kappl, H.-J. Butt, and A. Nikoubashman, Segregation in drying binary colloidal droplets, ACS Nano 13, 4972 (2019).