Thermodynamic Stability of Volatile Droplets and Thin Films Governed by Disjoining Pressure in Open and Closed Containers

Distributed thin films of water and their coexistence with droplets are investigated using a capillary description based on a thermodynamic fundamental relation for the film Helmholtz energy, derived from disjoining pressure isotherms and an accurate equation of state. Gas–film and film–solid interfacial tensions are derived, and the latter has a dependence on film thickness. The resulting energy functionals from the capillary description are discretized, and stationary states are identified. The thermodynamic stability of configurations with thin films in systems that are closed (canonical ensemble) or connected to a particle reservoir (grand canonical ensemble) is evaluated by considering the eigenvalues of the corresponding Hessian matrices. The conventional stability criterion from the literature states that thin flat films are stable when the derivative of the disjoining pressure with respect to the film thickness is negative. This criterion is found to apply only in open systems. A closer inspection of the eigenvectors of the negative eigenvalues shows that condensation/evaporation destabilizes the film in an open system. In closed systems, thin films can be stable even though the disjoining pressure derivative is positive, and their stability is governed by mechanical instabilities of a similar kind to those responsible for spinodal dewetting in nonvolatile systems. The films are stabilized when their thickness and disjoining pressure derivative are such that the minimum unstable wavelength is larger than the container diameter. Droplets in coexistence with thin films are found to be unstable for all considered examples in open systems. In closed systems, they are found to be stable under certain conditions. The unstable droplets in both open and closed systems are saddle points in their respective energy landscapes. In the closed system, they represent the activation barrier for the transition between a stable film and a stable droplet. In the open system, the unstable droplets represent the activation barrier for the transition from a film into a bulk liquid phase. Thin films are found to be the equilibrium configuration up to a certain value of the total density in a closed system. Beyond this value, there is a morphological phase transition to stable droplet configurations.


■ INTRODUCTION
In the literature, a distinction is made between thick (β-films) and thin films (α-films). The thermodynamic properties of thin films deviate from those of a bulk liquid phase at the same temperature and chemical potential. This deviation can be modeled by the disjoining pressure, a concept first introduced by Derjaguin in the 1930s. 1 The disjoining pressure describes the interaction between two interfaces in close proximity, such as the top and the bottom of a thin liquid film residing on a solid substrate. 2, 3 The formation of films and droplets at solid interfaces is of importance in many applications. Inside porous media, thin films provide an important mode for fluid transport 4 that is often neglected or under-resolved in flow modeling. 5−7 Both droplet formation and films impact the efficiency of water removal in fuel cells 8 and are important for atmospheric water collection. 9 Thin, nanometer-thick films are key components in paint, coatings, and different lubricants. 1 They are also important in thin-film evaporation 10,11 and boiling heat transfer. 12,13 With an emerging interest in micro- 14 and nanofluidics, 15 condensation on nanostructured surfaces, 9,16 and self-organization and pattern formation, 17−19 it becomes increasingly important to understand how thin films, often in coexistence with droplets, are influenced by confinement. 20 Previous works 21−25 have shown that the stability of heterogeneous structures is strongly affected by the choice of ensemble. For the simplest type of films, thick films with a negligible disjoining pressure, we recently showed that their thermodynamic stability was profoundly different in a closed system (canonical ensemble) and an open system (grand canonical ensemble). 25 In this work, we will study how thin films and their possible coexistence with droplets are influenced by confinement.
In the literature, thin films are often proclaimed to be unstable when dΠ/dh > 0, where Π is the disjoining pressure and h is the film thickness. 1,26,27 It has, however, been pointed out that this criterion is not necessarily valid for confined systems. 20,28 In fact, this has been exploited in computer simulations to calculate the disjoining pressure where it has a positive slope in h. 29 The aim of this work is to compare, consistently and on equal terms, the thermodynamic stability of thin films in combination with droplets in open and closed systems. To this end, we derive a thermodynamic fundamental relation for the film phase from the disjoining pressure combined with an equation of state (EOS) that represents the bulk-phase properties. Pure water is used as an example, but the methodology is equally applicable to other fluids and can straightforwardly be extended to mixtures.
By comparing with previous results from the literature, we will explain how different assumptions in the modeling of thin films make them appear as if in an open or a closed system. 3,20,26,27,30−33 For instance, in his study of volatile films, Sharma 34 used kinetic theory to model the gas−liquid mass transfer rate as a function of the difference between vapor pressure and coexistence pressure. The coexistence pressure was modeled by an extended Kelvin equation (including effects of the disjoining pressure). Subsequently, the mass transfer rate was incorporated into a partial differential equation for the film profile to study time evolution and stationary states of films and droplets on a flat substrate. He found that flat films are unstable only when dΠ/dh > 0, in agreement with the conventional stability criterion.
Dorfler et al., 20 on the other hand, used an effective interface Hamiltonian to describe the mechanical energy of a film. They showed that the mechanical instabilities are suppressed in systems that are smaller than the critical wavelength, in agreement with earlier findings. Although their effective interface Hamiltonian did not account for particle exchange, we will show that many of their findings can be recovered by considering a closed system, since the behavior of the system is then determined by mechanical instabilities of the film−gas interface.

