Microphysics of liquid water in sub-10 nm ultrafine aerosol particles

Ultrafine aerosol particles with sizes smaller than 50 nm have been shown in recent studies to serve as a large source of cloud condensation nuclei that can promote additional cloud droplet formation under supersaturation conditions. Knowledge of the microphysics of liquid water in these droplets remains limited, particularly in the sub-10 nm particle size range, due to experimental and theoretical challenges associated with the complexity of aerosol components and the small length scales of interest (e.g., difficulty of precisely sampling the liquid–air interface, questionable validity of mean-field theoretical representations). Here, we carried out molecular dynamics (MD) simulations of aerosol particles with diameters between 1 and 10 nm and characterized atomistic-level structure and water dynamics in well-mixed and phase-separated systems with different particle sizes, NaCl salinities, and pimelic acid (PML) organic surface loadings as a function of distance from the time-averaged Gibbs dividing interface or instantaneous water–air interface. We define a sphericity factor (ϕ) that can shed light on the phase-mixing state of nanodroplets, and we reveal an unexpected dependence of mixing state on droplet size. Our results also evidence an ion concentration enhancement in ultrafine aerosols, which should modulate salt nucleation kinetics in sub-10 nm droplets, and provide detailed characterization of the influence of droplet size on surface tension and on water self-diffusivity near the interface. Analysis of water evaporation free energy and water activity demonstrates the validity of the Kelvin equation and Köhler theory at droplet sizes larger than 4 nm under moderate salinities and organic loadings and the need for further extension to account for ion concentration enhancement in sub-10 nm aerosols, droplet-size-dependent phase separation effects, and a sharp decrease in the cohesiveness of liquid water in sub-4 nm droplets. Finally, we show that an idealized fractional surface coating factor (f_s) can be used to categorize and reconcile water accommodation coefficients (α^*) observed in MD simulations and experimental results in the presence of organic coatings, and we resolve the droplet size dependence of α^*.

1 Introduction

Water is a predominant component of atmospheric aerosol, with mass fraction typically higher than 70 % at moderate to high relative humidities (RH > 40 %), and plays crucial roles in the formation and growth of aerosol particles (Pierce et al., 2012; Karlsson et al., 2022). For example, aerosol liquid water facilitates the partitioning of gas phase water-soluble organic matter into the particle phase and provides a medium for aqueous reactions to form secondary aerosol particles (Faust et al., 2017; Sareen et al., 2017; Frank et al., 2020). In addition, water vapor condensation on nanoparticles can alter the phase state of aerosols and drive aerosol growth to climatically active sizes of about 80 nm and higher (Wu et al., 2018).

Because of water's abundance in atmospheric aerosols, knowledge of the microphysics of aerosol water is essential to resolving key unknowns in aerosol growth and evolution that underlie uncertainties in aerosol–cloud interactions in current weather and climate forecasts (Sato and Suzuki, 2019; Benjamin et al., 2019; Kreidenweis et al., 2019; Lawler et al., 2021). In particular, recent studies have shown that ultrafine aerosol particles with sizes smaller than 50 nm, which are abundant in the troposphere but are conventionally considered too small to affect cloud formation, can serve as an important source of cloud condensation nuclei (CCN) and promote additional cloud droplets under supersaturation conditions (Fan et al., 2018; Williamson et al., 2019). A key challenge in understanding this phenomenon is that the properties of liquid water in such ultrafine aerosol particles remain incompletely understood due to the difficulty of experimentally characterizing bulk and interfacial water in suspended nanodroplets (Ault and Axson, 2017; Bzdek and Reid, 2017; Ohno et al., 2021).

In the absence of detailed observations of aqueous chemistry in nanodroplets (Bzdek et al., 2020a), extant models rely on a simplified conceptual representation of these particles as a bulk aqueous phase, eventually with an accumulation or depletion of solutes in a negligibly thin interfacial region. In short, water in the core of aqueous aerosol particles is treated as bulk water, and the air–particle interface is represented as the surface of bulk water, even in particles that contain only a few hundred water molecules. This simplification underlies widely used models of aerosol growth thermodynamics and kinetics, organic chemistry, water freezing, salt efflorescence, and organic phase separation (McDonald, 1953; Laaksonen and Malila, 2016; Kulmala et al., 1997; Petters and Petters, 2016; Petters and Kreidenweis, 2007, 2008, 2013; Shiraiwa and Pöschl, 2021; Shiraiwa et al., 2012; Liu et al., 2019; Kostenidou et al., 2018; Semeniuk and Dastoor, 2020). For example, the hygroscopic growth of nanoaerosol particles is generally described within the framework of Köhler theory (McDonald, 1953; Laaksonen and Malila, 2016; Kulmala et al., 1997; Petters and Petters, 2016; Petters and Kreidenweis, 2007, 2008, 2013; Estillore et al., 2017), which accounts for osmotic and capillary pressure effects based solely on the average droplet salinity and surface curvature, neglecting any potential impacts of droplet size on water–water and water–solute molecular-level interactions. Similarly, descriptions of the mass transfer kinetics of water at the particle–gas interface largely rely on the Fuchs–Sutugin approximation with a fixed empirical accommodation coefficient (Shiraiwa and Pöschl, 2021; Shiraiwa et al., 2012), which is assumed invariant with droplet phase state, droplet size, and the diffusivity of the investigated compounds despite evidence that its value can vary by several orders of magnitude depending on experimental conditions (Liu et al., 2019; Kostenidou et al., 2018). Despite its widespread application, the conceptual simplification described above has not been rigorously evaluated in the case of aerosol particles. Analogous simplifications are known to be invalid in small pores and in reverse micelles, where liquid water adopts distinct properties when confined to spaces narrower than 10 nm (Bocquet and Charlaix, 2010; Bourg and Steefel, 2012).

Molecular dynamics (MD) simulations have been extensively used to gain insight into the properties of liquid water near interfaces at the atomistic level. In the case of water–air interfaces, most of these studies have focused on flat surfaces relevant to large aerosol particles (Brown et al., 2005; Thomas et al., 2007; D'Auria and Tobias, 2009; Tang et al., 2020; Luo et al., 2020). Comparatively few have examined the properties of liquid water in sub-10 nm aerosol particles (Chowdhary and Ladanyi, 2009; Li et al., 2011; Ma et al., 2011; Lovrić et al., 2016; Sun et al., 2013; Vardanega and Picaud, 2014; Radola et al., 2017; Karadima et al., 2019; Li et al., 2013; Radola et al., 2015, 2021). These simulations have shown that interfacial water adopts distinct properties from those of bulk liquid water and that the surface tension of water droplets decreases with increasing interfacial curvature at the nanometer scale (Li et al., 2010). In addition, interfacial water molecules have been shown to have highly distinct dynamics compared to bulk liquid water, including slower hydrogen bond dynamics and faster diffusion dynamics (Ni et al., 2013; Liu et al., 2004; Wick and Dang, 2005). Finally, a complicating factor revealed by these simulations and by experimental observations is that the shape of nanodroplets often exhibits significant fluctuations and is rarely a perfect sphere, which complicates the determination of interface thickness and water uptake properties (Li et al., 2016; Garrett et al., 2006).

Here, we use MD simulations to characterize nanodroplet structure, water diffusion, and water evaporation energetics in well-mixed and phase-separated systems with different droplet sizes, NaCl salinities, and pimelic acid (PML) organic surface loadings. Our approach builds on previous MD simulation studies of nanoaerosol droplets but includes two distinct features that have the potential to help resolve the challenges identified in previous studies. First, we use two methods of interface characterization – the time-averaged Gibbs dividing interface (Lau et al., 2015) and the instantaneous interface (Willard and Chandler, 2010) – to compare different perspectives on the properties of the interface. Second, we characterize a broad range of structural, dynamic, and energetic properties to help reveal links between these features. We re-examine the validity of the Kelvin equation and Köhler theory and, more broadly, of treating nanodroplet water as a bulk fluid phase. We find that the instantaneous interface enables a novel characterization of the droplet phase state and helps resolve apparent contradictions between experimental observation and MD simulation results on interfacial water dynamics. We also find that the Kelvin equation and Köhler theory accurately describe water microphysics in nanometer-scale droplets under moderate NaCl salinities or PML organic loadings. However, these theories underestimate water activity at droplet diameters below 4 nm and cannot account for the impact of droplet size on ion concentration enhancement in NaCl-bearing droplets and on phase separation in PML-bearing droplets. Overall, our results indicate that the properties of nanoaerosol water deviate strongly from those of bulk water in particles smaller than 4 nm. One implication of this finding is that molecular-scale insight is required to accurately represent the initial stages of droplet nucleation and growth.

2 Methods

2.1 Simulated systems

Molecular dynamics (MD) simulations of droplets containing 100, 500, 1000, 3000, or 5000 water molecules (N_w) with NaCl molality of 0, 1, 2, or 4 m were used to examine the effect of droplet size and salinity on water microphysics in nanoaerosol particles. In addition, droplets of N_w=100, 500, 1000, 3000, or 5000 with low or high organic surface loadings (0.7 or 3.3 pimelic acid molecules nm⁻², corresponding to surface loadings of 0.3 or 1.5 monolayers) were used to examine the impact of organic compounds. Sodium chloride is used as a generic representative inorganic salt in aerosol particles for two reasons: first, it accounts for a significant fraction of sea spray and mineral dust aerosols, which constitute the largest mass of particulate matter emitted to the atmosphere (Lewis and Schwartz, 2004; Wise et al., 2007), and, second, relative to other inorganic ions present in aerosol particles, such as ${\mathrm{NH}}_{\mathrm{4}}^{+}$ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$, alkali-halide salts such as NaCl have comparatively well-tested and accurate interatomic interaction potential models, the main input controlling the accuracy of MD simulation predictions (Smith and Dang, 1994; Joung and Cheatham, 2008, 2009). Pimelic acid (PML, C₇H₁₂O₄) was used as a generic representative organic component because dicarboxylic acids contribute a large fraction of the total identifiable resolved organic mass in fine aerosols (based on mass spectrometry results that are not strictly quantitative due to differences in ionization efficiencies) and have been identified and quantified in both rural and urban air (Hyder et al., 2012; Kawamura and Yasui, 2005; Tran et al., 2000). Furthermore, the $\mathrm{O}/\mathrm{C}$ ratio of PML (0.57) lies near the midpoint of the range commonly observed for aerosol organic materials (0.2 to 1.0) (Zhang et al., 2007; Hallquist et al., 2009; Song et al., 2018) and near the value below which liquid–liquid phase separation is commonly observed in aerosol particles (∼0.7 to 0.8 for organic–salt–water aerosols, ∼0.6 for organic–water aerosols) (Ciobanu et al., 2009; Song et al., 2012; You et al., 2014; Renbaum-Wolff et al., 2016; Rastak et al., 2017; Song et al., 2018). These observations suggest that NaCl and PML should have the potential to mimic key impacts of salts and organic substances in nanoaerosol droplets. We emphasize, however, that our focus is on identifying and characterizing deviations from ideality in the presence of a generic salt (NaCl) or organic compound (PML) as a function of droplet size and solute loading. Additional studies would be needed to establish the sensitivity of our results to the type of salt or organic compound. The sizes of droplets studied in this work range from 2 to 8 nm (Table A1).

2.2 System preparation and MD simulations

Molecular dynamics simulations were performed with the GROMACS program (Abraham et al., 2015) using cubic simulation cells subject to periodic boundary conditions in all three directions. The edge length of the simulation cell is set to be at least 6 nm larger than the droplet diameter to avoid interactions between the droplet and its periodic images (Fig. A1). Interatomic interactions were represented using the SPC/E model for water molecules (Berendsen et al., 1987), the Joung–Cheatham model for Na⁺ and Cl⁻ ions (Joung and Cheatham, 2008), and the OPLS-AA model for pimelic acid molecules (Jorgensen et al., 1996). In all simulations, bonds involving H atoms were restrained using the LINCS algorithm (Hess, 2008). Simulations were carried out in the canonical (a.k.a. NVT) ensemble with temperature maintained at 298.15 K using the Nosé–Hoover thermostat (Evans and Holian, 1985) and with the overall translational and rotational momentum of the simulated system set to zero. Coulomb and van der Waals interactions were resolved using a real space cutoff of 1.2 nm with a particle mesh Ewald sum treatment of long-range Coulomb interactions. Each system was simulated with a time step of 1 fs for a total simulation time of 105 ns, comprised of 5 ns for equilibration and 100 ns for data analysis. We note that the Ewald sum treatment of long-range Coulomb interactions used in our study is commonly used in MD simulations, including in simulations of aerosol particles, but care should be taken as it introduces artifacts that remain incompletely examined, particularly in the dynamics of charged and dipolar species (Frenkel and Smit, 2001; Hub et al., 2014).

2.3 Droplet structure

To calculate the radius of the water droplets in most systems (except at high organic loading), the time-averaged water number density profile ρ(r) was obtained for each droplet as a function of distance from the droplet center of mass. Results were fitted to the relation (Thompson et al., 1984)

$$\begin{array}{}\text{(1)}& \mathit{\rho }\left(r\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\mathit{\rho }}_{\mathrm{l}}+{\mathit{\rho }}_{\mathrm{v}}\right)-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\mathit{\rho }}_{\mathrm{l}}-{\mathit{\rho }}_{\mathrm{v}}\right)\tanh \left[{\displaystyle \frac{\mathrm{2}}{d_{\mathrm{i}}}}\left(r-r_{\mathrm{0}}\right)\right],\end{array}$$

