Surface tension of supercooled water nanodroplets from computer simulations

TL;DR

This study uses computer simulations to analyze the surface tension of supercooled water nanodroplets, revealing temperature-dependent behavior linked to structural changes and a potential liquid-liquid phase transition.

Contribution

It provides a detailed comparison of mechanical and thermodynamic methods for calculating surface tension and introduces the use of surface tension as a probe for phase transition behavior in supercooled water.

Findings

Different results from the two calculation routes for surface tension and Tolman length.

Surface tension increases more rapidly below 230 K, indicating structural changes.

Surface tension correlates with the formation of a tetrahedral network, suggesting a liquid-liquid phase transition.

Abstract

We estimate the liquid-vapour surface tension from simulations of TIP4P/2005 water nanodroplets of size $N$=100 to 2880 molecules over a temperature $T$ range of 180 K to 300 K. We compute the planar surface tension $\gamma_p$, the curvature-dependent surface tension $\gamma_s$, and the Tolman length $\delta$, via two approaches, one based on the pressure tensor (the "mechanical route") and the other on the Laplace pressure (the "thermodynamic route"). We find that these two routes give different results for $\gamma_p$, $\gamma_s$ and $\delta$, although in all cases we find that $\delta\ge 0$ and is independent of $T$. Nonetheless, the $T$ dependence of $\gamma_p$ is consistent between the two routes and with that of Vega and de Miguel [J. Chem. Phys. 126, 154707 (2007)] down to the crossing of the Widom line at 230 K for ambient pressure. Below 230 K, $\gamma_p$ rises more rapidly on cooling than predicted from behavior for $T\ge 300$ K. We show that the increase in $\gamma_p$ at low $T$ is correlated to the emergence of a well-structured random tetrahedral network in our nanodroplet cores, and thus that the surface tension can be used as a probe to detect behavior associated with the proposed liquid-liquid phase transition in supercooled water.

Figures

Equations

$$\Delta P=\frac{2{\gamma }_s}{R_s}.$$

$${\gamma }_s=\frac{{\gamma }_p}{\left(1+2\delta /R_s\right)}.$$

$$\delta =R_e-R_s,$$

$$\Delta P=\frac{2{\gamma }_p}{R_e}\left(\frac{1}{1+\delta /R_e}\right),$$

$$\frac{2}{\Delta P{R}_e}=\frac{1}{{\gamma }_p}\left(1+\delta /R_e\right),$$

$${\gamma }_s={\gamma }_p\frac{R_e-\delta }{R_e+\delta }.$$

$${\gamma }_s$$

$$=$$

$$R_s=\frac{{\int }_0^{\infty }{r}^2\left[P^{\alpha ,\beta }(r;R_s)-P_T(r)\right]dr}{{\int }_0^{\infty }r\left[P^{\alpha ,\beta }(r;R_s)-P_T(r)\right]dr}.$$

$${\int }_0^{\infty }{r}^2\left[P^{\alpha ,\beta }(r;R_s)-P_N(r)\right]dr=0,$$

$${\gamma }_s$$

$$R_s=\frac{{\int }_0^{\infty }{r}^2\left[P_N(r)-P_T(r)\right]dr}{{\int }_0^{\infty }r\left[P_N(r)-P_T(r)\right]dr}.$$

$${\gamma }_s^3$$

$$q_i$$

Full text

I Introduction

Microscopic and nanoscopic water droplets are of interest in many important research areas, such as the Earth’s climate baker ; wil , biological applications Ohno , interstellar space Tachibana , and numerous other systems klemp . In all of these areas, the surface tension of the liquid-vapour interface of the water droplet is a central physical property for understanding and predicting droplet behavior. For example, the surface tension is crucial for estimating the nucleation rate of liquid from the vapour using classical nucleation theory kulmala ; Debenedetti-book .

The surface tension is also the origin of the pressure difference that arises between the interior and exterior of a liquid droplet, as quantified by the Young-Laplace equation young1805 ; laplace1805 ,

[TABLE]

Here, $\Delta P=P_{l}-P_{v}$, where $P_{l}$ and $P_{v}$ are the respective pressures of the liquid interior and vapour exterior, and $\gamma_{s}$ is the surface tension of the curved interface. $R_{s}$ is the radius of the so-called “surface of tension” Rowlinson1982 . For macroscopic droplets, the width of the molecular interface is negligible compared to the droplet dimensions, and $R_{s}$ is simply the radius of the droplet. However, for nanoscale droplets, the interfacial width is significant compared to the size of the droplet itself, and various definitions for the radius of the droplet are possible.

It has long been understood that the surface tension of a curved interface deviates from that of a planar interface. For a curved surface, such as that of a droplet, the Tolman length $\delta$ quantifies how $\gamma_{s}$ deviates from the planar surface tension $\gamma_{p}$ as a function of $R_{s}$, via the expression Tolman1949 ,

