Water droplet excess free energy determined by cluster mitosis using guided molecular dynamics

Atmospheric aerosols play a vital role in a ff ecting climate by influencing the properties and lifetimes of clouds and precipitation. Understanding the underlying microscopic mechanisms involved in the nucleation of aerosol droplets from the vapour phase is therefore of great interest. One key thermodynamic quantity in nucleation is the excess free energy of cluster formation relative to that of the saturated vapour. In our current study, the excess free energy is extracted for clusters of pure water modelled with the TIP4P / 2005 intermolecular potential using a method based on nonequilibrium molecular dynamics and the Jarzynski relation. The change in free energy associated with the “mitosis” or division of a cluster of N water molecules into two N / 2 sub-clusters is evaluated. This methodology is an extension of the disassembly procedure used recently to calculate the excess free energy of argon clusters [H. Y. Tang and I. J. Ford, Phys. Rev. E 91 , 023308 (2015)]. Our findings are compared to the corresponding excess free energies obtained from classical nucleation theory (CNT) as well as internally consistent classical theory (ICCT). The values of the excess free energy that we obtain with the mitosis method are consistent with CNT for large cluster sizes but for the smallest clusters, the results tend towards ICCT; for intermediate sized clusters, we obtain values between the ICCT and CNT predictions. Furthermore, the curvature-dependent surface tension which can be obtained by regarding the clusters as spherical droplets of bulk density is found to be a monotonically increasing function of cluster size for the studied range. The data are compared to other values reported in the literature, agreeing qualitatively with some but disagreeing with the values determined by Joswiak et al. [J. Phys. Chem. Lett. 4 , 4267 (2013)] using a biased mitosis approach; an assessment of the di ff erences is the main motivation for our current study. C 2015 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution 3.0 Unported License. [http: // dx.doi.org / 10.1063 / 1.4935198]


I. INTRODUCTION
Atmospheric aerosols play a key role in climate, weather, pollution, and human health. These substances can influence the optical properties of clouds as well as their lifetimes and consequently have an effect on the distribution of precipitation and the global radiation budget (the difference between solar energy accumulated by the earth and the energy radiated into space). Aerosols can give rise to urban pollution episodes 1,2 and are thought to be responsible for significant morbidity in vulnerable populations. [3][4][5][6] There is therefore strong motivation to understand the underlying molecular mechanisms involved in atmospheric nucleation, namely, the formation of condensed phase molecular clusters from the ambient metastable vapour phase. 7-10 Significant progress has been made in recent years in identifying the species responsible for the initial stages of aerosol formation in the atmosphere [11][12][13][14][15] but the task of properly characterising the formation mechanism in detail still remains, particularly understanding the population dynamics and thermodynamics of the molecular clusters in question.
The most important quantity in nucleation modelling within a traditional kinetic/thermodynamic framework is the excess free energy which is often interpreted as the cost in free energy associated with the formation of the interface between the two phases. In classical nucleation theory (CNT) of a liquid cluster forming from a supersaturated vapour, the excess free energy F s is one of the two terms that comprise the work of cluster formation φ, where N is the number of molecules in the cluster and ∆µ is the difference in chemical potentials (which is related to the vapour supersaturation) between the condensed and vapour phases. The second term is the "bulk" contribution corresponding to the thermodynamic driving force ∆µ. In traditional CNT, the excess free energy is represented by the capillarity approximation, namely, F CNT s = γ ∞ A, where γ ∞ is the vapour-liquid interfacial tension and A is the surface area of a cluster regarded as a spherical liquid drop with a constant density matching that of the bulk condensed phase. A more general framework extending CNT can be devised from underlying microscopic models, such as the introduction of a size-dependent surface tension, 16 allowing for models which do not assume the vapour to be ideal, 17,18 and mean-field kinetic theories. 19 The CNT framework leads to a relationship for the nucleation rate J which has a sensitive dependence on the surface tension; where K is the kinetic prefactor that includes the aggregation rate of particles to the cluster, k B is the Boltzmann constant, T is the temperature, and v l is the volume per particle in the condensed phase.
In fact, the surface tension of curved interfaces is known not to be constant, and its size-dependence has been the focus of recent attention. 20,21 A variety of different methods have been employed for this purpose, ranging from an evaluation of the components of the pressure tensor [22][23][24] and thermodynamic perturbation methods [25][26][27] to techniques involving a finitesize scaling analysis in the grand canonical ensemble. [28][29][30] Consistency in the representation of the interfacial properties has not yet emerged in all cases. In some studies of water droplets, for example, the tension is found to be a monotonically increasing function of the droplet size until the planar limit is reached, but even in cases where the qualitative trend is the same, the cluster size where the tension essentially reaches the planar value differs from study to study. [31][32][33][34] It should be noted, however, that different models of water have been employed, and furthermore, the various studies were conducted under slightly different conditions, so that a direct comparison is not always possible. In stark contrast, Joswiak et al. 35 have predicted that the specific surface free energy of clusters of water molecules for the popular TIP4P/2005 model 36 will increase monotonically from the planar limit with decreasing cluster size. These authors employ the "mitosis method" using a hybrid molecular dynamics (MD)/Monte Carlo (MC) technique involving umbrella sampling to evaluate the free energy associated with splitting a water cluster into two equal sized sub-clusters. This free energy is then related to the vapour-liquid interfacial tension of the clusters. It is this recent investigation which motivates our current study.
Our aim is to re-assess the unusual nature of the cluster-size dependence obtained by Joswiak et al. 35 using the same TIP4P/2005 molecular model, a similar mitosis procedure, but a different fundamental approach to computing the thermodynamic properties. In our current work, we use a combination of nonequilibrium molecular dynamics (NEMD) and the Jarzynski equality 37 to extract the excess free energies. The Jarzynski equality is a simple but powerful formula relating an equilibrium change in the free energy to the mechanical work performed during the nonequilibrium simulations. We carry out the mitosis by dynamically driving the separation process and relating the extracted difference in free energy to the excess free energy of the initial and final clusters. The vapour-liquid interfacial tension can then be estimated by dividing the excess free energy by an area commensurate with the representation of the clusters as spherical drops with a constant density equal to that of the bulk condensed phase.
In Sections II A-II C, we describe the mitosis method and how the corresponding work calculated from the simulations can be related to the free energy difference through the Jarzynski equality. An analysis relating the free energy of mitosis extracted from the simulations to the excess free energy of the free cluster taking into account the perturbation due to the presence of the mechanical driving forces is made in Sections II D and II E. Further details of the analysis are provided in Appendix A. In Section III, we outline the procedure employed for the NEMD simulations. The results and discussion are presented in Section IV, and the conclusions of our work are given in Section V.

A. Guided molecular dynamics
The method we employ for evaluating the cluster free energy involves the disassembly of a cluster via nonequilibrium MD simulations, originally developed for studies of argon clusters. 38 The technique involves tethering every cluster molecule to an artificial "guide" particle using a harmonic potential. The guide particles are then arranged to move apart at constant speeds, driving cluster disassembly until a desired final state is reached. In our current study, the cluster is separated into two sub-clusters of equal size, corresponding to a mitosis process (cf. Figure 1). The strength of the harmonic potential can be varied throughout the course of the simulation. This allows the cluster to be initially weakly affected by their presence, ensuring minimal perturbation of the initial cluster structure. Towards the end stages of the simulation, the springs can be strengthened to ensure that the mitosis is completed. The work done during the mitosis process can then be related to the change in free energy.
It is important to mention that the computation of the change in free energy associated with the separation of a bulk FIG. 1. The cluster mitosis process. Guide particles and the cluster molecules are indicated in grey and blue, respectively. Initially, the cluster molecules are weakly tethered to the guide particles placed at the origin. The guide particles are then made to drift apart until the cluster is separated into two distinct sub-clusters. liquid into two (or more) liquid slabs separated by vapour has a long history. In the early days of the development of molecular simulation techniques, Mijazaki et al. 39 determined the excess surface free energy of a Lennard-Jones (LJ) fluid with a thermodynamic procedure involving the separation of a homogeneous liquid into two slabs separated by two vapour-liquid interfaces; a step-wise procedure was employed by introducing an external hard-wall potential to decouple the liquid stabs which were then allowed to relax to form the free interfaces. The vapour-liquid surface tension could then be estimated from the cost in free energy associated with the formation of the interfaces.

B. Work of mitosis
In general, when the Hamiltonian H is a function of an externally controlled parameter λ, the work done W on a system over some sampling time interval τ is For a cluster of N molecules tethered to their corresponding guide particles, the Hamiltonian consists of intermolecular interactions between the cluster molecules as well as harmonic interactions between the molecules and their guides. The total Hamiltonian is therefore a function of the time-dependent spring constant κ(t), the molecular positions {x i (t)}, and the guide positions {X i (t)}. The total work can be expressed as where dX i /dt is the velocity of the guide particle tethered to molecule i. The first term corresponds to the work done in tightening or loosening the spring constant over the course of the simulation, and the second term is the work done on the molecules by the moving guide particles through the springs.

C. Jarzynski equality
The work performed in the simulations needs to be related to the change in free energy between the initial and final states. According to the first and second laws of thermodynamics, the work done on a system as its Hamiltonian evolves is greater than or equal to the corresponding change in Helmholtz free energy: W ≥ ∆F. Only in the case of an infinitesimally slow and reversible (quasistatic) process will the work be equal to the free energy difference between the initial and final states. However, running simulations at the quasistatic limit in molecular dynamics is not practicable. On the other hand, one can perform many fast process simulations and estimate the change in free energy from the average of the work performed. Nevertheless, as a consequence of the second law, the outcome can only be an upper bound to the free energy difference: ⟨W ⟩ ≥ ∆F.
One solution to this problem is provided by the approach of Jarzynski, 37 who derived a nonequilibrium relationship, which is commonly known as the Jarzynski equality, where the difference in free energy between the initial and final states can be related exactly to an ensemble average of the Boltzmann factor of the work. The initial state needs to be represented in the canonical ensemble, but no conditions are placed on the statistics of the final state after the work process.
Since the Jarzynski equality can be used to relate the change in free energy to the work performed in nonequilibrium processes (far from the quasistatic limit), it is particularly suitable for use in conjunction with MD simulations. It is important to note, however, that because of limited sampling, the computed ∆F may depend on the rate of processing: slower processes will give rise to a narrower work distribution and hence an estimate of the change in free energy with lower statistical error than that derived from faster processes with a broader distribution of work.

D. Contributions to the excess free energy
The difference in free energy extracted from the MD simulations corresponds to that associated with the mitosis of tethered clusters. This quantity can be related to the excess free energy for the mitosis of a free (undistorted) cluster using a statistical thermodynamic analysis (details are given in Appendix A). We find that the difference in excess free energy between two N/2 sub-clusters and the parent N cluster (with N even) can be written as a sum of four terms, with where ρ vs is the saturated vapour density, κ i is the initial spring constant, v HO = (k B T/κ i ) 3/2 , and v l is the volume per particle in the condensed phase. The first term f 1 s corresponds to the change in free energy associated with the mitosis of a tethered cluster, as obtained from MD simulations. The second term f 2 s accounts for the correction that relates calculations derived from a MD simulation with distinguishable particles to a system of indistinguishable particles. The third term f 3 s is an entropic contribution associated with correcting for the confinement of the centre of mass of the cluster in the initial state, and the centres of mass of the two clusters in the final state, through the tethering process. The final term f 4 s corrects for the perturbation in the cluster energy, before and after mitosis, caused by the presence of the tethers.

E. Surface tension from mitosis
The excess free energy obtained from the mitosis procedure can be related to the surface tension. A similar analysis has been previously performed by Joswiak et al. 35 For a system with a planar interface, the thermodynamic surface tension is defined as the excess free energy per unit area γ = F s /A for an appropriate choice of dividing surface. For small liquid clusters, linking the excess free energy associated with the formation of the surface to the tension requires a similar analysis but requires the specification of an appropriate density of the condensed phase. In this study, we introduce a size-dependent specific surface free energy (tension) γ(N) while making the assumption that clusters are spherical (with a well defined radius and constant density corresponding to that of the bulk liquid phase), such that we can write the work of formation of a cluster of N particles from the vapour as the sum of a surface term and a volume term, where R(N) is the radius of the cluster and ρ l = 1/v l is the bulk number density of the condensed phase. The density corresponding to a bulk liquid phase at the same chemical potential as a coexisting vapour phase would be a more rigorous choice but it is more difficult to implement such a procedure; this would be consistent with the formal thermodynamic approaches of Gibbs 40 and Tolman. 16 We employ the aforementioned simpler approach to make qualitative comparisons of the size-dependent tension reported in the literature but are aware of the difficulties in assigning the macroscopic concept of a surface tension to small systems. To calculate the surface tension, Joswiak et al. used the expression of Tolman 16 relating the surface free energy to the surface curvature of large droplets. The Tolman expression can be useful for larger clusters but is expected to break down at smaller length scales and we opt not to employ it in our current work. With the excess free energy of an N molecule cluster written in the form of a surface term F s = 4π[R(N)] 2 γ(N) and employing the incompressibility condition R(N) = 2 1/3 R(N/2), the excess free energy extracted from the mitosis procedure can be expressed as In order to progress, we assume that the planar limit is valid for the largest cluster studied (of size N max ), allowing us to insert the value for γ(N max ) in Eq. (9). We check the validity of this assumption by ensuring that the calculated surface tensions of the N max cluster and the N max /2 cluster differ negligibly. The surface tensions of the smaller clusters can then be calculated using Eq. (9) from the excess free energy obtained from the simulation.

III. METHODOLOGY
Molecular dynamics simulations are carried out for the mitosis of clusters for N = 6, 12, 24, 48, 76, 96, and 128 TIP4P/2005 36 water molecules. Simulations are performed in the canonical ensemble in a nonperiodic box at a temperature of 300 K maintained using a Langevin thermostat 41 with a friction constant of 1 ps −1 . For each cluster, 1000 initial configurations are sampled from an equilibrium trajectory of duration 10 ns with the guide particles fixed at the origin. An initial spring constant κ i = 0.05 kJ mol −1 Å −2 is used, deemed to be minimally perturbing since the mean tether energy of a molecule for all cluster sizes is then less than the typical thermal energy k B T. The radii of the clusters are calculated based on the bulk liquid density (of the TIP4P/2005 model) at the prevailing temperature.
Starting from the initial configurations, half of the guide particles are made to drift at a constant speed in the positive direction parallel to one of the Cartesian axes while the other half drift in the negative direction. The spring constant is time-dependent, smoothly varying from the initial value where t i is the time when the spring constant starts changing and t f is the time when the spring constant reaches its final value. In our simulations, the values t i and t f are pre-set as the following fractions of the simulation time: t i = t sim /5 and t f = 4t sim /5 and we use a total simulation time of t sim = 10.24 ns. Further details on the justification of our particular choices of κ i and t sim are given in Appendix B. The guide particles drift until the two sets of guide particles are separated by a distance of 40 Å as illustrated in Figure 2. At this point, the guide particle velocities are reduced to zero and the sub-clusters are allowed to relax while the tether spring constants are loosened back to their initial value κ i over a further period of t sim .
The differences in free energy obtained from the simulations using the Jarzynski relation are converted to excess free energy following the analysis presented in Section II. The

IV. RESULTS
An example of the distribution of work obtained for the mitosis of the N = 48 cluster is depicted in Figure 3. The distributions for all clusters exhibit a smooth, near-Gaussian form suggesting that the statistics are good enough to extract accurate free energy differences. Values for the excess free energy for different cluster sizes obtained from a sum of the four contributions f i s are shown in Table I. The vapour-liquid surface tensions derived from the procedure described in Section II are also presented in Table I. The surface tension for the N = 96 cluster is set to that for the planar vapour-liquid interface (γ ∞ = 69.3 mN m −1 for this model including longrange corrections 43 ) and serves as an initial input in Eq. (9) for the surface tensions of the smaller clusters. The results for N = 76 and N = 128 are obtained using interpolated values for γ at N = 38 and N = 64 (derived from the results for N in the range 24-96) and using them as inputs in Eq. (9) to obtain γ at N = 76 and N = 128, since these sizes were those studied  by Joswiak et al. 35 The interpolated values for the tension are estimated from a cubic spline of the data points at N = 24, 48, and 96. A plot of the surface tension normalised by γ(N = 96) is depicted in Figure 4. It is apparent that this ratio is close to unity for N = 48 (there is a small difference of ∼2% in the value between N = 48 and N = 96), suggesting that the tension has almost reached the planar limit for this cluster size. This justifies our choice of setting the tension at N = 96 to the planar value although we are aware that, in principle, the planar value can only be reached in the limit N → ∞. The excess free energy F s as a function of cluster size N is also illustrated in Figure 4. In addition, we present curves of the excess free energy obtained from CNT and internally consistent classical theory (ICCT), 44 (1)], respectively. In these cases, the tension is assumed to be that of a planar vapour-liquid interface of the TIP4P/2005 system at 300 K (i.e., 69.3 mN m −1 ). Our values for the excess free energy are therefore consistent with CNT for the larger cluster sizes. For clusters of N = 3 and N = 6 particles, our results tend towards the ICCT curve. CNT is known to break down at small cluster sizes since it fails to predict a zero excess free energy at N = 1. ICCT addresses this inconsistency by shifting the CNT free energy curve down by a constant so that F s = 0 at N = 1. For N = 12 and N = 24, our values lie between the CNT and ICCT curves, indicating that there is an intermediate length scale corresponding to a transition from ICCT to CNT.
The vapour-liquid surface tension estimated for the drops of TIP4P/2005 water as a function of cluster size is plotted again in Figure 5 for comparison with results obtained in other studies. We observe the tension to be a monotonically increasing function of size that approaches the planar limit. The capillarity approximation, namely, the assumption that the clusters all have a tension equal to the planar tension, is surprisingly accurate at this temperature, at least down to cluster sizes of about N = 40 molecules. After this point, the tension drops off rapidly to about 60% of the planar value for our smallest system (three water molecules). Our results do not agree with those obtained by Joswiak et al. 35 who find that the tension increases from the planar limit upon decreasing the size of the cluster. However, we find that our results are in near-quantitative agreement with the values obtained independently by Samsonov et al., 32 by Ghoufi and Malfreyt, 31 and by Factorovich et al. 33 Samsonov et al. used the Stockmayer potential to represent water and equated the change in free energy of formation of a drop to the change in internal energy due to excision from the bulk liquid phase at 300 K. Ghoufi and Malfreyt performed mesoscale simulations of water at 298 K using a dissipative particle dynamics (DPD) model and employed a method to calculate the local normal and transverse components of the pressure tensor using volume perturbation; the curvature dependent surface tension was then computed from the pressure components. Factorovich et al. took a different approach by analysing the saturated vapour pressure over different sized drops of water represented with the mW coarse-grained model at 298 K. The latter authors found that the macroscopic Kelvin relation holds down to clusters of about 40 molecules, suggesting that the tension has little dependence on curvature until the dimensions become very small; the values of the tension reported for Factorovich et al. in Figure 5 are deduced from a rearrangement of the Kelvin equation.
In a recent study, 34 we applied the test-area (TA) method to TIP4P/2005 water clusters consisting of between N = 16 and N = 1000 molecules at 293 K. In contrast to the other studies discussed thus far, including our current work, the tension obtained with the test-area approach was found to depart from the planar limit at significantly larger cluster sizes. Nevertheless, the trend of a decreasing tension at small N is still in qualitative agreement with the other studies. At least at a temperature of T ≈ 300 K, there is qualitative agreement between most of the existing literature and our current work. The general consensus appears to indicate that the capillarity approximation only starts to break down for clusters with fewer than about 40 particles. Care should however be employed in the analysis of macroscopic thermodynamic properties such as the surface tension (or for that matter the average liquid density) for such small clusters.
The differences between the size dependence of the vapour-liquid interfacial tension for water obtained using the TA and mitosis methodologies may appear to be rather puzzling. We should, however, acknowledge that the two approaches characterise the free energy associated with a (nominal) surface area of the cluster in rather different ways.
In a macroscopic thermodynamic treatment, 46 the interfacial tension for a system in the canonical ensemble can be expressed as the derivative of the Helmholtz free energy with respect to the surface area, (∂F/∂ A) N,V,T . In the case of a system with a curved interface, the relation still holds but one has to be careful with the definition of the dividing surface. The surface of tension, useful for a mechanical analysis employing the differential surface tension, is distinct from the equimolar Gibbs dividing surface that most conveniently characterises the surface tension as a specific surface free energy per unit area. If the clusters are large enough, the radii defining these surfaces are similar (the Tolman length is small 21 ), and the definition of the surface area is not problematic with the differential and specific surface tensions not differing significantly.
However, for small clusters, an interfacial tension defined using the excess Helmholtz free energy F s (derived from a microscopic statistical mechanical analysis) will generally differ from the differential surface tension. In the macroscopic thermodynamic analysis, 47 the counterpart to F s is the surface grand potential Ω σ = γ E A E , where A E is the equimolar surface area and γ E is a specific surface free energy defined with respect to this area. Such an analysis tells us that the statistical mechanical approach should produce results consistent with the macroscopic surface tension for large enough clusters, though it unfortunately has little to say about the behaviour of the excess free energy at small sizes.
In the case of the TA method, we determine a quantity associated with deforming the area of a mechanically stable (long lived) cluster and relate it to the differential definition of the surface tension. It is a statistical mechanical approach based on an evaluation of (∂F/∂ A) N,V,T using a distortion of the metric of space as a proxy for differentiation with respect to a surface area. 34 We again expect the results to tend towards the macroscopic surface tension for large clusters, but for smaller sizes, we expect differences to arise. Furthermore, we do not expect coincidence between the resulting differential surface tension and the excess free energy divided by the area of equimolar surface. We properly regard these as distinct cluster properties.
The methods are most appropriate for different ranges of cluster size. Difficulties associated with the TA approach applied to very small drops include the lack of mechanical stability, stabilization effects due to keeping particles within a finite-size simulation cell (potentially distorting the statistics), uncertainties arising from the proportionality of the free energy of deformation to second order in the distortion parameter (large numerical fluctuations are seen in Figure 11 of Ref. 34), and also the non-sphericity of small clusters which can make the method impracticable. It may therefore not be advisable to employ the methodology for very small clusters.
In the case of the mitosis approach, we determine a well defined excess Helmholtz free energy of cluster formation F s , free of any continuum assumptions, even for small clusters of arbitrary geometry. The only concession to a continuum treatment is that in relating F s to a surface free energy per unit area, we impose that F s (N) = 4π[R(N)] 2 γ(N) and make an assumption of incompressibility of the condensed phase in order to define the equimolar radius R(N). Such a conversion is not necessary in an assessment of cluster stability and is only done for convenience of presentation. It is noted that the density at the interior of small LJ clusters can be found to be considerably larger than that of the corresponding bulk system, 20 indicating some degree of compressibility, though it is not clear that this will be an issue for water at the conditions investigated.
So while we can reasonably assume that the TA and mitosis methods will provide an equivalent estimate of the thermodynamic surface tension for large clusters, for small clusters we are clearly computing different thermodynamic quantities. For example, it is straightforward to show that the TA "surface tension" of a harmonically bound dimer is zero, while the excess free energy of the same cluster is not. Differences will persist for more complex interactions or larger clusters. We regard the TA surface tension as the quantity most appropriate for comparison with the thermodynamic surface tension and the excess free energy per unit area as the quantity most relevant for an accurate characterisation of cluster stability.
The results for the size dependence of various formulations of the interfacial tension, normalised by the planar value, are collected together in Figure 6 for both LJ and TIP4P/2005 systems. We acknowledge that the different studies are carried out at a variety of temperatures and cutoffs of the potential. Results obtained with the TA method for stable drops of LJ 27 and TIP4P/2005 34 water particles appear to follow the same trend (although they are not expected in general); the water drops are also found to be stable for smaller clusters as might be expected for a system with stronger intermolecular forces. The behaviour obtained with the disassembly and mitosis methods for the LJ 38 and TIP4P/2005 clusters, respectively, also displays close correspondence. There are clear differences between the surface tensions derived by the two approaches for both systems when cluster size is decreased.
Interestingly, results obtained from calculations with a fundamental measure theory (FMT), a non-local density functional theory (DFT), 20  for LJ drops, involving a finite-size scaling analysis within the grand canonical ensemble, is somewhere between those obtained with the TA and mitosis approaches. The various methodologies tend to agree for large systems. Unfortunately, we are not able to rationalise the apparent increase in the interfacial tension observed by Joswiak et al. although their results also tend towards the planar value in the large R limit. 35

