Kinetic Monte Carlo and hydrodynamic modelling of droplet dynamics on surfaces, including evaporation and condensation

TL;DR

This paper introduces a combined kinetic Monte Carlo and hydrodynamic model to simulate droplet dynamics on surfaces, capturing evaporation, condensation, and fluid motion at microscopic and macroscopic scales.

Contribution

It develops a novel lattice-gas model integrating microscopic particle interactions with macroscopic fluid dynamics for droplet behavior on surfaces.

Findings

Droplets can join, spread, and slide under gravity.

Evaporation and condensation significantly affect droplet dynamics.

The model links microscopic parameters to macroscopic quantities like contact angle.

Abstract

We present a lattice-gas (generalised Ising) model for liquid droplets on solid surfaces. The time evolution in the model involves two processes: (i) Single-particle 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) Larger-scale collective moves, modelling advective hydrodynamic fluid motion, determined by considering the dynamics predicted by a thin-film 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.

Equations

$${\gamma }_{lv}\cos \theta ={\gamma }_{sv}-{\gamma }_{sl}.$$

$$E=-\epsilon \sum _{i,j}{c}_{i,j}{l}_i{l}_j+\sum _i{V}_i{l}_i,$$

$$c_{i,j}=\{\begin{array}{ll}1& \text{if }j\in NNi,\\ \frac{1}{2}& \text{if }j\in NNNi,\\ 0& \text{otherwise,}\end{array}$$

$$V_i=\{\begin{array}{ll}\infty & \text{if }i=(i,0),\\ -{\epsilon }_w& \text{if }i=(i,1),\\ 0& \text{if }i=\text{otherwise},\end{array}$$

$$P_n=\frac{1}{Z}{e}^{-\beta E_n},$$

$$\frac{\partial P_n(t)}{\partial t}=-\sum _{m\ne n}[P_n(t)W_{n\to m}-P_m(t)W_{m\to n}],$$

$$P_n(t)W_{n\to m}=P_m(t)W_{m\to n}.$$

$$\frac{P_n(t)}{P_m(t)}=\frac{W_{m\to n}}{W_{n\to m}}=e^{-\beta \Delta E},$$

$$\frac{\partial h}{\partial t}=-\nabla \cdot \left[\frac{h^3}{3\eta }\left(\nabla [{\gamma }_{lv}{\nabla }^{2}h+\Pi (h)]+\mathcal{G}\right)\right],$$

$$\Pi (h)=-\frac{\partial g}{\partial h}.$$

$$\theta ={\cos }^{-1}\left(1+\frac{g(h_{min})}{{\gamma }_{lv}}\right).$$

$$g(h)=a_1\sin (2\pi a_{2}h+a_3)+a_4\exp (-h/b)+a_5\exp (-2h/b)+a_6\exp (-3h/b),$$

$$\lambda =\frac{M_1\tau }{M_{2}N\Delta t}.$$

$$\tau =3\eta \beta {\sigma }^3.$$

$${\rho }_i=⟨l_i⟩,$$

$$F(\{{\rho }_i\})=k_{B}T\sum _i[{\rho }_i\ln {\rho }_i+(1-{\rho }_i)\ln (1-{\rho }_i)]-\epsilon \sum _{i,j}{c}_{i,j}{\rho }_i{\rho }_j+\sum _i{V}_i{\rho }_i.$$

$$\Omega (\{{\rho }_i\})=F(\{{\rho }_i\})-\mu \sum _i{\rho }_i,$$

$${\rho }_l=1-{\rho }_v.$$

$$f=\frac{F}{A}=k_{B}T[\rho \ln \rho +(1-\rho )\ln (1-\rho )]-3\epsilon {\rho }^2,$$

$$p$$

$$=-k_{B}T\ln (1-\rho )-3\epsilon {\rho }^2.$$

$$\frac{k_{B}T}{\epsilon }=\frac{3(2\rho -1)}{\ln (\rho /(1-\rho ))}.$$

$$\mu$$

$$\displaystyle=k_{B}T\ln\bigg{(}\dfrac{\rho}{1-\rho}\bigg{)}-6\varepsilon\rho.$$

$${\mu }_{coex}=-3\epsilon .$$

$$\frac{{\partial }^{2}f}{\partial {\rho }^2}=0.$$

$$\frac{k_{B}T}{\epsilon }=6\rho (1-\rho ),$$

$${\gamma }_{lv}=\frac{\Omega -{\Omega }_0}{2L},$$

$${\Omega }_0=-pA,$$

$${\omega }_{ex}(h)\equiv \frac{\Omega -{\Omega }_0}{L}={\gamma }_{sl}+{\gamma }_{lv}+g(h)+\Gamma \delta \mu ,$$

$$\Gamma =\frac{1}{L}\sum _i({\rho }_i-{\rho }_v),$$

Full text

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 Ala-Nissila 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 single-particle 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 lattice-gas (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 non-equilibrium droplets and liquid films on surfaces have been developed Rabani et al. (2003); Sztrum et al. (2005); Pauliac-Vaujour 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 non-equilibrium 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 large-scale collective particle moves by considering the dynamics that is predicted by a thin-film equation. This is a time evolution equation for the liquid film thickness profile $h$, i.e. the height of the liquid-vapour 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 Navier-Stokes equations for a film of liquid on a surface by making a long-wave approximation, which greatly simplifies the analysis. When generating these ‘thin-film moves’, we first determine the liquid film height profile from the configuration of occupied lattice sites, we then evolve the thin-film equation for a short amount of time, and then map back to the lattice. In our model, an important quantity is the parameter $\lambda$, defined below in Eq. (13), which is proportional to the ratio of the number of KMC attempted particle moves to the number of thin-film moves. When $\lambda=0$, our model reduces to just evolving the thin-film equation, whilst in the limit $\lambda\to\infty$, our model is a pure KMC model. Therefore, this parameter $\lambda\propto D/\eta$, where $D$ is a single-particle diffusion coefficient and $\eta$ is the viscosity.

Key quantities that characterises how a liquid wets a surface are the spreading parameter $S=\gamma_{sv}-(\gamma_{sl}+\gamma_{lv})$, where $\gamma_{sv}$, $\gamma_{sl}$, and $\gamma_{lv}$ are the solid-vapour, solid-liquid and liquid-vapour interfacial tensions, respectively De Gennes et al. (2013), and also the contact angle $\theta$, which is given by Young’s equation De Gennes et al. (2013):

[TABLE]

The contact angle is the angle that the height profile for equilibrium droplets makes with the surface. $\gamma_{lv}$ is a parameter that must be input into the thin-film equation, whilst $\theta$ is determined by the disjoining pressure $\Pi(h)$, which is a function that must also be input into the thin-film 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 $g(h)$, from which one obtains $\Pi(h)=-\partial g(h)/\partial h$, and also the interfacial tensions. Thin-film 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 thin-film equation with the quantities $g(h)$ and $\gamma_{lv}$ determined from DFT to generate collective multi-particle 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 coarse-grained mesoscopic theories such as those described in Refs. Thiele et al. (2009); Robbins et al. (2011); Chalmers et al. (2017b). Thin-film 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 thin-film equation models Grün et al. (2006); Durán-Olivencia et al. (2019). Macroscopic scale approaches include just directly solving the Navier-Stokes equations or via methods such as lattice-Boltzmann 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 single-particle dynamics algorithm and the thin-film 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 liquid-vapour 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 $(x,y,z)$ with the $x$- and $y$-axis being parallel to the solid surface and the $z$-axis being perpendicular to the solid surface. The system is discretised onto a regular grid. To simplify, here we treat the system as being two-dimensional (2D), so that there are no variations in the $y$ direction, but the model can easily be generalised to three dimensions (3D). We define $\sigma$ to be the lattice spacing, which we set to $\sigma=1$, throughout. We define $l_{\mathbb{i}}$ to be the occupation number of lattice site $\mathbb{i}$, where $\mathbb{i}=(i,j)$ is the 2D discrete position vector. We assume that $1\leq i\leq N_{x}$ and $0\leq j\leq N_{y}$, with $j=0$ corresponding to the solid surface. When $l_{\mathbb{i}}=1$ then site $\mathbb{i}$ is occupied by a particle and when it is unoccupied, then $l_{\mathbb{i}}=0$. We typically initiate the liquid in one or two semi-circular 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 $\sigma$ 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 $E$ is modelled by the Hamiltonian

[TABLE]

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 $\varepsilon$. The interaction energy between pairs of particles at lattice sites $\mathbb{i}$ and $\mathbb{j}$ depends on the distance between them and is given by $c_{\mathbb{i},\mathbb{j}}$. Here, we use

[TABLE]

where $NN\mathbb{i}$ and $NNN\mathbb{i}$ denote the nearest neighbours of $\mathbb{i}$ and the next nearest neighbours of $\mathbb{i}$, 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 $c_{\mathbb{i},\mathbb{j}}=0$, for $\mathbb{j}\in{NNN\mathbb{i}}$, then rectangular shaped droplets are liable to be formed, particularly at low temperatures. The choice of $c_{\mathbb{i},\mathbb{j}}$ in Eq. (3) minimises the dependence of the vapour-liquid surface tension on the interface-orientation. The values of $c_{\mathbb{i},\mathbb{j}}$ 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 $V_{\mathbb{i}}$ due to the surface on which liquid is deposited. We model it as follows:

[TABLE]

where $\varepsilon_{w}$ 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 straight-forward: we set $l_{\mathbb{i}}=0$ for all lattice sites with ${j}<1$ 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 ${\cal G}$ on the particles parallel to the surface, for example due to gravity. In this case, a term $j{\cal G}$ is added to the expression for the external potential $V_{\mathbb{i}}$ 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 many-particle moves, generated by considering the dynamics from a thin-film 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:

[TABLE]

where the integer index $n$ labels the sequence of configurations $\{l_{\mathbb{i}}\}$ generated by the algorithm, each having energy $E_{n}=E\left(\{l_{\mathbb{i}}\}\right)$, given by the Hamiltonian (2). The normalisation factor $Z$ is the partition function and $\beta=1/k_{B}T$, where $k_{B}$ is Boltzmann’s constant and $T$ is the temperature. The master equation governing the time evolution of the non-equilibrium probability $P_{n}(t)$ is Landau and Binder (2014):

[TABLE]

where $W_{n\to m}$ is the transition probability for the system to go from state $n$ to state $m$. Eq. (6) has the form of a continuity equation, so the total probability remains normalised $\left(\sum_{n}P_{n}(t)\equiv 1\right)$ for all time. At equilibrium, where $\partial P_{n}/\partial t=0$, Eq. (6) yields

[TABLE]

Since the equilibrium probability distribution function for state $n$ is given by Eq. (5), and similarly for state $m$, the ratio of the probabilities also gives the transition rate ratio

[TABLE]

where $\Delta E=E_{n}-E_{m}$. Below is outlined the Metropolis Monte Carlo algorithm that we implement, which samples with the above ratio of probabilities:

Pick a random particle. 2. 2.

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

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

If $\Delta E<0$, then move the particle to its new position. 5. 5.

If $\Delta E>0$, allow the move with probability $e^{-\beta\Delta E}$.

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 thin-film equation

While the above KMC algorithm correctly samples the equilibrium states with a diffusive dynamics, this is not in fact what non-equilibrium 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 thin-film equation from hydrodynamic theory Thiele (2007). The effects of collective advective droplet hydrodynamic motion is illustrated in Fig. 1. The film height profile $h(x,t_{0})$ for the non-equilibrium droplet at time $t_{0}$ is given by the red dashed line, while the black solid line is the new height profile $h(x,t_{1})$ that the droplet assumes at the later time $t_{1}$. In 3D, the thin-film equation is:

[TABLE]

where $h(x,y,t)$ is the film thickness profile, $\gamma_{lv}$ is the liquid-vapour interfacial tension, $\eta$ is the viscosity and $\nabla=({\partial}/{\partial x},{\partial}/{\partial y})$. The term $\gamma_{lv}\nabla^{2}h$ represents the Laplace or curvature pressure. Here, where we assume the system is 2D, we simply replace $\nabla\to{\partial}/{\partial x}$. $\mathcal{G}$ represents an additional force that may be present (e.g. gravity) acting in the positive $x$ direction. The term $\Pi(h)$ is the disjoining pressure, which describes the effective molecular forces between the substrate and the liquid-vapour 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, $g(h)$, as follows:

[TABLE]

The disjoining pressure is in fact a constrained excess free energy. The equilibrium film thickness is given by the minimum value of $g(h_{min})$, where $h_{min}$ is the value of $h$ that minimises $g(h)$. 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 $h\to\infty$, then the liquid is said to be wetting. Additionally, the equilibrium contact angle $\theta$ is determined by the value of $g(h_{min})$, as follows Churaev (1995); Rauscher and Dietrich (2008); Hughes et al. (2015):

[TABLE]

In the literature, there are many different approximations used for the disjoining pressure $\Pi(h)$, 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 $g(h)$ depends on the state point, i.e. on the particular values of $\varepsilon$, $\varepsilon_{w}$ and the temperature $T$. 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):

[TABLE]

where $a_{1}$, $a_{2}$, $a_{3}$, $a_{4}$, $a_{5}$, $a_{6}$, and $b$ are constants that depend on the state point. We also use the DFT to calculate the liquid-vapour interfacial tension $\gamma_{lv}$, 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 thin-film steps is as follows. After completing a prescribed number of KMC steps, $M_{1}$, we first determine the discretised liquid film height profile from the configuration of occupied lattice sites. This is done by computing for each $i$ the number of occupied lattice sites above the lattice site $(i,0)$ 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 $h(x_{i},0)$, where $x_{i}=i$. If for a certain $x_{i}$ it turns out that $h(x_{i},0)=0$, we set it equal to the precursor film thickness $h_{min}$. This profile is then used as an initial condition for the thin-film equation. We solve the thin-film equation numerically using finite-difference approximations to represent spatial derivatives, and evolve the solution in time using e.g. Euler’s method. We choose a time step, $\Delta t$, and solve the thin-film equation for a certain number of time steps, $M_{2}$. Next, we map the computed thin-film 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 thin-film 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 thin-film steps as described above for a certain number of times, $n$.

III.3 Diffusion coefficient and time scales

An important quantity in the procedure described above is the parameter $\lambda$ proportional to the ratio of the number of KMC attempted particle moves per particle $M_{1}/N$ to the number of thin-film moves $M_{2}$, or more strictly the total time the thin film equation is evolved for at each step, $M_{2}\Delta t$. In dimensionless time units, this ratio becomes:

[TABLE]

Here, $N$ is the total number of the particles in the system and $\tau$ is proportional to the viscosity of the liquid, $\eta$:

[TABLE]

When $\lambda=0$, the system is evolved by using only the thin-film equation, whilst the limit $\lambda\rightarrow\infty$ corresponds to evolving the system using only KMC particle moves. Since the single-particle KMC moves generate diffusive type motion, a large $M_{1}/N$ leads to a higher diffusion coefficient. Therefore, $\lambda\propto D/\eta$, where $D$ is a single-particle 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 $D\approx 2.6\times 10^{-4}\sigma^{2}$ KMC steps*-1*, where $x$ KMC steps means $x$ attempted moves per lattice site. In the present 2D system, one should expect $D$ to be of a similar magnitude in the bulk liquid. As regards the diffusion coefficient for isolated particles on the surface, in this regime $D\sim\sigma^{2}\lambda\propto M_{1}$. Recalling the Einstein relation $D=(6\pi\sigma\eta\beta)^{-1}$, this may be used if seeking to match results from the present model to particular experimental systems. Moreover, this allows to see that the matched value of $\lambda$ is largely determined by the viscosity and the temperature.

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 $g(h)$ and also the liquid-vapour interfacial tension $\gamma_{lv}$. These quantities are then input into the thin-film 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 one-body average density profile, which on a lattice is the quantity:

[TABLE]

where $\langle\,\boldsymbol{\cdot}\,\rangle$ denotes an ensemble average. The approximation for the Helmholtz free energy that we use is

[TABLE]

A derivation of this is given in Ref. Hughes et al. (2014). There is an obvious (mean-field) 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 $\{\rho_{\mathbb{i}}\}$ is obtained by minimising the grand potential

[TABLE]

where $\mu$ 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 $T$ is less than the critical temperature $T_{c}$. 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 $T_{c}$ 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, $T_{c}\approx 1.1\varepsilon/k_{B}$. Recall that the Onsager result is for the standard 2D Ising which just has nearest neighbour interactions, i.e. setting $c_{\mathbb{i},\mathbb{j}}=0$ for $\mathbb{j}\in{NNN\mathbb{i}}$. We see that in the region near $T_{c}$, the mean-field DFT fails, but for temperatures ${k_{B}}T/\varepsilon<1$ 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 lattice-gas model proves to be rather useful. Note that this symmetry is absent in continuum models of fluids. If we replace $l_{\mathbb{i}}=1-o_{\mathbb{i}}$ in the Hamiltonian (2), then we see that this leaves the form of it unchanged. We can think of $o_{\mathbb{i}}$ as being a ‘hole’ occupation number. As a consequence of this particle-hole symmetry, we therefore have the following relation between the density of the liquid $\rho_{l}$ and the vapour $\rho_{v}$ at coexistence:

[TABLE]

From Eq. (16), we obtain the Helmholtz free energy per unit area in the absence of any external field:

[TABLE]

where $A$ 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 $F/V$ in 3D, where $V$ is the volume. Note also that we have used the result $\varepsilon\sum_{\mathbb{i},\mathbb{j}}c_{\mathbb{i},\mathbb{j}}=6\varepsilon/2$, to obtain Eq. (19). This is the sum over the interactions of a particle with all of its neighbours [see Eq. (3)]. The factor $1/2$ is present to avoid double counting the interactions between each of the pairs. The pressure in the system is therefore Hansen and McDonald (2013)

[TABLE]

Invoking Eq. (18) and solving $p(\rho)=p(1-\rho)$ leads to the binodal curve

[TABLE]

This is the curve displayed in Fig. 2. The maximum on the binodal corresponds to the critical point temperature, above which there is no vapour-liquid phase separation. Substituting the critical density $\rho=1/2$ into Eq. (21) gives the critical temperature predicted by the DFT to be $T_{c}=1.5\varepsilon/k_{B}$, while from our KMC simulations we find that in reality it is $T_{c}\approx 1.1\varepsilon/k_{B}$. We can also calculate the chemical potential

[TABLE]

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

[TABLE]

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

[TABLE]

Using this together with Eq. (19) we obtain the following result for the spinodal

[TABLE]

which is the curve plotted in Fig. 2. Note that the particle-hole symmetry discussed above is manifested in all the curves in Fig. 2. Of course, real continuum liquids do not have this symmetry. This symmetry can be broken by introducing a 3-body term $\propto\sum_{\mathbb{i},\mathbb{j},\mathbb{k}}\tilde{c}_{\mathbb{i},\mathbb{j},\mathbb{k}}l_{\mathbb{i}}l_{\mathbb{j}}l_{\mathbb{k}}$ in the Hamiltonian (2) – see Ref. MB85 for more details. We do not pursue this possible extension to the model here.

IV.2 Properties of the vapour-liquid 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 vapour-liquid interface, calculated over a range of different temperatures. The density profiles are calculated by minimising the grand potential (17) with the Helmholtz free energy $F$ given by Eq. (16) and the chemical potential $\mu$ set to be the value at coexistence, given in Eq. (23). To calculate these we set the external potential $V_{\mathbb{i}}=0$ and assume that the profiles only vary in one direction (along the $z$-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 vapour-liquid interfaces, but only one of these is displayed in Fig. 3.

Having calculated these profiles, we can then calculate the liquid-vapour interfacial tension $\gamma_{lv}$. 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 $2L\gamma_{lv}$, where $2L$ 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 $\Omega$ for the system containing the interface, the interfacial tension is then

[TABLE]

where

[TABLE]

is the grand potential for a 2D box having the same area $A$, but only containing one of the two phases, and where $p$ is pressure. For example, when $\beta\varepsilon=2.5$, this approach gives the vapour-liquid interfacial tension to be $\beta\gamma_{lv}=2.48$.

In Fig. 4 we display the surface tension of the planar liquid-vapour interface plotted as a function of temperature. We see that as the temperature $T$ approaches from below the critical temperature $T_{c}=1.5\varepsilon/k_{B}$, then $\gamma_{lv}\rightarrow 0$, due to the fact that the difference in the density of the two coexisting phases goes to zero as $T\to T_{c}$. In addition to the full DFT results for $\gamma_{lv}$, in Fig. 4 we also display the results from a so-called ‘sharp-kink’ 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 $\rho_{l}$, to the vapour with density $\rho_{v}$, 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 sharp-kink type approximation that assumes that in the liquid $l_{\mathbb{i}}=1$ 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 $l_{\mathbb{i}}=0$, with a step-like 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 sharp-kink approximation for all temperatures and also with the second sharp-kink approximation for temperatures $k_{B}T/\varepsilon<0.8$, indicating that in this regime the major contribution to the interfacial tension is energetic, rather than entropic.

The calculations of the solid-vapour and solid-liquid interfacial tensions, $\gamma_{sv}$ and $\gamma_{sl}$, respectively, are done in a similar manner to the calculation of $\gamma_{lv}$ 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)

[TABLE]

where $L$ is the length of the interface with the wall, $\Omega_{0}$ is given by Eq. (27), $\delta\mu=(\mu_{coex}-\mu)$, $g(h)$ is the binding potential and $\Gamma$ is the adsorption, given by

[TABLE]

which is the excess (over bulk) number of particles adsorbed on the interface. Note that strictly we should consider $g$ to be a function of the adsorption $\Gamma$, not of $h$, since the latter is not a well defined quantity Hughes et al. (2015). However, when $h$ is large Eq. (29) gives

[TABLE]

Therefore, throughout this paper we discuss $h$ and $\Gamma$ almost interchangeably, and treat Eq. (30) as an equality which defines the value of $h$.

When the system is at vapour-liquid phase coexistence we have $\delta\mu=0$ and so the last term in Eq. (28) is zero. We also see that when $h\to\infty$ we must have $g(h)\to 0$, since when $\delta\mu=0$ the only contributions to $\omega_{ex}$ must be those from the wall-liquid and the liquid-vapour interfaces, i.e. $\gamma_{sl}$ and $\gamma_{lv}$, respectively. These considerations also enables us to see what $g(h)$ really is: it is the additional contribution to $\omega_{ex}$ 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 $g(h)$ has a minimum at a finite value of $h$, then $g(h)$ ‘binds’ these two interfaces together, hence the name given to $g(h)$. 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 $\delta\mu=0$, this is obtained by solving $\Pi(h)=-\frac{\partial g}{\partial h}=0$. For all other values of $h$, the system is out of equilibrium and therefore to determine $g(h)$ for other values of $h$, one must either move away from coexistence or apply a constraint, the constraint being that there is a specified thickness $h$ of liquid at the wall.

The approach used to calculate $g(h)$ 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 $\Gamma$ is a specified value, which via Eq. (30) is equivalent to minimising subject to the constrain that the film thickness $h$ 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 $\Gamma$, enabling us to calculate $g(\Gamma)$ over the whole range of values of $\Gamma$. This data set is subsequently fitted to the form given in Eq. (12) and then input into the thin-film 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 $\beta\varepsilon=0.85$, which corresponds to a temperature that is not that far below the critical temperature $T_{c}$. In Fig. 5(a) are results for a sequence of different values of the wall attraction strength $\varepsilon_{w}$. We see that for small values of $\beta\varepsilon_{w}$, which corresponds to a weakly attracting solvophobic wall, the global minimum of $g(\Gamma)$ is at a small value of $\Gamma$, i.e. the liquid does not wet the wall. As $\beta\varepsilon_{w}$ is increased, there is a first order wetting transition when $\beta\varepsilon_{w}\approx 1.36$ and for $\beta\varepsilon_{w}>1.36$, the global minimum in $g(\Gamma)$ is at $\Gamma\rightarrow\infty$; i.e. for these values of $\beta\varepsilon_{w}$ 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 $h\to\infty$. In Fig. 5(b) we show the result for $\beta\varepsilon_{w}=1.02$ to demonstrate the excellent agreement between the numerical results for $g(\Gamma)$ 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 $\Gamma$ Hughes et al. (2015, 2017); Tewes et al. (2017).

In Fig. 6 we display results for the binding potential for $\beta\varepsilon=2.5$, 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 $\varepsilon_{w}$ and we see that for $\beta\varepsilon_{w}\approx 5$ the minimum in $g(\Gamma)$ close to the wall dissapears, so that for $\beta\varepsilon_{w}>5$ the global minimum in $g(\Gamma)$ is at a larger value of $\Gamma$. However, we note that in this case, due to the presence of oscillations in $g(\Gamma)$ that do not decay in amplitude, there is an infinite sequence of minima in $g(\Gamma)$ as $\Gamma\to\infty$. 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 $g(h)$ corresponds to the vapour-liquid 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 $\beta\varepsilon_{w}=2.5$.

V Contact angle and static droplet profiles from KMC

As described in the previous section, it is straightforward to calculate the three interfacial tensions $\gamma_{lv}$, $\gamma_{sl}$ and $\gamma_{sv}$ using our lattice DFT (16). These can then be inserted into Young’s equation (1) to determine the equilibrium contact angle $\theta$. This is also the contact angle that equilibrium droplet solutions of the thin-film equation have, since $\theta$ is determined by $g(h)$, as Eq. (11) shows. It is important that the value of $\theta$ determined by $g(h)$ in the thin-film 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 thin-film 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 $T_{c}$ 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. 7 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 $\beta{{\varepsilon}_{w}}$ and for the temperature $k_{B}T/\varepsilon=1/2.5$. For these simulations we do not include any of the thin-film collective dynamics, since at this stage we are solely interested in equilibrium properties. In each case there are $N=897$ particles in the system, which has a total domain size of $100\times 100$ lattice sites and with periodic boundary conditions between the left and right sides. We initiate the system by setting $l_{\mathbb{i}}=1$ for all the lattice sites within a semi-circular region in the middle upon the substrate and we set $l_{\mathbb{i}}=0$ outside the semi-circle. We see that for the largest value of the wall attraction $\beta{\varepsilon}_{w}=5$, the liquid spreads out over the substrate; i.e. it wets the substrate. For smaller values of $\beta{{\varepsilon}_{w}}$, the liquid forms a droplet on the substrate with an increasingly large contact angle as $\beta{{\varepsilon}_{w}}$ is decreased; i.e. for these cases the liquid is partially wetting. This behaviour and the observed wetting transition at $\beta{{\varepsilon}_{w}}\approx 5$ is in very good agreement with what we expect based on Fig. 6. Note too that the droplet profiles in Fig. 7 are similar to the droplet profiles calculated in Refs. Hughes et al. (2014, 2015) using DFT for a very similar lattice-gas 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. 7, following e.g. Ref. Chalmers et al. (2017a), we can calculate the contact angle by fitting a circle to the location of the liquid-vapour interface, which is defined as the points determined by linear interpolation between pairs of neighbouring lattice sites, where the density $\rho=0.5$. Results for the contact angle determined via this method are plotted in Fig. 8 for the temperature $\beta\varepsilon=2.5$ and for a range of different values of the wall attraction $\beta\varepsilon_{w}$. 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. 8 show that the contact angle $\theta=0$ for $\beta\varepsilon_{w}>5$, where the liquid wets the wall. As $\varepsilon_{w}$ is decreased, decreasing the strength of the attraction between the particles and the surface so that the surface becomes more solvophobic, the contact angle $\theta$ increases and the liquid increasingly beads-up upon the surface.

VI Results for non-equilibrium droplets on surfaces

The results in the preceding sections essentially constitute a series of checks that the model has the desired, correct and self-consistent thermodynamics. In this section, we present results for the non-equilibrium 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 non-attracting ($\varepsilon_{w}=0$) wall on the top.

We first present results for the behaviour of two droplets using only KMC simulations. The total domain is of size $120\sigma\times 60\sigma$, with periodic boundary conditions between the left and right sides. The system is initiated by setting $l_{\mathbb{i}}=1$ for all the lattice sites within two semi-circular regions of radii $12\sigma$ with the centres of the corresponding circles at the sites $(40,1)$ and $(80,1)$, and we set $l_{\mathbb{i}}=0$ outside of these two semi-circles. This gives $N=368$ particles in the system.

The results in Fig. 9 are for the temperature $k_{B}T/\varepsilon=1/1.5$ and we choose the wall attraction parameter $\beta\varepsilon_{w}=2.1$. The contact angle for these parameter values is $\theta\approx 60^{\circ}$. The different panels in Fig. 9 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 $\theta=90^{\circ}$ closer to $\theta\approx 60^{\circ}$, the equilibrium value (c.f. Fig. 7). 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. 9 (i.e. for the temperature $k_{B}T/\varepsilon=1/1.5$), 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. 9 it is clearly the Ostwald mode that is dominant. For more background on Ostwald ripening, see e.g. Ref. Ratke and Voorhees (2013).

Figure 10 correspond to a lower temperature $k_{B}T/\varepsilon=0.4$ and we choose the wall attraction parameter $\beta\varepsilon_{w}=3.5$ so that the equilibrium contact angle is still $\theta\approx 60^{\circ}$. 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 thin-film equation dynamics

We now consider behaviour of the system when the fluid dynamics governed by a combination of the single-particle KMC moves in conjunction with the collective particle moves generated by the thin-film equation. We use the same domain and the same initial configuration as in the previous subsection. Figure 11 corresponds to the same temperature and wall attraction strength as in Fig. 9. For this simulation, we used $n=600$ cycles of $M_{1}=10^{4}$ KMC steps followed by $M_{2}=10^{4}$ thin-film steps with the time step $\Delta t=10^{-5}$, which gives $\lambda=272$, where $\lambda$ is the parameter defined by Eq. (13). As in Fig. 9, the system coarsens into one droplet. However, the speed of coarsening is increased when the thin-film steps are included, which is not too surprising, since the total number of possible moves has been increased by the addition of the thin-film 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 left one in this simulation) keeps decreasing, but the size of the other droplet (the right one in this simulation) starts to increase. As in the case in Fig. 9, the centres of the droplets in Fig. 11 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. 11 how the numbers of particles that cross the vertical line in the middle of the domain at $i=60$ varies in time. These are split into two populations: (i) those that migrate over the surface at height $j\leq h_{i}$ ($N_{s}$, see the blue line) and (ii) those that travel through the vapour at height $j>h_{i}$ ($N_{v}$, see the red line). Increasing/decreasing $N_{s}$ or $N_{v}$ means that the particle moves to the left/right, respectively, crossing the $i=60$ 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 thin-film equation would not allow us to capture this, and this justifies the importance of our hybrid model which combines both the thin-film equation governed dynamics with the KMC diffusive dynamics.

Figure 12 corresponds to the same values of the temperature and the wall attraction strength parameter as in Fig. 10. For this simulation, we used $n=3000$ cycles of $M_{1}=10^{4}$ KMC steps followed by $M_{2}=10^{4}$ thin-film steps with the time step $\Delta t=10^{-5}$, which gives the same value of $\lambda$ as for the results in Fig. 11. As in the higher temperature case in Fig. 11, we find that the speed of the coarsening is increased when the thin-film steps are included. As in the pure KMC case in Fig. 10, the droplets initially spread over the surface to reach a contact angle close to the equilibrium value $\theta\approx 60^{\circ}$, but with only a small amount of evaporation due to the low temperature. As time proceeds, the size of the right droplet starts to decrease and the left 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. 12 the time evolution of the numbers of the particles that cross the vertical line in the middle of the domain, at $i=60$, that migrate over the surface, $N_{s}$, and through the vapour $N_{v}$. 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 ${\cal G}=0.2$ to the right on each particle. The domain is a box of size $100\sigma\times 50\sigma$, 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 $l_{\mathbb{i}}=1$ for all the lattice sites within a semi-circular region of radius $14\sigma$ with the centre at the site $(25,1)$, and we set $l_{\mathbb{i}}=0$ outside of this semi-circle. This gives $N=305$ 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 pile-up of particles on the right-hand boundary wall.

Figure 13 shows results for when the dynamics is solely generated by the KMC algorithm, with no thin-film dynamics. The initial configuration consists of the semi-circular droplet displayed in the top left panel, which has contact angle $\theta\approx 60^{\circ}$. We set $\beta\varepsilon=1.5$ and the wall attraction strength $\beta\varepsilon_{w}=2.1$. We see from the results in Fig. 13 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 ${\cal G}\neq 0$. 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 ${\cal G}=0.2$, 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 14 shows snapshots from a simulation where all the model parameters are the same as for the simulation in Fig. 13 except we now include some collective multi-particle moves generated by considering the thin-film equation (9) with ${\cal G}=0.2$. We set $M_{1}=1.5\times 10^{3}$, $M_{2}=10^{3}$ and $\Delta t=10^{-6}$, which gives $\lambda=4918$. We see that the thin-film moves enable the droplet to migrate across the surface under the lateral driving. By comparing with the previous result in Fig. 13, it is clear that this dynamics must be a collective multi-particle (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 thin-film moves.

In Fig. 15 we show results for a system with an even higher proportion of particle moves generated by the thin-film equation. In this case we have $M_{1}=1.5\times 10^{3}$, $M_{2}=10^{5}$ and $\Delta t=10^{-6}$, which gives $\lambda=0.49$. We again observe that the thin-film 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 rain-droplets sliding down a window. The appearance and bifurcations of such tails is discussed in Ref. Wilczek et al. (2017) for the thin-film equation.

VII Concluding remarks

We have developed a hybrid model for droplet dynamics on surfaces, that combines both single-particle diffusive dynamics via KMC and advective hydrodynamic multi-particle moves via a thin-film 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 coarse-graining that we performed in order to bridge from the microscopic Hamiltonian (2) to the mesoscopic thin-film 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 Lenard-Jones 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 thin-film 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 straight-forward 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 thin-film 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 lattice-DFT models of this type. There has previously been some testing of the mean-field 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 Ostwald-type 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 thin-film 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 many-particle hydrodynamic motion is essential for droplets to slide over surfaces. The thin-film 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.

For the driven droplets in Fig. 13, where the dynamics is solely due to the KMC ($\lambda=\infty$), we see that the advancing and receding contact angles appear to be the same as the equilibrium contact angles. However, when $\lambda$ is finite (Figs. 14 and 15), clearly the advancing and receding contact angles are different. It would be of interest to determine how these angles depend on $\lambda$ and on the strength of the driving ${\cal G}$ and to compare with predictions in the literature, such as Refs. BH69 ; V76 . However, to get the same amount of data in such non-equilibrium situations as was used to determine the equilibrium contact angle values with the statistical precision displayed in Fig. 8 would require numerous repeated simulations. It is a much more complex scenario than the equilibrium case and so would require a very careful analysis of the results to match together all the different simulations, for robust conclusions to be obtained. This would enable the calculation of the non-equilibrium ensemble-average density profiles (the analogue of Fig. 7), which would yield further information about the complex relationship between the contact line structure (density profile) and the wetting dynamics Wang . Therefore, we leave this for future work.

One limitation of our model, inherited from using the thin-film equation (9), is that for contact angles $\theta\geq 90^{\circ}$ the model becomes inaccurate. The KMC part of the model works fine in this regime, it is only the thin-film equation part that fails. Thus, for example, the model could not be applied in its present form to study droplets on superhydrophobic surfaces. A further limitation of the model is that in certain situations it does not correctly capture the fluid dynamics within the droplet. For example, sometimes there can be a circulation of the liquid within a droplet as it slides Thampi2013 ; Karapetsas2016 . Some aspects of circulatory motion is described by the thin film equation, see e.g. Dmitri . However, in our model if one were to tag a particular lattice particle, one would not observe circulatory motion, since our collective particle moves do not couple to the liquid velocity field, and are determined from the variations in the hight profile.

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. We are also grateful for constructive comments from two anonymous referees.