[TABLE]

The magnitude of $\delta$ is generally found to be 10-20% of the molecular diameter.

However, the sign of $\delta$ is a subject of continuing debate Blokhuis2009 . While modeling on the basis of classical density functional theory has predicted negative values of $\delta$ for liquid Lennard-Jones (LJ) droplets Blokhuis2013 ; reguera2015 , simulations of droplets have estimated both negative and positive values of $\delta$. For example, Yan, et al. Yan2016 performed MD simulations of liquid argon nanodroplets with sizes ranging from 800 to 2000 atoms at 78 K, as modelled using the LJ potential. They evaluated the pressure tensor, and using the Young-Laplace equation they concluded that $\delta$ is positive for LJ nanodroplets. However, Giessen and Blokhuis Blokhuis2009 estimated a negative value of $\delta$ for LJ nanodroplets.

A similar disagreement regarding the magnitude and sign of $\delta$ appears in water simulations. Leong and Wang Leong2018 performed MD simulations using the BLYPSP-4F water potential BLYP on nanoscale droplets with radii varying between 2 and 8 nm at temperature $T=298$ K. Using an empirical correlation between the pressure and density, they estimated $\delta=-0.048$ nm. A similar value for $\delta$ was obtained by measuring the free energy of droplet mitosis in a study by Joswiak, et al. Joswiak for the mW model of water mW . On the other hand, Lau, et al. Lau2015 used a test-area method and obtained a positive value of $\delta$ for the TIP4P/2005 model of water vega2005 . Simulation studies of cavitation for TIP4P/2005 find relatively large positive values of $\delta$ for vapour bubbles with magnitudes in the range of 0.12 to 0.195 nm menzl2016 ; min2019 . This result implies that for a TIP4P/2005 water droplet of the same size, $\delta$ should be of similar magnitude, but negative. It is evident from this recent work that disagreement exists on both the magnitude and sign of $\delta$, even when the same water model is used.

The variation of the surface tension with $T$ for deeply supercooled water has also been investigated, in particular as a way to test for evidence of a possible liquid-liquid phase transition (LLPT) in supercooled water Poole1992-Phase . Theoretical studies have shown that if a LLPT occurs, then at low $T$ the surface tension should increase faster with decreasing $T$ than is expected otherwise feeney ; hruby2004 ; hruby2005 . Some computer simulations studies are consistent with this behavior lu2006a ; lu2006b while others are not chen ; vier . Recent careful experiments by Hruby and coworkers do not find evidence for a change in the $T$ dependence of the surface tension for $T$ as low as $-26$ ∘C hruby2014 ; hruby2015 ; hruby2017 . However, it is possible that the anomalous increase in the surface tension will only be observed for $T$ below the Widom line that is associated with the LLPT, a range of $T$ that has only recently begun to be probed in experiments kim2017 .

The sign of $\delta$ determines whether $\gamma_{s}$ decreases or increases with $R_{s}$. For a positive $\delta$, $\gamma_{s}$ decreases as $R_{s}$ decreases. Moreover, $\delta$ relates the equimolar radius $R_{e}$ and $R_{s}$ Tolman1949 ,

[TABLE]

where $R_{e}$ is the radius of a sphere that has a uniform density equal to that of the interior part of the droplet and that has the same number of molecules as the droplet. Since determining $R_{e}$ is more straightforward than determining $R_{s}$, we can rewrite Eqs. 1 and 2 in terms of $R_{e}$,

[TABLE]

or in the form,

[TABLE]

and

[TABLE]

The above equations provide the basis for a procedure to find $\gamma_{p}$, $\gamma_{s}$ and $\delta$, which following past practise we refer to here as the “thermodynamic route” Thompson1984 . As we will see below, computer simulations of water nanodroplets allow us to directly estimate $\Delta P$ and $R_{e}$. If we obtain $\Delta P$ and $R_{e}$ for a range of droplet sizes at fixed $T$, we can use Eq. 5 to estimate $\gamma_{p}$ and $\delta$ by curve fitting. From $\delta$ an estimate of $R_{s}$ is obtained from Eq. 3, and so an estimate of $\gamma_{s}$ can be computed using Eq. 1.

Aside from the Laplace equation, Rowlinson and Widom proposed a model to derive $\gamma_{s}$ from the tangential and normal components of the pressure tensor as functions of the radial distance $r$ from the centre of mass of a droplet, $P_{T}(r)$ and $P_{N}(r)$ Rowlinson1982 . The model assumes two homogeneous fluid phases, with homogeneous pressures $P^{\alpha}$ and $P^{\beta}$ far from the interface, and an inhomogeneous interface between them. Under the model assumption that the surface tension acts at a single value of $r=R_{s}$, the mechanical requirements for static equilibrium, i.e. force and torque balance, yield,

[TABLE]

