Kinetic Monte Carlo and hydrodynamic modelling of droplet dynamics on surfaces, including evaporation and condensation
Abstract
We present a latticegas (generalised Ising) model for liquid droplets on solid surfaces. The time evolution in the model involves two processes: (i) Singleparticle moves which are determined by a kinetic Monte Carlo algorithm. These incorporate into the model particle diffusion over the surface and within the droplets and also evaporation and condensation, i.e. the exchange of particles between droplets and the surrounding vapour. (ii) Largerscale collective moves, modelling advective hydrodynamic fluid motion, determined by considering the dynamics predicted by a thinfilm equation. The model enables us to relate how macroscopic quantities such as the contact angle and the surface tension depend on the microscopic interaction parameters between the particles and with the solid surface. We present results for droplets joining, spreading, sliding under gravity, dewetting, the effects of evaporation, the interplay of diffusive and advective dynamics, and how all this behaviour depends on the temperature and other parameters.
I Introduction
Understanding and predicting the collective dynamics of fluid particles (atoms or molecules) over substrates is a crucial aspect of surface science. For instance, the rate of diffusion of the particles is a key controlling factor in many dynamical processes occurring on surfaces, such as spreading, chemical reactions and the growth of islands and epitaxial layers AlaNissila et al. (2002). Factors such as the temperature, the shape of the particles, the degree to which they adhere to the surface, how they adhere to each other and the particle density on the surface, greatly influence the rate of diffusion over the surface. Whether the motion over the surface is largely a singleparticle phenomenon or a collective process is fundamental to determining the structures that are formed on surfaces.
The situations of particular interest here concern the behaviour of liquid droplets on solid surfaces and the influence of evaporation and wettability. This is a common situation arising in a vast range of applications, such as in coating, printing and lubrication, making such systems important to study and understand De Gennes et al. (2013). During the drying process droplets can interact, joining and fusing in a coarsening process that involves hard to observe phenomena such as the evaporation and condensation of particles and the diffusion of particles over the surface. This behaviour is even more rich if the liquid has a tendency to dewet from the substrate. To describe situations like this, here we develop a coupled kinetic Monte Carlo (KMC) and hydrodynamic model for liquids on solid surfaces that incorporates the thermodynamics of evaporation and condensation, as well as the hydrodynamics of fluid flow. We assume a simple latticegas (generalised Ising) model for the underlying Hamiltonian that incorporates the microscopic interactions between the fluid particles and the surface. The benefit of modelling in this way is that although it is a simplified (coarse grained) model, it does correctly include much of the relevant fluid structure and thermodynamics and allows for the study of the interfacial properties of the model system. For example, it is fairly straightforward to find the wetting or dewetting behaviour of the liquid and the connection to the microscopic properties of the system such as the particle interaction energies.
The use of lattice models to study equilibrium fluid wetting behaviour goes back to the 1980s Pandit and Wortis (1982); Pandit et al. (1982); Binder and Landau (1988) and they have been extremely useful for building our understanding of the physics. More recently, dynamical lattice models for nonequilibrium droplets and liquid films on surfaces have been developed Rabani et al. (2003); Sztrum et al. (2005); PauliacVaujour et al. (2008); Vancea et al. (2008); Blunt et al. (2010); Stannard (2011); Vancea (2011); Kim et al. (2011); Jung et al. (2014); Crivoi and Duan (2014); Tewes et al. (2017); Chalmers et al. (2017a), with the focus largely being on modelling and understanding the patterns formed by the drying of particle suspensions, i.e. the coffee ring stain effect, etc. Some of these models solely implement a dynamics governed by the Metropolis Monte Carlo algorithm (of course, based on the model Hamiltonian), whilst other models add other kinds of particle moves, with the aim being to incorporate the influence of the fluid hydrodynamics, which leads to particle advection. For example, the model in Kim et al. (2011) biases particle moves depending on the distance to the nearest contact line and the model in Crivoi and Duan (2014) assumes a certain form for the fluid velocity field and uses this to bias particle moves. Such moves do not satisfy detailed balance, which is required if one is to use Monte Carlo methods to simulate the dynamic properties of Hamiltonian systems Fichthorn and Weinberg (1991). However, it is not clear to us that one should enforce detailed balance when constructing coarse grained models of such nonequilibrium systems, since droplets on surfaces are highly dissipative systems, due to the evaporation, the strong coupling to the substrate, etc. Of course, at equilibrium, there should be detailed balance in the model.
To model collective hydrodynamic motion in our model, we generate largescale collective particle moves by considering the dynamics that is predicted by a thinfilm equation. This is a time evolution equation for the liquid film thickness profile , i.e. the height of the liquidvapour interface above the surface De Gennes et al. (2013); Oron et al. (1997); Kalliadasis and Thiele (2007); Thiele (2007); Craster and Matar (2009). This equation is derived from the NavierStokes equations for a film of liquid on a surface by making a longwave approximation, which greatly simplifies the analysis. When generating these ‘thinfilm moves’, we first determine the liquid film height profile from the configuration of occupied lattice sites, we then evolve the thinfilm equation for a short amount of time, and then map back to the lattice. In our model, an important quantity is the parameter , defined below in Eq. (13), which is proportional to the ratio of the number of KMC attempted particle moves to the number of thinfilm moves. When , our model reduces to just evolving the thinfilm equation, whilst in the limit , our model is a pure KMC model. Therefore, this parameter , where is a singleparticle diffusion coefficient and is the viscosity.
Key quantities that characterises how a liquid wets a surface are the spreading parameter , where , , and are the solidvapour, solidliquid and liquidvapour interfacial tensions, respectively De Gennes et al. (2013), and also the contact angle , which is given by Young’s equation De Gennes et al. (2013):
(1) 
The contact angle is the angle that the height profile for equilibrium droplets makes with the surface. is a parameter that must be input into the thinfilm equation, whilst is determined by the disjoining pressure , which is a function that must also be input into the thinfilm equation. All these quantities depend on the fluid state point (i.e. the temperature and the pressure) and also on the nature of the interactions between the fluid particles and with the substrate. Here, we use the approach of Refs. Archer and Evans (2011); Hughes et al. (2015) (see also Hughes et al. (2017); Buller et al. (2017)) that uses classical density functional theory (DFT) Evans (1979, 1992); Hansen and McDonald (2013) to calculate the binding potential , from which one obtains , and also the interfacial tensions. Thinfilm equations together with binding potentials from DFT have been used previously to elucidated the influence of the fluid microscopic structure on droplet spreading Yin et al. (2017). Here, we instead use the thinfilm equation with the quantities and determined from DFT to generate collective multiparticle moves for our lattice model.
There are, of course, other approaches that one can take to model the behaviour of liquid droplets on surfaces: Fully microscopic (particle resolved) approaches include (i) applying molecular dynamics computer simulations to solve Newton’s equations of motion for all the individual atoms or molecules in the fluid – see e.g. Ref. Tretyakov et al. (2013), or (ii) Monte Carlo simulations – see e.g. MacDowell et al. (2002). As mentioned above, (iii) microscopic DFT can be applied, such as in Ref. Hughes et al. (2015), or DFT can also be applied to develop coarsegrained mesoscopic theories such as those described in Refs. Thiele et al. (2009); Robbins et al. (2011); Chalmers et al. (2017b). Thinfilm equation based models are also widely used De Gennes et al. (2013); Oron et al. (1997); Kalliadasis and Thiele (2007); Thiele (2007); Craster and Matar (2009) and if thermal fluctuations are important, these can also be incorporated, creating stochastic thinfilm equation models Grün et al. (2006); DuránOlivencia et al. (2019). Macroscopic scale approaches include just directly solving the NavierStokes equations or via methods such as latticeBoltzmann Krüger et al. (2017). However, when implementing such methods it is hard to connect to microscopic system properties (i.e. the underlying Hamiltonian), as we do here.
This paper is structured as follows: In Sec. II, we introduce the lattice Hamiltonian that defines the energetics the model. In Sec. III, the Monte Carlo singleparticle dynamics algorithm and the thinfilm evolution equation used to describe the effects of collective advective droplet dynamics over solid substrates are described. In Sec. IV, the DFT based statistical mechanical theory used for calculating the binding potential and other thermodynamic quantities is presented. We present results for the bulk fluid phase diagram, the liquidvapour interfacial tension and the binding potential, showing how all of these depend on the temperature, on the fluid interaction parameters and the parameter governing the strength of the attraction of the fluid to the wall. In Sec. V, we present equilibrium droplet properties, including calculating droplet density profiles using KMC and comparing the contact angle of droplets with the predictions from DFT. In Sec. VI, we present results from simulating nonequilibrium droplets on surfaces. This includes droplets joining via evaporation and condensation process, droplets sliding together and joining and also droplets sliding under lateral driving, e.g. due to gravity. Finally, in Sec. VII, we draw our conclusions.
Ii Lattice model and Hamiltonian
In Fig. 1 is illustrated the system of interest and how it is modelled here. It consists of a planar solid surface upon which are a collection of interacting particles, with the majority gathered together forming a liquid droplet. We consider a Cartesian coordinate system with the  and axis being parallel to the solid surface and the axis being perpendicular to the solid surface. The system is discretised onto a regular grid. To simplify, here we treat the system as being twodimensional (2D), so that there are no variations in the direction, but the model can easily be generalised to three dimensions (3D). We define to be the lattice spacing, which we set to , throughout. We define to be the occupation number of lattice site , where is the 2D discrete position vector. We assume that and , with corresponding to the solid surface. When then site is occupied by a particle and when it is unoccupied, then . We typically initiate the liquid in one or two semicircular droplet configurations on the surface.
Note that one need not necessarily assume that each lattice site is of the same size as the individual molecules. One can also consider the model to be a coarse grained model, with being much larger than the diameter of the molecules. In this case, when we refer to a lattice site as being ‘occupied by a particle’ we mean that it contains mostly liquid and when a lattice site is referred to as ‘unoccupied’ or ‘not containing a particle’, then we mean that it mostly contains vapour.
The total energy of the system is modelled by the Hamiltonian
(2) 
where the first term corresponds to a sum over pairs of lattice sites, and the overall strength of the interaction between particles is defined by the parameter . The interaction energy between pairs of particles at lattice sites and depends on the distance between them and is given by . Here, we use
(3) 
where and denote the nearest neighbours of and the next nearest neighbours of , respectively. The choice of values in Eq. (3) is important, since with these any equilibrium droplets that form on the surface tend to have the shape in the form of a segment of a circle Chalmers et al. (2017a); Archer et al. (2010); Robbins et al. (2011). If, for example, one were to assume only nearest neighbour interactions, i.e. with , for , then rectangular shaped droplets are liable to be formed, particularly at low temperatures. The choice of in Eq. (3) minimises the dependence of the vapourliquid surface tension on the interfaceorientation. The values of in Eq. (3) come from considering how to discretise the Laplacian whilst minimising errors due to the discretisation Archer et al. (2010); Robbins et al. (2011).
In the second term in Eq. (2), is a sum over all the lattice sites and is the contribution to the potential energy from the external potential due to the surface on which liquid is deposited. We model it as follows:
(4) 
where is the parameter which determines the interaction strength between the particles and the solid surface (the wall). With this potential, the boundary condition for the particles at the wall is straightforward: we set for all lattice sites with in Eq. (4). As regards the other boundaries, for some of the simulations we use periodic boundary conditions on the left and right boundaries that are perpendicular to the wall and for other simulations we make these hard purely repulsive walls. For the top boundary, we either assume that the system is closed, i.e. so that this is also a hard purely repulsive wall, or that the system is open, so that it is an absorbing boundary.
We also consider cases where there is a constant force on the particles parallel to the surface, for example due to gravity. In this case, a term is added to the expression for the external potential in Eq. (4).
Iii Dynamics of the liquid
Having defined the Hamiltonian, we must now specify the dynamics of the system. We use a combination of two types of particle moves to evolve the system forward in time: (i) Single particle moves, generated via the Metropolis Monte Carlo algorithm and (ii) collective manyparticle moves, generated by considering the dynamics from a thinfilm equation. We describe the moves of type (i) first.
iii.1 Monte Carlo dynamics
The Metropolis Monte Carlo algorithm Landau and Binder (2014) is commonly used to evolve systems in time, since it generates configurations with the correct (Boltzmann) equilibrium distribution:
(5) 
where the integer index labels the sequence of configurations generated by the algorithm, each having energy , given by the Hamiltonian (2). The normalisation factor is the partition function and , where is Boltzmann’s constant and is the temperature. The master equation governing the time evolution of the nonequilibrium probability is Landau and Binder (2014):
(6) 
where is the transition probability for the system to go from state to state . Eq. (6) has the form of a continuity equation, so the total probability remains normalised for all time. At equilibrium, where , Eq. (6) yields
(7) 
Since the equilibrium probability distribution function for state is given by Eq. (5), and similarly for state , the ratio of the probabilities also gives the transition rate ratio
(8) 
where . Below is outlined the Metropolis Monte Carlo algorithm that we implement, which samples with the above ratio of probabilities:

Pick a random particle.

Pick randomly a nearest neighbour lattice site to the selected particle.

Check whether the neighbouring site is empty. If it is, calculate the change in energy using Eq. (2) corresponding to moving the particle into the empty site.

If , then move the particle to its new position.

If , allow the move with probability .
The assumption we make (the standard assumption in KMC algorithms) is that Eq. (8) still holds out of equilibrium and so we can use the above algorithm to evolve the system forward in time. The KMC algorithm outlined above leads to the particles performing a diffusive dynamics that at equilibrium generates configurations with the (correct) probability distribution (5).
iii.2 Hydrodynamics via the thinfilm equation
While the above KMC algorithm correctly samples the equilibrium states with a diffusive dynamics, this is not in fact what nonequilibrium liquids on surfaces actually do. We therefore introduce additional particle moves to better model the hydrodynamics of the liquid on the surface. The additional particle moves are collective moves, transferring several particles simultaneously. These are governed by considering a thinfilm equation from hydrodynamic theory Thiele (2007). The effects of collective advective droplet hydrodynamic motion is illustrated in Fig. 1. The film height profile for the nonequilibrium droplet at time is given by the red dashed line, while the black solid line is the new height profile that the droplet assumes at the later time . In 3D, the thinfilm equation is:
(9) 
where is the film thickness profile, is the liquidvapour interfacial tension, is the viscosity and . The term represents the Laplace or curvature pressure. Here, where we assume the system is 2D, we simply replace . represents an additional force that may be present (e.g. gravity) acting in the positive direction. The term is the disjoining pressure, which describes the effective molecular forces between the substrate and the liquidvapour interface. The disjoining pressure can be written as the derivative of an excess free energy De Gennes et al. (2013), referred to as the binding potential, , as follows:
(10) 
The disjoining pressure is in fact a constrained excess free energy. The equilibrium film thickness is given by the minimum value of , where is the value of that minimises . When the minimum is at a small and finite value, then the liquid is said to be partially wetting. In contrast, when the minimum is at , then the liquid is said to be wetting. Additionally, the equilibrium contact angle is determined by the value of , as follows Churaev (1995); Rauscher and Dietrich (2008); Hughes et al. (2015):
(11) 
In the literature, there are many different approximations used for the disjoining pressure , depending on the nature of the fluid and the interactions with the surface De Gennes et al. (2013); Davis (1996); Israelachvili (2011); Yin et al. (2017). Here, we use the binding potential that we calculate from DFT using the approach of Ref. Hughes et al. (2015). The form of depends on the state point, i.e. on the particular values of , and the temperature . The DFT we use is also described below in Sec. IV. The DFT data is fitted to the form Hughes et al. (2015, 2017); Yin et al. (2017):
(12) 
where , , , , , , and are constants that depend on the state point. We also use the DFT to calculate the liquidvapour interfacial tension , as described below in Sec. IV. These are then used together with Eqs. (9) and (10) to determine the liquid droplet dynamics.
The relation between the KMC and thinfilm steps is as follows. After completing a prescribed number of KMC steps, , we first determine the discretised liquid film height profile from the configuration of occupied lattice sites. This is done by computing for each the number of occupied lattice sites above the lattice site up to a point where we find a vacancy of a certain number (we normally choose two) of consecutive lattice sites. This determines the height profile , where . If for a certain it turns out that , we set it equal to the precursor film thickness . This profile is then used as an initial condition for the thinfilm equation. We solve the thinfilm equation numerically using finitedifference approximations to represent spatial derivatives, and evolve the solution in time using e.g. Euler’s method. We choose a time step, , and solve the thinfilm equation for a certain number of time steps, . Next, we map the computed thinfilm profile back to the lattice by appropriately removing/adding particles where the film profile became thicker/thinner, as illustrated in Fig. 1. Note that this does not change the position of any particles in the vapour phase. To ensure the conservation of the number of the particles in the system, we then count the total number of the particles and add or remove a few particles, if needed, at random positions above the computed thinfilm profile. The resulting lattice site configuration can then be used again for the KMC computations. The evolution of the system is obtained by successively repeating the KMC and thinfilm steps as described above for a certain number of times, .
iii.3 Diffusion coefficient and time scales
An important quantity in the procedure described above is the parameter proportional to the ratio of the number of KMC attempted particle moves per particle to the number of thinfilm moves , or more strictly the total time the thin film equation is evolved for at each step, . In dimensionless time units, this ratio becomes:
(13) 
Here, is the total number of the particles in the system and is proportional to the viscosity of the liquid, :
(14) 
When , the system is evolved by using only the thinfilm equation, whilst the limit corresponds to evolving the system using only KMC particle moves. Therefore, , where is a singleparticle diffusion coefficient. The bulk liquid diffusion coeficient for a 3D model very similar to that considered here, with dynamics governed by the same KMC algorithm, was determined in Ref. Chalmers et al. (2017a). For temperatures in a similar range to those considered here, they obtained the result that KMC steps, where KMC steps means attempted moves per lattice site. In the present 2D system, one should expect to be of a similar magnitude in the bulk liquid. As regards the diffusion coefficient for isolated particles on the surface, in this regime .
Iv Lattice DFT for the fluid
We use classical DFT Evans (1979); Hansen and McDonald (2013) adapted to lattice models Hughes et al. (2014), in order to calculate the binding potential and also the liquidvapour interfacial tension . These quantities are then input into the thinfilm equation to create the collective hydrodynamic moves, as described in the previous section. We also use the DFT to calculate the bulk fluid phase diagram.
DFT is a statistical mechanical theory for the fluid onebody average density profile, which on a lattice is the quantity:
(15) 
where denotes an ensemble average. The approximation for the Helmholtz free energy that we use is
(16) 
A derivation of this is given in Ref. Hughes et al. (2014). There is an obvious (meanfield) connection between the Hamiltonian (2) and the terms on the second line in Eq. (16), and the first term is essentially entropic in origin Hughes et al. (2014). The equilibrium fluid density profile is obtained by minimising the grand potential
(17) 
where is the chemical potential Evans (1979); Hansen and McDonald (2013).
iv.1 Bulk fluid phase diagram from DFT
Determining the phase diagram of the bulk fluid, away from the influence of any interfaces, is a prerequisite to embarking on any study of the behaviour of the fluid at walls. Phase coexistence between a high density liquid phase and a low density vapour phase is exhibited when the temperature is less than the critical temperature . Thermodynamic coexistence of two phases occurs when the temperature, chemical potential and pressure are all equal in the two phases. The locus in the phase diagram of these coexisting states is referred to as the binodal and for our model is displayed in Fig. 2, plotted in the density versus temperature plane. We display the result from the DFT (16) (see the following paragraph for the details of how this is calculated) together with an estimate for the true location of the binodal. The latter is made by determining using our Monte Carlo simulations, based on Eq. (2), and then by scaling the exact Onsager 2D Ising result for the binodal Baxter (2012), so that the critical temperature matches that from our simulations, . Recall that the Onsager result is for the standard 2D Ising which just has nearest neighbour interactions, i.e. setting for . We see that in the region near , the meanfield DFT fails, but for temperatures the agreement between the DFT and the Monte Carlo results is fairly good, as we also see below where we present results for droplets on surfaces and for the wetting behaviour.
To calculate the binodal predicted by the free energy (16), a symmetry in the latticegas model proves to be rather useful. Note that this symmetry is absent in continuum models of fluids. If we replace in the Hamiltonian (2), then we see that this leaves the form of it unchanged. We can think of as being a ‘hole’ occupation number. As a consequence of this particlehole symmetry, we therefore have the following relation between the density of the liquid and the vapour at coexistence:
(18) 
From Eq. (16), we obtain the Helmholtz free energy per unit area in the absence of any external field:
(19) 
where is the area of the system. Recall that we are considering a 2D system; one should, of course, consider the free energy per unit volume in 3D, where is the volume. Note also that we have used the result , to obtain Eq. (19). This is the sum over the interactions of a particle with all of its neighbours [see Eq. (3)]. The factor is present to avoid double counting the interactions between each of the pairs. The pressure in the system is therefore Hansen and McDonald (2013)
(20)  
Invoking Eq. (18) and solving leads to the binodal curve
(21) 
This is the curve displayed in Fig. 2. The maximum on the binodal corresponds to the critical point temperature, above which there is no vapourliquid phase separation. Substituting the critical density into Eq. (21) gives the critical temperature predicted by the DFT to be , while from our KMC simulations we find that in reality it is . We can also calculate the chemical potential
(22)  
The value of the chemical potential at bulk phase coexistence (i.e. on the binodal) is found by substituting Eq. (21) into Eq. (22), which gives
(23) 
In Fig. 2 we also display the spinodal predicted by the DFT. Spontaneous phase separation occurs within this region of the phase diagram, since the uniform fluid is unstable for state points within the spinodal curve. The condition to determine the spinodal is
(24) 
Using this together with Eq. (19) we obtain the following result for the spinodal
(25) 
which is the curve plotted in Fig. 2.
iv.2 Properties of the vapourliquid interface
Having determined the bulk fluid phase diagram, one can then determine the properties of the interface between the coexisting liquid and vapour phases. In Fig. 3 we display a sequence of density profiles for the vapourliquid interface, calculated over a range of different temperatures. The density profiles are calculated by minimising the grand potential (17) with the Helmholtz free energy given by Eq. (16) and the chemical potential set to be the value at coexistence, given in Eq. (23). To calculate these we set the external potential and assume that the profiles only vary in one direction (along the axis). We also assume the system has periodic boundary conditions and initiate the system with half of it containing the liquid and the other half containing the vapour. The final equilibrium therefore contains two vapourliquid interfaces, but only one of these is displayed in Fig. 3.
Having calculated these profiles, we can then calculate the liquidvapour interfacial tension . Recall that a planar interface between the vapour and the liquid in a 3D system adds to the free energy of the system an excess (over bulk) contribution equal to the interfacial tension multiplied by the area of the interface Evans (1979); Hansen and McDonald (2013). Similarly, for the present 2D system the interfacial contribution is equal to , where is the length of the interface. The factor 2 arises because there are two interfaces in our system with periodic boundary conditions. Having calculate the grand potential for the system containing the interface, the interfacial tension is then
(26) 
where
(27) 
is the grand potential for a 2D box having the same area , but only containing one of the two phases, and where is pressure. For example, when , this approach gives the vapourliquid interfacial tension to be .
In Fig. 4 we display the surface tension of the planar liquidvapour interface plotted as a function of temperature. We see that as the temperature approaches from below the critical temperature , then , due to the fact that the difference in the density of the two coexisting phases goes to zero as . In addition to the full DFT results for , in Fig. 4 we also display the results from a socalled ‘sharpkink’ approximation (c.f. Dietrich (1988); Stewart and Evans (2005)), where we assume that the interface consists of a sharp step from the liquid with density , to the vapour with density , both being the values at phase coexistence; see Eqs. (18) and (21). The excess free energy for this density profile is then evaluated via Eq. (16). We also show in Fig. 4 the results from an even cruder sharpkink type approximation that assumes that in the liquid and then use the Hamiltonian to evaluate the energy change when such a system is cleaved in two to create an interface with a vacuum, i.e. we effectively assume for this calculation that in the vapour , with a steplike planar interface between the two. We see from Fig. 4 that there is good agreement between the value for the surface tension calculated from DFT and from the first sharpkink approximation for all temperatures and also with the second sharpkink approximation for temperatures , indicating that in this regime the major contribution to the interfacial tension is energetic, rather than entropic.
The calculations of the solidvapour and solidliquid interfacial tensions, and , respectively, are done in a similar manner to the calculation of described above, but of course the wall potential is retained and the system is initialised with just either the vapour or the liquid density values, respectively. For additional details on how these calculations are done, see e.g. the discussion in Ref. Hughes et al. (2014).
iv.3 Calculation of the binding potential using DFT
The excess grand potential per unit length (per unit area, in 3D) for a film of liquid adsorbed between a wall and a bulk vapour can be written as Dietrich (1988); Hughes et al. (2014, 2015, 2017)
(28) 
where is the length of the interface with the wall, is given by Eq. (27), , is the binding potential and is the adsorption, given by
(29) 
which is the excess (over bulk) number of particles adsorbed on the interface. Note that strictly we should consider to be a function of the adsorption , not of , since the latter is not a well defined quantity Hughes et al. (2015). However, when is large Eq. (29) gives
(30) 
Therefore, throughout this paper we discuss and almost interchangeably, and treat Eq. (30) as an equality which defines the value of .
When the system is at vapourliquid phase coexistence we have and so the last term in Eq. (28) is zero. We also see that when we must have , since when the only contributions to must be those from the wallliquid and the liquidvapour interfaces, i.e. and , respectively. These considerations also enables us to see what really is: it is the additional contribution to from the interaction between these two interfaces when they become close to one another. For the case when the liquid is partially wetting, i.e. when has a minimum at a finite value of , then ‘binds’ these two interfaces together, hence the name given to . This interaction between the interfaces stems from the molecular interactions between fluid particles and with the wall and can be thought of as arising in a similar manner to the surface tension.
The equilibrium film thickness is that minimises Eq. (28). When , this is obtained by solving . For all other values of , the system is out of equilibrium and therefore to determine for other values of , one must either move away from coexistence or apply a constraint, the constraint being that there is a specified thickness of liquid at the wall.
The approach used to calculate here is that developed in Refs. Archer and Evans (2011); Hughes et al. (2015), which consists of a constrained free energy minimisation. This method uses an iterative algorithm to minimise the grand potential (17) subject to the constraint that the adsorption is a specified value, which via Eq. (30) is equivalent to minimising subject to the constrain that the film thickness is a specified value. It turns out that this constraint is equivalent to applying a fictitious additional external potential that stabilises the liquid film of the specified thickness Archer and Evans (2011); Hughes et al. (2015). The calculation is repeated for a sequence of different values of , enabling us to calculate over the whole range of values of . This data set is subsequently fitted to the form given in Eq. (12) and then input into the thinfilm equation (9).
In Figs. 5 and 6 we display results for the binding potential determined via this method for a range of different state points. In Fig. 5 we display results for , which corresponds to a temperature that is not that far below the critical temperature . In Fig. 5(a) are results for a sequence of different values of the wall attraction strength . We see that for small values of , which corresponds to a weakly attracting solvophobic wall, the global minimum of is at a small value of , i.e. the liquid does not wet the wall. As is increased, there is a first order wetting transition when and for , the global minimum in is at ; i.e. for these values of the liquid wets the wall. The inset of Fig. 5(a) shows the same results with the vertical axis plotted on a logarithmic scale to show more clearly the monotonic decay behaviour as . In Fig. 5(b) we show the result for to demonstrate the excellent agreement between the numerical results for and the fit function in Eq. (12). This shows clearly that this simple form is highly appropriate for fitting the binding potential over the whole range of values of Hughes et al. (2015, 2017); Tewes et al. (2017).
In Fig. 6 we display results for the binding potential for , which corresponds to a much lower temperature than the results in Fig. 5. The results in Fig. 6(a) are for a range of different values of and we see that for the minimum in close to the wall dissapears, so that for the global minimum in is at a larger value of . However, we note that in this case, due to the presence of oscillations in that do not decay in amplitude, there is an infinite sequence of minima in as . Oscillations that decay in amplitude can occur in the binding potential as a result of the nature of the molecular interactions Hughes et al. (2017); Yin et al. (2017). However, in the present system the oscillations are due the low temperature and the discrete (lattice) nature of the model Hughes et al. (2015). Each subsequent minimum in corresponds to the vapourliquid interface moving by a distance of one lattice site further from the wall. The inset of Fig. 6(a) shows the same results with the vertical axis plotted on a logarithmic scale to show more clearly the oscillatory decay behaviour. We see in Fig. 6(b) that, as for the previous case, the fit function (12) gives an excellent representation of the numerical DFT results, in this case for .
V Contact angle and static droplet profiles from KMC
As described in the previous section, it is straightforward to calculate the three interfacial tensions , and using our lattice DFT (16). These can then be inserted into Young’s equation (1) to determine the equilibrium contact angle . This is also the contact angle that equilibrium droplet solutions of the thinfilm equation have, since is determined by , as Eq. (11) shows. It is important that the value of determined by in the thinfilm equation is close to the value of the contact angle that the KMC simulations drive the system towards. If there is too big a difference between the contact angle from DFT and from the KMC simulations, then when we have both the KMC and thinfilm dynamics in our model, these will work against each other, which would be unsatisfactory. Thus, the contact angle from the DFT and from the KMC must be close in value. This is also a check of the accuracy of the DFT. Of course, at higher temperatures near the DFT is not accurate; that is clear from the phase diagram in Fig. 2. However, at lower temperatures we do find that the DFT becomes rather accurate, as we now show.
In Fig. 9 we display a sequence of density profiles, each obtained by averaging over a time series from our KMC simulations, for various different values of the wall attraction parameter and for the temperature . For these simulations we do not include any of the thinfilm collective dynamics, since at this stage we are solely interested in equilibrium properties. In each case there are particles in the system, which has a total domain size of lattice sites and with periodic boundary conditions between the left and right sides. We initiate the system by setting for all the lattice sites within a semicircular region in the middle upon the substrate and we set outside the semicircle. We see that for the largest value of the wall attraction , the liquid spreads out over the substrate; i.e. it wets the substrate. For smaller values of , the liquid forms a droplet on the substrate with an increasingly large contact angle as is decreased; i.e. for these cases the liquid is partially wetting. This behaviour and the observed wetting transition at is in very good agreement with what we expect based on Fig. 6. Note too that the droplet profiles in Fig. 9 are similar to the droplet profiles calculated in Refs. Hughes et al. (2014, 2015) using DFT for a very similar latticegas model; see also the 3D droplet density profile calculated using Monte Carlo simulations in Ref. Chalmers et al. (2017a).
From droplet density profiles such as those in Fig. 9, following e.g. Ref. Chalmers et al. (2017a), we can calculate the contact angle by fitting a circle to the location of the liquidvapour interface, which is defined as the points determined by linear interpolation between pairs of neighbouring lattice sites, where the density . Results for the contact angle determined via this method are plotted in Fig. 10 for the temperature and for a range of different values of the wall attraction . We also display the results from DFT and via Young’s equation (1). We see excellent agreement between the results from the two approaches, demonstrating that at such temperatures (away from the critical point) the DFT is accurate and so the binding potential obtained from the DFT should also be reliable.
The results in Fig. 10 show that the contact angle for , where the liquid wets the wall. As is decreased, decreasing the strength of the attraction between the particles and the surface so that the surface becomes more solvophobic, the contact angle increases and the liquid increasingly beadsup upon the surface.
Vi Results for nonequilibrium droplets on surfaces
The results in the preceding sections essentially constitute a series of checks that the model has the desired, correct and selfconsistent thermodynamics. In this section, we present results for the nonequilibrium behaviour predicted by the model.
vi.1 Pure KMC at two different temperatures
For droplets to be able to equilibrate with the vapour phase, the system must be enclosed (i.e. treated in the canonical ensemble), with a fixed number of particles in the system. This is achieved by assuming a hard but nonattracting () wall on the top.
We first present results for the behaviour of two droplets using only KMC simulations. The total domain is of size , with periodic boundary conditions between the left and right sides. The system is initiated by setting for all the lattice sites within two semicircular regions of radii with the centres of the corresponding circles at the sites and , and we set outside of these two semicircles. This gives particles in the system.
The results in Fig. 11 are for the temperature and we choose the wall attraction parameter . The contact angle for these parameter values is . The different panels in Fig. 11 correspond to a sequence in time, i.e. for increasing numbers of the KMC steps, as indicated in the top right of each. The droplets initially evaporate and somewhat decrease in size, equilibrating with a low density vapour phase. They also adjust their average contact angle from the initial value of closer to , the equilibrium value (c.f. Fig. 9). However, as time proceeds, due to the random (thermal) fluctuations, the size of one of the droplets (the left one in this KMC simulation) keeps decreasing, but the size of the other droplet (the right one in this KMC simulation) starts to increase. The centres of the droplets remain at approximately the initial positions. Eventually, the left droplet completely disappears, leaving the system containing only a single droplet, which of course is the global free energy minimum state. In the case in Fig. 11 (i.e. for the temperature ), this coarsening process occurs mainly via the diffusive migration of the particles through the vapour, evaporating from one droplet and condensing onto the other, although there is additionally a small particle flux over the surface. In Ref. Pototsky et al. (2014), when clusters (equivalent to our droplets) aggregate in a similar manner, it is referred to as joining via the ‘Ostwald mode’, in contrast to joining by the droplets moving towards one another, which is referred to as joining via the ‘translation mode’. For the model system considered in Ref. Pototsky et al. (2014), it is possible to show that there is an eigenfunction associated with each of these two possible modes and that it is the mode with the largest eigenvalue that dominates the observed dynamics, although in principle it is possible for the two eigenvalues to be similar in value. However, in the case in Fig. 11 it is clearly the Ostwald mode that is dominant. For more background on Ostwald ripening, see e.g. Ref. Ratke and Voorhees (2013).
Figure 12 correspond to a lower temperature and we choose the wall attraction parameter so that the equilibrium contact angle is still . The droplets in this case initially spread a little over the surface because the KMC time evolution drives the system towards achieving the equilibrium contact angle for both droplets. There is a much smaller rate of evaporation because of the lower temperature, or, equivalently, due to the stronger attraction between the particles in this case. Then over time the droplets eventually merge and form a single droplet in the centre of the domain. The coarsening process in this case appears to mainly be due to overall translation of the droplets over the surface, i.e. coarsening due to the translational mode Pototsky et al. (2014); Ratke and Voorhees (2013).
vi.2 Combination of KMC and thinfilm equation dynamics
We now consider behaviour of the system when the fluid dynamics governed by a combination of the singleparticle KMC moves in conjunction with the collective particle moves generated by the thinfilm equation. We use the same domain and the same initial configuration as in the previous subsection. Figure 13 corresponds to the same temperature and wall attraction strength as in Fig. 11. For this simulation, we used cycles of KMC steps followed by thinfilm steps with the time step , which gives , where is the parameter defined by Eq. (13). As in Fig. 11, the system coarsens into one droplet. However, the speed of coarsening is increased when the thinfilm steps are included, which is not too surprising, since the total number of possible moves has been increased by the addition of the thinfilm moves. As before, in the initial stages the droplets partially evaporate and spread, to reach a contact angle close to the equilibrium value. In the later stages, it can be observed that the size of one of the droplets (the right one in this simulation) keeps decreasing, but the size of the other droplet (the left one in this simulation) starts to increase. As in the case in Fig. 11, the centres of the droplets in Fig. 13 remain more or less fixed. The coarsening in this case is mainly due to a flux of particles evaporating from one droplet, travelling though the vapour and condensing onto the other, although there is also a smaller flux due to a diffusive migration of the particles over the surface. To illustrate this in more detail, we show in the bottom right panel of Fig. 13 how the numbers of particles that cross the vertical line in the middle of the domain at varies in time. These are split into two populations: (i) those that migrate over the surface at height (, see the blue line) and (ii) those that travel through the vapour at height (, see the red line). Increasing/decreasing or means that the particle moves to the left/right, respectively, crossing the line. In this particular simulation, it can be observed that there is relatively low migration of the particles over the surface, whilst there is a significant net transfer of particles from left to right through the vapour. Thus, we can conclude that in this case the main coarsening mechanism is due to particles evaporating from one droplet, diffusing through the vapour, and then condensing onto the other droplet. Using only the thinfilm equation would not allow us to capture this, and this justifies the importance of our hybrid model which combines both the thinfilm equation governed dynamics with the KMC diffusive dynamics.
Figure 14 corresponds to the same values of the temperature and the wall attraction strength parameter as in Fig. 12. For this simulation, we used cycles of KMC steps followed by thinfilm steps with the time step , which gives the same value of as for the results in Fig. 13. As in the higher temperature case in Fig. 13, we find that the speed of the coarsening is increased when the thinfilm steps are included. As in the pure KMC case in Fig. 12, the droplets initially spread over the surface to reach a contact angle close to the equilibrium value , but with only a small amount of evaporation due to the low temperature. As time proceeds, the size of the left droplet starts to decrease and the right droplet starts to grow, and apparently this happens mainly due to motion of particles from one to the other over the surface. To understand the coarsening mechanisms in more detail, we present in the last panel of Fig. 14 the time evolution of the numbers of the particles that cross the vertical line in the middle of the domain, at , that migrate over the surface, , and through the vapour . In this particular simulation, it can be observed that there is a relatively low rate of migration of the particles through the vapour, whilst the majority of the mass transfer is via the surface. Thus, we can conclude that in this case of lower temperature, the main coarsening mechanism is due to the transfer of the particles over the surface.
vi.3 Droplet evolving under lateral driving
In this subsection, we consider the dynamics predicted by the model when the fluid is under the influence of a constant lateral driving force parallel to the surface, i.e. we have a constant force to the right on each particle. The domain is a box of size , with hard walls on both the left and right boundaries, as well as the attracting surface on the bottom and the hard wall at the top considered previously. The system is initiated by setting for all the lattice sites within a semicircular region of radius with the centre at the site , and we set outside of this semicircle. This gives particles in the system. As a consequence of the lateral driving and the walls at the left and right boundaries, in all of the results displayed in this section in the final stages of the dynamics we see a pileup of particles on the righthand boundary wall.
Figure 15 shows results for when the dynamics is solely generated by the KMC algorithm, with no thinfilm dynamics. The initial configuration consists of the semicircular droplet displayed in the top left panel, which has contact angle . We set and the wall attraction strength . We see from the results in Fig. 15 that the droplet centre of mass hardly moves over time. However, since there is a constant (but slow) evaporation rate of particles from the surface, once these leave the droplet they are then carried away to the right by the driving force . Over time, these then build up on the wall on the right side. Thus, although on every particle in the system there is a force to the right , due to the KMC dynamics only considering one particle at a time, we see that there is little or no collective motion resulting from this force. As we see next, this is not the case when we also allow collective motion by considering the dynamics from the thin film equation.
Figure 16 shows snapshots from a simulation where all the model parameters are the same as for the simulation in Fig. 15 except we now include some collective multiparticle moves generated by considering the thinfilm equation (9) with . We set , and , which gives . We see that the thinfilm moves enable the droplet to migrate across the surface under the lateral driving. By comparing with the previous result in Fig. 15, it is clear that this dynamics must be a collective multiparticle (hydrodynamic) phenomenon. Though the majority of the particle motion to the right is now via the droplet sliding over the surface, there is still a small flux through the vapour. Note too that the overall timescale for the dynamics is reduced, due to the additional thinfilm moves.
In Fig. 17 we show results for a system with an even higher proportion of particle moves generated by the thinfilm equation. In this case we have , and , which gives . We again observe that the thinfilm moves enable the droplet to slide from left to right under the driving. Additionally, for this choice of parameter values we see that the droplet takes a characteristic shape with a ‘tail’ trailing behind, typical sometimes of randroplets sliding down a window. The appearance and bifurcations of such tails is discussed in Ref. Wilczek et al. (2017) for the thinfilm equation.
Vii Concluding remarks
We have developed a hybrid model for droplet dynamics on surfaces, that combines both singleparticle diffusive dynamics via KMC and advective hydrodynamic multiparticle moves via a thinfilm equation. We have shown that the model incorporates the correct thermodynamics and hydrodynamics, microscopic structure, evaporation, condensation and the diffusion of particles over a surface. The model contains parameters which can be determined by relating to experiments. For example, from measurements of the equilibrium contact angle of drops of the liquid on the relevant surfaces and the surface tension, the interaction parameters can be estimated and the mobility coefficients estimated from the bulk liquid viscosity and diffusion coefficient – see also the discussion in Ref. Chalmers et al. (2017a) about determining these coefficients.
An important part of the coarsegraining that we performed in order to bridge from the microscopic Hamiltonian (2) to the mesoscopic thinfilm equation (9), is the use of DFT. Here we considered a lattice model, to make this bridging more easily. However, we see no reason why this could not also be done for a continuum model. The main requirement is to have an accurate DFT, that is able to describe the structure of the liquid of interest. For example, for the LenardJones model liquid, much of the work has already been done in Ref. Hughes et al. (2017). The one thing that would need to be adapted and developed further is the algorithm for determining the liquid droplet height profile and the algorithm for moving the particles as determined by the thinfilm equation. For the former, we believe that calculating a local adsorption would be sufficient. However, the present algorithm for inserting the particles would surely need some testing and refining, since in the present work the underlying lattice makes this process much more straightforward to implement. We should also emphasise how valuable it is to make the connection with DFT, since this gives access to thermodynamic quantities such as the surface tensions, equilibrium contact angles, pressure and the bulk fluid phase behaviour. These are important things to know and are much harder to obtain via other methods; for example, to obtain these from KMC simulations is more complicated and requires more lengthy computations.
For the purposes of the present study the DFT results presented here in Sec. IV were performed in order to make sure the thermodynamic underpinnings of the thinfilm equation used in the time evolution is correct. However, we believe the present work is also a valuable contribution to the development and testing of latticeDFT models of this type. There has previously been some testing of the meanfield DFT (16) by comparing to simulations – see e.g. Ref. Chalmers et al. (2017b) – but the results presented here showing that the DFT is also accurate at some temperatures for describing the contact angle and the other related interfacial properties have not been performed before for a 2D model, to the best of our knowledge.
We have presented striking results showing the process of how a pair of droplets on a surface join together depends on the temperature and also the other system parameters: at higher temperatures the coarsening process largely occurs via evaporation of particles from the smaller droplets, diffusion through the vapour and then condensation onto the larger droplets, i.e. via the Ostwald mode Pototsky et al. (2014). In contrast, at lower temperatures drops coalesce via the translation mode, moving towards one another over the surface and joining. In the present lattice system it is not possible to calculate the eigenfunctions and eigenvalues associated with each of these modes, as it is for the simpler continuum model in Ref. Pototsky et al. (2014). However, it is clear that these processes are present here. In fact, there are tantalising hints in some of our results that there may in fact be two different surface modes: an Ostwaldtype mode, determined by the diffusion of particles over the surface (rather than via the vapour), which is somewhat different in character to the collective surface motion due to capillary forces that is described by the thinfilm equation. However, to confirm this conjecture and really distinguish that there are two distinct surface modes would require a much simpler continuum model that is amenable to the analysis approach used in Ref. Pototsky et al. (2014).
The results presented in Sec. VI.3, where the system evolves under lateral driving, show that collective manyparticle hydrodynamic motion is essential for droplets to slide over surfaces. The thinfilm equation dynamics is also crucial to whether the droplet takes a sliding shape with a tail or not. If the particles are assumed to move one at a time, then the droplets do not slide. The knowledge that the overall droplet shape is linked to the underlying particle motion gives valuable insight into the behaviour of these types of systems.
Viii Acknowledgements
The authors would like to thank Tyler Shendruk and Uwe Thiele for valuable comments. MA gratefully acknowledges the KSA Ministry of Higher Education and University of Tabuk for financial support.
References
 AlaNissila et al. (2002) T. AlaNissila, R. Ferrando, and S. C. Ying, Adv. Phys. 51, 949 (2002).
 De Gennes et al. (2013) P.G. De Gennes, F. BrochardWyart, and D. Quéré, Capillarity and wetting phenomena: drops, bubbles, pearls, waves (Springer, 2013).
 Pandit and Wortis (1982) R. Pandit and M. Wortis, Phys. Rev. B 25, 3226 (1982).
 Pandit et al. (1982) R. Pandit, M. Schick, and M. Wortis, Phys. Rev. B 26, 5112 (1982).
 Binder and Landau (1988) K. Binder and D. P. Landau, Phys. Rev. B 37, 1745 (1988).
 Rabani et al. (2003) E. Rabani, D. R. Reichman, P. L. Geissler, and L. E. Brus, Nature 426, 271 (2003).
 Sztrum et al. (2005) C. G. Sztrum, O. Hod, and E. Rabani, J. Phys. Chem. B 109, 6741 (2005).
 PauliacVaujour et al. (2008) E. PauliacVaujour, A. Stannard, C. Martin, M. O. Blunt, I. Notingher, P. Moriarty, I. Vancea, and U. Thiele, Phys. Rev. Lett. 100, 176102 (2008).
 Vancea et al. (2008) I. Vancea, U. Thiele, E. PauliacVaujour, A. Stannard, C. P. Martin, M. O. Blunt, and P. J. Moriarty, Phys. Rev. E 78, 041601 (2008).
 Blunt et al. (2010) M. O. Blunt, A. Stannard, E. PauliacVaujour, C. P. Martin, I. Vancea, M. Suvakov, U. Thiele, B. Tadic, and P. Moriarty, in Oxford Handbook of Nanoscience and Technology (2010).
 Stannard (2011) A. Stannard, J. Phys.: Condens. Matter 23, 083001 (2011).
 Vancea (2011) I. Vancea, Ph.D. thesis, Loughborough University (2011).
 Kim et al. (2011) H.S. Kim, S. S. Park, and F. Hagelberg, J. Nanoparticle Research 13, 59 (2011).
 Jung et al. (2014) N. Jung, C. S. Yoo, and P. H. Leo, J. Phys. Chem. B 118, 2535 (2014).
 Crivoi and Duan (2014) A. Crivoi and F. Duan, Scientific Reports 4 (2014).
 Tewes et al. (2017) W. Tewes, O. Buller, A. Heuer, U. Thiele, and S. V. Gurevich, J. Chem. Phys. 146, 094704 (2017).
 Chalmers et al. (2017a) C. Chalmers, R. Smith, and A. J. Archer, J. Phys.: Condens. Matter 29, 295102 (2017a).
 Fichthorn and Weinberg (1991) K. A. Fichthorn and W. H. Weinberg, J. Chem. Phys. 95, 1090 (1991).
 Oron et al. (1997) A. Oron, S. H. Davis, and S. G. Bankoff, Rev. Mod. Phys. 69, 931 (1997).
 Kalliadasis and Thiele (2007) S. Kalliadasis and U. Thiele, Thin Films of Soft Matter (Springer, 2007).
 Thiele (2007) U. Thiele, in Thin films of soft matter (Springer, 2007), pp. 25–93.
 Craster and Matar (2009) R. V. Craster and O. K. Matar, Rev. Mod. Phys. 81, 1131 (2009).
 Archer and Evans (2011) A. J. Archer and R. Evans, Mol. Phys. 109, 2711 (2011).
 Hughes et al. (2015) A. P. Hughes, U. Thiele, and A. J. Archer, J Chem. Phys. 142, 074702 (2015).
 Hughes et al. (2017) A. P. Hughes, U. Thiele, and A. J. Archer, J. Chem. Phys. 146, 064705 (2017).
 Buller et al. (2017) O. Buller, W. Tewes, A. J. Archer, A. Heuer, U. Thiele, and S. V. Gurevich, J. Chem. Phys. 147, 024701 (2017).
 Evans (1979) R. Evans, Adv. Phys. 28, 143 (1979).
 Evans (1992) R. Evans, Fundamentals of Inhomogeneous Fluids (Dekker, New York, 1992).
 Hansen and McDonald (2013) J.P. Hansen and I. R. McDonald, Theory of simple liquids: with applications to soft matter (Academic Press, 2013).
 Yin et al. (2017) H. Yin, D. N. Sibley, U. Thiele, and A. J. Archer, Phys. Rev. E 95, 023104 (2017).
 Tretyakov et al. (2013) N. Tretyakov, M. Müller, D. Todorova, and U. Thiele, J. Chem. Phys. 138, 064905 (2013).
 MacDowell et al. (2002) L. G. MacDowell, M. Müller, and K. Binder, Colloids Surf. A 206, 277 (2002).
 Thiele et al. (2009) U. Thiele, I. Vancea, A. J. Archer, M. J. Robbins, L. Frastia, A. Stannard, E. PauliacVaujour, C. P. Martin, M. O. Blunt, and P. J. Moriarty, J. Phys.: Condens. Matter 21, 264016 (2009).
 Robbins et al. (2011) M. J. Robbins, A. J. Archer, and U. Thiele, J. Phys.: Condens. Matter 23, 415102 (2011).
 Chalmers et al. (2017b) C. Chalmers, R. Smith, and A. J. Archer, Langmuir 33, 14490 (2017b).
 Grün et al. (2006) G. Grün, K. Mecke, and M. Rauscher, J. Stat. Phys. 122, 1261 (2006).
 DuránOlivencia et al. (2019) M. A. DuránOlivencia, R. S. Gvalani, S. Kalliadasis, and G. A. Pavliotis, J. Stat. Phys. 174, 579 (2019).
 Krüger et al. (2017) T. Krüger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva, and E. M. Viggen, The lattice Boltzmann method (Springer, 2017).
 Archer et al. (2010) A. J. Archer, M. J. Robbins, and U. Thiele, Phys. Rev. E 81, 021602 (2010).
 Landau and Binder (2014) D. P. Landau and K. Binder, A guide to Monte Carlo simulations in statistical physics (Cambridge university press, 2014).
 Churaev (1995) N. V. Churaev, Adv. Colloid Interface Sci. 58, 87 (1995).
 Rauscher and Dietrich (2008) M. Rauscher and S. Dietrich, Annu. Rev. Mater. Res. 38, 143 (2008).
 Davis (1996) H. T. Davis, Statistical Mechanics of Phases, Interfaces and Thin Films (WileyVCH, 1996).
 Israelachvili (2011) J. N. Israelachvili, Intermolecular and surface forces (Academic press, 2011).
 Hughes et al. (2014) A. P. Hughes, U. Thiele, and A. J. Archer, Am. J. Phys. 82, 1119 (2014).
 Baxter (2012) R. J. Baxter, J. Stat. Phys. 149, 1164 (2012).
 Dietrich (1988) S. Dietrich, in Phase Transitions and Critical Phenemona, edited by C. Domb and J. Lebowitz (1988), vol. 12.
 Stewart and Evans (2005) M. C. Stewart and R. Evans, Phys. Rev. E 71, 011602 (2005).
 Pototsky et al. (2014) A. Pototsky, U. Thiele, and A. J. Archer, Phys. Rev. E 89, 032144 (2014).
 Ratke and Voorhees (2013) L. Ratke and P. W. Voorhees, Growth and coarsening: Ostwald ripening in material processing (Springer Science & Business Media, 2013).
 Wilczek et al. (2017) M. Wilczek, W. Tewes, S. Engelnkemper, S. V. Gurevich, and U. Thiele, Phys. Rev. Lett. 119, 204501 (2017).