Nonlinear dynamics of phase separation in thin films


We present a long-wavelength approximation to the Navier–Stokes Cahn–Hilliard equations to describe phase separation in thin films. The equations we derive underscore the coupled behaviour of free-surface variations and phase separation. We introduce a repulsive substrate-film interaction potential and analyse the resulting fourth-order equations by constructing a Lyapunov functional, which, combined with the regularizing repulsive potential, gives rise to a positive lower bound for the free-surface height. The value of this lower bound depends on the parameters of the problem, a result which we compare with numerical simulations. While the theoretical lower bound is an obstacle to the rupture of a film that initially is everywhere of finite height, it is not sufficiently sharp to represent accurately the parametric dependence of the observed dips or ‘valleys’ in free-surface height. We observe these valleys across zones where the concentration of the binary mixture changes sharply, indicating the formation of bubbles. Finally, we carry out numerical simulations without the repulsive interaction, and find that the film ruptures in finite time, while the gradient of the Cahn–Hilliard concentration develops a singularity.

I Introduction

Below a certain critical temperature, a well-mixed binary fluid spontaneously separates into its component parts, forming domains of pure liquid. This process can be characterised by the Cahn–Hilliard equation, and numerous studies describe the physics and mathematics of phase separation (1); (2); (3); (4); (5). In this paper we study phase separation in a thin layer, in which the varying free-surface and concentration fields are coupled through a pair of nonlinear evolution equations.

Cahn and Hilliard introduced their eponymous equation in (1) to model phase separation in a binary alloy. Since then, the model has been used in diverse applications: to describe polymeric fluids (6), fluids with interfacial tension (7); (8), and self-segregating populations in biology (9). An analysis of the Cahn–Hilliard (CH) equation was given by Elliott and Zheng (10), where they obtain existence, uniqueness, and regularity results. Several authors have developed generalisations of the CH equation: a variable-mobility model was introduced by Elliott and Garcke (11), while nonlocal effects were considered by Gajweski and Zacharias (12). These additional features do not qualitatively change the phase separation, and we therefore turn to one mechanism that does: the coupling of a flow field to the Cahn–Hilliard equation (2). In this case, the Cahn–Hilliard concentration equation is modified by an advection term, and the flow field is either prescribed or evolves according to some fluid equation. Ding and co-workers (7) provide a derivation of coupled Navier–Stokes Cahn–Hilliard (NSCH) equations in which the velocity advects the phase-separating concentration field, while concentration gradients modify the velocity through an additional stress term in the momentum equation. A similar model has been produced by Lowengrub and Truskinowsky (8). Such models have formed the basis of numerical studies of binary fluids (13), while other studies without this feedback term highlight different regimes of phase separation under flow (14); (5); (15). Here, the NSCH equations form the starting point for our asymptotic analysis.

As in other applications involving the Navier–Stokes equations, the complexity of the problem is reduced when the fluid is spread thinly on a substrate, and the upper vertical boundary forms a free surface (16); (17). Then, provided lateral gradients are small compared to vertical gradients, a long-wavelength approximation is possible, in which the full equations with a moving boundary at the free surface are reduced to a single equation for the free-surface height. In the present case, the reduction yields two equations: one for the free surface, and one for the Cahn–Hilliard concentration. The resulting thin-film Stokes Cahn–Hilliard equations have already been introduced by the authors in (18), although the focus there was on control of phase separation and numerical simulations in three dimensions. Here we confine ourselves to the two-dimensional case: we derive the thin-film equations from first principles, present analysis of the resulting equations, and highlight the impossibility of film rupture, once a regularizing potential is prescribed. We use numerical simulations to show that in the absence of this regularizing potential, the film does indeed rupture, an event that coincides with the development of a singularity in the concentration gradient.

Along with the simplification of the problem that thin-film theory provides, there are many practical reasons for studying phase separation in thin layers. Thin polymer films are used in the fabrication of semiconductor devices, for which detailed knowledge of film morphology is required (19). Other industrial applications of polymer films include paints and coatings, which are typically mixtures of polymers. One potential application of the thin-film Cahn–Hilliard theory is in self-assembly (23); (22); (24); (20); (21). Here molecules (usually residing in a thin layer) respond to an energy-minimisation requirement by spontaneously forming large-scale structures. Equations of Cahn–Hilliard type have been proposed to explain the qualitative features of self-assembly (20); (25), and knowledge of variations in the film height could enhance these models. Indeed in (18) the authors use the present thin-film Cahn–Hilliard model in three dimensions to control phase separation, a useful tool in applications where it is necessary for the molecules in the film to form a given structure.