where $P^{\alpha,\beta}(r;R_{s})$ is $P^{\alpha}$ for $r<R_{s}$ and $P^{\beta}$ for $r>R_{s}$. These equations in turn give an expression for $R_{s}$,

[TABLE]

With the assumption that the two phases are homogeneous, we can assume that $P^{\alpha}=P_{l}$ and $P^{\beta}=P_{v}$. Since $P^{\alpha,\beta}(r;R_{s})$ depends on $R_{s}$, Eq. 9 must be evaluated numerically.

From the condition of mechanical stability, $\nabla\cdot{\bf P}=0$, it can be shown that,

[TABLE]

and hence $\gamma_{s}$ can be obtained using the $P_{N}(r)$ component of the pressure in Eqs. 7 and 8, yielding,

[TABLE]

and Rowlinson1982 ,

[TABLE]

We can also find $\gamma_{s}$ while avoiding the need for $R_{s}$ by combining Eqs. 1 and 11,

[TABLE]

Eqs. 7-13 provide an alternative pathway, referred to as the “mechanical route” Thompson1984 , to compute $\gamma_{s}$ and $R_{s}$, as well as $\gamma_{p}$ and $\delta$. First, $P_{N}(r)$ and $P_{T}(r)$ are calculated from simulations of nanodroplets, with Eqs. 7–13 yielding values for $\gamma_{s}$ and $R_{s}$. Estimates for $\gamma_{p}$ and $\delta$ are then obtained through Eqs. 2 and 3.

In this study, we use the TIP4P/2005 model to simulate water nanodroplets over a wide range of temperatures and sizes and determine $\gamma_{s}$, $\gamma_{p}$ and $\delta$ using both the thermodynamic and mechanical routes. In Section II, we provide details of our simulations. In Section III-V, we show that the thermodynamic and mechanical routes give different results for the surface tension and Tolman length for water. Despite these differences, we find that all methods demonstrate that $\gamma_{p}$ increases more rapidly upon cooling through the Widom line temperature for TIP4P/2005. In Section VI, we show that the rapid increase in $\gamma_{p}$ at low $T$ is consistent with a crossover within the core of our nanodroplets from the high density liquid phase (HDL) to the low density liquid phase (LDL) of the LLPT. We present a discussion and our conclusions in Section VII.

II Simulations

We recently studied the thermodynamic and structural properties of simulated water nanodroplets ranging in size from $N=$ 100 to 2880 molecules, over a $T$ range of 180 to 300 K NATURECOMM , with molecules interacting through the TIP4P/2005 model vega2005 . The same data set is used in the present study. We summarize the simulation details below for the reader’s convenience.

We carry out the simulations in the canonical ensemble – constant $N$, volume $V$, and $T$. The droplets are located in a periodic cubic box of side length $L$ that increases with $N$ and ranges from 10 to 20 nm. We ensure the box is large enough to avoid any direct interaction between the water droplet and its periodic images, and small enough to ensure that at most only a few molecules are in the vapour phase. We use a potential cutoff of $L$/2, ensuring that all molecules in the droplet interact without truncation of the potential. We use Gromacs v4.6.1 GROMACS to carry out our molecular dynamics (MD) simulations. We hold the temperature constant using the Nosé-Hoover thermostat with time constant 0.1 ps. The equations of motion are integrated with the leap-frog algorithm with a time step of 2 fs.

The data set is generated from two kinds of MD runs: conventional “single long runs” (SLR), and using a “swarm relaxation” method (SWRM) swarm . For droplet sizes $N=100$, 200, 360, 776, 1100, 1440, and 2880, we use SLRs. For $N=1440$ and 2880, we start our simulations by placing $N$ molecules randomly within the simulation box, and run long enough for the molecules to condense into a single droplet. We harvest an equilibrated $N=1440$ configuration, and progressively remove molecules from the droplet surface to obtain starting configurations for the other droplet sizes. The slowest relaxation times are approximately 12 ns, and our longest post-equilibration simulations last 2.8 $\mu$s.

For droplet sizes $N=205$, 301, 405, 512, 614, and 729, we use SWRM. To generate initial configurations for each of these droplet sizes, we first remove molecules from the surface of an equilibrated $N=2880$ configuration to obtain the desired size. We first conduct SLRs for each size at $T=200$ K for not less than 350 ns. We then take the last configuration of each run and randomize the velocities using the Maxwell-Boltzmann distribution at $T=220$ K to generate $M$ different configurations, which are used to initiate our swarm relaxation runs. We determine the relaxation time $\tau_{s}$ for a swarm ensemble from the potential energy autocorrelation function of the system. See Ref. swarm for details. The final equilibrated $M$ configurations of the ensemble are then used to initiate an ensemble of runs at $T=200$ K. Similarly, we take the final equilibrated $M$ configurations of the ensemble at $T=200$ K to start a swarm ensemble at $T=180$ K.