■ FUNDAMENTAL THEORY OF THIN FILMS
In the following, we present the fundamental theory that will be used to describe thin films. Starting with the original definition of the disjoining pressure by Derjaguin, 1 we derive a thermodynamic fundamental relation for the liquid film phase. Next, we show that the fundamental relation is consistent with Derjaguin's well-known relation for the macroscopic contact angle. 2 Results from the literature on the mechanical stability of thin films are briefly reviewed, which will be used in the subsequent stability analysis.
Disjoining Pressure. Since there appear to be somewhat differing interpretations of the disjoining pressure in the literature, 1,2, 35 we start by describing the definition used in this work. When two interfacial regions are brought in close proximity to form a thin film, they experience either attractive or repulsive forces. These forces can be described by the disjoining pressure Π. They lead to anisotropic stresses in the film, manifested by a pressure p ⊥ normal to the interfaces that differs from the pressure p ∥ parallel to them.
To define the disjoining pressure, we will use the example illustrated in Figure 1. In the left container, there is a thin liquid film and a bulk gas phase with pressure p g , while in the right, there is a bulk liquid phase with pressure p b . The film phase in the left container is connected to, and in chemical equilibrium with, the bulk liquid phase in the right through a tube. The system is also in mechanical force balance, but the pressures p g and p b are, in general, different.
Inside the region enclosed by dashed lines in Figure 1, the gas−film interface has negligible curvature. Due to the force balance, the normal pressure in the film is equal to the gas pressure The disjoining pressure is then defined, in accordance with Churaev et al. 1 (p 36), as the difference between the normal pressure in the film and the pressure in the bulk liquid phase with which the film is in chemical equilibrium. This may be expressed as (2) Fundamental Relation. We will next derive a fundamental relation for the liquid film phase, i.e., the Helmholtz energy of an infinitesimally small film section. The fundamental relation may be divided by the infinitesimal substrate area covered by the film section to obtain an expression that can be integrated over the entire substrate to calculate the total Helmholtz energy of a distributed film with varying thicknesses.
Consider a section of the film covering a small area A fs of a flat solid surface, as illustrated by the dotted white lines in Figure 2.
Since the section is small, any variation in the film thickness across it may be considered small with respect to (w.r.t.) the thickness h in the middle and can be approximated as linear. The Helmholtz energy differential for the film section may then be expressed as Figure 1. Two connected containers, where the left contains a thin liquid film (blue) and a bulk gas phase (white) with pressure p g and the right contains a bulk liquid phase with pressure p b . As indicated by the curved liquid menisci along the walls of the left container, p g ≠ p b in general. The disjoining pressure is Π = p g − p b . The region enclosed by dashed lines is drawn in Figure 2. Any effects of gravity have been neglected. Since the section is small, any variations in the film thickness across it may be considered small and linear. The symbol h refers to the thickness in the middle of the film section. A somewhat exaggerated slope in the gas−film interface is used to illustrate that, in general, A gf ≠ A fs .
Langmuir pubs.acs.org/Langmuir Article where S f is the entropy, T is the temperature, μ f is the chemical potential, N f is the number of particles, V f is the volume of the film, and A gf is the gas−film interfacial area. The interfacial tensions of the gas−film and film−solid interfaces are γ gf and γ fs , respectively. The reason for the appearance of p ⊥ in eq 3 is that the only way to change V f at constant A gf and A fs is to change the film thickness h. The work required to change h must be performed against the normal pressure in the film. The Helmholtz energy can be expressed as a function of T, V f , A gf , A fs , and N f by integrating from a thick film of volume V ∞ f , which has the desired areas A gf and A fs and is unaffected by the disjoining pressure, to a thin film with volume V f Herein, γ ∞ gl and γ ∞ ls are the standard macroscopic gas−liquid and liquid−solid interfacial tensions, respectively, and F b is a bulkphase description of the fluid Helmholtz energy as given, e.g., by an EOS. Since the integration is performed with constant areas A gf and A fs , we may replace dV f with A fs dh. Replacing the integrand in eq 4 by p ⊥ = p b + Π from eq 2 and absorbing the resulting integral over p b into F b , we obtain where, implicitly, h = V f /A fs . For convenience, we introduce the shorthand notation that will be used in further derivations. Equation 5 is a fundamental relation for the film section, and many other thermodynamic quantities may be derived from it by differentiation. 36 In particular, the chemical potential is meaning that the chemical potential of the film is the same as in a bulk liquid at the same temperature and density. This is a consequence of the disjoining pressure being a function of the film thickness only. Hence, μ f is a function of T, V f , and N f . Furthermore, the gas−film interfacial tension is the same as the macroscopic gas−liquid interfacial tension The film−solid interfacial tension, on the other hand, becomes a function of the film thickness through the action of the disjoining pressure ls f gf f (9) From eq 5, we observe that the Helmholtz energy is first-order Euler homogeneous in V f , A gf , A fs , and N f f where β is a variable that describes a scaling of the system size. Differentiating w.r.t. β and setting β = 1 yields the Euler relation Choosing instead β = 1/A fs yields the Helmholtz energy per area of solid substrate where h = V f /A fs , Γ = N f /A fs , and α = A gf /A fs . From the Euler relation (eq 11), we get Macroscopic Wetting Properties. We previously showed that the film−solid interfacial tension is a function of the film thickness; see eq 9. This has implications for the wetting properties of a macroscopic liquid droplet on a solid surface covered by a thin film.
The spreading coefficient associated with the spreading of a thick liquid layer, whose surface tensions are unaffected by the disjoining pressure, onto a solid surface covered by a thin film is  (14) since γ γ = ∞ gf gl according to eq 8. For a given disjoining pressure isotherm, ψ is a function of h only. Young's equation for the contact angle θ results from a balance of interfacial forces. In terms of the spreading coefficient, it can be expressed as Replacing γ fs using eqs 9 and eq 3 in the equations above gives This result is identical to Derjaguin's equation for the macroscopic contact angle. 1−3 Since the droplet is macroscopic, and thus largely unaffected by the disjoining pressure and the Young−Laplace pressure difference, it is reasonable to approximate the film thickness by h 0 , which is such that Π(h 0 ) = 0. This simplifies eq 16 to as discussed in detail in the insightful review by Boinovich and Emelyanenko. 2 An alternative derivation of eq 16 based on a simplified droplet model similar to that of Dorfler et al. 20 is provided in Appendix A.
Langmuir pubs.acs.org/Langmuir Article Linear Stability Analysis. The mechanical stability of thin films, not taking into account condensation and evaporation, has been studied extensively in the literature. 18,20,[30][31][32][33]37 One approach has been to use the film profile equation, see, e.g., refs 18 and 31, to perform a linear stability analysis and determine the growth rate of mechanical disturbances of a thin flat film with different wavelengths λ. To represent an unstable disturbance, λ must satisfy the well-known criterion A discussion and derivation of this equation can be found in ref 18.
The criterion (eq 18) cannot be satisfied for any wavelength when dΠ/dh < 0, and films corresponding to a negative slope in the disjoining pressure isotherm are therefore mechanically stable. When dΠ/dh > 0, on the other hand, the interfacial tension acts to suppress disturbances with short wavelengths, while disturbances with wavelengths longer than λ 0 can grow. Films forming droplets by succumbing to such instabilities are said to undergo spinodal dewetting. 20 Locally stable films forming stable, energetically favorable droplets by overcoming some energy barrier do so through nucleation. If the substrate size is smaller than the shortest unstable wavelength, flat films may be stable even though dΠ/dh > 0. 20,28 ■ MODELS We will next develop models for the two composite systems illustrated in Figure 3: a distributed flat film and a droplet system. Details will be provided for the representation of the disjoining pressure and the equation of state. Water at 293.15 K will be used as the example fluid, for which γ ∞ gl = 0.073 N m −1 . Since γ ∞ ls only adds a constant to the Helmholtz energy that has no qualitative effect on the results, we have set γ ∞ ls = 0. In the section Fundamental Relation, V f , A gf , A fs , and N f refer to the volume, interfacial areas, and particle number of a small section of a film. From this point on, these symbols will refer to the entire distributed film.
Disjoining Pressure Model. There are many different models for the disjoining pressure isotherms, depending on the nature of the interface interactions. These vary from the classical van der Waals curves where Π ∝ 1/h 3 (ref 1, p 99) to more complex curves 38 that exhibit one or more local extrema; see, e.g., Figure 3 in ref 27. Disjoining pressure isotherms are usually modeled by adding terms that account for different types of forces acting between two interfaces, 39 resulting in a plethora of possibilities for combinations of terms. We will restrict the attention to a type of model that has been used to describe water films on a solid graphite surface. 27,40 This model has two terms VdW sr (19) accounting for van der Waals forces and short-ranged forces, respectively. Any effects of substrate size on the disjoining pressure are neglected. 28,41 Other surfaces and fluids are likely to require different terms and parameter values. The van der Waals term is proportional to 1/h 3 VdW 3 (20) where A is the effective Hamaker constant. This is modeled by applying a mixing rule to the liquid and solid Hamaker constants Here, we will use A ll = 4.4 × 10 −20 J and A ss = 4.7 × 10 −19 J. 27,40 The short-ranged contribution is, in line with previous works, 27,40 modeled by an exponential where K sr and L sr are the strength and correlation length of the interactions, respectively. We will use L sr = 0.6 nm 27,40 and consider K sr = 0 and three different negative values of K sr . The four resulting disjoining pressure curves are illustrated in Figure  4. The effect of the negative K sr values is to create a minimum in  Langmuir pubs.acs.org/Langmuir Article the disjoining pressure curve that makes the liquid partially wetting and allows the existence of droplets in equilibrium with thin films. For each isotherm in the figure, the legends state the contact angle θ 0 predicted by eq 17 for a macroscopic droplet in mechanical force balance with a thin film. Equation of State for Bulk-Phase Properties. A necessary component of the models developed here is a thermodynamic description of the bulk phases. For this purpose, any EOS capable of predicting Helmholtz energies of both gas and liquid phases is applicable. We will use the cubic-plusassociation modification of the Soave−Redlich−Kwong EOS. 42 This EOS is implemented in our in-house thermodynamics library. 43 Distributed Film in a Cylindrical Container. We consider a cylindrical container with radius R, height H, and volume V = πR 2 H. Furthermore, we assume for simplicity that the film thickness h is a cylindrically symmetric function of the radial coordinate r. The film−solid interfacial area is then a constant The film volume may be obtained by integrating h over the area (24) and, similarly A homogeneous bulk gas phase fills the container volume that is not occupied by the film.
The container has a fixed volume V, and we will only consider the case when it is connected to a thermal reservoir with constant temperature T. If the container is closed, i.e., it is in the canonical ensemble and contains a fixed number of particles N, equilibrium is defined by the state that maximizes the total entropy of the system in the container and the reservoir subject to these constraints on T, V, and N. This is equivalently described as a minimum in the Helmholtz energy for the system. 36 Equilibrium in an open container, which also has a fixed volume V and temperature T, but where the gas phase can exchange particles with an external particle reservoir at constant chemical potential, corresponds to a minimum in the grand canonical energy Ω.
The total Helmholtz energy for the combined film-and-gas system is The Helmholtz energy of the gas is Inserting for f f using eq 13, and subsequently combining with eqs 8, eq 4, and eq 2, we get (30) where α = +ḣ 1 2 with ḣ= dh/dr. We restrict our attention to cases where there is chemical equilibrium within the film phase, i.e., μ f is uniform, and thus Adding together eqs 28 and 5, we obtain the total Helmholtz energy This expression is a functional of h and a function of N f . Assuming chemical equilibrium, i.e., μ g = μ f , at a given value of where F μ°i s a constant and the integrand is We observe that F μ closely resembles an effective interface Hamiltonian as presented, e.g., by MacDowell 28 and Dorfler et al. 20 A film profile at mechanical equilibrium must have a vanishing first variation of F μ and satisfy the Euler−Lagrange equation 44 We differentiate L and insert into the Euler−Lagrange equation to get are the curvatures of the film−gas interface. Equation 36 is the augmented Young−Laplace equation for the film. It is a second-order ordinary differential equation that may be solved to yield a film profile that satisfies the mechanical force balance. Solving eq 36 requires two boundary conditions. The first is a homogeneous Neumann boundary condition at r = 0̇= which follows from the cylindrical symmetry of the problem. The second is a Dirichlet condition at r = R where h R is some specified constant value.
Langmuir pubs.acs.org/Langmuir Article For a given value of h R , the Euler−Lagrange equation may have several solutions. We will consider two qualitatively different types of solutions. The first type is the trivial solution where h is a constant. This corresponds to a flat film (Figure 3a). In this case, the curvature terms in eq 36 are zero and Δp = p g − p b is equal to Π, which is constant along the film profile. The other type of solution is a droplet and a film, as illustrated in Figure 3b. In this case, the disjoining pressure and curvature terms vary along the film profile, but their sum remains constant. Crucially, as the gas−liquid interface of the droplet approaches the solid surface, the disjoining pressure term balances the curvature terms in such a way that the gas−liquid interface flattens into a flat film.
The grand canonical energy can be calculated from the Helmholtz energy by The criteria for a stationary state in Ω are the same as for F. However, the stability of these stationary points, i.e., whether they are minima, maxima, or saddle points, is likely to differ.