The mathematical analysis of thin-film equations was given great impetus by Bernis and Friedman in (26). They focus on the basic thin-film equation,


with no-flux boundary conditions on a line segment, and smooth nonnegative initial conditions. For this equation describes a thin bridge between two masses of fluid in a Hele–Shaw cell, for it is used in slip models as  (27), while for it gives the evolution of the free surface of a thin film experiencing capillary forces (16). Using a decaying free-energy functional, they analyzed Eq. (1) and obtained results concerning the existence, uniqueness, and regularity of solutions, as a function of the exponent . Only for does a classical, smooth solution exist; this fact is established by construction of an entropy functional. This paper (26) has inspired other work on the subject (28); (29); (30), in which the effect of a Van der Waals term on Eq. (1) is investigated. These works provide results concerning regularity, long-time behaviour, and film rupture in the presence of an attractive Van der Waals force. More relevant to the present work is the paper by Wieland and Garcke (31), in which a pair of partial differential equations describes the coupled evolution of free-surface variations and surfactant concentration. The authors derive the relevant equations using the long-wavelength theory, obtain a decaying energy functional, and prove results concerning the existence and non-negativity of solutions.

When the binary fluid forms a thin film on a substrate, we shall show in Sec. II that a long-wave approximation simplifies the density- and viscosity-matched Navier–Stokes Cahn–Hilliard equations, which reduce to a pair of coupled evolution equations for the free surface and concentration. If is the scaled free-surface height, and is the binary fluid concentration, then the dimensionless equations take the form


Here is the capillary number, measures the strength of coupling between the concentration and free-surface variations (backreaction), and is the scaled interfacial thickness — sometimes called the Cahn number. Additionally, is the dimensionless, spatially-varying surface tension, and is the body-force potential acting on the film. In this paper we focus on a class of potentials that models a repulsive interaction between the film and the substrate, thus preventing rupture. This enables us to focus on late-time phase separation, which is synonymous with a tendency to equilibrium. However, rupture is in itself an important feature in thin-film equations (16); (28); (29): we therefore use numerical methods to highlight the possibility of rupture in the absence of a repulsive film-substrate interaction.

The paper is organised as follows. In Sec. II we discuss the Navier–Stokes Cahn–Hilliard equation and the scaling laws that facilitate the passage to the long-wavelength equations, and we derive Eq. (2). In Sec. III we perform linear and non-linear analyses of the long-wavelength equations. Our non-linear analysis centres on finding a Lyapunov functional for a given class of potentials. We derive a priori bounds for a given positive () solution, and estimate the minimum value of the free-surface height. In Sec. IV we outline a series of numerical studies, with and without a regularizing potential, and we discuss the dependence of the minimum free-surface height on the problem parameters. Finally, in Sec. V we present our conclusions.

Ii The model equations: derivation

In this section we introduce the two-dimensional Navier–Stokes Cahn–Hilliard (NSCH) equation set. We focus on the so-called matched case, wherein both components of the binary mixture have the same density and viscosity. We discuss the assumptions underlying the long-wavelength approximation. We enumerate the scaling rules necessary to obtain the simplified equations. Finally, we arrive at a set of equations that describe phase separation in a thin film subject to arbitrary body forces.

The full NSCH equations describe the coupled effects of phase separation and flow in a binary fluid. If the fluids are density- and viscosity matched, then the models described in the references (7); (8) agree; this is the case we study. If is the fluid velocity and is the concentration of the mixture, where indicates total segregation, then these fields evolve as




is the stress tensor, is the fluid pressure, is the body potential and is the constant density. The constant is the kinematic viscosity, , where is the dynamic viscosity. Additionally, is a constant with units of , is a constant that gives the typical width of interdomain transitions, and a diffusion coefficient with dimensions .

We impose the following boundary conditions (BCs). On the lateral boundaries (-direction), the velocity must satisfy the no-slip condition, while and must vanish too. Alternatively, we may enforce periodicity, and demand that the velocity, , , and be periodic functions of the lateral coordinate. Finally, if the system has a free surface in the vertical or -direction, then the vertical BCs are the following:

while on the free surface they are

where is the unit normal to the surface, is the unit vector tangent to the surface, is the surface coordinate, is the surface tension, and is the mean curvature,

these free-surface conditions are standard and are discussed in the review papers (16); (17). This choice of BCs guarantees the conservation of the total mass and volume,


