Interfacial thermodynamics of spherical nanodroplets: Molecular understanding of surface tension via hydrogen bond network

Nucleation of nanoscale water droplets in the atmospheric clouds is crucial for the control of earth climate Baker1997 ; Kulmala2004 . When homogeneous nucleation initiates with the formation of water nanodroplets from the supersaturated vapor, there exists the thermodynamic energy barrier formed by the difference between the favorable volume- and unfavorable surface-energy. The barrier height, which determines the nucleation rate of the condensed phase, is governed by the surface tension of the critical nucleus Kalikmanov2005 . The classical nucleation theory has employed the value of surface tension of planar interface in estimating the nucleation rate and has suffered from poor predictability, which is even orders of magnitude off from the experimental results of liquid droplets such as the simple liquid and water Oxtoby1988 ; Holten2005 ; Bruot2016 . Moreover, the surface energy model with the planar tension produces (i) inaccurate estimation of the solvation free energy of hydrophobic molecules Sedlmeier2012 and (ii) improper description of the heterogeneous nucleation processes of the concave nanodroplets formed between two hydrophilic surfaces Kwon2018 ; Kim2018 . Despite numerous theoretical Factorovich2014 ; Sampayo2010 ; Koga1998 ; Malijevsky2012 ; Rehner2018 and experimental Bruot2016 ; Kim2018 ; Kwon2018 studies performed for decades, determination of surface tension of the curved interface Bruot2016 ; Rowlinson1982 ; Joswiak2013 ; Azouzi2013 ; Wilhelmsen2015 ; Kwon2018 ; Kim2018-2 of the nucleus is still a topic of ongoing controversies. Therefore, quantitative and consistent study of surface tension of the highly curved liquid-vapor interface should be carried out for the fundamental understanding of the nucleation phenomena.

The first attempt to describe the curvature dependence of surface tension was made by the Tolman’s equation in 1949 Tolman1949 ,

where the constants $\gamma_{0}$ and $\delta_{0}$ are the surface tension and Tolman length of the planar interface, respectively. Notice that $\delta_{0}$ is defined as the difference between two possible definitions of the nanodroplet radius, equimolar radius $R_{\textrm{e}}$ and the radius of surface of tension $R_{\textrm{s}}$, when the nanodroplets become infinitely large; $\delta_{0}\equiv\lim_{R_{s}\rightarrow\infty}\delta(R_{\textrm{s}})=\lim_{R_{s}\rightarrow\infty}\left(R_{\textrm{e}}-R_{\textrm{s}}\right)$. Therefore, Tolman’s equation provides a general description of surface tension of nanodroplets in terms of only a ‘single’ parameter $\delta(R_{\textrm{s}})$ for a given $R_{\textrm{s}}$. Theoretical studies suggest that the magnitude of the constant $\delta_{0}$ is comparable to the molecular diameter Rowlinson1982 , but there are still many debates concerning the magnitude and sign of $\delta_{0}$. Moreover, despite substantial development of experimental methods Bruot2016 ; Kim2018 , there does not exist a consensus on either the value of $\delta_{0}$ nor the change of surface tension at the molecular scale because quantitative information of molecules in the nanodroplets still lacks, which may behave very differently from the planar case. Therefore, for a unified understanding of the molecular-scale surface tension and the nanoscale liquid-vapor interface, single-molecule study is on high demand, including the investigation of the free energy as well as hydrogen bond (HB) networks of individual molecules.

Here, we employ the molecular dynamics (MD) and density functional theory (DFT) to calculate the surface tension of liquid nanodroplet from the pressure profile and the thermodynamic free energy. We find curvature dependence of surface tension below 1 nm radius and derive analytic expression from the classic Tolman’s equation combined with the radius-dependent Tolman length Tolman1949 . The derived analytic expression predicts surface tension within the error below 1$\%$. Moreover, to describe the origin of the change of surface tension in the nanoscale droplet, we use the statistical analysis of individual molecules that provides the free energy of nanodroplets and also use the HB network analysis with graph theory approach. They show that surface tension is governed by the energy change of interfacial molecules, associated with the disruption of the HB network.