V. CONCLUSIONS
We have developed a method based on mitosis by guided molecular dynamics and the use of the Jarzynski equality to extract the excess free energy and estimate the vapourliquid surface tension of water clusters. Simulations have been carried out for clusters of N = 6 to N = 128 TIP4P/2005 water molecules at a temperature of 300 K. The size-dependence of the excess free energy has been extracted and compared to the CNT and ICCT models. Our results are consistent with CNT for larger clusters (N ≥ 40) and with ICCT for the smallest clusters (N = 3 and N = 6). For clusters of intermediate size (N = 12 and N = 24), our values transition between the ICCT and CNT curves. The surface tension is observed to be a monotonically increasing function of cluster size, converging to the planar limit for clusters of N ≥ 40 molecules but reducing significantly in magnitude at the lower end of the size scale. Our results are not found to agree with those obtained by Joswiak et al., 35 who carried out a similar mitosis procedure using umbrella sampling; by contrast, these authors found that the tension increased from the planar limit with decreasing cluster size. On the other hand, our results are in qualitative agreement with other studies. [31][32][33][34] We conclude that mitosis simulation in nonequilibrium molecular dynamics appears to be a powerful and efficient method for the determination of excess free energy of small clusters.

APPENDIX A: THERMODYNAMICS OF TETHERED CLUSTER MITOSIS
The change in free energy extracted from the MD simulations is associated with the mitosis of tethered clusters. This quantity needs to be related to the excess free energy between an untethered (free) cluster of N molecules and two untethered sub-clusters of size N/2. We follow a procedure similar to that used recently for cluster disassembly. 38 We start by expressing the free energy of an untethered cluster of N indistinguishable particles as where h is Planck's constant, and H({x j }) is the Hamiltonian of the system, a function of the molecular positions x j and momenta p j . We can write Z F = V Z c F , where Z c F is the partition function of a free cluster with its centre of mass fixed at the origin, and V is the system volume. This is and recasting the equation in terms of particle coordinates with respect to the centre of mass x c . The partition function becomes where x ′ j = x j − x c are the transformed particle coordinates. The free energy of a tethered cluster is In a similar fashion to the case of the free cluster, the partition function of the tethered cluster can be written as where Z c T is the partition function of a cluster with the particles tethered to the origin as well as being constrained to have its centre of mass at the origin, Assuming that the mean tethering energy is less than the typical thermal energy scale k B T, the difference in free energy between the free and tethered clusters can be expressed as 38 where ρ N c (0) = (N κ i /2πk B T) 3/2 is the probability density function, evaluated at the origin, of the centre of mass of the N tethered particles in the cluster. ρ N 0 (y) represents the spatial density profile of a single particle in a free cluster of N particles with the centre of mass constrained to the origin.
The cluster excess free energy can be defined as where φ(N) is the work of formation of an N-cluster, ∆µ is the chemical potential difference between the condensed and vapour phases at a given vapour pressure, F F (N) is the Helmholtz free energy of an N-cluster, and µ s is the coexistence chemical potential. Eq. (A6) is subtly different from Eq. (1) in that the work of formation of a cluster minus that of a monomer is the appropriate quantity that emerges from careful development of nucleation theory. 45 Eq. (1) is commonly taken as an approximation in CNT. If the vapour phase is considered ideal, then the monomer Helmholtz free energy is F(1) = −k B T ln(V/Λ 3 ) and the coexistence chemical potential is µ s = k B T ln(ρ vs Λ 3 ), where ρ vs is the saturated vapour density, and Λ is the thermal de Broglie wavelength. Combining Eqs. (A6) and (A5), the excess free energy of an N-cluster can now be written as It should be recognized that the difference in free energy obtained from guided molecular dynamics simulations corresponds to a process involving distinguishable molecules. However, we require the free energy difference between a cluster consisting of N indistinguishable molecules and two sub-clusters which are distinguishable from each other but each consisting of N/2 indistinguishable molecules. To relate the change in free energy obtained from simulation to the desired quantity, we need to apply the usual classical corrections for indistinguishability. The partition function corresponding to an initial cluster of N indistinguishable molecules is Z N = 1/(N)!Z dist N = exp(−F F (N)/k B T), where the superscript "dist" denotes a system of distinguishable particles. The partition function corresponding to a final state of two distinguishable clusters of N/2 indistinguishable molecules is 1/((N/2)!) 2 (Z dist N /2 ) 2 . The free energy difference is thus ∆F }. However, the difference extracted from the molecular dynamics simulations is ∆F mit = −k B T ln{(Z dist N /2 ) 2 /Z dist N } and thus ∆F = ∆F mit + k B T ln[(N/2)! 2 /N!].
We now turn our attention to the excess free energy between an N cluster and two N/2 subclusters, which after grouping terms reduces to the following four contributions: We note that there is a cancellation of the terms involving the system volume V , as required in an expression representing excess free energies. The third term can be written as where v HO = (k B T/κ i ) 3/2 is the volume scale characterising the confinement of particles in the final harmonic potential. For the fourth term, we represent the single particle density profile as a step function ρ N 0 (y) = ρ l /N for 0 < y < r max . In this approximation, the cluster is imagined to be spherical with a constant molecular density ρ l and a radius r max . Inserting this into the fourth term gives where v l is the volume per molecule in the bulk condensed phase. Inserting Eqs. ( 3v l 4π which is the form given in Eq. (7) in the main text.

APPENDIX B: CONVERGENCE OF THE FREE ENERGY
An analysis of the convergence of the change in free energy associated with the mitosis of a tethered cluster of TIP4P/2005 water obtained from our MD simulations for the 6, 12, and 24 particle systems as a function of the guided separation time t sim is depicted in Figure 7. A value of t sim = 10.24 ns is seen to be sufficient to obtain a converged change in free energy without being too computationally expensive. The time scales are in line with what was found for the disassembly of argon clusters, 38 though one should note that in such a process longer times are expected than for the mitosis method.
Another possible issue with the calculation of the excess free energy with the mitosis methodology is sensitivity of the results to the value of the spring constant employed to tether the clusters. The excess free energy obtained to split a cluster of six TIP4P/2005 water molecules (into two clusters of 3 molecules) as a function of the initial value of the spring constant is illustrated in Figure 8; the largest variation between the values is found to be significantly smaller than the FIG. 7. The change in free energy associated with the mitosis of N = 6, 12, and 24 tethered water clusters obtained from MD simulations for a range of separation times.   .1k B T). The corresponding value of the excess free energy per unit area for the N = 3 cluster system is depicted in Figure 9 and found to be quite insensitive to the value of the initial spring constant with small deviations of ∼0.2 mN m −1 . These results suggest that κ i = 0.05 kJ mol −1 Å −2 is a suitable (and sufficient) choice for the value of initial spring constant. It is also noted that for lower spring constants, the clusters are not sufficiently long lived in some trajectories within the duration of the simulation, so those are not included in the computation of the change in free energy as shown in Table II.