where ρ_l is the liquid number density, ρ_v is the vapor number density, and d_i is the interfacial width. The density profile was used to locate the equimolar (Gibbs) dividing surface, which corresponds to a vanishing adsorption in a one-component system. The radius R_e of the equimolar dividing surface was calculated, using numerical quadrature, based on the relation

$$\begin{array}{}\text{(2)}& R_{\mathrm{e}}^{\mathrm{3}}={\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{v}}-{\mathit{\rho }}_{\mathrm{l}}}}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{r}^{\mathrm{3}}{\displaystyle \frac{\mathrm{d}\mathit{\rho }\left(r\right)}{\mathrm{d}r}}\mathrm{d}r.\end{array}$$

For droplets with high organic loading, the radius was calculated instead based on the molecular volumes of water and pimelic acid:

$$\begin{array}{}\text{(3)}& R_{\mathrm{e}}^{*}={\left[{\displaystyle \frac{\mathrm{3}}{\mathrm{4}\mathit{\pi }}}\left({\displaystyle \frac{N_{\mathrm{w}}}{{\mathit{\rho }}_{\mathrm{w}}}}+{\displaystyle \frac{N_{\mathrm{org}}}{{\mathit{\rho }}_{\mathrm{org}}}}\right)\right]}^{\mathrm{1}/\mathrm{3}},\end{array}$$

where N_w is the number of water molecules, N_org is the number of PML molecules, ρ_w is the bulk water density (34.02±0.19 molec. nm⁻³), and ρ_org is the bulk PML density (4.373±0.003 molec. nm⁻³) calculated in the conditions of our simulations.

2.4 Diffusion coefficient

The self-diffusion coefficient of water molecules was calculated by analyzing water trajectories at 1 ps intervals during 25 ns simulations for all droplets. To analyze the dependence of diffusion on position relative to the interface, we binned water molecules into different 3 Å thick layers based on their distance relative to the interface. Any molecule initially present in a specified layer that passed through the layer boundaries no longer contributed to the calculation. The self-diffusion of water molecules was investigated in directions parallel and normal to the interface.

The parallel diffusion coefficient (D_∥) was calculated using the Einstein relation (Liu et al., 2004):