Additionally, we carry out simulations for bulk liquid TIP4P/2005 with $T$ varying from $300$ to $180$ K. We simulate 360 molecules with density varying approximately between $0.96$ and $1.12$ g/cm3 using the protocols described in Ref. saika2013 .

The mechanical route to finding the surface tension of a droplet requires the determination of both $P_{T}(r)$ and $P_{N}(r)$. We compute kinetic and configurational contributions to the pressure inside our droplets; see Ref. malek2 for details on applying to TIP4P/2005 a coarse-grained method Ikeshoji2011 based on the Irving-Kirkwood IrvingKirkwood choice of contour in defining the microscopic pressure. Fig. 1 shows all contributions to the pressure for two example cases. We define $R_{L}$ such that the configurational contributions to the normal and tangential pressures, ${P}_{c,N}$ and ${P}_{c,T}$ respectively, are equal to each other within error for $r<R_{L}$ (dashed line in Fig. 1), noting that they differ near the surface. To define the pressure in the interior of the droplets $P_{L}$, we average the total (isotropic) pressure $P_{tot}(r)={P}_{c,N}(r)/3+2{P}_{c,T}(r)/3+\rho_{\circ}(r)k_{\rm B}T$ over the spherical volume of radius $R_{L}$, where $\rho_{\circ}(r)$ is the local number density.

All error bars reported in this work indicate one standard deviation in the mean.

III Thermodynamic route

In the thermodynamic route, we use the Young-Laplace equation in the form of Eq. 4 to determine $\gamma_{p}$ and $\delta$. We equate $P_{L}$ with $\Delta P$, since the vapour pressure is negligible, and plot isotherms of $P_{L}$ as a function of $R_{e}^{-1}$ in Fig. 2. The isotherms show that there is a significant pressure that naturally builds up in the interior of the droplets, and it can reach more than 200 MPa for $R_{e}^{-1}\simeq 1.2$ nm*-1* ($R_{e}\simeq 0.83$ nm).

In Fig. 2 we study the curvature correction to $P_{L}$ as a function of $R_{e}$. Assuming $\delta=0$, the fits in Fig. 2a using Eq. 4 (fitting only for $\gamma_{p}$) show that there is no obvious curvature correction to the Young-Laplace equation. To see how small $\delta$ is in our range of droplet sizes, we fit $P_{L}$ as a function of $R_{e}^{-1}$ at each $T$ with Eq. 4 (fitting for both $\gamma_{p}$ and $\delta$), as shown in Fig. 2b. We report the value of $\delta$ as a function of $T$ in Fig. 3, and can discern no clear dependence of $\delta$ on $T$. The average small (positive) value of the Tolman length $\mkern 1.5mu\overline{\mkern-1.5mu\delta\mkern-1.0mu}\mkern 3.0mu=0.055$ nm explains the absence of strong curvature in the isotherms of Fig. 2. As an alternative way of obtaining $\gamma_{p}$ and $\delta$, we plot isotherms of $2/(P_{L}R_{e})$ as functions of $R^{-1}_{e}$ in Fig. 4. Since $\delta$ does not have an apparent dependence on $T$, we fit the isotherms in Fig. 4 to Eq. 5 assuming a single common value of the fitting parameter $\delta$ for all $T$. As shown in Fig. 4, this global fit reasonably describes all the isotherms, and gives a value of $\delta=0.056$ nm that is similar to $\mkern 1.5mu\overline{\mkern-1.5mu\delta\mkern-1.0mu}\mkern 3.0mu$. The intercepts in Fig. 4 yield $1/\gamma_{p}$ for each $T$. As shown in Fig. 5a, the values of $\gamma_{p}$ obtained in this way increase as $T$ decreases.

In Fig. 5a, we compare estimates for $\gamma_{p}$ assuming $\delta$=0 (obtained from the fits in Fig. 2a) and $\delta\neq 0$ (obtained from the fits in Fig. 2b). At $T=180$ K the discrepancy in $\gamma_{p}$ between assuming $\delta=0$ and $\delta\neq 0$ appears to be outside of error, with the curvature-corrected result yielding a value of $\gamma_{p}$ approximately 10% higher (blue squares versus red circles in Fig. 5a). For $T\geq 220$ K our estimates of $\gamma_{p}$ are also consistent with the extrapolation down to low $T$ of $\gamma_{p}$ obtained using the test-area method, taken from Eq. 6 in the work of Vega and de Miguel Vega2007 .

To obtain $\gamma_{s}$ from the thermodynamics route, we combine Eqs. 1 and 3 to obtain $\gamma_{s}=P_{L}\,(R_{e}-\delta)/2$, and use the values of $\delta$ for each $T$ obtained from fits shown in Fig. 2b. The resulting values of $\gamma_{s}$ are shown in Figs. 6a and 7.

IV Mechanical route