■ NUMERICAL METHODS
In this section, we will describe the numerical methods used to solve the models presented. This amounts to identifying stationary states of the two energy functionals, F and Ω, and determining their stability by characterizing them as maxima, minima, or saddle points.
To identify stationary states, we adopt the same strategy as described in ref 25. We first find states that are in mechanical force balance by solving eq 36, the Euler−Lagrange equation. The solution specifies the geometry of the film configuration, i.e., the thickness h as a function or r, and corresponds to a specific value of Δp = p g − p b , the difference between the gas phase pressure and the pressure of the hypothetical bulk liquid phase with which the film would be in chemical equilibrium. A phase equilibrium calculation is performed next, where equality of the chemical potentials is used to determine the number of particles in the film and gas phases. The next step after determining a stationary state is to evaluate its stability. To this end, we discretize the Helmholtz and grand canonical energy functionals and examine the eigenvalues of their Hessian matrices.
Phase Equilibrium Calculations. Any film profile in mechanical force balance satisfies eq 36 for a specific value of Δp = p g − p b . Through eq 24, the film profile also specifies the volumes V f and V g . Phase equilibrium for a particular film profile can therefore be determined by finding the particle numbers N f and N g that give the necessary pressure difference and equality of the film and gas chemical potentials. This is accomplished by solving the nonlinear system of equations°°= where and μ b = μ f according to eq 7. The pressures, chemical potentials, and the derivatives required to compute F and its Jacobian are provided by the EOS. The scaling quantities are°= where R is the universal gas constant.
The nonlinear system defined by eq 42 was solved using Newton's method. Initial guess values for N f and N g were obtained by a standard phase equilibrium calculation 45,46 at the specified temperature and saturation pressure.
Solving the Euler−Lagrange Equation.
The Euler−Lagrange equation for the film eq 36 was solved to obtain film profiles in mechanical force balance. This equation is a second-order ordinary differential equation that has two boundary conditions, one at r = 0 and the other at r = R. Together, these boundary conditions and eq 36 constitute a two-point boundary value problem, which was solved using solve_bvp from Scipy's integrate module. 47 Discrete Description of the Distributed Film. Following the procedures described above results in stationary states of the Helmholtz and grand canonical energy functionals. These are states that have a vanishing first variation for any perturbation of the film profile. However, these procedures do not give any information about whether or not the stationary states are stable. This information is contained in the second variation or higher-order variations if the second variation is zero. A stationary state is a minimum if the second variation is positive for any perturbation of the film profile. To evaluate this numerically, we use a discrete approach similar to Gjennestad and Wilhelmsen. 25 The idea is to represent the continuous function h(r) by its values at discrete points and then use a quadrature rule to approximate the functional integrals. Stability of the stationary states can then be determined by considering the eigenvalues of a Hessian matrix, as will be described in the section Stability Analysis.
An equivalent way of expressing the Helmholtz energy of the distributed film model (eq 32) is The functional A gf may then be approximated by the sum where the vector of unknowns x includes N f in addition to h The Jacobian vector of the discretized Helmholtz energy can be shown to be and for i ∈ {2, M + 1}. The Hessian matrix may be found by further differentiation. Details of the procedure for calculating the derivatives of A gf , W, and V f w.r.t. x i and a validation of the discretization procedure are given in ref 25 and its accompanying Supporting Information. An analogous procedure is applied to obtain a discrete approximation of the grand canonical energy and its derivatives. Stability Analysis. Stationary states x* obtained with the above procedures will have Jacobian vectors equal to zero for both F and Ω. The change in, e.g., F due to a small perturbation dx of the stationary state x* is thus determined by the Hessian matrix of F The Hessian matrix is symmetric and can therefore be decomposed as where the matrix Λ is a diagonal matrix of eigenvalues and Q is a matrix whose column i is the eigenvector q (i) corresponding to the eigenvalue Λ (i) . We will use the convention that all eigenvectors q (i) have length 1 in the L 2 -norm. Since the Hessian is symmetric, the eigenvectors are orthogonal. 25 The stationary state x* is locally stable in a closed container and represents a local minimum in F if all of the eigenvalues of the Hessian matrix Λ (i) are positive. On the other hand, if one or more eigenvalues are negative, it is possible to choose the perturbation dx along one of the associated eigenvectors q − (i) . The subscripted minus sign indicates the negative eigenvalue. In this case, we choose dx = dx − ∝ q − (i) and this will give a negative dF. The stationary state x* is then not a local minimum and is considered locally unstable. Stability in open containers is evaluated in a similar manner. Unstable states with both positive and negative eigenvalues are called saddle points, while only negative eigenvalues characterize a maximum.
The eigenvectors associated with negative eigenvalues give information about the direction in configuration space that the system can go to reduce its energy and hence what causes an instability to the stationary state.
Before the calculation of eigenvalues and eigenvectors, it was ensured that the discrete Jacobian vectors of the stationary states obtained by the methods in sections Solving the Euler−Lagrange Equation and Phase Equilibrium Calculations were indeed zero. This was accomplished through a small number of iterations with Newton's method using the Jacobian vectors and Hessian matrices defined by the discrete approach.   Eigenvalues and eigenvectors were calculated with eigh from the linalg module in Numpy. 48 This function again uses the *syevd routines from LAPACK. 49