Figure 1(a) summarizes the main results of our work. The MD snapshots of nanodroplets for $N=32$, 96 and 256 water molecules (or equivalently $R_{\textrm{s}}=0.45$, 0.8 and 1.2 $\textrm{nm})$ are provided in the top row. The bottom row shows the network representations of the snapshots. The points (nodes) and lines (edges) of the HB network correspond to the water molecules and HBs, respectively. The strength of the interaction is qualitatively represented by the number of HB, $N_{\textrm{HB}}$, which is reflected by the color of the nodes. The color distributions in the networks and the average values of the HB number, $N_{\textrm{HB,avg}}$, demonstrate that the overall interaction strength, which is directly related to the surface tension of nanodroplets, decreases considerably when $R_{\textrm{s}}<1\ \textrm{nm}$.

To study nanodroplet systems, we have performed the MD simulation of the TIP4P/2005 water model and the DFT calculation of the Lennard-Jones liquid (refer to Supplementary Information for detailed information about MD and DFT). In the MD simulation, the number of water molecules ranges from $N=32$ to 512 at $T$ = 290 K. In DFT, we deal with nanodroplets with a radius of $R^{*}=1.77$ to 30 at the temperature $T^{*}=0.53$ to 0.91. We chose the smallest radius of droplet in DFT as the lower size limit without thermal instability. In our simulations, the reduced units are employed: $r^{*}=r/\sigma$, $R^{*}=R/\sigma$, $\delta^{*}=\delta/\sigma$, $\rho^{*}=\rho\sigma^{3}$ and $T^{*}=T/T_{c}$, where $\sigma$ is the diameter of the liquid molecule and $T_{c}$ is the critical temperature of the liquid. For the TIP4P/2005 water model Vega2006 , we used $T_{c}=641.4$ K and $\sigma=0.3159$ nm.

Let us first calculate the theoretical values of the Tolman length, $\delta^{*}(R_{\textrm{s}}^{*})=R_{\textrm{e}}^{*}-R_{\textrm{s}}^{*}$, which are presented in Fig. 1(b) as a function of the inverse radius of the surface of tension, $1/R_{\textrm{s}}^{*}$ (definitions of $R_{\textrm{e}}^{*}$ and $R_{\textrm{s}}^{*}$ are summarized in Supplementary Information). For the curvature dependence, we have adopted (and justified later) the ansatz in the linear form, $\delta^{*}\left(R_{\textrm{s}}^{*}\right)=\delta_{0}^{*}+\delta_{1}^{*}/R_{\textrm{s}}^{*}$, which fits well with the MD and DFT resultsBruot2016 ; Block2010 . The resulting numerical values of the coefficients $\delta_{0}^{*}$ and $\delta_{1}^{*}$ are listed in Table 1. Interestingly, at $T^{*}$ = 0.45 (water droplet), the sign of the Tolman length changes from negative to positive at $R_{\textrm{s}}\sim 1.2\ \textrm{nm}$, which marks the region where surface tension begins to decrease. Notice that the Tolman length we find for the planar interface of water is $\delta_{0}=-0.59\ \mathring{\textrm{A}}$, which is consistent with the previous theoretical and experimental resultsWilhelmsen2015 ; Azouzi2013 ; Joswiak2013 .

Using the simulation results of $\delta_{0}$ and $\delta_{1}$ presented in Fig. 1(b) and Table 1, we now can derive the analytic governing equation of surface tension from the rigorous Gibbs-Tolman-Koenig-Buff (GTKB) equation Gibbs1906 ; Koenig1950 ; Buff1951 ; Rowlinson1982 ,

Here we have used the series expansion using the ansatz $\delta(R_{\textrm{s}})=\delta_{0}+\delta_{1}/R_{\textrm{s}}$ and thus the resulting equation can be written as, up to the second order of $1/R_{\textrm{s}}$,