As discussed in Sec. I, $\gamma_{p}$ and $\delta$ can also be obtained using the mechanical route. To find $\gamma_{p}$ and $\delta$, we first evaluate $\gamma_{s}$ using Eq. 8, where we set $P^{\alpha}=P_{L}$ and $P^{\beta}=0$ since the vapour pressure in our simulations is negligible. In Fig. 8a we show isotherms of $\gamma_{s}$ as a function of $R_{s}^{-1}$, where $R_{s}$ is obtained from Eq. 9. We see that $\gamma_{s}$ decreases with increasing $R_{s}^{-1}$ along isotherms, indicating that $\delta$ is positive. Fitting these isotherms with Eq. 2 yields curves from which $\gamma_{p}$ is estimated from the intercept at $R_{e}^{-1}=0$. Our results for $\gamma_{p}$ obtained in this way are shown in Fig. 5b.

Another way of evaluating $\gamma_{s}$ is through Eq. 11. Isotherms of $\gamma_{s}$ from Eq. 11 as a function of $R_{s}^{-1}$, where $R_{s}^{-1}$ is estimated using Eq. 12, are shown in Fig. 8b. Although the trend seems to indicate that $\gamma_{s}$ decreases with $R_{s}^{-1}$ using Eqs. 11 and 12, the noise resulting from subtracting $P_{N}(r)$ and $P_{T}(r)$ prevents useful fitting of $\gamma_{s}$. The solid curves shown in Fig. 8b are simply the fits taken from Fig. 8a, and show a general consistency between using Eqs. 8 and 9, and using Eqs. 11 and 12, with the former set suffering from less statistical scatter.

Unlike both Eqs. 8 and 11, which require the determination of $R_{s}$ to evaluate $\gamma_{s}$, Eq. 13 does not involve calculating $R_{s}$. We plot $\gamma_{s}$ obtained from Eq. 13 as a function of $R_{s}^{-1}$ in Fig. 9a. We choose $R_{s}$ from Eq. 12 because $\gamma_{s}$ in Eq. 13 is derived from Eq. 11. The absence of $R_{s}$ in Eq. 13 seems to suppress the noise from $P_{N}(r)$. We fit the isotherms in Fig. 9a to Eq. 2, and we obtain values of $\gamma_{p}$ and $\delta$ similar to those obtained from the isotherms in Fig. 8a, as shown in Fig. 5b and Fig. 3.

To avoid any difficulty inherent in calculating $R_{s}$, another way of representing $\gamma_{s}$ is as a function of $R_{e}^{-1}$, as shown in Fig. 9b. Regardless of which variant of the mechanical route is taken, we observe that $\gamma_{s}$ decreases as $R_{e}$ and $R_{s}$ decrease, $\delta$ is positive with little evidence for a dependence on $T$, and $\gamma_{p}$ increases as $T$ decreases. The values of $\delta$ obtained from each variant of the mechanical route are shown for each $T$ in Fig. 3.

V Comparison of thermodynamic and mechanical routes

Fig. 5b shows that the mechanical route yields values of $\gamma_{p}$ approximately 10% larger than does the thermodynamic route, though the trends with $T$ are similar. The mechanical route values for $\gamma_{p}$ are also generally higher than those from Vega and de Miguel Vega2007 . However, the best agreement with Ref. Vega2007 between $T=240$ K and 300 K comes from using Eq. 13 for calculating $\gamma_{s}$, Eq. 12 for calculating $R_{s}$, and fitting the resulting isotherms (shown in Fig. 9a) with Eq. 2 to obtain $\gamma_{p}$ and $\delta$.

For all variants of both routes, $\gamma_{p}$ shows a striking departure from the extrapolated low $T$ behaviour presented in Ref. Vega2007 . The sharper than expected increase in $\gamma_{p}$ with decreasing $T$ occurs between 220 and 240 K, and is therefore consistent with crossing the Widom line at 230 K for bulk TIP4P/2005 at ambient pressure vega2010 . In order for bulk properties of the liquid to influence $\gamma_{p}$ as obtained from nanodroplets, it is reasonable to expect that nanodroplet interiors are structurally similar to the bulk, an expectation for which we provide evidence in Section VI.

For an independent comparison of $\gamma_{s}$, we show in Fig. 7 our results for $\gamma_{s}$ as a function of $R_{e}$ as obtained from both the thermodynamic and mechanical routes at $T=300$ K along with the values from Lau et al. Lau2015 obtained using the test-area method at $T=293$ K. We see that our results for $\gamma_{s}$ from the thermodynamic route are consistent with Lau, et al. However, the mechanical route gives significantly smaller values of $\gamma_{s}$. Smaller values of $\gamma_{s}$ and $R_{s}$ and larger values of $\delta$ for the mechanical route are also observed in nanodroplets interacting through the Lennard-Jones potential studied by Thompson, et al Thompson1984 .