■ RESULTS
We first discuss the thermodynamic stability of flat films and droplet configurations. Next, we compare the energies of the two types of configuration and study the morphological phase transition between them.
Except where explicitly stated otherwise, we consider a container with R = 20 nm and H = 10 nm. The number of grid points used in the discrete approximation is M = 400.
Flat Films. Figure 5 compares the local stability of flat films as a function of film thickness for four different disjoining pressure isotherms. Stability in a closed container is illustrated in Figure  5a In the open container, films are stable when dΠ/dh < 0. This is in agreement with the conventional stability condition presented in the literature; see, e.g., the review paper by Boinovich and Emelyanenko 2 (or ref 1, p 56). A similar result was also obtained through dynamic considerations by Sharma, 34 assuming a constant gas pressure and a model for the evaporation rate. In the closed container, on the other hand, some films are found to be stable even if dΠ/dh > 0.
Eigenvectors q − (i) of Hessian matrices that are associated with negative eigenvalues correspond to perturbations dx − ∝ q − (i) that the film is unstable against. An unstable stationary state may have one or more such negative eigenvalues. The vectors dx − are composed of two parts: dN − f , the component describing the perturbation of the number of particles in the film; and dh − , a vector representing the perturbation of the film profile. The vectors dh − are here referred to as instability modes. Figure 6 compares the instability modes for the θ 0 = 60°i sotherm and two different film thicknesses in an open container. Specifically, Figure 6 displays the instability modes obtained when h = 1.06 nm. This film thickness is well into the unstable region of the disjoining pressure isotherm with dΠ/dh > 0. Each mode has resemblance to a sinusoidal function with some wavelength and corresponds to a specific number of internal extrema, indicated by their color (black, 0; blue, 1; yellow, 2; etc.). Instabilities in the closed system are similar, except that the zero-extrema mode (black) is absent. By moving along this mode, the open system can reduce its energy by condensing or evaporating the film while the film retains its flat profile. We call this mode a condensation−evaporation instability and find that it is present in the open system whenever the derivative of the disjoining pressure curve is positive. This type of mode thus causes the open system to be unstable whenever dΠ/dh > 0. The other instability modes involve some degree of rearrangement of the film profile, and we therefore call these mechanical instability modes. The stability of the closed system is determined by this kind of instability. We emphasize that the mechanical instabilities also involve some exchange of particles. The distinction between condensation/evaporation and mechanical instabilities is based on whether or not the alteration of the gas− film interface shape occurs. We further observe that the Langmuir pubs.acs.org/Langmuir Article exchange of particles is larger in the open system than in the closed system, which is due to the constrained total number of particles in the closed system. As the film thickness is increased, the number of unstable mechanical modes is gradually reduced. They disappear in order of decreasing number of internal maxima (i.e., cyan first, then magenta, etc.) until there is only one left. As the film thickness is increased further and approaches the end of the unstable interval for the closed container (Figure 5a), the maximum of the mode with a single internal extremum moves toward the center of the container. Figure 6 shows the instability modes for a film of thickness h = 2.68 nm in an open container. This thickness corresponds to the rightmost red dot on the θ 0 = 60°isotherm of the closed container displayed in Figure 5a. The mechanical instability mode resembles in this case a sinusoidal function with a wavelength close to the diameter of the container. For thicker films, there are no unstable modes in the closed system and it is stable. In the open system however, the condensation/ evaporation instability persists until the disjoining pressure derivative again becomes negative. The existence of condensation/evaporation inabilities in open systems and their absence in closed systems have recently been reported also for thick films in pores. 25 The mechanical instabilities seen here are of a similar kind as those that cause spinodal dewetting. As shown, e.g., by Dorfler et al., 20 the mechanical instabilities must have wavelengths λ that are smaller than the diameter of the container. The finite extent of the container makes instabilities with longer wavelengths impossible. Furthermore, the wavelengths of unstable perturbations must be bounded from below by the surface tension. Perturbations with very short wavelengths will require the creation of a large amount of surface area, relative to the amount of energy than can be gained from the disjoining pressure by shifting the liquid around. In a translationally invariant system of infinite extent, this lower bound is given by eq 18. A consequence of eq 18 is that dΠ/dh > 0 is a necessary condition for mechanical instabilities. This results in specific intervals of film thicknesses where mechanical instabilities may occur.
To show that these two effects indeed constrain the mechanical instabilities observed here, we estimate for each unstable state the wavelengths of the mechanical instability modes. For the modes with many extrema, the horizontal crestto-trough distance changes slightly along the r-axis. We therefore estimate the wavelength as the distance from the container sidewall to the nearest internal extremum, multiplied by a factor 2. Figure 7 displays the resulting wavelengths as functions of film thickness. Results are presented for disjoining pressure isotherms with θ 0 = 40° (Figure 7a) and θ 0 = 60° (Figure 7b). Figure 7c displays results for the 40°isotherm and a wider container. The figures also show the container diameter 2R (solid horizontal line) and the shortest unstable wavelength λ 0 as predicted by eq 18 (dashed line). The dotted vertical lines indicate where dΠ/dh = 0 and λ 0 → ∞. The disjoining pressure derivative is positive between them and, according to the linear stability analysis, mechanical instabilities should thus occur only for film thicknesses in between these lines.
All measured wavelengths correspond to film thicknesses between the dotted lines, i.e., where dΠ/dh > 0. We observe no mechanical instabilities with wavelengths longer than 2R. As a concrete illustration, compare Figure 7a where R = 20 nm with Figure 7c where the container size is increased to R = 30 nm. In the latter case, the increased container size gives room for instabilities with longer wavelengths. This results in an increased interval of film thicknesses where films are mechanically unstable. The larger container also gives room for a larger number of extrema in instabilities with a given wavelength. For instance, the shortest-wavelength instabilities in Figure 7a have two internal extrema, whereas instabilities with approximately the same wavelength in Figure 7c have four. For both container sizes, it is clear that the mechanical instabilities with the longest wavelength disappear when they reach the 2R-bound.
In addition to being smaller than 2R, the measured unstable wavelengths are longer than λ 0 . Although derived with an assumption of translational invariance that does not apply in the present example, eq 18 appears to provide an accurate estimate for the lower bound for the unstable wavelengths.
From the form of eq 18, we may expect that instabilities of shorter wavelengths are possible when the disjoining pressure derivatives are larger. This is indeed what we observe, e.g., when comparing results from the θ 0 = 40°isotherm (Figure 7a) to results from the θ 0 = 60°isotherm (Figure 7b). The latter disjoining pressure curve extends to larger negative values for the disjoining pressure. This results in unstable modes with shorter wavelengths and a larger number of internal extrema.
In large closed containers, more mechanical instability modes will be present as the upper bound of 2R increases. The films will then become mechanically unstable for a larger part of the interval where dΠ/dh > 0. However, the dotted lines in Figure 7 display the film thicknesses for which the estimated lower bound for unstable wavelengths λ 0 diverges. The divergence of λ 0 means that there may be an interval on the h-axis where a finite container is not large enough to support sufficiently long wavelengths for the film to be unstable no matter how large the container is. However, this region quickly becomes narrow as the container diameter is increased. As an example, consider the θ 0 = 40°isotherm and a container with 2R = 1 μm. The lower bound, λ 0 , is larger than 1 μm only for film thicknesses between 5.41 nm and the thickness for which λ 0 diverges, 5.56 nm. The effect of mechanical stabilization due to finite container size is therefore expected to be small for large containers.
In summary, we find that the criterion for thermodynamic stability of films provided in the literature, dΠ/dh < 0, applies only to open systems due to a condensation/evaporation instability that is present whenever this criterion is not satisfied. The stability of films in closed systems is governed by mechanical instabilities of a similar kind as those responsible for spinodal dewetting in nonvolatile systems. Similar to nonvolatile films, 20 we find that films in small closed containers may be stable even though dΠ/dh > 0 due to the finite size of the container.
Droplets and Films. We shall next compare the thermodynamic stability of a thin film in coexistence with a droplet in an open and a closed system. The stability of different droplet configurations in a closed container is shown in Figure 8 for the θ 0 = 60°isotherm. Here, the configurations with the largest droplets are stable. As the droplet size is reduced, the contact angle decreases and the film thickness increases. Thus, the droplet configurations gradually converge to a flat film as the film thickness approaches the minimum on the disjoining pressure curve where dΠ/dh| h=h R = 0. At some point before the disjoining pressure in the thin-film part of the configurations reaches this minimum, however, the droplets become unstable. As we shall see in the next section, this point is associated with a minimum in the total density in the container. A similar behavior was obtained also for the θ 0 = 40°isotherm. Langmuir pubs.acs.org/Langmuir Article Dorfler et al. 20 found that stable nonvolatile droplets eventually became unstable if the substrate size was increased by making the film-covered substrate area larger while keeping the size of the droplet fixed. We observe the same for the droplets in a closed system studied here.
If the container is open instead of closed, then all of the droplet configurations in Figure 8 are unstable.
The unstable droplet states have one negative Hessian eigenvalue, in both the open and the closed container. The unstable droplets thus represent saddle points in both free energy landscapes, and the systems can reduce their energies by moving along (or against) one particular direction given by the associated eigenvector. Examples of instability modes dh − are shown in Figure 9. These modes are similar for the open and closed systems. However, as discussed for the flat film, the unstable perturbations dx − ∝ q − (i) for the open system involve a larger degree of particle exchange with the gas phase.
By initially perturbing the unstable droplets in the open system along the unstable eigenvector dx − ∝ q − (i) , we were able to form a path in the configuration space, with monotonically decreasing grand canonical energy, to a homogeneous liquid phase filling the entire container. By perturbing the system in the opposite direction, a path to a stable flat film was obtained. Sharma 34 observed a similar behavior by use of dynamic considerations. He assumed a constant gas pressure and a model for the evaporation rate. Time-stationary droplet solutions were then unstable, and droplets slightly smaller than the timestationary state evaporated, shrunk, and eventually became a flat film, while slightly bigger droplets condensed further and grew in size. The unstable droplets in the open system thus represent activation barriers in the energy landscape, which correspond to stationary states with gas pressures above the saturation pressure.
By perturbing the unstable droplets in the closed system, we were able to form paths with monotonically decreasing Helmholtz energy to either a larger stable droplet or a stable flat film. The unstable droplets in the closed system therefore represent the activation barrier for wetting/dewetting of the solid surface.
A notable difference between the paths in the open and closed systems was that in the closed system, the number of particles in the film phase N f changed very little. As an example, the decay of the unstable droplet in Figure 9 to a flat film resulted in a 5.5% reduction in N f in the open system and only a 0.074% reduction in the closed system. A consequence of this is that the assumption of a nonvolatile film phase in the effective interface Hamiltonian approach used, e.g., by Dorfler et al. 20 appears reasonable for water in a closed system at the current temperature. However, this might not be the case for fluids closer to their critical temperature.
Film−Droplet Transition. By comparing the Helmholtz energy of different locally stable flat film and droplet configurations with the same total density ρ = N/V, we determine the equilibrium configuration in a closed system. This comparison is shown in Figure 10 for the θ 0 = 60°isotherm. The film thickness h R at the container sidewall is plotted as a function of the total density for both types of configurations. Unstable configurations are indicated by red markers, locally stable configurations are indicated in green, and stable configurations that represent the lowest Helmholtz energy for a particular value of ρ are shown in blue. The gas spinodal density at the considered temperature is 0.35 mol L −1 , and the liquid spinodal   Langmuir pubs.acs.org/Langmuir Article density is 43 mol L −1 . Since the total densities of the configurations lie between the gas and liquid spinodals, a homogeneous fluid phase can be ruled out as a possible alternative equilibrium state since it is thermodynamically unstable. 46 Figure 10 shows that the point on the dotted line where droplet configurations go from stable to unstable occurs at a minimum in the total density. Below this minimum, no stationary droplet states can be found. A similar effect has been observed for free bubbles and droplets in a closed system. 22, 23 The unstable branch of droplet configurations converges toward the flat film configurations. The two types of configurations merge at the film thickness corresponding to a minimum in the disjoining pressure curves.
The flat films are stable at low densities. They are also the equilibrium configuration up to a certain value of the total density, indicated by the vertical solid line in Figure 10, where there is a morphological phase transition and the equilibrium state becomes a stable droplet configuration. The transition density is 3.03 mol L −1 , a value that depends both on the size and on the shape of the container.
Tracing the blue markers from left to right in Figure 10 shows that in a closed container at equilibrium, the system will initially consist of a thin flat film that grows in thickness as the number of particles is increased. Eventually, a droplet will form when the transition density is exceeded. The formation of a droplet will deplete liquid from the thin-film region such that the film thickness is reduced. This transition, from a stable film state to a stable droplet state, occurs through dewetting by nucleation.
The nucleation regime corresponds to the interval of densities in Figure 10 for which there are three possible stationary states: one stable flat film state and two droplet states. The two droplet states correspond to a small droplet, which is unstable, and a larger droplet, which is stable. An example of three such states at ρ = 3.1 mol L −1 is shown in Figure 11b. A transition from the stable flat film to the stable droplet state will pass through the unstable droplet, as it represents a saddle point in the thermodynamic energy landscape.
In Figure 11a, the Helmholtz energy difference ΔF between the droplet configurations and the corresponding film configuration with the same total density is plotted. Like in Figure 10, the stationary droplet states fork out into a stable and an unstable branch. The Helmholtz energy of the unstable branch is always larger than the stable branch. Furthermore, the Helmholtz energy of the unstable branch approaches that of the flat film as the total density in increased.
The Helmholtz energy differences ΔF of the three states in Figure 11b are indicated by markers of corresponding colors in Figure 11a. It is clear that the large droplet has the lowest Helmholtz energy and is the equilibrium configuration. However, a transition from the stable flat configuration, which passes through the unstable small-droplet configuration, must overcome an energy barrier through an activated nucleation process. The probability of the transition to occur increases as the energy barrier decreases and eventually goes to zero as the total density is increased. The opposite transition, from a large droplet to a flat film, may also occur as an activated nucleation process below the transition density. A similar analysis for the θ 0 = 40°isotherm gives qualitatively similar results but a somewhat higher transition density of 3.88 mol L −1 . Since the activation barriers for the film−droplet transition displayed in Figure 11a are rather small, we expect for this example that the morphological change will occur close to the transition density.