$$\begin{array}{}\text{(4)}& D_{\parallel }={\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\begin{array}{c}\lim \\ \mathit{\tau }\to \mathrm{\infty }\end{array}{\displaystyle \frac{\mathrm{d}\langle l_{\parallel ,i}^{\mathrm{2}}\rangle }{\mathrm{d}\mathit{\tau }}},\end{array}$$

where $l_{\parallel ,i}^{\mathrm{2}}$ is the mean square displacement (MSD) of water molecules in directions parallel to the interface, averaged over all molecules of interest remaining in the ith layer during the entire time interval τ. The infinite-time limit in Eq. (4) was approximated by evaluating the slope of MSD versus τ for τ=2 to 6 ps based on the linear relationship between MSD and τ on this timescale; use of longer timescales yielded less precise results as fewer molecules remained consistently located in the ith layer (Fig. A6).

The diffusion coefficient normal to the interface (D_⟂) was evaluated using the anisotropic Smoluchowski equation in the r direction:

$$\begin{array}{}\text{(5)}& {\displaystyle \frac{\partial p}{\partial t}}={\displaystyle \frac{\mathrm{1}}{r^{\mathrm{2}}}}{\displaystyle \frac{\partial }{\partial r}}\left[r^{\mathrm{2}}{D}_{⟂}\left({\displaystyle \frac{\partial p}{\partial r}}+\mathit{\beta }{f}_{\mathrm{r}}p\right)\right],\end{array}$$

where $p=p(r,t|r_{\mathrm{0}},t_{\mathrm{0}})$ is the conditional probability distribution function, $\mathit{\beta }=\frac{\mathrm{1}}{k_{\mathrm{B}}T}$, k_B is Boltzmann's constant, T is absolute temperature, $f_{\mathrm{r}}=-\frac{\partial F\left(r\right)}{\partial r}$ is the mean force acting on the molecule, and $F\left(r\right)=-k_{\mathrm{B}}T\ln \frac{\mathit{\rho }\left(r\right)}{{\mathit{\rho }}_{\mathrm{l}}}$ is the potential of mean force (PMF) calculated from the water density profile. Values of D_⟂ in different water layers were determined by matching the water survival probabilities calculated from MD simulation and from the numerical solution to the Smoluchowski equation as described by Wick and Dang (2005).

2.5 Surface tension

Surface tension was calculated by analyzing the normal component of the Irving–Kirkwood pressure tensor (Irving and Kirkwood, 1950) at 1 ps intervals during 25 ns simulations for all droplets. According to the formula given by Thompson et al. (1984), the normal pressure tensor P_N(r) can be written as

$$\begin{array}{}\text{(6)}& P_{\mathrm{N}}\left(r\right)=P_{\mathrm{K}}\left(r\right)+P_{\mathrm{U}}\left(r\right)=k_{\mathrm{B}}T\mathit{\rho }\left(r\right)+S^{-\mathrm{1}}\sum _k{f}_k,\end{array}$$

where P_K(r)=k_BTρ(r) is the kinetic term, $P_{\mathrm{U}}\left(r\right)=S^{-\mathrm{1}}\sum _k{f}_k$ is the configuration term, S is the surface area of the spherical surface of radius r, and f_k is the normal component of the force between a pair of molecules acting across the surface S. The surface tension is evaluated as

$$\begin{array}{}\text{(7)}& \mathit{\sigma }\approx {\displaystyle \frac{\mathrm{3}W}{\mathrm{4}\mathit{\pi }{R}_{\mathrm{e}}^{\mathrm{2}}}},\end{array}$$

where W is the work of formation as described by Zakharov et al. (1997):

$$\begin{array}{}\text{(8)}& W=\mathrm{2}\mathit{\pi }\underset{\mathrm{0}}{\overset{R_{\mathit{\beta }}}{\int }}{P}_{\mathrm{N}}\left(r\right)r^{\mathrm{2}}\mathrm{d}r-\mathrm{2}\mathit{\pi }{P}^{\mathit{\beta }}{R}_{\mathit{\beta }}^{\mathrm{3}}/\mathrm{3}.\end{array}$$

In Eq. (8), P^β is the vapor pressure of the droplet, and R^β is the radius of a sphere in the vapor region, which can be approximated as the equimolar radius R_e. We note that the surface tension should be rigorously calculated at the radius R_s, where the surface tension applies (McGraw and Laaksonen, 1997). Our use of R_e in Eq. (7) introduces a systematic error of less than 8 % in the σ calculation for droplets with N_w≥500 (see Appendix A for more details). This potential source of systematic error is commensurate with the statistical error of our predicted σ values as shown by our reported error bars of $\sim \pm \mathrm{8}$ %.

2.6 Free energy of evaporation

The potential of mean force (PMF) associated with the evaporation of water molecules as a function of distance r from the droplet center F(r) was calculated using two methods. The first method consisted in calculating the PMF based on the water density profile during 100 ns equilibrium simulations by

$$\begin{array}{}\text{(9)}& F\left(r\right)=-k_{\mathrm{B}}T\ln {\displaystyle \frac{\mathit{\rho }\left(r\right)}{{\mathit{\rho }}_{\mathrm{l}}}}.\end{array}$$

The second method is the well-known umbrella sampling technique (Torrie and Valleau, 1977), whereby F(r) is calculated by successively tethering a single water molecule at different distances from the droplet center of mass using a harmonic potential, here with a force constant of 1000 kJ mol⁻¹ nm⁻¹ and 100 ns of simulation time at each tethering distance. The evaporation free energy of water molecules was calculated as the free energy difference between a water molecule in the vapor phase and in the liquid phase as

$$\begin{array}{}\text{(10)}& \mathit{\delta }F=F_{\mathrm{v}}-F_{\mathrm{l}}.\end{array}$$

2.7 Accommodation coefficient

To calculate the accommodation coefficient for a single droplet, the phase state of every water molecule in the system was recorded during a 100 ns trajectory with a time step of 1 fs. The phase state of individual water molecules was divided into three categories (Varilly and Chandler, 2013): (1) liquid phase (state −1) – any water molecule that is hydrogen-bonded to at least one water molecule and is at most 4 Å away from the nearest water molecules (where the “position” of a water molecule refers to the position of its center of mass), (2) vapor phase (state 1) – any water molecule that has no hydrogen-bonded water neighbors and is more than 8 Å away from its nearest water neighbor, and (3) intermediate phase (state 0) – any water molecule that is neither in state 1 nor in state −1. An accommodation event is identified if the state of a water molecule changes from 1 at time step t₁ to −1 at time step t₂, without adopting a non-zero state value during the time interval (t₁, t₂). A reflection event is identified if the state of a single water molecule changes from 1 at t₁ to 1 at t₂ and stays in state 0 during the time interval (t₁, t₂). The accommodation coefficient α was calculated as

$$\begin{array}{}\text{(11)}& \mathit{\alpha }={\displaystyle \frac{N_{\mathrm{acc}}}{N_{\mathrm{acc}}+N_{\mathrm{ref}}}},\end{array}$$

where N_acc is the number of accommodation events and N_ref is the number of reflection events. The reported normalized accommodation coefficient α^* was calculated as ${\mathit{\alpha }}^{*}=\mathit{\alpha }/{\mathit{\alpha }}_{\mathrm{0},\mathrm{w}}$, where α_0,w (the value for the flat surface of pure water) was approximated as the α value of the largest simulated pure-water droplet of N_w=5000.

3 Results and discussions

3.1 Morphology and phase state characterization

Figure 1a–e show representative snapshots of simulated droplets that illustrate the variety of predicted morphologies and phase-mixing states. Additional droplet snapshots are presented in Fig. A3. For NaCl droplets, Na⁺ and Cl⁻ ions are repelled from the water–air interface and distributed inside the droplets in all simulations. These ions are crystallized in the center of the droplet at a 4 m concentration, forming a phase state illustrated in Fig. 1e, whereas they are well mixed with water molecules at 1 and 2 m concentrations with the phase state illustrated in Fig. 1a. The NaCl solubility implied by our results (between 2 and 4 m) is consistent with that observed in previous MD studies that examined bulk NaCl solutions using the same interatomic potential (Mester and Panagiotopoulos, 2015a, b). For droplets with PML, the arrangement of organic molecules varies with droplet size and organic surface loading. At low organic loading, PML molecules are predominantly adsorbed at the water–air interface with minimal clustering, as illustrated in Fig. 1b. At high organic loading, when droplet size is sufficiently large (diameter D≥4.1 nm), PML molecules retain their affinity for the water–air interface and cluster at the droplet surface, forming the “engulfed” morphology shown in Fig. 1c. In smaller droplets with high organic loading, however, the organic molecules cluster in the center of the droplet and displace water molecules towards the surface as shown in Fig. 1d. This transition suggests the existence of a threshold, at a droplet diameter of 3 to 4 nm, in the competition between water hydrogen bonding, which promotes water clustering, and hydrophobic attraction between PML molecules, which promotes organic clustering. Our results are consistent with previous findings that the hydrogen bond network formed by bulk liquid water involves cooperativity over length scales of 1 to 2 nm and becomes significantly weaker upon attenuation of long-range water–water orientational correlations (Alper and Levy, 1989; Ohmine and Tanaka, 1993).

Figure 1 Sphericity factor (ϕ) of nanodroplets with different compositions as a function of droplet diameter. Simulation snapshots are shown for the droplets with water (blue color), Na⁺ (light-orange spheres), Cl⁻ (green spheres), and pimelic acid (colored chains): (a) 1 m NaCl with N_w=1000, (b) low PML surface loading with N_w=1000, (c) high PML surface loading with N_w=5000, (d) high PML surface loading with N_w=100, and (e) 4 m NaCl with N_w=3000.

To characterize the various phase states identified above, we propose a quantitative variable, referred to hereafter as the sphericity factor (ϕ), defined as the ratio of the surface area of a perfect spherical water droplet S_w,sp to the instantaneous interface surface area formed by water molecules with the same amount of water S_w,inst:

$$\begin{array}{}\text{(12)}& \mathit{\varphi }={\displaystyle \frac{S_{\mathrm{w},\mathrm{sp}}}{S_{\mathrm{w},\mathrm{inst}}}}.\end{array}$$

The ϕ value is in the range of (0, 1] and reflects the tendency towards water clustering. In particular, ϕ=1 indicates that water forms a perfect spherical droplet with no capillary waves or disruption due to phase separation.

Figure 1 presents the ϕ values of all simulated systems as a function of droplet diameter D. By visual inspection of our simulation trajectories, we identify droplets as being in a well-mixed state if no single phase is formed except for the water cluster and in a phase-separated state if an additional phase exits. For the conditions examined in this study, a threshold at ϕ≈0.84 (illustrated by the horizontal gray dashed line in Fig. 1) distinguishes between well-mixed (above the threshold) and phase-separated droplets (below the threshold). For well-mixed systems, at a fixed ion concentration or organic surface excess, ϕ decreases slightly with increasing droplet size, indicating that larger droplets have a more fluctuating water–air interface. This observation is consistent with our results showing that the water–air interface width observed in time-averaged density profiles increases with droplet size (Table A1) and agrees with the expectation that the capillary wave contribution to interfacial fluctuation increases with interfacial area (Ismail et al., 2006; Lau et al., 2015). We also find that ϕ increases with ion concentration but decreases with organics at the interface, which indicates that NaCl reduces but PML enhances interfacial fluctuations. This observation is consistent with the variations in interfacial width with salinity and organic surface loading (Table A1) and with the expectation that NaCl and PML should, respectively, increase or decrease surface tension because of their negative (NaCl) or positive (PML) adsorption at the water–air interface, as confirmed below.

Phase-separated systems (i.e., droplets with either 4 m NaCl or high PML loading) also show a small decrease in ϕ with increasing droplet size in most cases. A significant exception is presented by the smallest droplet with high organic loading, which shows a sharp decrease in ϕ with decreasing droplet size between D=4.1 and 2.7 nm. This decrease in ϕ coincides with a change in droplet phase-mixing state from a water-core structure at D≥4.1 nm to an organic-core structure at D=2.7 nm and shows that the phase-mixing state of nanodroplets can abruptly vary with particle size. Organic particles with sizes below 5 nm and different PML–water ratios, which lie in the abovementioned phase-mixing transition region and may shed light on the early stages of hygroscopic growth of secondary organic aerosols (Qin et al., 2021), will be examined in more detail in an upcoming study.

3.2 Interfacial structure of nanodroplets

To gain additional insight into the phase state behaviors outlined above, the radial density profiles of water, ions, and PML molecules in our nanodroplets were calculated relative to either the time-averaged interface (Gibbs dividing surface) (Lau et al., 2015) or the instantaneous interface (Willard and Chandler, 2010). As shown in previous studies of flat water–air interfaces (Willard and Chandler, 2010) and in Fig. A4, atomic density profiles obtained under the instantaneous interface scheme provide more detailed information about molecular structure at the interface (e.g., water and solute density layering). Figure 2a shows ion density relative to the instantaneous interface in 1 m NaCl droplets with different numbers of water molecules. Our results show that both Na⁺ and Cl⁻ ions are negatively adsorbed at the water–air interface, but Na⁺ is more strongly excluded from the interface than Cl⁻ such that the interface is negatively charged, in agreement with previous studies of flat water–air interfaces (Jungwirth and Tobias, 2000; Underwood and Greenwell, 2018) (see Fig. A4 for more density profiles). Although all droplets highlighted in Fig. 2a were initialized as 1 m uniformly mixed NaCl solutions, ion concentrations are strongly enhanced (up to >2 times the initial concentration) in the core of the droplets at equilibrium due to the negative adsorption of ions at the water–air interface. This ion concentration enhancement increases with decreasing droplet size, particularly in droplets smaller than 4 nm (N_w<1000). We note that the interface propensity of halide anions varies with anion size and polarizability and that the simple interatomic potential models used here may underestimate the affinity of halide anions for the interface (Levin et al., 2009; Jungwirth and Tobias, 2006; Caleman et al., 2011). However, the negative adsorption of alkali cations predicted here is well established (Underwood and Greenwell, 2018; Levin et al., 2009; Jungwirth and Tobias, 2006; Caleman et al., 2011), and the overall negative adsorption of NaCl salt at the water–air interface is correctly predicted as shown by previous MD simulation results on the salinity dependence of water–air interfacial tension (Jungwirth and Tobias, 2001, 2006; D'Auria and Tobias, 2009; Horinek et al., 2009; Sun et al., 2012).

Figure 2 (a) Density profiles of Na⁺ and Cl⁻ ions (ρ_ion) relative to the instantaneous interface in 1 m NaCl droplets with different numbers of water molecules. (b) Salt concentration enhancement factor ϵ as a function of droplet diameter D for droplets with an average salinity of 1 m (blue circles: MD simulation results; lines: predictions for different salts based on a simplified model, whereby the salt exclusion distance at the water–air interface equals 1.2 times the first hydration radius of the cation; horizontal dashed lines: ϵ values where different salts would be expected to precipitate, based on their solubility in bulk liquid water, in a droplet with an average salinity of 1 m). (c) Deliquescence relative humidity (DRH) of KCl and NaCl nanoparticles as a function of diameter (D): shaded area marks model predictions obtained in this study; blue and orange diamonds are previous experimental results for KCl and NaCl (Biskos et al., 2006a; Giamarelou et al., 2018); orange triangles are previous MD simulation results for NaCl (Bahadur and Russell, 2008); and the orange dashed line is the NaCl DRH prediction based on a thermodynamic model (Russell and Ming, 2002). (d) Atomic density profiles of water molecules (ρ_w) and PML molecules (ρ_org) relative to the Gibbs dividing surface under the time-averaged interface scheme at low organic loading.

To illustrate the expected effect of droplet size on ion concentration, we define a salt concentration enhancement factor (ϵ) as the ratio between salt concentration in the core of the droplet and average salt concentration in the entire droplet. Predicted ϵ values for different alkali-halide salts are shown in Fig. 2b based on the simplifying assumption that the salt is excluded from the vicinity of the water surface, in a region with a thickness equal to 1.2 times the first hydration radius of alkali cations (see Appendix B for sensitivity analysis), quantified as the first minimum of the ion–water radial distribution function predicted using the SPC/E water model and compatible ion models (the first hydration shell radii of alkali cations Li⁺, Na⁺, and K⁺ are 0.28, 0.32, and 0.36 nm, respectively) (Joung and Cheatham, 2009). Our simplistic model, shown as solid curves in Fig. 2b for different salts, is consistent with our MD simulation predictions (blue circles, shown for our systems with 1 m NaCl), indicating that the length scale associated with salt exclusion at the water–air interface is invariant with surface curvature. Horizontal lines in Fig. 2b show the ϵ values necessary to achieve crystallization of different salts based on their experimental solubility (Joung and Cheatham, 2009) in droplets where the average concentration is 1 m. Overall, both our simplistic model and our detailed atomistic simulation predictions indicate that the concentration enhancement of simple alkali-halide salts is well above 1.25 for sub-10 nm droplets and may be more than 10 for droplets with a diameter below 1.5 nm. In other words, in droplets with common sea-salt ionic species, crystallization should occur more readily than expected based on the average salt concentration, particularly at small droplet sizes. Our results suggest that the tendency for salt deliquescence should be diminished and salt efflorescence should be enhanced in ultrafine nanoparticles. The ion concentration enhancement discussed here is consistent with experimental findings that both the deliquescence and efflorescence relative humidity (DRH, ERH) of NaCl and KCl particles increase as droplet size decreases from 60 nm to less 10 nm, with a sharper decrease at smaller droplet sizes (Biskos et al., 2006a, b; Giamarelou et al., 2018; Hämeri et al., 2001). Shaded areas in Fig. 2c show the KCl and NaCl DRH as a function of nanoparticle diameter predicted from the ion concentration enhancement model discussed above (see Appendix C for calculation details). Our predictions are consistent with both the size dependence of NaCl DRH reported in previous studies (Biskos et al., 2006a; Bahadur and Russell, 2008; Russell and Ming, 2002) and previous observations that KCl has a larger DRH than NaCl (Giamarelou et al., 2018). Overall, the concentration enhancement discussed above may help explain experimental observations of a size dependence of aerosol shape, in which nanoaerosol particles generated from well-controlled sea spray chambers are more cubic at smaller sizes but more spherical at larger sizes (Zieger et al., 2017). Our observations indicate that the ϵ factor becomes particularly significant for sub-10 nm particles and should be taken into consideration in predictions of salt microphyics in ultrafine particles.

Figure 2d shows the density profiles of water and PML molecules based on the location of their center of mass under the time-averaged interface scheme at low organic surface loading. The PML density profiles show a sharp peak near the water–air interface, confirming the high affinity of the organic compound for the interface. The location of the PML density peak shifts into the liquid water region with decreasing droplet size, particularly at sizes below 3 nm (N_w<500). This observation suggests that the hydrophilic character of the water phase decreases with decreasing number of water molecules in the droplet, as already noted above in the case of the smallest droplet at high organic loading. The PML density curves at high organic loading (Fig. A5d) reveal a more complex and variable pattern of organic–water phase-mixing states, in agreement with the transition between water- and organic-core morphologies identified above.

3.3 Water diffusion kinetics in nanodroplets

Water diffusion in nanodroplets is a fundamental phenomenon that influences interfacial mass fluxes and chemical transformation rates (Chan et al., 2014; Slade and Knopf, 2014; Julin et al., 2013; Evoy et al., 2020). Details of this phenomenon are complicated by the inherently anisotropic and non-uniform nature of water diffusion near the interface. In this study, the water self-diffusion coefficient is calculated in discrete 3 Å thick water layers in directions parallel (D_∥) or normal (D_⟂) to the interface. As explained in the Methods section, we calculated D_∥ using the Einstein relation and D_⟂ using the anisotropic Smoluchowski equation (Liu et al., 2004). Diffusion coefficients were calculated relative either to the center of mass of the droplet (“com” scheme) or to the instantaneous interface (“inst” scheme). As expected, both schemes yield identical results for D_∥, while for D_⟂ the com scheme yields larger values due to the presence of capillary waves at the interface (Appendix D2.1). The following discussion focuses primarily on our predicted D_∥ values. Our D_⟂ values are somewhat less precise and are discussed in Appendix D2.2.

Figure 3a shows D_∥ in different simulated systems relative to the distance from the Gibbs dividing surface. Our results show that D_∥ values are bulk-liquid-like at distances > 9 Å (i.e., >3 water layers) from the interface for all well-mixed systems. The D_∥ values at the interface are 2 to 3 times larger than in the core of the droplet, with potential implications for mass transfer and reaction kinetics at the interface. In addition, D_∥ in the droplet core decreases from 0.25 to 0.15 Ų ps⁻¹ as ion concentration increases from 0 to 2 m, while D_∥ at the interface decreases with increasing PML surface loading. The D_∥ difference between interfacial and bulk liquid regions and the D_∥ change with salinity and organic addition observed in our results are consistent with trends reported in previous studies (Kim et al., 2012; Paul and Chandra, 2005; Tsimpanogiannis et al., 2019) as shown in Fig. 3a. In all well-mixed NaCl droplets (1 or 2 m), we consistently observe a droplet size dependence of water diffusion near the interface as illustrated in Fig. 3b. More precisely, D_∥ at the interface decreases with droplet size, while D_∥ in the few monolayers below the interface increases with droplet size. This observation is consistent with our results showing that for droplets containing water and NaCl, smaller droplets have lower ion density at the interface and greater ion density in the droplet core. The minor offset in Fig. 3a, between our D_∥ values in the bulk-liquid-like region and those reported by Kim et al. (2012) at the same salinity (calculated in bulk liquid NaCl solutions), is also consistent with the ion density enhancement in the droplet core noted above. Finally, the faster water diffusivity at the interface observed in the case of D_∥ is not observed in the case of D_⟂, as discussed in Appendix D2. In fact, results obtained under the inst scheme yield $D_{⟂}\approx \mathrm{0.15}$ to 0.17 Ų ps⁻¹ near the interface; i.e., water diffusivity normal to the interface is smaller than in the droplet core, in agreement with studies showing slower reorientational dynamics of water near hydrophobic surfaces (Laage et al., 2009; Fayer and Levinger, 2010; Fayer, 2012).

Figure 3 Water diffusion coefficient in the direction parallel to the interface (D_∥) as a function of distance from the time-averaged interface (com scheme) (a) in 5.5 nm droplets (N_w=3000) with different aqueous chemistries and (b) in 1 m NaCl droplets with different water contents (symbol sizes are scaled with droplet diameter). Literature data shown in panel (a) were collected from previous MD simulation studies using the SPC/E water model at 0 m (Kim et al., 2012; Paul and Chandra, 2005; Tsimpanogiannis et al., 2019) and 1 and 2 m (Kim et al., 2012) and for a methanol–water solution at low organic concentrations (Paul and Chandra, 2005). All literature D_∥ values were obtained for systems with a flat water–air interface.

3.4 Surface tension of nanodroplets

Surface tension is an important factor that governs the hygroscopic growth and the activation of aerosol particles and is incorporated in climate models through Köhler theory. However, the surface tension of nanoaerosol particles, especially at the sub-10 nm scale, is not well known and remains arduous to characterize experimentally (Noziere et al., 2014; Bzdek et al., 2016, 2020b; Lee and Tivanski, 2021). Figure A12 shows surface tensions (σ) calculated for our different droplets based on the Irving–Kirkwood pressure tensor (Irving and Kirkwood, 1950; Thompson et al., 1984). A representative example of the instantaneous normal pressure tensor is shown in Fig. A11. All calculated σ values are also tabulated in Table A1. We note that σ values reported in different MD simulation studies differ significantly, even for a pure-liquid-water system (by up to ≈20 mN m⁻¹), depending on the interatomic potential models and surface tension determination methods employed (Bahadur et al., 2007; Bahadur and Russell, 2008; Li et al., 2010; Lau et al., 2015). For well-mixed NaCl droplets with 0, 1, and 2 m NaCl concentrations, our results show that surface tension increases with droplet size and ion concentration and approaches 55, 56, and 58 mN m⁻¹, respectively, when N_w≥3000 (i.e., D>5 nm). The predicted increase in surface tension with NaCl salinity is consistent with experimental values reported for flat water–air interfaces (1.3 to 1.8 mN m⁻¹ per 1 m of salinity) (Hård and Johansson, 1977; Washburn et al., 1926; Abramzon and Gaukhberg, 1993; Tuckermann, 2007). To highlight the impacts of curvature, salinity, and organic coatings on σ and compare our results with those of other studies, we report normalized surface tension values $\mathit{\sigma }/{\mathit{\sigma }}_{\mathrm{0}}$ in Fig. 4, where σ₀ is the surface tension reported for the flat pure-liquid-water–air interface in the same study. Our σ₀ value is 58 mN m⁻¹, in agreement with previous estimates of the surface tension of SPC/E water (Ismail et al., 2006). The solid light-blue line in Fig. 4 shows the fit of our $\mathit{\sigma }/{\mathit{\sigma }}_{\mathrm{0}}$ results for 0 m droplets to the second-order Helfrich expansion, $\mathit{\sigma }/{\mathit{\sigma }}_{\mathrm{0}}=\mathrm{1}-\mathrm{2}{\mathit{\delta }}_{\mathrm{0}}/r+c/r^{\mathrm{2}}$, where δ₀ is the Tolman length, r is the droplet radius, and c is the rigidity coefficient (Wilhelmsen et al., 2015). The best fit is achieved using δ₀=0.06 nm (95 % confidence interval [0.032, 0.094] nm) and $c=-\mathrm{0.025}$ nm² (95 % confidence interval [−0.096, 0.045] nm²). The positive value of δ₀ in this study, which indicates that σ increases with droplet size, agrees with the trends observed in other studies as shown in Fig. 4. In addition, although the absolute σ values reported in different studies differ depending on the choice of simulation methods, $\mathit{\sigma }/{\mathit{\sigma }}_{\mathrm{0}}$ values from previous studies (Bahadur and Russell, 2008; Li et al., 2010; Lau et al., 2015) agree well with our simulation results and collapse near the light-blue fitting curve. This observation shows that, for pure-water droplets, σ decreases sharply as droplet size decreases below 4 nm and can even be less than half of σ₀ at sizes below 1 nm. For droplets with PML organic surface loading, surface tension not only shows the same trend with particle size as for pure-water droplets but also decreases with organic surface loading for all systems where the droplet has a water core (N_w≥500), in agreement with experimental observations that water surface tension decreases with organic adsorption at the water–air interface (Tuckermann, 2007; Kiss et al., 2005). An exception is observed in Fig. 4 in the case of the smallest droplets with low and high PML loadings (leftmost yellow and orange circles in Fig. 4), which both contain the same amount of water (N_w=100), yet the droplet with high organic loading has a higher surface tension. The surface tension of the organic-core configuration (i.e., the smallest droplet with high organic loading) is close to that observed for droplets of pure PML (on the order of 30 mN m⁻¹ at D≈3 nm), suggesting that the unexpected increase in surface tension with organic loading is associated with the observed morphology reversal from a water core to an organic core. In short, the common assumption in aerosol microphysics studies that surface tension invariably decreases with increasing organic loading in aqueous aerosols appears invalid, at least in some cases, for droplets smaller than 4 nm. This finding may have important implications, in particular, in understanding the initial nucleation and growth of secondary organic aerosols (Enghoff and Svensmark, 2008; Kürten et al., 2018; Guo et al., 2020).

Figure 4 Normalized surface tension ($\mathit{\sigma }/{\mathit{\sigma }}_{\mathrm{0}}$) of droplets as a function of droplet diameter (D). Colored circles show results obtained in this study using the SPC/E water model and the Irving–Kirkwood (I–K) pressure tensor method; symbols with black edges show previous MD simulation results obtained for pure-water droplets at ambient temperature using the TIP4P/2005 or SPC/E water models and the test-area (TA), energy difference (ED), or Irving–Kirkwood (I–K) methods (Bahadur and Russell, 2008; Li et al., 2010; Lau et al., 2015); the light-blue solid line shows a fit of our 0 m surface tension data to the second-order Helfrich expansion, yielding the Tolman length δ₀=0.06 nm.

3.5 Validity of the Kelvin equation

The evaporation free energy δF, i.e., the free energy difference associated with transferring a single water molecule from the droplet phase to the vapor phase, is an important factor that governs water evaporation kinetics through the relation (Seinfeld and Pandis, 2016)

$$\begin{array}{}\text{(13)}& j\approx c\exp \left(-{\displaystyle \frac{\mathit{\delta }F}{RT}}\right),\end{array}$$

where j is the evaporation rate and c is a constant factor. From the classical thermodynamic Kelvin theory, the water evaporation free energy for a droplet of radius r (δF(r)) relates to the value for a bulk liquid system with a flat interface and with the same aqueous chemistry δF_b through the relation

$$\begin{array}{}\text{(14)}& \mathit{\delta }F\left(r\right)-\mathit{\delta }{F}_{\mathrm{b}}=-RT\ln {\displaystyle \frac{P_{\mathrm{v}}\left(r\right)}{P_{\mathrm{v},\mathrm{b}}}}=-{\displaystyle \frac{\mathrm{2}\mathit{\sigma }{V}_{\mathrm{m}}}{r}},\end{array}$$

where P_v(r) is the saturated vapor pressure of the droplet, P_v,b is the saturated vapor pressure of the bulk liquid system with a flat interface and the same aqueous chemistry, σ is the surface tension of the droplet, and V_m is the molar volume of the liquid. Therefore, if the Kelvin equation is valid and if solute concentration is invariant with r, a plot of δF as a function of $\mathrm{2}\mathit{\sigma }/D$ (where D=2r) should yield a straight line with a slope of −2V_m. Such a plot is presented in Fig. 5a.

Figure 5 (a) Water evaporation free energy calculated from the density profile in equilibrium MD simulations (“Dens”) or from the umbrella sampling method (“Um”): the black line shows the reference slope of δF versus $\mathrm{2}\mathit{\sigma }/D$ from the classical Kelvin equation prediction; the blue line shows the best linear fit obtained for well-mixed NaCl systems with droplet sizes larger than 3.5 nm (N_w≥1000). Marker size is scaled with droplet size. Literature data were collected from previous MD simulations carried out using the SPC/E water model (Garrett et al., 2006; Li et al., 2018; Varilly and Chandler, 2013). (b) Water activity a_w (circles) versus $\mathrm{1}/D$ for different droplets: colored lines show best fits obtained using Köhler theory (Eq. 15).

For each simulated system, we calculated δF using two different approaches. The first method consisted in determining the free energy of water evaporation using the umbrella sampling method (Torrie and Valleau, 1977), a biased MD simulation method in which a single water molecule is progressively steered from the particle to the vapor phase. The second method consisted in determining the equilibrium water densities in the liquid and vapor phases, ρ_l and ρ_v, from which the evaporation free energy was obtained as $\mathit{\delta }F=-RT\ln {\mathit{\rho }}_{\mathrm{v}}/{\mathit{\rho }}_{\mathrm{l}}$. The second method was applied only for systems where a well-defined aqueous phase (and corresponding ρ_l value) could be identified. From Fig. 5a, we observe that the umbrella sampling and density profile methods yield consistent δF values for all well-mixed systems (0, 1, and 2 m NaCl or low PML loading). We note that there are limited δF literature data, even for a pure-water system. Even for MD simulations of the flat water–air interface carried out with the same water model as in the present study (SPC/E), δF values reported in previous studies vary significantly, as shown in Fig. 5a (Garrett et al., 2006; Li et al., 2018; Varilly and Chandler, 2013). The comprehensive δF dataset reported in our study, with longer simulation times, multiple analysis approaches, and a range of droplet sizes and water chemistry conditions, should provide a useful guide for future examinations of the evaporation free energy of water molecules from aerosol particles.

Calculated δF values in our study increase with droplet size (i.e., water evaporates more readily from smaller droplets) as expected from Eq. (14). Although the impact of ion concentration on evaporation energy cannot be precisely resolved, a linear relationship between δF and $\mathrm{2}\mathit{\sigma }/D$ is clearly observed in the case of well-mixed droplets, in agreement with the Kelvin equation. The theoretical slope for pure-water droplets according to the Kelvin equation prediction ($k=-\mathrm{36}\times {\mathrm{10}}^{-\mathrm{6}}$ m³ mol⁻¹, calculated using the molar volume of liquid water predicted by the SPC/E water model) is indicated by the black line in Fig. 5a. This theoretical slope lies within the 95 % confidence interval of the slope obtained by linear regression of the data for well-mixed NaCl droplets with sizes larger than 3.5 nm (N_w≥1000) (blue solid line in Fig. 5a with slope $k_{\mathrm{f}}=-\mathrm{41}\pm \mathrm{13}\times {\mathrm{10}}^{-\mathrm{6}}$ m³ mol⁻¹). However, if the evaporation data obtained for droplets smaller than 3 nm (N_w≤500) are included in the regression, the best fitting line is deflected downward and k falls outside of the 95 % confidence interval of k_f. Our observation that the Kelvin equation is valid at droplet sizes larger than 4 nm under moderate salinities has also been revealed by experimental studies (Winkler et al., 2012).

For droplets with low PML loading, the trend shown in Fig. 5a is generally consistent with that observed for well-mixed water and NaCl droplets. However, the smallest low-organic-loading droplet represents a clear outlier in Fig. 5a, as also observed in Fig. 4, suggesting a possible link between its low surface tension and its low δF value. For the 4 m NaCl droplet system, the trend in δF (dark-purple circles) is similar to that observed for the well-mixed NaCl systems, suggesting that although the ions in the center of the droplet form a solid, the water surface is similar to that of the well-mixed ion–water systems. In contrast, for the phase-separated high-PML-loading system (orange circles), predicted δF values are consistently much smaller than those of the well-mixed systems. This indicates that at high PML loading, the organic coating strongly enhances water evaporation and invalidates Eq. (14). For the smallest droplet with high PML loading, δF is only 24 kJ mol⁻¹, deviating significantly even from the values obtained for larger droplets with high PML surface loading, an observation that may be related to the reversal from water-core to organic-core morphology identified above.

To conclude, the Kelvin equation accurately predicts water evaporation energetics for salty NaCl droplets with sizes larger than 3.5 nm, even at salinities up to 4 m where salt precipitation is observed. However, it overestimates the stability of droplets smaller than 3 nm and of droplets with high PML organic surface loading, in agreement with the expectation that water becomes a somewhat distinct (less hydrophilic) fluid when the collective nature of the water hydrogen bond network is disrupted.

3.6 Validity of Köhler theory

Water activity a_w, the ratio of the partial pressure of water in solution and the saturated water vapor pressure of bulk pure water, is an important factor that governs the hygroscopic growth and mass concentration of aerosol particles at a given relative humidity (Seinfeld and Pandis, 2016). According to Köhler theory (Seinfeld and Pandis, 2016), the water activity of an aerosol droplet can be expressed as

$$\begin{array}{}\text{(15)}& a_{\mathrm{w}}=P_{\mathrm{v}}\left(r\right)/P_{\mathrm{v},\mathrm{0}}=x_{\mathrm{w}}\exp {\displaystyle \frac{\mathrm{2}\mathit{\sigma }{V}_{\mathrm{m}}}{RTr}},\end{array}$$

where x_w is the effective mole fraction of water molecules in the solution, which accounts for the effect of solutes as represented using Raoult's law, and P_v,0 is the saturated vapor pressure of a pure-water system with a flat interface.

To evaluate Köhler theory at the nanometer scale, we calculated a_w in all simulated droplet systems and compared our results with Eq. (15) as shown in Fig. 5b. For each system, we calculated the water activity as $a_{\mathrm{w}}={\mathit{\rho }}_{\mathrm{v}}/{\mathit{\rho }}_{\mathrm{v},\mathrm{0}}$, where ρ_v,0 is the density of water in the vapor phase in a system with pure water and a flat water–air interface. The theoretical prediction was determined by fitting a_w values obtained for droplets with sizes larger than 3.5 nm using Eq. (15), with σ and V_m obtained from simulations of systems with bulk-liquid-like water and a flat water–air interface (with the same salinity or PML surface loading) and x_w fitted to minimize the root mean square deviation from simulation results. From Fig. 5b, we see that a_w increases with decreasing droplet size. For pure-water droplets, within the precision of our results, the simulation data match the theoretical prediction well for droplets larger than 3.5 nm in diameter. For droplets smaller than 3 nm, water activity is underestimated by Köhler theory, indicating that small nanodroplets are less stable than expected (i.e., they exhibit higher saturated vapor pressure than expected, as noted above). For 2 m NaCl droplets, within the precision of our results, the simulation results are consistent with Eq. (15). However, we note that under the ideal-solution assumption that underlies Köhler theory, a 2 m NaCl droplet should have x_w=0.93, whereas our results are consistent with Eq. (15) with an effective water mole fraction of $x_{\mathrm{w}}=\mathrm{0.86}\pm \mathrm{0.02}$, which indicates that the impact of NaCl on the activity of water is roughly twice as large as predicted by Raoult's law. Since the interatomic potential models used in this study accurately represent experimental water activity in bulk liquid water up to 3 m NaCl salinity (Mester and Panagiotopoulos, 2015a), this observation likely relates to the presence of the water–air interface. The magnitude of the observed deviation is consistent with the ion concentration enhancement effect identified above. For organic aerosols, as the mixing state of water and organic matter is often poorly defined, the effective water mole fraction is often fitted to match measured values of a_w (Seinfeld and Pandis, 2016). From our fitting results, the effective x_w value in our low-PML-loading system equals 0.996±0.006, indicating that at low surface coverage PML has essentially no impact on water activity. In contrast, in our high-PML-loading systems, the effective x_w value equals 0.890±0.011, indicating that PML can lower water activity by 10 % in engulfed structure droplets. Overall, our results are consistent with Köhler theory for droplets with sizes larger than 3.5 nm. However, Köhler theory overestimates the stability of water in droplets smaller than 3 nm (a phenomenon that may be related to the lower surface tension and water cohesion in small droplets identified above) and underestimates the impact of ions on water activity (a phenomenon that may be related to the ion concentration enhancement noted above).

3.7 Accommodation coefficient of water and condensation kinetics

The mass accommodation coefficient of water molecules on nanoaerosol particles, defined as the probability that a water vapor molecule that impinges on the particle surface becomes incorporated in the particle, is a fundamental factor that governs the hygroscopic growth kinetics of atmospheric aerosols (Diveky et al., 2021; Seinfeld and Pandis, 2016). Its precise value in atmospheric conditions remains highly uncertain (within a range of 0.01 to 1) despite decades of study because of its complex dependence on particle composition, phase-mixing state, and surface curvature (Diveky et al., 2021; Shiraiwa and Pöschl, 2021; Barclay and Lukes, 2019; Davis, 2006; Winkler et al., 2004). In our study, the water accommodation coefficient (α) was calculated by tracking the trajectory and phase state of water molecules for 100 ns for each simulated system. We note that even based on atomistic trajectories determined using MD simulations, there are several possible ways of defining whether an impinging water molecule is absorbed or reflected, and, as a result, different α values (ranging from 0.9 to 1.0) are reported even for the flat surface of pure liquid water (Julin et al., 2013; Barclay and Lukes, 2019; Skarbalius et al., 2021; Tsuruta and Nagayama, 2004; Ohashi et al., 2020). To highlight the impacts of curvature, salinity, and organic coatings on α and compare our results with those of other studies, we report normalized values as ${\mathit{\alpha }}^{*}=\mathit{\alpha }/{\mathit{\alpha }}_{\mathrm{0},\mathrm{w}}$, where α_0,w is the value reported for the flat surface of pure liquid water in the same study. For studies that did not report α_0,w, we used ${\mathit{\alpha }}_{\mathrm{0},\mathrm{w}}=\mathrm{1}$ for simplicity. We convert organic concentrations into a surface coating factor f_s, defined as the ratio of the surface area occupied by an idealized monolayer of organic molecules (with the thickness of a single carbon monolayer, 0.35 nm) to the surface area of an idealized spherical water droplet. For consistency, we restrict our comparison to studies that considered semi-volatile or non-volatile organic molecules with chain-like structures and more than five carbon atoms.

Figure 6 shows α^* as a function of $\mathrm{1}/D$ based on our results (circles) and previous studies of aerosols containing liquid water and organic matter (stars) (Ma et al., 2011; Chakraborty and Zachariah, 2008; Li et al., 2019; Takahama and Russell, 2011; Diveky et al., 2019; Duffey et al., 2013; Diveky et al., 2021). The α^* values of the droplets containing pure or salty water (f_s=0) show no apparent dependence on salinity but increase with droplet size. This size dependence of α^* is consistent with previous observations that the mass accommodation coefficient of small water droplets depends on their surface curvature, particularly at droplet diameters below 3 nm (Barclay and Lukes, 2019). Although the origin of this size dependence is unclear, we note that previous studies have demonstrated that the angle of incidence of impinging molecules can significantly impact α^*, with more tangential collisions resulting in smaller α^* values (Nagayama and Tsuruta, 2003; Garrett et al., 2006). For an idealized geometry (with point-sized water molecules impacting a perfectly spherical droplet), the distribution of incidence angles of impinging water vapor molecules should be invariant with droplet size. The size dependence of α^* suggests that deviations from ideality, such as particle surface roughness or the finite size of incident water molecule, may impact the outcome of particle–water collision events.

Figure 6 Normalized water accommodation coefficient α^* versus $\mathrm{1}/D$ for different droplets. Circles show results calculated in this study. Stars show previous MD simulation (Ma et al., 2011; Chakraborty and Zachariah, 2008; Li et al., 2019; Takahama and Russell, 2011) and experimental results (Diveky et al., 2019; Duffey et al., 2013; Diveky et al., 2021) obtained at ambient temperature categorized by their f_s values.

For our low-PML-loading systems (f_s=0.30), shown as yellow circles, calculated α^* values show similar size dependence as for NaCl droplets and converge to 0.70 in droplets larger than 3 nm. This observation is consistent with the expectation that for droplets with a sub-monolayer organic coating, the accommodation coefficient can be roughly approximated as the fraction of the droplet surface that is not coated by organic matter, i.e., ${\mathit{\alpha }}^{*}\simeq \mathrm{1}-f_{\mathrm{s}}$, as expected based on the simplistic assumption that ${\mathit{\alpha }}^{*}\simeq \mathrm{1}$ on the surface of liquid water and ${\mathit{\alpha }}^{*}\simeq \mathrm{0}$ on the organic-coated surface. This relationship between f_s and α^* is also consistent with previous MD simulation results obtained with water droplets containing azelaic acid, a C9 dicarboxylic acid (f_s=0.35, ${\mathit{\alpha }}^{*}=\mathrm{0.60}$) (Ma et al., 2011). For our high-PML-loading systems (f_s=1.50), shown as orange circles, predicted α^* values are significantly smaller than for our low-PML-loading systems, as expected, with ${\mathit{\alpha }}^{*}\simeq \mathrm{0.23}$ at the largest droplet sizes. Again, our α^* results for high-organic-loading systems are consistent with values obtained in previous MD simulation and experimental studies for systems with f_s>1, which range from 0 to 0.3 as shown in Fig. 6. Interestingly, α^* now decreases with increasing droplet size, displaying a size dependence opposite to that observed for the other systems discussed above. This unexpected behavior appears to be caused by the distinct size dependence of the phase-mixing state observed for organic-rich droplets, with smaller droplets displaying an organic-core, water-coated morphology, whereas larger droplets have an engulfed structure that exposes less water at its surface.

4 Conclusions

To the best of our knowledge, this work presents an unprecedentedly broad MD simulation examination of the properties of liquid water in sub-10 nm aerosol particles. Specifically, we characterize droplet shape, water and solute distribution, water diffusivity parallel and perpendicular to the interface, surface tension, water evaporation free energy, water activity, and the water accommodation coefficient in aerosol particles with six different aqueous chemistries for a range of particle sizes. We build upon previous MD simulation studies that have carried out more piecemeal characterizations of water structure, dynamics, and energetics in aerosol particles (Karadima et al., 2019; Chowdhary and Ladanyi, 2009; Li et al., 2011, 2010). We propose a sphericity factor (ϕ) that enables a quantitative characterization of the mixing structure of water and organic- or salt-rich phases and helps reveal the links between aerosol particle morphology and the distinct energetic and kinetic properties of water molecules at the nanometer scale.

A key finding of this study is that in aerosol particles larger than 4 nm, the size dependence of water microphysics in ultrafine aerosol particles is relatively mild and can be represented using continuum-scale models. In contrast, in sub-4 nm diameter particles, multiple behaviors are observed that require fundamental understanding of molecular-level interactions. For particles consisting of a water-rich phase, smaller droplets display less water cohesion, stronger ion concentration enhancement in the droplet core, faster interfacial water diffusion, smaller surface tension, a smaller water accommodation coefficient, and higher water activity. For organic-rich aerosols with the dicarboxylic acid PML as the organic compound, a micelle structure (with an organic core and water coating) can form at small particle sizes, which results in distinct water microphysics compared to water-core aerosols at similar sizes.

Finally, we reiterate that our simulations use highly simplified chemistries, with a single type of salt (NaCl) or organic solute (PML), in order to focus on the impacts of droplet size and solute loading. The organic compound, PML, was selected for its similarity to the compounds most abundantly observed in organic aerosols in terms of molecular weight (Thornton et al., 2020), number of C atoms (Zheng et al., 2021), O:C ratio (Zheng et al., 2021; Song et al., 2018), and functional groups (Zhao et al., 2018; Qi et al., 2022). Although the use of such an idealized proxy facilitates our effort to evaluate the sensitivity of aerosol properties to salt and organic loading, additional studies are required to evaluate the sensitivity of our predictions to the type of salt or organic solutes (or solute mixtures). In particular, future studies should extend this work to other organic compounds abundantly found in organic aerosols, such as alcohols or other dicarboxylic acids.

Appendix A: Discussion of the use of the equimolar radius R_e in surface tension calculations

As noted in the main text, the surface tension should be rigorously calculated at the radius R_s, where the surface tension applies, instead of R_e. Our use of R_e (instead of R_s) in Eq. (7) is an approximation and can be justified for the following reasonings. (1) For theoretical analysis, the impact of the difference between R_e and R_s on the σ calculation is small in most of our simulated systems. The Tolman length δ can be approximated as a measurement of the difference of (R_e−R_s). Previous studies, including MD simulations, theoretical analysis, and experimental measurement, have reported δ values ranging from −0.1 to 0.2 nm for liquid water systems (e.g., water droplets, bubbles, water cavitation in minerals, and water with a hydrophobic solute) (Min and Berkowitz, 2019; Wilhelmsen et al., 2015; Kim and Jhe, 2020; Burian et al., 2017; Menzl et al., 2016; Azouzi et al., 2013; Joswiak et al., 2013; Sedlmeier and Netz, 2012). Among these studies, Kim and Jhe (2020) directly calculated (R_e−R_s) from MD simulations of pure-water droplets with radii ranging from 0.5 to 1.5 nm. Their results show that for nanodroplets with radii larger than 0.9 nm, (R_e−R_s) is in the range of [−0.05, 0.05] nm. Our results in Sect. 3.4 similarly show that δ is very small (∼0.06 nm). The sensitivity of Eq. (7) to the choice of R value can be evaluated using the expression

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}{f}_{\mathrm{err}}& ={\displaystyle \frac{\mathrm{d}\mathit{\sigma }/\mathrm{d}R}{\mathit{\sigma }}}\left(R_{\mathrm{s}}-R_{\mathrm{e}}\right)={\displaystyle \frac{-\mathrm{2}{R}^{-\mathrm{3}}}{R^{-\mathrm{2}}}}\left(R_{\mathrm{s}}-R_{\mathrm{e}}\right)\\ & ={\displaystyle \frac{-\mathrm{2}}{R}}\left(R_{\mathrm{s}}-R_{\mathrm{e}}\right).\end{array}\end{array}$$

For $(R_{\mathrm{s}}-R_{\mathrm{e}})<\mathrm{0.06}$ nm, which is the case in our study, the maximum error introduced by the use of R_e in σ calculation is less than 8 % for droplets with N_w≥500. This systematic error is commensurate with the statistical error of our predicted σ values (as shown by our reported error bars of $\sim \pm \mathrm{8}$ %). (2) In MD simulation practice, as noted in our article, the droplets simulated here do not have a perfectly spherical shape. As a result, it is impossible to accurately determine the surface tension surface R_s in a traditional manner. The surface tension calculation scheme used in our study (Eqs. 6–8), which uses R_e for approximation, has been widely used in previous MD simulation studies for nanodroplet surface tension calculation (Li et al., 2010; Zhao et al., 2019; Li et al., 2013; Zakharov et al., 1997).

Table A1 Values of equimolar radius (R_e), interfacial width (d_i), inverse of the sphericity factor ($\mathrm{1}/\mathit{\varphi }$), and surface tension (σ) for water–NaCl mixtures with different ion concentrations (C) and for water–pimelic acid droplets with different values of surface organic loading (C_s) for all simulated systems. R_e values were determined from the location of water's Gibbs dividing surface except where denoted with a ^*, where they were calculated as $\left[\mathrm{3}\right(N_{\mathrm{w}}/{\mathit{\rho }}_{\mathrm{w}}+N_{\mathrm{org}}/{\mathit{\rho }}_{\mathrm{org}})/\mathrm{4}\mathit{\pi }{]}^{\mathrm{1}/\mathrm{3}}$ using the bulk water density 34.02±0.19 molec. nm⁻² and bulk organic density 4.373±0.003 molec. nm⁻² calculated from separate MD simulations of each fluid.

Figure A1 A 1 m NaCl droplet with 5000 water molecules in a simulation cell of edge length of 21 nm. Simulation cell boundaries are represented as blue lines.

Figure A2 Predicted ϵ value of NaCl droplets as a function droplet diameter D with different d_exc values.

Figure A3 MD simulation snapshots of droplets with water (gray color), Na⁺ (dark-blue spheres), Cl⁻ (light-blue spheres), and pimelic acid (colored chains with O atoms in red, C atoms in light blue, and H atoms in white): (a) N_w=3000; (b–d) 1, 2, and 4 m NaCl solutions with N_w=3000; (e, f) low and high organic loadings with N_w=3000; and (g, h) low and high organic loadings with N_w=100. Panel (i) shows the same system as panel (h) (high-organic-loading cluster with N_w=100) but with water represented as a dark-blue surface.

Figure A4 Atomic density profiles of different components relative to distance from the Gibbs dividing surface (a, c, g, i, k, m) and the instantaneous interface (b, d, f, h, j, l, n). The top row (a, b) shows water density profiles in pure-water droplets. Rows 2–3 (c–f) show water and ion density profiles in droplets with 1 m NaCl. Rows 4–7 (g–n) show results at other aqueous chemistries.

Figure A5 Atomic density profiles of different components relative to distance from the Gibbs dividing surface for the phase-separated systems (4 m NaCl or high organic loading).

Figure A6 Mean square displacement (MSD) of water molecules as a function of time evaluated in different interfacial layers in the direction parallel to the time-averaged Gibbs dividing surface in simulations carried out with different thermostats: (a) NVT, (b) NVE (number of particles–volume–energy), (c) NVT with fixed angular momentum, and (d) Langevin dynamics (LD) with a friction constant of γ=0.5 ps⁻¹ for a pure-water droplet with N_w=1000.

Figure A7 Diffusion coefficient D_∥ of water molecules (i.e., in the direction parallel to the interface) as a function of distance from the interface in simulations of a pure-water droplet with N_w=1000 carried out using different thermostats (calculated from the slope, from 2 to 6 ps, of the MSD data shown in Fig. A4).

Figure A8 Water diffusion coefficient in the direction parallel to the interface, as a function of distance from the interface: (a) under the center of mass (com) and the instantaneous interface (inst) schemes for droplets containing 3000 water molecules with salinity of 0, 1, and 2 m NaCl and low organic surface loading (“lorg”) and (b) in different simulated systems under the com scheme.

Figure A9 Validation of our numerical solution of the anisotropic Smoluchowski equation: (a) survival probability of water molecules as a function of time in a bulk-liquid-like water layer (with a distance from the Gibbs dividing surface of [1, 1.5] nm) in 5000 water droplets under the com scheme and (b) the same calculation for a bulk-liquid-like water layer (with a distance from the instantaneous interface of [1.5, 1.8] nm) under the inst scheme. The water diffusion coefficient perpendicular to the interface (D_⟂) is calculated by optimizing the agreement between the water survival probability in MD simulations (blue circle) and the Smoluchowski prediction obtained with different D_⟂ values (black symbols). The results show that both schemes show consistent results for D_⟂ in bulk liquid droplet regions, with $D_{⟂,\mathrm{bulk}}=(\mathrm{0.22},\mathrm{0.24})$ Ų ps⁻¹.

Figure A10 Survival probability of water molecules in a 0.3 nm thick interfacial layer in 5000 water droplets: (a) under the com scheme and (b) under the inst scheme. The results show that $D_{⟂,\mathrm{interface}}=(\mathrm{0.20},\mathrm{0.22})$ Ų ps⁻¹ under the com scheme and $D_{⟂,\mathrm{interface}}=(\mathrm{0.15},\mathrm{0.17})$ Ų ps⁻¹ under the inst scheme.

Figure A11 (a) Instantaneous normal pressure tensor in 1000 water droplets at the equatorial plane perpendicular to the z axis. (b) Instantaneous value of r²P_N(r) versus the radial coordinate in the droplet containing 1000 water molecules.

Figure A12 Surface tension (σ) of droplets with different compositions as a function of droplet diameter.

Appendix B: Sensitivity of ion concentration enhancement to ion interfacial exclusion width

The first hydration shell radii (r_hyd) of alkali cations Li⁺, Na⁺, and K⁺ are 0.279, 0.318, and 0.359 nm, respectively, as reported by Joung and Cheatham (2009). The interfacial width parameter (exclusion width d_exc) was determined following two steps: (1) from the ion density profile relative to the Gibbs dividing interface (Fig. A4e), we can see that the distance where Na⁺ density reaches the well-mixed density of 1 m ranges from 0.3 to 0.5 nm for all simulated droplets, which provides estimates of the upper and lower bounds of d_exc; (2) we tested different values of the ion exclusion width and calculated the corresponding ion concentration enhancement factor and, then, compared the calculated results with our MD simulation results. Our results show that setting d_exc to 1.2r_hyd yielded the best match to our MD simulation results. The sensitivity of ϵ values of NaCl droplets to d_exc is illustrated in Fig. A2.

Appendix C: Deliquescence relative humidity as a function of nanoparticle size

To evaluate the DRH based on our simulation results, we assume that for a nanoparticle at the DRH, ion concentration in the droplet core equals the solubility of the salt C_sol. We denote the wet particle diameter (including water molecules) as D_wet and the cation hydration radius as r_hyd. If we follow the same assumptions as in Sect. 3.2 (i.e., that the salt is excluded from the vicinity of the water surface in a region with a thickness equal to 1.2r_hyd), the ion concentration corresponding to a well-mixed system C can be calculated from

$$\begin{array}{}\text{(C1)}& \mathrm{1}/\mathrm{6}\mathit{\pi }{D}_{\mathrm{wet}}^{\mathrm{3}}C=\mathrm{1}/\mathrm{6}\mathit{\pi }{\left(D_{\mathrm{wet}}-\mathrm{2.4}{r}_{\mathrm{hyd}}\right)}^{\mathrm{3}}{C}_{\mathrm{sol}}.\end{array}$$

The dry diameter of this deliquescent particle D can be approximately evaluated as falling within the range [D_salt, D_wet], where D_salt is the diameter of a pure salt crystal approximated as a sphere with density equal to the bulk solid density. The DRH as a function of D (DRH(D)) can be calculated from

$$\begin{array}{}\text{(C2)}& \mathrm{DRH}\left(D\right)=a_{\mathrm{w}}\left(C\right),\end{array}$$

where the expression for a_w(C) was taken from experimental studies of bulk aqueous solutions (El Guendouzi and Dinane, 2000).

Appendix D: Diffusion coefficient calculation

D1 Sensitivity to the simulation thermostat

We tested the impact of different thermostat settings – NVT, NVE, NVT with fixed system angular momentum, and Langevin dynamics (LD) with friction constant γ=0.5 ps⁻¹ (this value results in a friction that is lower than the internal friction of water, while still providing efficient thermostatting) – on the value of the water diffusion coefficient D_∥ using the Einstein relation for a water droplet with N_w=1000. Resulting plots of water mean square displacement (MSD) versus time are shown in Fig. A6. Different symbols in Fig. A6 represent the MSD values of water molecules in different 3 Å layers categorized based on distance from the time-averaged Gibbs dividing interface. As shown in Fig. A6, the slope of MSD versus time decreases from the interfacial region to the core of the droplet, indicating that D_∥ decreases from the interfacial region to the core. As discussed in the main text, the time interval $\mathit{\tau }=[\mathrm{2},\mathrm{6}]$ ps was chosen for the calculation of D_∥, that is $D_{\parallel }=\mathrm{1}/\mathrm{4}\times \frac{\mathrm{\Delta }\mathrm{MSD}}{\mathrm{\Delta }\mathit{\tau }}=\mathrm{1}/\mathrm{4}\times \frac{\mathrm{MSD}\left(\mathrm{6}\mathrm{ps}\right)-\mathrm{MSD}\left(\mathrm{2}\mathrm{ps}\right)}{\mathrm{4}\mathrm{ps}}$.

The corresponding D_∥ values determined for different water layers in simulations carried out with different thermostats are shown in Fig. A7. From Fig. A7, we can see that the NVT thermostat provides the largest D_∥ values for all the water layers analyzed, followed by the NVE thermostat. Results obtained using the NVT with a fixed-angular-momentum thermostat are consistent with the results obtained using the LD thermostat and approach the correct bulk liquid water diffusion coefficient value D=0.25 Ų ps⁻¹ in the droplet core. The difference in D_∥ values obtained using different thermostats is consistent with our visual observation of droplet configuration during the MD simulations. In simulations with the NVT thermostat, the simulated isolated droplet shows consistent rotational motion around the droplet's center of mass (com) such that the standard Einstein relation overestimates D_∥, particularly near the interface. This effect is inhibited in simulations carried out using the NVT with a fixed-angular-momentum thermostat or the LD thermostat. Given that the bulk-liquid-like region of individual droplets should show D_∥ values equal to the bulk liquid water diffusion coefficient, D=0.2 Ų ps⁻¹, diffusion coefficients reported in the main text are based on simulations carried out using the NVT thermostat with fixed system angular momentum.

D2 Sensitivity to the frame of reference (com and inst scheme)

D2.1 Diffusion coefficient parallel to the interface (D_∥)

The Einstein relation was used to calculate D_∥ under both the com and inst schemes, i.e., relative to the average or instantaneous interface. Both schemes are expected to provide mutually consistent results for D_∥, irrespective of distance from the interface because capillary waves exist only at the interface and are mostly transverse waves that should have negligible impact on water displacement parallel to the interface. A comparison of D_∥ values calculated under both schemes is presented in Fig. A8a for droplets of N_w=3000. As expected, the two schemes provide mutually consistent values. D_∥ values calculated for different simulated systems under the com scheme are shown in Fig. A8b.

D2.2 Diffusion coefficient normal to the interface (D_⟂)

As discussed in the main text, D_⟂ values were determined as a function of distance from the interface by comparing simulation trajectories to the solution of the anisotropic Smoluchowski equation. The com and inst schemes are expected to yield different D_⟂ values near the interface (larger in the case of the com scheme) but similar results in the bulk-liquid-like region because the com scheme intrinsically includes motions due to capillary waves as part of the interfacial water dynamics, whereas the inst scheme excludes the interface fluctuation and only quantifies water displacement relative to the instantaneous interface. Since the magnitude of capillary waves increases with droplet size (as observed based on both ϕ values and interfacial width as noted above), the difference between the D_⟂ values at the interface calculated with the two schemes should increase with droplet size.

As shown in Fig. A9, we obtained mutually consistent D_⟂ values near 0.25 Ų ps⁻¹ under both com and inst schemes for the bulk liquid region in the pure-liquid-water droplets of N_w=5000 (D=6.6 nm), which validates our calculation method. For the same systems, for the interfacial water layer with a thickness of 3 Å, we obtained a D_⟂ value in the range of [0.20, 0.22] Ų ps⁻¹ under the com scheme and a smaller value in the range of [0.15, 0.17] under the inst scheme as shown in Fig. A10. We note that although the observation of smaller interfacial D_⟂ values under the inst scheme compared with the com scheme is consistent with our expectation, our interfacial water D_⟂ values under both schemes are smaller than the values reported in previous simulation studies for the flat water–air interface (Liu et al., 2004; Wick and Dang, 2005). This discrepancy might be due to different water models, different interfacial layer definitions, or differences in the curve-fitting procedure required to determine D_⟂. As noted in the main text, our observations of relatively low D_⟂ values near the interface are consistent with reports of slower water reorientational kinetics near hydrophobic surfaces. The discrepancies between different reported values of D_⟂ for interfacial water should be carefully examined in future studies.