Here represents the time-dependent domain of integration, owing to the variability of the free-surface height. Note that in view of the concentration BC (5d), the stress BC (5b) and does not contain or its derivatives.

These equations simplify considerably if the fluid forms a thin layer of mean thickness , for then the scale of lateral variations is large compared with the scale of vertical variations . Specifically, the parameter is small, and after nondimensionalisation of Eq. (3) we expand its solution in terms of this parameter, keeping only the lowest-order terms. For a review of this method and its applications, see (16); (17). For simplicity, we shall work in two dimensions, but the generalisation to three dimensions is easily effected (18).

In terms of the small parameter , the equations nondimensionalise as follows. The diffusion time scale is and we choose this to be the unit of time. Then the unit of horizontal velocity is so that , where variables in upper case denote dimensionless quantities. Similarly, the vertical velocity is , and the free-surface height is . The dimensionless coordinates are introduced through the equations , . Finally, for the equations of motion to be half-Poiseuille at (in the absence of the backreaction) we choose and . We also stipulate the following form for the surface tension:

where the exponent is yet to be determined. Using these scaling rules, the dimensionless momentum equations are




The choice of ordering for the Reynolds number allows us to recover half-Poiseuille flow at . We delay choosing the ordering of the dimensionless group until we have examined the concentration equation, which in nondimensional form is


where . By switching off the backreaction in the momentum equations (corresponding to ), we find the trivial solution to the momentum equations, , . The concentration boundary conditions are then on , which forces so that the Cahn–Hilliard equation is simply

To make the lubrication approximation consistent, we take


We now carry out a long-wavelength approximation to Eq. (11), writing , , . We examine the boundary conditions on first. They are on ; on these conditions are simply , while on the surface derivatives are determined by the relations

Thus, the BCs on are simply on , which forces . Similarly, we find , and

In the same manner, we derive the results . Using these facts, Eq. (11) becomes

We now integrate this equation from to and use the boundary conditions

After rearrangement, the concentration equation becomes


is the vertically-averaged velocity. Introducing

the thin-film Cahn–Hilliard equation becomes


We are now able to perform the long-wavelength approximation to Eqs. (7) and (8). At lowest order, Eq. (8) is , since , and hence

We introduce a dimensionless group to measure the strength of the interaction between the concentration and velocity fields; we also specify its order of magnitude:


Later on we refer to this quantity as the ‘backraction strength’, since it is a measure of the extent to which concentration gradients feed back into the flow field. Using this dimensionless group, Eq. (7) becomes

Using this becomes


At lowest order, the BC (5b) reduces to


which combined with Eq. (15) yields the relation

Here is the dimensionless, spatially-varying component of the surface tension. Making use of the BC on and integrating again, we obtain the result


The vertically-averaged velocity is therefore


where we used the standard Laplace–Young free-surface boundary condition to eliminate the pressure, and


Finally, by integrating the continuity equation in the -direction, we obtain, in a standard manner, an equation for free-surface variations,


The inclusion of the surface-tension terms requires further elucidation. In the long-wave limit, the normal-stress condition reduces to , which in dimensionless form is

By promoting the constant to , and by taking , this equation reduces to , as in Eq. (18). Similarly, the normal-stress condition reduces to , which in dimensionless form reads

Taking gives a contribution to the shear-stress balance at lowest order, , , as in Eq. (16). Going to higher exponents suppresses this contribution.

Let us assemble our results, restoring the lower-case fonts and omitting ornamentation over the constants. The height equation (20) becomes

while the concentration equation (13) becomes

and where we have the nondimensional constants


The boundary conditions are inherited from the full Navier–Stokes Cahn–Hilliard equations: Since contains a depth-averaged velocity, we write it as . Thus, the concentration , the chemical potential , and the flux are either periodic functions in the lateral direction, or satisfy the no-flux conditions on the lateral boundaries. Equations (21), together with the boundary conditions described, are the thin-film NSCH equations. The integral quantities defined in Eq. (6) are manifestly conserved, while the free surface and concentration are coupled. The term in the flux represents a driving force, which can be externally prescribed, or a function of the concentration . In either case, the inclusion of this term can have a substantial effect on the behaviour of the system. For the rest of this study, this term is set to zero; its inclusion is discussed elsewhere by the present authors (18), and by others (32).