■ CONCLUSIONS
The thermodynamic properties of thin films deviate from those of a bulk liquid phase at the same temperature and chemical potential. This deviation can be modeled by the disjoining pressure. Thin films occur in a variety of applications, such as in flow through porous media, fuel cells, in evaporation as well as in micro-and nanofluidics.
This work is a study of thin films of spatially varying thicknesses and their coexistence with droplets, in open systems and under confinement. Based on Derjaguin's concept of a disjoining pressure Π, which is dependent on film thickness h, and the existing bulk-phase equations of state, we derived a thermodynamic fundamental relation for a thin liquid film phase. From this fundamental relation, the film−gas and film− solid interfacial tensions were obtained, the latter of which was dependent on the film thickness. We verified that Derjaguin's equation for the contact angle of a macroscopic droplet was reproduced by the fundamental relation, and it was next employed to derive a capillary model for a distributed film of varying thicknesses, with a homogeneous gas phase above it.
The model was used to study the thermodynamic stability of stationary states that represent flat films and droplets of water in closed systems (canonical ensemble) or connected to a particle reservoir (grand canonical ensemble). The stationary states were found by solving the Euler−Lagrange equations derived Langmuir pubs.acs.org/Langmuir Article from the Helmholtz energy functional (closed system) or the grand canonical energy functional (open system). By discretizing the energy functionals, the thermodynamic stability of the configurations could be inferred from the eigenvalues of the corresponding Hessian matrices. In closed systems, the flat films were occasionally stable even when the conventional stability criterion often presented in the literature, dΠ/dh < 0, was not satisfied. A closer inspection of the eigenvectors associated with the negative eigenvalues revealed that the stability was governed by mechanical instabilities of a similar kind as those responsible for spinodal dewetting in nonvolatile systems. In line with earlier works, 20,28 the films were stabilized when their thickness and disjoining pressure derivatives were such that the minimum unstable wavelength became larger than the container diameter.
Flat films in open systems, on the other hand, were found to follow the conventional stability criterion. The reason for this was the presence of an additional type of instability whenever the criterion was not satisfied. This type of instability corresponded to condensation/evaporation of the film while the film profile remained flat. This instability was not present in the closed systems.
Droplets in coexistence with thin films were found to be unstable for all considered examples in open systems. In closed systems, they were stable under certain conditions. When two different droplet states were possible at the same density, the large droplet was stable and the small droplet was unstable. Smaller and smaller droplets exhibited smaller and smaller contact angles and eventually converged to flat films.
The unstable droplets in both open and closed systems were found to be saddle points in their respective energy landscapes. In the closed system, they represent the activation barrier for the transition between a stable flat film and a stable droplet. This activation barrier was quantified by evaluating the Helmholtz energies of both the stable and unstable stationary states. In the open system, the unstable droplets represent the activation barrier for the transition from a flat film to a bulk liquid phase through condensation. In a closed system, flat films were found to be the equilibrium configuration up to a certain value for the total density. Beyond this value, there was a morphological phase transition to stable droplet configurations.
The framework presented can readily be extended to study thin-film configurations in multicomponent systems and in different geometries such as around fibers and within porous media.

■ APPENDIX A. MACROSCOPIC DROPLET MODEL
In this section, we derive a simplified model for a macroscopic droplet configuration. This model serves two purposes. (1) Analysis of the macroscopic model will provide an alternative way to derive Derjaguin's relation for the macroscopic contact angle (eq 16), using the fundamental relation for the film phase, that does not presume a balance of interfacial forces. (2) Since the macroscopic model is consistent with a known result, Derjaguin's contact angle relation, it provides an opportunity to validate the distributed model and the discrete approach used to solve it. We will show that the two give the same interface profiles in the large-droplet limit.
We consider a cylindrical container of base radius R, base area A = πR 2 , and height H, giving a total volume V = AH. It contains a spherical droplet, large enough for the pressure inside it to be unaffected by the disjoining pressure, covering part of A, while a thin film covers the rest. A gas phase occupies the remaining container volume. The gas−liquid interface of the droplet is a spherical cap, and the droplet volume can be expressed as where h is the film thickness. The gas volume is In a closed container in contact with a thermal reservoir, equilibrium is a stationary state of the Helmholtz energy. Using the capillary model approach, the Helmholtz energy differential for the system is This work was partly supported by the Research Council of Norway through its Centres of Excellence funding scheme, project number 262644.