Equation (3) is our analytic result and is expected to be more exact compared to Eq. (1) since the curvature dependent Tolman length is considered. To demonstrate the validity of Eq. (3), we compare it with the surface tension directly calculated from MD and DFT results without the ansatz (Fig. 1(c)). The surface tension of the TIP4P/2005 water droplet is computed from the difference between the normal and transverse pressure tensorsGhoufi2011 . Surface tension of Lennard-Jones liquid droplet is calculated from the grand potential energy of the system (refer to Supplementary Information for detailed information). As shown in Fig. 1(c), Eq. (3) predicts excellently the surface tension of nanodroplets down to the smallest radius we consider, $R_{\textrm{s}}=0.45\ \textrm{nm}$ for water. We estimate the accuracy of Eq. (3), for example, at $T^{*}=0.53$ of Lennard-Jones droplet; Eq. (3) shows an error within 0.5 $\%$ until the lower size limit of nanodroplet, while the error jumps over 10 $\%$ without the second order term (Supplementary Fig. 1). The observed Tolman length and surface tension equation with second order correction agrees well with the recent results of Rehner and GrossRehner2018 . Our study shows that the Tolman length approach can be applied to the hydrogen-bond liquid with strong directional bond, and also the second-order correction is good for describing surface tension of 1 nm nanodroplet.

Interestingly, the second-order term in Eq. (3) is closely related to the bending rigidityHelfrich1973 , which is typically used as an empirical parameter to describe the free energy of the curved membrane. The bending rigidity, $k_{s}$, can be related to the Tolman length as $k_{\textrm{s}}=\gamma_{0}(3\delta_{0}^{2}-\delta_{1})$Helfrich1973 . For water nanodroplets, we find the negative rigidity as shown in Table 1, $k_{\textrm{s}}=-1.34k_{\textrm{B}}T$, which is comparable to the MD result Sedlmeier2012 , $k_{\textrm{s}}\approx-3k_{\textrm{B}}T$, that considered the interface between hydrophobic solute and water. Here, the negative sign of bending rigidity is the general feature of the liquid nanodroplet of water, Lennard-Jones liquid, and polymeric liquid Rehner2018 .

Let us now elucidate the physical origin of the smaller surface tension with respect to the planar interface from three perspectives; (i) pressure tensor, (ii) free energy, and (iii) HB network. First, we calculate the radial pressure tensor of molecules, $P(r)$, which comes from the sum of kinetic pressure, $\rho(r)k_{\textrm{B}}T$, obtained from radial density $\rho(r)$ and configurational pressure, $-\left<\Delta U(r)/\Delta V(r)\right>$, obtained from potential energy difference induced by volume perturbation (refer to Supplementary Information for detailed derivation). We decompose $P(r)$ into the transverse ($P_{\textrm{T}}(r)$) and normal ($P_{\textrm{N}}(r)$) profiles under the hydrostatic equilibrium condition ($\nabla\cdot\vec{P}(r)=0$) and the boundary condition ($P(r)=1/3P_{\textrm{N}}(r)+2/3P_{\textrm{T}}(r)$). At the interface of $R_{\textrm{s}}$ = 1.5 nm droplet (Fig. 2(a)), $P_{\textrm{T}}(r)$ deviates from $P_{\textrm{N}}(r)$ and shows a negative minimum, which indicates the origin of surface tension from the imbalanced interaction of the interfacial molecules. Notice that the values of negative minimum (Fig. 2(b)) and pressure imbalance (Fig. 2(c)) decrease as the droplet radius decreases, which shows that the interaction between interfacial molecules weakens accordingly and is responsible for the curvature dependence of surface tension. As expected, a direct comparison between surface tension and pressure imbalance (Fig. 2(d)) evidences the relationship between the two.