While Fig. 7 shows consistency in the value of $\gamma_{s}$ and its $R_{e}$ dependence between our thermodynamic route and the test-area method employed by Lau et al. Lau2015 , other studies have found that $\gamma_{s}\approx\gamma_{p}$ for droplets as small as approximately 40 molecules ($R_{e}\approx 0.6$ nm). These studies employed techniques including excision of spherical portions from a bulk liquid Samsonov2003 , a volume perturbation method allowing for a thermodynamic determination of the pressure tensor components Ghoufi2011 , and a mitosis method by Joswiak et al. Joswiak . Interestingly, Lau et al. also carried out a mitosis method in a different study Lau2015-2 (also finding that $\gamma_{s}$ depends more weakly on $R_{e}$), and offered some discussion on the disparity between the mitosis and test-area methods. All these other studies point to the validity of approximating $\gamma_{s}$ with $\gamma_{p}$ in estimating the Laplace pressure for very small nanodroplets.

To compare the difference in $\gamma_{s}$ as obtained from the mechanical and thermodynamic routes, we plot in Fig. 6 isotherms of $\gamma_{s}$ as a function of $R_{e}$. Fig. 6b shows a significant change in $\gamma_{s}$ obtained from Eq. 13 as droplet size varies. For a change in the nanodroplet radius from 1 to 3 nm, there is a 50% increase in $\gamma_{s}$ at $T=180$ K, and 44% at $T=300$ K. However, if we compare this with $\gamma_{s}$ estimated from the thermodynamic route shown in Fig. 6a, we see that the isotherms are almost flat for $T\geq 220$ K, while there is only a 15% difference in $\gamma_{s}$ across the droplet size range for $T\leq 200$ K. We also can see that $\gamma_{s}$ from the thermodynamic route is systematically larger than the mechanical route, which is once again consistent with Thompson, et al Thompson1984 .

As shown by Sampayo et al. Sampayo2010 , disparity between thermodynamic and mechanical routes arises when energy fluctuations (and not merely the average change in energy) become important in determining the free energy change when surface area is increased. While for planar interfaces such fluctuations do not contribute significantly to the surface tension, the contribution for small droplets can be significant. Thermodynamic routes include these fluctuations, while mechanical routes would require the addition of hypervirial terms to include the effects of fluctuations. Ignoring these fluctuations in the mechanical route (i.e., when using pressure tensor components only) can lead, for example, to a change in sign in the determination of $\delta$ for Lennard-Jones droplets. See also Ref. Malijevsky2012 for discussion.

VI Local structure ordering

To quantify the structure of the interior of our water nanodroplets, we calculate the distance $d_{5}(r)$ between a molecule located at a distance $r$ from the centre of the droplet and its fifth-nearest-neighbour molecule (using distances between centres of mass). A large value of $d_{5}(r)$ indicates that molecules tend to be four-coordinated, i.e. that the local tetrahedral network is well formed ssp2001 .

In Fig. 10a we show $d_{5}(r)$ over a wide range of $N$ and $T$. We observe that $d_{5}(r)$ for droplet size $N=100$ is small and stays rather constant with $T$. The low value of $d_{5}$ indicates a collapse of the second neighbour shell around each molecule. This collapse is characteristic of the HDL form of water. The absence of any change in $d_{5}(r)$ with $r$ as we approach the surface indicates a disturbance in the tetrahedral network in the whole droplet. The overlap of the curves at different $T$ for $N=100$ suggests that droplets at this small size remain HDL-like both in the interior and at the surface regardless of how deeply we supercool them.

As we increase the droplet size to $N=360$, the profiles systematically shift to higher value of $d_{5}$ in the interior as we cool to 180 K. This change is a signature of a crossover from HDL at high $T$ to LDL at low $T$. However, for $T\leq 220$ K and $N=360$ in Fig. 10a, there is a decrease in $d_{5}$ going from interior to surface, which indicates a disturbance of the tetrahedral network and an increase in density at the surface. For larger droplets, such as $N=776$, we see similar behaviour as for $N=360$, but the transformation spans a wider range of $d_{5}$. Moreover, for $N=776$ at $T=180$ K we see a monotonic decrease in $d_{5}$ as we approach the surface. This may reflect the emergence of structural transformation within the droplet. The same scenario presents itself for $N=1440$. At $T=180$ K, as $N$ increases from 100 to 1440, $d_{5}$ in the interior monotonically increases with $N$. This indicates that as $N$ increases, a better LDL forms in the interior of the droplets.

To further probe the ordering inside the nanodroplets, we compute the local tetrahedral order parameter qt ,

[TABLE]

where $\psi_{jik}$ is the angle between an oxygen atom $i$ and its nearest neighbour oxygen atoms $j$ and $k$. Subsequently, we define $q_{T}(r)$ as the average value of $q_{i}$ for all molecules within a spherical shell bounded by radii $r\pm\Delta r/2$, where $\Delta r=0.05$ nm.