In view of the severe constraint , some discussion about the applicability of Eqs. (21) to real systems is warranted. This constraint is the condition that the mean thickness of the film be much smaller than the transition-layer thickness. In experiments involving the smallest film thicknesses attainable ( m) (33), this condition is naturally satisfied. Furthermore, in certain situations far from this limiting case, variations in the domain structure in the vertical direction are suppressed, and a system of equations with no vertical (-) dependence, such as Eqs. (21), is appropriate. This kind of situation arises when external effects such as the air-fluid and fluid-substrate interactions do not prefer one binary fluid component or another; hence, the dimensionality of the film is reduced, and the balance laws implied by Eqs. (21) are applicable.

Iii The model equations: analysis

The choice of potential determines the behaviour of solutions. In this section, we perform a linear-stability analysis on the model equations (21) and identify the pattern-formation mechanism. We also develop results for the non-linear regime using the theory of a priori bounds. These are bounds on norms of the solution that are obtained without assuming any prior knowledge of the solution. Throughout this section, the driving force resulting from surface-tension gradients is set to zero.

The first step in our analytical study is to find the circumstances under which the constant state is unstable to a small-amplitude, initial perturbation . This perturbation evolves in time to a state , which satisfies the linearized version of equations (21). By writing down a wave ansatz , we obtain an eigenvalue equation,

with eigenvalues


Thus, there are two routes to instability. The system can become unstable as a result of substrate-film interactions if . Such an interaction will often lead to rupture (16). If this route is suppressed, then film rupture may be prevented, but the second route to instability is also relevant. This is accessible when is in the spinodal range . Thus, even when the first route to instability is not accessible, a critical mixed state will phase separate in a manner similar to the classical Cahn–Hilliard fluid, as described in Sec. I.

While the linear analysis is helpful to describe early-stage evolution, it sheds no light on the behaviour at later times. We therefore turn to the non-linear analysis of the problem (21). The non-linear analysis centres on finding bounds for a given solution . To do this, it is necessary to construct a Lyapunov functional, that is, a non-negative, non-increasing functional of the solution . It is certainly the case that we can find a non-increasing functional based on the solution pair , which is a kind of energy for the problem:

Proposition 1 (Existence of a decreasing functional)

Given a smooth solution to the equations (21), positive in the sense that , and a continuous potential function , then the functional


is non-increasing, .

The proof of this claim is readily obtained by a straightforward time-differentiation of and application of the equations (21), together with the no-flux or periodic boundary conditions. We find,


We build on this result by focussing on the following class of potential whose anti-derivative is positive:


where is an arbitrary reference height. Using Prop. 1 and the condition (26), we obtain a Lyapunov functional for the positive solution , :

Proposition 2 (Existence of a positive Lyapunov functional)

For a smooth solution to the equations (21), positive in the sense that , and for a potential function with positive anti-derivative, there is an associated Lyapunov functional.

To verify this claim, it suffices to note that since for the class of potential functions under consideration, all terms in the functional (Eq. (24)) are positive, and thus is a positive, non-increasing function of time, i.e. a Lyapunov functional.

The boundedness result provides a regularity condition on the height , although this is valid only in a single spatial dimension.

Proposition 3 (Hölder continuity of )

If is a smooth, positive solution to the equations (21), in the sense that , and if the potential function has a positive anti-derivative, then is Hölder continuous, with time-independent Hölder constant .

Proposition 3 follows from the a priori bounds


and from the use of Hölder’s inequality on the following string of relations:

where is the time-independent Hölder constant. As an immediate corollary of this result, we obtain an upper bound on the height field:

Proposition 4 (An upper bound on the height field)

If is a smooth, positive solution to the equations (21), in the sense that , and if the potential function has a positive anti-derivative, then is bounded above.

Since the free energy contains a term in the -norm of , a similar result exists for the concentration field, provided everywhere. Loss of control over the minimum value of therefore implies loss of control over the concentration gradient. This suggests that blowup of gradients and film rupture are related, a claim which we demonstrate numerically in Sec. IV.2. Such extreme events are avoided when a repulsive film-substrate interaction is present, in which case a positive lower bound on exists; it is to that result that we now turn.

Given the form of Eqs. (21), regularity of a given solution is guaranteed only when a lower bound on is obtained, in addition to the upper bounds just provided. To derive such a result, we first specialize to the potential


from which a more general result will follow.

Proposition 5 (No-rupture condition for the potential in Eq. (28))

If is a smooth, positive solution to the equations (21), in the sense that , and if the potential function has the form given by Eq. (28), then there is an a priori, time-independent lower bound on .