Second, we investigate the physical mechanism with regard to interfacial thermodynamics at the molecular level. For this purpose, we compute the energy ($E$), entropy ($S$), and Helmholtz free energy ($F$) of nanodroplets (Fig. 3). Entropy is derived from the integration of the density of states of molecules, which are obtained via the Fourier transform of the velocity autocorrelation function. Formalism is established by the two-phase thermodynamics (2PT) model in calculating the entropy of liquid molecules Huang2011 (refer to Supplementary Information for detailed information). Helmholtz free energy is calculated by the relation $F=E-TS$. Figures 3(a) - (f) show the probability distributions of $E$, $S$, and $F$ for the droplets of $R_{\textrm{s}}$ = 0.45 nm and 1.5 nm. By dividing two regions across $R_{\textrm{s}}$, each nanodroplet is well separated into the ‘inside’ region having low energy and low entropy, and the ‘surface’ region with high energy and high entropy. Figures 3(g) and (h) compile all the thermodynamic quantities of water molecules in the inside and surface region as a function of nanodroplet radius, respectively. For example, the average Helmholtz free energy of the surface molecules for $R_{\textrm{s}}$ = 0.45 nm is about 2 kJ/mol higher than that of $R_{\textrm{s}}$ = 1.5 nm, which evidences that the change of surface tension is associated with the decrease of intermolecular interaction strength of nanodroplets. Notice that the free energy change is dominated by the energy, which compensates for the entropy change that behaves oppositely.

Lastly, we address the molecular description in terms of the relation between the configurational energy change of nanodroplets and the change of HB network characteristics (refer to Supplementary Information for detailed information). Two water molecules are considered to be connected via HB when (i) the distance between two oxygen atoms is smaller than 0.35 nm and (ii) the angle between hydrogen and donor-acceptor oxygen atoms is smaller than $30^{\circ}$ Kumar2007 . We represent the HB network of nanodroplets by the adjacency matrix $\bf{A}$, which is a square symmetric matrix with elements of 1 if the molecular pair $(i,j)$ has an HB, and 0 if otherwise. Figure 4(a) shows the change of average number of HB, $N_{\textrm{HB,avg}}$, for the inside as well as the surface molecules in comparison to $N_{\textrm{HB,avg}}$ of $R_{\textrm{s}}$ = 1.5 nm nanodroplet. We find that the smaller nanodroplets have relatively smaller $N_{\textrm{HB,avg}}$, which is already expected by Figs. 3(g) and (h). The disruption of the HB network shown in Fig. 4(a) can be explained by the decreasing coordination number (Fig. 4(b)). Notice that these results are independent of the choice of radius, which is still observed even if $R_{\textrm{e}}$ is used instead of $R_{\textrm{s}}$ (Supplementary Fig. 2). The overall change of the HB network characteristics can be observed by the clustering coefficient and cyclic number distribution of the network presented in Figs. 4(c) and (d), respectively. The clustering coefficient of the HB network is defined as the probability that the two hydrogen-bonded neighboring molecules are connected via the HB and form a triangular network. The cyclic number is defined as the number of edges connected in a closed HB loop in a network. As shown, the nanodroplet network has a larger average clustering coefficient and a smaller average cyclic number than the bulk water network, which demonstrates that the nanodroplets are more structured than the bulk water, and their structure is similar to that of the high-density water or Ice-VII Shin2019 . Moreover, from the eigenvalue spectra of $\bf{A}$ (Figs. 4(e) and (f)), one can observe that the change of network structure is concentrated in the local network topology with large eigenvalue of $\bf{A}$ (or large number of HB), which shows that the disruptions of HB networks are focused on the four or five HBs. Because of the large energy content of water HB, disruption of the HB network directly leads to the increase of free energy and the resulting decrease of surface tension.

As a conclusion, we emphasize that the above three perspectives provide the unified understanding of the decreasing surface tension in terms of the decrease of interaction strength (Fig. 2), an increase of energy (Fig. 3) and the decrease of HB number as well as the change of HB network structure (Fig. 4) in the molecular-scale nanodroplets. In particular, the large HB energy of water (i) compensates the entropy change and (ii) dominates the change of Helmholtz free energy and the change of pressure imbalance.

In summary, we have demonstrated the curvature dependence of surface tension from the perspectives of interfacial thermodynamics using the complementary methods of MD and DFT simulations. We show that the surface tension is a decreasing function of the inverse radius of nanodroplets with the positive Tolman length. Statistical analysis of the Helmholtz free energy of individual water molecules shows that the change of surface tension is governed by the surface molecules having higher free energy than that of the planar interface, which is closely related to the disruption of the HB network at the interface. Our results may trigger further theoretical development for an accurate understanding of the nucleation of nanoclusters at the molecular level. Moreover, the thermodynamic approach presented here can be used to study the chemically heterogenous interface where the local energy varies at the nanoscale.