Stability of nanoparticle laden aerosol liquid droplets

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 m in less than a few seconds, while smaller droplets evaporate rapidly, as long as the temperature and pressure conditions are such that the vapor of the volatile liquid is in 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,¹ showed that the average sedimentation time scales with $1/R02$⁠, where R₀ is the initial droplet radius, while the evaporation time scales with $R02$⁠. Therefore, a droplet of pure water always either sediments or evaporates sufficiently quickly, preventing it from traveling a significant horizontal distance from the source. In contrast to the droplets described by these simple estimates, aerosols produced by people when they breathe, talk, cough, or sneeze can stay airborne and, therefore, are potentially dangerous for much longer periods, 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 airborne. The droplet is stabilized by the sufficiently large concentration of nanoparticles within. In the equilibrium state, the rate at which liquid evaporates from this droplet is balanced by the condensation of the vapor phase onto the droplet. Clearly, the temperature and relative humidity of the vapor atmosphere have to be taken into account.^3,4

Therefore, studying such droplets allows us 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–6 The fact that this non-volatile material can stabilize the droplets goes some way toward explaining why it was observed that influenza viruses remain stable and infectious in aerosols across a wide range of relative humidities.⁷ Note that the air humidity is directly related to the vapor pressure or, equivalently, is determined by the chemical potential of the vapor. 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 how long clouds persist in the sky before the water they contain returns to the ground as rain.⁸ They are also often the host locations for chemical reactions between airborne species,⁹ 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.¹⁰ 

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 stabilized by the particles within. The recent study by Netz³ 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–15 thereby demonstrating the stability of the droplets and elucidating the 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 the diffusive exchange of solvent molecules from the liquid to the vapor 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 vapor phase and do not explicitly treat 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 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 successfully applied to the study of various properties of liquids condensing and adsorbed in pores and porous media^16–18 and to liquid droplets on surfaces.^19–22 The recent study in Ref. 23 is particularly noteworthy because comparisons with experimental measurements on water droplets containing nanoparticles drying on various surfaces were 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 provide a good description of the thermodynamics and dynamics of aerosol droplets. Additionally, in a future study, we will support our lattice DFT results with additional off-lattice DFT calculations.

This paper proceeds as follows: In Sec. II, we develop a generalized 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 it 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.

In Fig. 1, we illustrate a cross-section through an aerosol liquid droplet containing nanoparticles. To model this, we discretize the system onto a lattice, as also illustrated in Fig. 1. We choose the lattice spacing to be the diameter of the nanoparticles, σ. Therefore, 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 liquid molecules within. As an example of a concrete physical system consisting of the case where the solvent is water and for σ ≈ 100 nm (roughly the diameter of a COVID virion), we would find $≈3.3×107$ 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 (generalized Ising model) and thereby to investigate the thermodynamics of the system. The resulting model can also be thought of as a discretized partial differential equation for the density distributions of the liquid and nanoparticles.²⁴ 

As mentioned earlier, 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 ρ^l = 0.034 and ρ^l = 0.966. For the interaction strengths with the nanoparticles, we initially set βɛⁿⁿ = 0.9 and βɛ^nl = 1.5. This choice has ɛⁿⁿ < ɛ^ll < ɛ^nl, which ensures that the nanoparticles prefer to stay well-dispersed within the liquid and do not aggregate. If ɛⁿⁿ were larger and/or if ɛ^nl were smaller, this could potentially lead to the nanoparticles and the liquid demixing.³⁰ 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 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 to a moist environment, a typical comfortable day-to-day level, and to a fairly dry level, where one may need extra moisturizing 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. (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 dry air, which of course is to be expected.

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.

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

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 that 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. Therefore, we conclude that the simple CA is indeed able to capture the effects of varying attraction strengths on determining the equilibrium droplet size.

We assume that the non-equilibrium dynamics is governed by dynamical density functional theory (DDFT), as described in detail in Ref. 21. We could include hydrodynamics, e.g., following Ref. 22, but here, we 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 ensemble-averaged densities at site i, $ρil$ and $ρin$ 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 Ref. 21, in particular, alternating the direction of the spatial finite difference. Here, we apply this method in both spatial directions, in contrast to Ref. 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.³³ 

In Figs. 6 and 7, we show the results for two DDFT simulations. Both start from the DFT equilibrium droplet for βɛ^ll = 1.2, βɛⁿⁿ = 0.9, and βɛ^nl = 1.5 with βμ = −4.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 with 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 and 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 ℓ₁ norm) between the final dynamic profiles and the corresponding equilibria are below 1% in all cases after time 4 × 10⁴; this can be further reduced by running the DDFT simulations for longer times.

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 $Φil$ and $Φin$ due to a surface, we could model the slow impact of droplets on 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 and 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 cases when βɛ^ll = 1.2, βɛⁿⁿ = 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 center of a single equilibrium droplet for this set of parameter values, while the density outside the circles is set to be that of the corresponding vapor in the single-droplet DFT calculation. The radii of the two circles (i.e., the initial radii of the two droplets) are 5 and 10, with centers 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. Therefore, the only difference between the two simulations is that the smaller droplet is moved slightly closer to the bigger droplet for the right hand set of results. We see, however, that this very small change makes a big difference in the dynamics. In the case on the left, the small droplet shrinks and joins the larger droplet via diffusion through the vapor, 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 and 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 linearizing 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.³⁴ 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.³⁷ 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.

In Fig. 9, we present results corresponding to an initial state where the density of the liquid is uniform throughout the system while the nanoparticles are gathered within a square region in the center of the box. We then perform both the DFT Picard minimization 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. In 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 the mixing parameter α = 0.01.^19,38 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.

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 vapor temperature, humidity, the 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 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–vapor 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,³⁹ crop spraying, and the myriad of different aerosols playing important roles in the world’s atmosphere.^8–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–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, in future work, one could replace the lattice DFT used here with a theory based on the lattice DFT of Refs. 44 and 45, which gives an improved approach for dealing with the nearest-neighbor attractions between particles on a lattice. Alternatively, one could improve the model by using an accurate continuum DFT,¹² 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–vapor interface. For example, size-selectivity can occur in the drying of colloidal films containing two sizes of nanoparticles,^47–49 so such effects may occur in the drying of aerosol droplets containing particle mixtures, potentially resulting in highly structured final states.⁵⁰ 

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

The authors have no conflicts to disclose.

A. J. Archer: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Project administration (equal); Writing – original draft (equal); Writing – review & editing (equal). B. D. Goddard: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Project administration (equal); Writing – original draft (equal); Writing – review & editing (equal). R. Roth: Conceptualization (equal); Funding acquisition (equal); Project administration (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.