Note first of all that the potential (28) has a positive anti-derivative, , where the reference height in Eq. (26) is set to . Thus, there is a Lyapunov functional for the solution , and hence,


Using the Hölder continuity of ,


where is the Hölder constant for . Using results (29) and (30),


By integrating this equation, we arrive at the relation


Let us examine the properties of the function . It tends to infinity as , and tends to zero as . It is also monotone-decreasing over . Thus, the equation has precisely one positive root for , which we call . To satisfy the inequality (32), it must be the case that

For the potential (Eq. (28)), the root can be obtained explicitly:

This completes the proof of the no-rupture condition for the potential (28).

To arrive at a no-rupture condition for a general potential, we introduce the function


Based on Prop. 5, we write down sufficient conditions on and for the existence of a positive lower bound on :

Proposition 6 (A sufficient condition to avoid rupture)

If is a smooth, positive solution to the equations (21), in the sense that , if the anti-derivative of the potential function is a positive, non-increasing function, and moreover, if satisfies the conditions


then a positive a priori lower bound on exists, independent of time.

The proof of Prop. 6 is in the same spirit as that of Prop. 5. Given the positivity of the anti-derivative , a Lyapunov exponent exists (Prop. 1), and thus we have the bound . Using the Hölder continuity of (Prop. 3), and the condition that the should be a non-increasing function,

Using the bound ,

Given the conditions (34), there exists at least one solution to the equation , for . Indeed, since is non-increasing, so too is , and thus this equation has precisely one solution, which we call . Then, for the condition to be satisfied on the interval , we must take

Finally, using this theory, we investigate the potential


and calculate the -values for which a no-rupture condition can be found.

Proposition 7 (Conditions on the potential (35) to avoid rupture)

If is a smooth, positive solution to the equations (21), in the sense that , and if the potential function is given by Eq. (35), then a no-rupture condition is guaranteed to hold for .

The proof of Prop. 7 follows by a straightforward evaluation of the integral (33). The case is covered by Prop. 5. Thus, we focus on the case , where the integral has the value


and where we have set . For the conditions (34) hold: , ; note also that is non-increasing. Thus, the equation has exactly one positive root , and this serves as a lower bound on , . Note

Figure 1: A plot of the integral as a function of . For the integral diverges as and tends to zero as , allowing for the existence of a positive root for the equation .

that the relation given in Eq. (36) fails to satisfy conditions (34) when . Thus, a sufficient condition for the solution to possess a no-rupture condition is for the potential to have the form given in Eq. (35), with . This analysis is also described schematically in Fig. 1, where a plot of the integral is shown as a function of . For the integral diverges as and tends to zero as , which provides for the existence of a positive root for the equation .

Two questions arise from our results. The first involves the form of the potential used in the derivation of a no-rupture condition: can we replace the repulsive power-law form with a more general function and still obtain a no-rupture condition? The answer to this question comes readily through the observation that our construction of a positive lower bound for relies only on the surface-tension and body-force components of the free energy, namely

Thus, Propositions 57 can be viewed in the context of the PDE theory for a single variable (the free-surface height), where a no-rupture condition exists (34), both for the power-law type potential considered here and for the Lennard–Jones potential , where and are positive constants and  (34). Hence, our regularity results are generalizable to this wider class of potential.

Having weakened the sufficient condition for the avoidance of rupture, it is reasonable to ask, is the presence of a suitable potential even a necessary condition? This question is motivated by the single-variable theory for the free-surface height, where under certain conditions an entropy functional facilitates the construction of an a priori lower bound on the height  (26). The entropy is obtained through the following steps, which we describe for the single-variable case :

  1. Identify the power of that is a factor in the flux ; this is the mobility, . For the single-variable case, .

  2. Obtain the function .

  3. The entropy is then defined as ; for the single-variable case, this is (we have omitted the unimportant constant of proportionality).

For the single-variable case,


and the entropy is a non-increasing functional of , . Setting in Eq. (37) gives the relation


For , Hölder continuity combined with the bound in Eq. (38) enables the construction of a pointwise lower bound on . When (the case considered in this paper), no such pointwise bound exists; then, estimates for the entropy based on the Hölder continuity of are non-singular in the minimum height. However, a more definitive obstacle than this exists when we seek to construct an entropy functional for the system (21), namely, that the entropy functional implied by Eqs. (21) fails to be a non-increasing function of time:

Proposition 8

(The time-derivative of the entropy associated with Eqs. (21) is not sign definite) Given a smooth, positive solution to the equations (21), in the sense that