We show how $q_{T}(r)$ changes in Fig. 10b. We see that $q_{T}$ is low for $N=100$ and it increases as we cool the droplet. Similar behaviour appears for $N=360$, 776, and 1440. However, for $T=200$ and 180 K, the increase in $q_{T}$ upon increasing $N$ becomes quite dramatic, supporting the suggestion that a better tetrahedral network forms as $N$ increases. For $N=1440$ at $T=180$ K, the core reaches 90% of perfect tetrahedral order. The monotonic decrease in $q_{T}$ with $r$ for $N=776$ and $1440$ at $T=180$ K is consistent with the decrease with $r$ that we observe in $d_{5}$.

Our results for $d_{5}$ and $q_{T}$ suggest the progressive formation of LDL-like structure in the interior of our droplets as $T$ decreases. The structural results presented here are also consistent with the evolution of the density profiles of our droplets presented in Ref. NATURECOMM . The transformation from HDL to LDL in the droplet interior can also explain the change in behaviour in $\gamma_{p}$ for $T<240$ K shown in Fig. 5. Since LDL is a more structured liquid than HDL, with a better formed hydrogen bond network, we expect that the interface between LDL and the vapour phase will have a higher surface tension than for the interface between HDL and the vapour.

To illustrate the structurally bulk-like character of our droplet interiors, we plot $d_{5}$ and $q_{T}$ as functions of density in Fig. 11. To compute the density, we define the density within the core of our droplets as $\rho_{c}=m\langle{\cal N}/{\cal V}\rangle$, where ${\cal N}$ is the number of O atoms within a defined core radius $r_{c}=0.5$ nm of the droplet centre, ${\cal V}$ is the total volume of the Voronoi cells for these atoms voronoi , and $m$ is the mass of a water molecule. Since in the smallest droplets surface effects extend closer to the centre of droplet, we use $r_{c}=0.25$ nm for $N\leq 205$. Similarly, we define $d_{5}$ and $q_{T}$ for droplet interiors by averaging the corresponding local quantities for particles within $r_{c}$ of the droplet center (molecules for $d_{5}$ and O atoms for $q_{T}$). Fig. 11 shows the agreement between $d_{5}$ and $q_{T}$ as functions of density for bulk systems and droplets. This correspondence demonstrates that the core of the droplets for our range of $N$ is bulk-like. These structurally bulk-like interiors are consistent with the possibility that our extrapolated values of $\gamma_{p}$, obtained from the behavior of nanoscale droplets, approximate those for bulk planar liquid-vapour interfaces, and hence that the anomalous increase in $\gamma_{p}$ we observe below 230 K reflects the bulk liquid anomalies associated with crossing the Widom line of the LLPT.

VII Discussion and conclusions

We estimate the surface tension of water nanodroplets using the TIP4P/2005 model over a wide range of $N$ and $T$. We do so from an evaluation of the components of the pressure tensor inside the droplets NATURECOMM using the method described in Ref. malek2 . From the pressure tensor components, we determine the isotropic pressure $P_{L}$ in the interior of the droplets. This allows us to calculate the surface tension with two approaches: using the Young-Laplace equation directly, and using the variation of the pressure tensor components with distance from the droplet center. The direct route, which we call the thermodynamic route, requires $P_{L}$ and $R_{e}$ to estimate $\gamma_{s}$, $\gamma_{p}$ and $\delta$ as fit parameters, and the mechanical route evaluates $\gamma_{s}$ and $R_{s}$ from the pressure tensor components, and yields $\gamma_{p}$ and $\delta$ from fitting. It should be noted, however, that the analysis carried out by Gibbs (see e.g. Ref. Rowlinson1982 ) and reiterated by Tolman Tolman1949 , imply that the interior pressure that should be used in Eq. 1 is that of the bulk fluid with the same chemical potential as the droplet interior. It would be interesting to quantify the differences in the calculated surface tension that arise from using this definition instead of the directly-calculated pressure, particularly for smaller droplets where differences may be significant.

Isotherms of $P_{L}$ plotted as a function of $R_{e}^{-1}$ on the assumption that the surface of tension acts at $R_{e}$ (i.e. that $\delta=0$) show a linear dependence between $P_{L}$ and $R_{e}^{-1}$ that is valid for droplets as small as 0.86 nm in radius. To validate this apparent linearity, we insert the Tolman length correction into the Young-Laplace equation and find that $\delta$ is positive and small with a value of $0.055\pm 0.021$ nm. Moreover, $\gamma_{p}$ values for $T\geq 220$ K from this thermodynamic route, regardless of whether we assume $\delta$ is zero or not, are consistent with the extrapolation of $\gamma_{p}$ obtained for TIP4P/2005 using the test-area method Vega2007 , a thermodynamic method, as shown in Fig. 5a.

