Stability of nanoparticle laden aerosol liquid droplets

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 2 0 , where R 0 is the initial droplet radius, while the evaporation time scales with R 2 0 .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 finitesize 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 [4][5][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 nanopar-FIG.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.
ticle 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 capillarityapproximation 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
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, σ.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 σ ≈ 100nm (roughly the diameter of a COVID viron), we would find ≈ 3.3 × 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 i.In three dimensions, we have 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 i = (i, j).We introduce two occupation numbers for each lattice site, p n i and p l i , for the nanoparticles and for the liquid, respectively.If a lattice site is empty, then both p n i = p l i = 0.If lattice site i is occupied by a nanoparticle, then p n i = 1, while p l i = 0.If instead lattice site i is full of liquid, then p n i = 0 and p l i = 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 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 Φ l i and Φ n i on the liquid and nanoparticles, respectively.Here, we assume Φ l i = Φ n i = 0. Note that i = Mx i=1 My j=1 , 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, i,j denotes a sum over all pairs of lattice sites.The pair interaction terms involved the discretised pair-potentials ε ll ij , ε nl ij and ε nn ij , which we assume are given by the matrices where the parameters ε ll , ε nl and ε nn determine the overall strength of the pair-interactions and the matrix where N N i denotes the lattice sites that are nearest neighbours to site i, while N N N 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 [21,24], 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 ij 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 [21,22,25].Note that with the sign convention used in Eq. ( 1) for the pair interaction terms, positive values of ε ll , ε nl and ε nn correspond to attractive interactions between neighbouring particles, while negative values (not considered here) correspond to repulsive interactions.The parameter ε 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.

A. 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 Thus, we apply an extension of DFT [11,12] to lattice-systems.We use the following approximation for the Helmholtz free energy of the system [21,29]: where k B is Boltzmann's constant and T is the temperature.To determine the equilibrium density profiles ρ l i and ρ n i corresponding to an aerosol droplet, we minimise the free energy F in the semigrand ensemble, i.e. with fixed total number of nanoparticles in the system but with the chemical potential of the liquid µ ≡ µ l (or equivalently the relative humidity) as an input, so that the total amount of liquid in the system 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 {ρ l i , ρ n i } are those that minimize the semi-grand free energy subject to the constraint that N n , given by ( 5), is equal to the desired value, Nn .In other words, we solve the coupled set of equations N n = Nn .
Differentiating Ω with respect to ρ l i and ρ n i and rearranging gives the conditions We solve this set of equations iteratively, starting from an initial guess for the profiles.A standard Picard solver would take the current state ρ m,old i and replace it with the result of evaluating the right hand sides of Eqs. ( 11) and ( 12) with these densities, denoted ρ m,rhs i .Here, and below, m ∈ {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 {ρ n i } at each step, so that Eq. ( 13) is satisfied.
Note also that it is often necessary to mix the results from the previous step ρ m,old i with the result from evaluating the right hand sides, ρ m,rhs i : where α is typically small, e.g., 0.01 ≤ α ≤ 0.1.This mixing increases the robustness of the scheme, in particular by preventing ρ m,new i 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 ρ l i = ρ l = N l /V and ρ n i = ρ n = N n /V , constants for all i, the Helmholtz free energy per unit volume, f = F/V , where V is the volume of the system, is given by: From this, we can obtain the following relation between the chemical potential of the liquid and the particle densities, The pressure (or equivalently minus the semi-grand potential density) is obtained from the relation For the case where no nanoparticles are present, i.e.where ρ n → 0, the system exhibits vapourliquid phase separation.The critical point occurs for βε ll = 2/3, where . Bulk liquid-vapour phase coexistence occurs for µ = −3ε 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
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 βε ll = 1.2, where the densities of the coexisting (pure) liquid and gas are Note that the droplet decreases in size with decreasing µ, i.e. with decreasing air humidity.Moreover, the liquid density within the droplet is also less.
ρ l = 0.034 and ρ l = 0.966.For the interaction strengths with the nanoparticles, we initially set βε nn = 0.9 and βε nl = 1.5.This choice has ε nn < ε ll < ε nl , which ensures that the nanoparticles prefer to stay well-dispersed within the liquid and do not aggregate.If ε nn were larger and/or if ε 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 ε 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σ × 55σ.In  15)-( 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.βε nn = 0.9, βε nl = 0.6 and liquid chemical potential βµ = −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 ε 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 ε 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 ε 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 where p vap is the bulk pressure in the vapour phase.In our model we assume this is given by Eq. ( 17) in the limit ρ n → 0. Note that the vapour density ρ l is determined by the chemical potential µ, 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 ρ n → 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 The first term is the bulk contribution due to the droplet which has volume (area in 2D) equal to π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; γ is the interfacial tension.The quantity of relevance in the calculations that follow is the difference between these, ∆Ω = Ω drop − Ω vap , which is given by In 3D, the above should be replaced by the following Here we assume that the interfacial tension γ is that of the pure vapour-liquid system at the given temperature.For our model, when βε ll = 1.2, we obtain βγσ = 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) or (in 3D) Substituting Eq. ( 22) into Eq.( 17) [together with Eq. ( 15)], we obtain the following expression for the pressure in the drop (in 2D) Substituting this into Eq.( 20), we obtain an expression for ∆Ω = ∆Ω(ρ l , R), that is a function of (ρ l , R).The equilibrium droplet radius R is the value that minimises the free energy, therefore we require Similarly, substituting Eq. ( 22) into Eq.( 16) yields the following second equation that is also a function of the two unknowns (ρ 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 µ.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, ρ l and R. From these, we can then easily obtain the total amount of liquid in the droplet as (in 2D) from the capillarity approximation (CA), where Γ l is calculated via Eq.( 27), with those from DFT, where Γ l is calculated via Eq.( 29).The results here are for βε ll = 1.2, βε nn = 0.9, and βε nl = 1.or (in 3D) In Fig. 4 we compare the results from this simple CA with results from the DFT for the case when βε ll = 1.2, βε nn = 0.9, βε nl = 1.5 and for various values of the chemical potential µ (i.e. for various values of the humidity).To compare, we plot Γ 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 where ρ l max = sup({ρ l 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 ρ l i > 0.99ρ l max is somewhat arbitrary.However, changing this to e.g.ρ l i > 0.9ρ l max 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 γ obtained for the pure liquid system.In reality, the presence of the nanoparticles changes the value of γ.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 βε 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: 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 i, ρ l i and ρ n i 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].
In Figs. 6 and 7 we show the results for two DDFT simulations.Both start from the DFT equilibrium droplet for βε ll = 1.2, βε nn = 0.9, and βε nl = 1.5 with βµ = −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 µ.This results in the DDFT simulation equilibrating to the corresponding DFT result.We show two cases, βµ = −3.8 and βµ = −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 ℓ 1 norm) between the final dynamic profiles and the corresponding equlibria are below 1% in all cases after time 4 × 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 i.e. the amount of liquid decreases exponentially over time with the rate constant λ = 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 λ 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 ≈ 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.
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 Φ l i and Φ n i 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 βε ll = 1.2, βε nn = 0.9, βε nl = 1.5 and βµ = −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, 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.
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 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 (15) used here with a more accurate equation of state.For example, the Mansoori-Carnahan-Starling-Leland equation of state for hard-sphere mixtures [38,40] 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 [38,46].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].

5 FIG. 2 .
FIG.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 βε ll = 1.2, βε nn = 0.9, βε nl = 1.5 and liquid chemical potential βµ = −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 βµ = −4.5 (comfortable moisture level, with H r = 38%), and all other parameters the same.Bottom row: the corresponding case for βµ = −5 (dry air, with H r = 23%), again with all other parameters the same.
the top row are displayed results for the case when βµ = −3.8, in the middle row for βµ = −4.5 and in the bottom row for βµ = −5.For the temperature considered here (i.e. for the value of βε ll = 1.2, used here), the chemical potential at gas-liquid coexistence is βµ 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 (µ) = 100 × p(µ)/p(µ coex ), where the corresponding pressures are calculated via Eqs.(

8 FIG. 3 .
FIG. 3. Density profiles (liquid on the left, nanoparticles on the right), for the case when βε ll = 1.2,

6 FIG. 4 .
FIG.4.Plots of the amount of liquid Γ l in the droplet as a function of the number of nanoparticles in the droplet N n , for various values of the chemical potential µ, i.e. for varying humidity.We compare results 5. Note that βµ = −3.8corresponds to to a relative humidity of H r = 80%; βµ = −4 corresponds to H r = 65%; βµ = −5 corresponds to H r = 23%; and βµ = −6 corresponds to H r = 8%.Note also that the CA lines for βµ = −5 and βµ = −6 are barely visible, since the corresponding DFT result is almost identical.

6 FIG. 5 .
FIG.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 βε 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 βµ = −5 and βµ = −6 are barely visible, since the corresponding DFT result is almost identical.

8 FIG. 6 .
FIG. 6. Snapshots of DDFT simulations, starting from the DFT equilibrium with βµ = −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 βµ = −3.8(βµ = −5).As expected from the equilibrium calculations, the droplet on the left (right) grows (shrinks) over time.

FIG. 7 .Fig. 6 .
FIG.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) βµ from the value of −5.5 which provides the initial condition causes the amount of liquid in the box to decrease (increase).