We compute $\gamma_{p}$ from the mechanical approach by first finding $\gamma_{s}$ and $R_{s}$ using Eqs. 8 and 9; then by using Eqs. 11 and 12, which produces consistent, but noisier results; and finally by using Eqs. 13 and 12. For our range of $T$ and $N$, we show that $\gamma_{s}$ decreases as $R_{s}$ decreases. Fitting these results with Eq. 2 results in positive and rather large values of $\delta=0.32\pm 0.02$ nm from Fig. 8a, and $\delta=0.21\pm 0.01$ nm from Fig. 9a. Although these two values do not overlap within error, they both suggest that $\delta$ from the mechanical route is significantly larger than the value from the thermodynamic route. Moreover, estimates of $\gamma_{p}$ obtained from fitting mechanical-route results tend to be higher than thermodynamic-route results, as apparent in Fig. 5b. However, if we consider $\gamma_{s}$ from Eq. 13 as a function of $R_{s}^{-1}$ as shown in Fig. 9a, the $\gamma_{p}$ values resulting from fitting with Eq. 2 are consistent with the thermodynamic route and with Vega and de Miguel’s extrapolation for $T\geq 240$ K.

We also conclude that $\gamma_{s}$ from the thermodynamic route remains relatively constant as we vary $R_{e}$ for $T\geq 220$ K, but shows larger variation at $T=200$ and 180 K, where it changes by 15% over the range of droplet sizes we use. In contrast, $\gamma_{s}$ from the mechanical route increases significantly with $R_{e}$, resulting in almost a 50% change in $\gamma_{s}$ at $T=180$ K. These results are equivalent to $\delta$ being small for the thermodynamic route and large for the mechanical route.

At 300 K, our thermodynamic results for $\gamma_{s}$ as a function of droplet size are consistent with those of Lau, et al. Lau2015 , while those from the mechanical route are not. One might conclude, therefore, that the mechanical route for determining $\gamma_{s}$ and $\delta$ lacks validity, and the relatively large value of $\delta=0.2$ – 0.3 nm should be rejected in favour of the smaller value of $\delta\approx 0.06$ determined from the thermodynamic route. However, as $\delta$ is the difference between $R_{e}$ and $R_{s}$, which is understood to be where the surface tension acts, values in the range of 0.2 to 0.3 nm are reasonable given the locations of $R_{e}$ and the negative pressure minima in Fig. 1. In sum, our work confirms the discrepancy between the mechanical and thermodynamic routes that has been previously noted in the literature, and so supports the need for a better theoretical understanding of the connection between the two.

The marked increase in $\gamma_{p}$ for $T<220$, as shown in Fig. 5, approximately coincides with the crossing of the Widom line at $T=230$ K for bulk TIP4P/2005 water at ambient pressure vega2010 , and hence, is correlated to the LLPT occurring in this water model. This increase in $\gamma_{p}$ is consistent across both the mechanical and thermodynamic routes. Our results thus confirm the scenario predicted theoretically in Refs. feeney ; hruby2004 ; hruby2005 , in which $\gamma_{p}$ increases more rapidly with decreasing $T$ when the system enters the $T$ regime below the Widom line where LDL-like properties begin to dominate the bulk behavior. We also note that Ref. hruby2004 predicts that the surface of a deeply supercooled water nanodroplet will exhibit a dense surface layer relative to the bulk-like density of the droplet interior. This prediction is confirmed by the density profiles presented in Ref. NATURECOMM , and is consistent with the radial variation of the structural properties presented in Section VI. The sudden increase in $\gamma_{p}$ at low $T$ that we infer from our droplet simulations was also observed in simulations of planar interfaces using the WAIL potential for water and was also interpreted as evidence for the LLPT scenario Rogers2016 .

Characterizing how local structure varies with radial distance from the center of the droplet with $d_{5}$ and $q_{T}$, we see behavior consistent with the formation of a well-ordered random tetrahedral network at low $T$ and large $N$ within droplet interiors. Furthermore, the dependence of these structural measures on local density match that of bulk TIP4P/2005 water. Hence, from a structural perspective, the interiors of our nanodroplets are characteristic of the bulk.

We conclude that $\gamma_{s}$ and $R_{s}$ determined from the mechanical route are smaller than the values evaluated in the thermodynamic route, leading to larger values of $\delta$ and $\gamma_{p}$. However, both routes give a positive value of $\delta$ for our range of $T$ and $N$, and suggest that $\delta$ is independendent of $T$. Moreover, assuming the validity of thermodynamic route, for $R_{e}\geq 1$ nm we can ignore the curvature correction and use the planar surface tension to estimate the Laplace pressure inside water nanodroplets to within 15% down to 180 K. This last point is of practical importance for the estimation of the interior pressure in real water nanodroplets, for which the Laplace pressure is not easily measured directly.

Acknowledgements.

ISV and PHP thank NSERC for support. PHP also acknowledges support from the Dr. W. F. James Research Chair Program. Computational resources were provided by ACENET and Compute Canada.