Nonlinear dynamics of phase separation in thin films
Abstract
We present a longwavelength approximation to the Navier–Stokes Cahn–Hilliard equations to describe phase separation in thin films. The equations we derive underscore the coupled behaviour of freesurface variations and phase separation. We introduce a repulsive substratefilm interaction potential and analyse the resulting fourthorder equations by constructing a Lyapunov functional, which, combined with the regularizing repulsive potential, gives rise to a positive lower bound for the freesurface 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 freesurface 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 wellmixed 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 freesurface 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 selfsegregating 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 variablemobility 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 coworkers (7) provide a derivation of coupled Navier–Stokes Cahn–Hilliard (NSCH) equations in which the velocity advects the phaseseparating 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 longwavelength approximation is possible, in which the full equations with a moving boundary at the free surface are reduced to a single equation for the freesurface height. In the present case, the reduction yields two equations: one for the free surface, and one for the Cahn–Hilliard concentration. The resulting thinfilm 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 twodimensional case: we derive the thinfilm 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 thinfilm 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 thinfilm Cahn–Hilliard theory is in selfassembly (23); (22); (24); (20); (21). Here molecules (usually residing in a thin layer) respond to an energyminimisation requirement by spontaneously forming largescale structures. Equations of Cahn–Hilliard type have been proposed to explain the qualitative features of selfassembly (20); (25), and knowledge of variations in the film height could enhance these models. Indeed in (18) the authors use the present thinfilm 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 thinfilm equations was given great impetus by Bernis and Friedman in (26). They focus on the basic thinfilm equation,
(1) 
with noflux 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 freeenergy 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, longtime 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 freesurface variations and surfactant concentration. The authors derive the relevant equations using the longwavelength theory, obtain a decaying energy functional, and prove results concerning the existence and nonnegativity of solutions.
When the binary fluid forms a thin film on a substrate, we shall show in Sec. II that a longwave approximation simplifies the density and viscositymatched Navier–Stokes Cahn–Hilliard equations, which reduce to a pair of coupled evolution equations for the free surface and concentration. If is the scaled freesurface height, and is the binary fluid concentration, then the dimensionless equations take the form
(2a)  
where  
(2b)  
(2c) 
Here is the capillary number, measures the strength of coupling between the concentration and freesurface variations (backreaction), and is the scaled interfacial thickness — sometimes called the Cahn number. Additionally, is the dimensionless, spatiallyvarying surface tension, and is the bodyforce 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 latetime phase separation, which is synonymous with a tendency to equilibrium. However, rupture is in itself an important feature in thinfilm equations (16); (28); (29): we therefore use numerical methods to highlight the possibility of rupture in the absence of a repulsive filmsubstrate 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 longwavelength equations, and we derive Eq. (2). In Sec. III we perform linear and nonlinear analyses of the longwavelength equations. Our nonlinear 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 freesurface 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 freesurface height on the problem parameters. Finally, in Sec. V we present our conclusions.
Ii The model equations: derivation
In this section we introduce the twodimensional Navier–Stokes Cahn–Hilliard (NSCH) equation set. We focus on the socalled matched case, wherein both components of the binary mixture have the same density and viscosity. We discuss the assumptions underlying the longwavelength 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
(3a)  
(3b)  
(3c) 
where
(4) 
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 noslip 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:
(5a)  
while on the free surface they are  
(5b)  
(5c)  
(5d) 
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 freesurface 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,
(6) 
Here represents the timedependent domain of integration, owing to the variability of the freesurface 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 lowestorder 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 freesurface height is . The dimensionless coordinates are introduced through the equations , . Finally, for the equations of motion to be halfPoiseuille 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
(7) 
(8) 
(9) 
where
(10) 
The choice of ordering for the Reynolds number allows us to recover halfPoiseuille flow at . We delay choosing the ordering of the dimensionless group until we have examined the concentration equation, which in nondimensional form is
(11) 
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
(12) 
We now carry out a longwavelength 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
where
is the verticallyaveraged velocity. Introducing
the thinfilm Cahn–Hilliard equation becomes
(13) 
We are now able to perform the longwavelength 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:
(14) 
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
(15) 
At lowest order, the BC (5b) reduces to
(16) 
which combined with Eq. (15) yields the relation
Here is the dimensionless, spatiallyvarying component of the surface tension. Making use of the BC on and integrating again, we obtain the result
(17) 
The verticallyaveraged velocity is therefore
(18) 
where we used the standard Laplace–Young freesurface boundary condition to eliminate the pressure, and
(19) 
Finally, by integrating the continuity equation in the direction, we obtain, in a standard manner, an equation for freesurface variations,
(20) 
The inclusion of the surfacetension terms requires further elucidation. In the longwave limit, the normalstress 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 normalstress condition reduces to , which in dimensionless form reads
Taking gives a contribution to the shearstress balance at lowest order, , , as in Eq. (16). Going to higher exponents suppresses this contribution.
Let us assemble our results, restoring the lowercase fonts and omitting ornamentation over the constants. The height equation (20) becomes
(21a)  
while the concentration equation (13) becomes  
(21b)  
where  
(21c)  
and  
(21d) 
and where we have the nondimensional constants
(22) 
The boundary conditions are inherited from the full Navier–Stokes Cahn–Hilliard equations: Since contains a depthaveraged 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 noflux conditions on the lateral boundaries. Equations (21), together with the boundary conditions described, are the thinfilm 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 transitionlayer 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 airfluid and fluidsubstrate 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 linearstability analysis on the model equations (21) and identify the patternformation mechanism. We also develop results for the nonlinear 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 surfacetension 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 smallamplitude, 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
(23) 
Thus, there are two routes to instability. The system can become unstable as a result of substratefilm 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 earlystage evolution, it sheds no light on the behaviour at later times. We therefore turn to the nonlinear analysis of the problem (21). The nonlinear analysis centres on finding bounds for a given solution . To do this, it is necessary to construct a Lyapunov functional, that is, a nonnegative, nonincreasing functional of the solution . It is certainly the case that we can find a nonincreasing 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
(24) 
is nonincreasing, .
The proof of this claim is readily obtained by a straightforward timedifferentiation of and application of the equations (21), together with the noflux or periodic boundary conditions. We find,
(25) 
We build on this result by focussing on the following class of potential whose antiderivative is positive:
(26) 
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 antiderivative, 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, nonincreasing 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 antiderivative, then is Hölder continuous, with timeindependent Hölder constant .
Proposition 3 follows from the a priori bounds
(27) 
and from the use of Hölder’s inequality on the following string of relations:
where is the timeindependent 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 antiderivative, 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 filmsubstrate 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
(28) 
from which a more general result will follow.
Proposition 5 (Norupture condition for the potential in Eq. (28))
Note first of all that the potential (28) has a positive antiderivative, , where the reference height in Eq. (26) is set to . Thus, there is a Lyapunov functional for the solution , and hence,
(29) 
Using the Hölder continuity of ,
(30) 
where is the Hölder constant for . Using results (29) and (30),
(31) 
By integrating this equation, we arrive at the relation
(32) 
Let us examine the properties of the function . It tends to infinity as , and tends to zero as . It is also monotonedecreasing 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 norupture condition for the potential (28).
To arrive at a norupture condition for a general potential, we introduce the function
(33) 
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 antiderivative of the potential function is a positive, nonincreasing function, and moreover, if satisfies the conditions
(34) 
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 antiderivative , 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 nonincreasing function,
Using the bound ,
Given the conditions (34), there exists at least one solution to the equation , for . Indeed, since is nonincreasing, 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
(35) 
and calculate the values for which a norupture condition can be found.
Proposition 7 (Conditions on the potential (35) to avoid rupture)
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
(36) 
and where we have set . For the conditions (34) hold: , ; note also that is nonincreasing. Thus, the equation has exactly one positive root , and this serves as a lower bound on , . Note
that the relation given in Eq. (36) fails to satisfy conditions (34) when . Thus, a sufficient condition for the solution to possess a norupture 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 norupture condition: can we replace the repulsive powerlaw form with a more general function and still obtain a norupture condition? The answer to this question comes readily through the observation that our construction of a positive lower bound for relies only on the surfacetension and bodyforce components of the free energy, namely
Thus, Propositions 5–7 can be viewed in the context of the PDE theory for a single variable (the freesurface height), where a norupture condition exists (34), both for the powerlaw 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 singlevariable theory for the freesurface 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 singlevariable case :

Identify the power of that is a factor in the flux ; this is the mobility, . For the singlevariable case, .

Obtain the function .

The entropy is then defined as ; for the singlevariable case, this is (we have omitted the unimportant constant of proportionality).
For the singlevariable case,
(37) 
and the entropy is a nonincreasing functional of , . Setting in Eq. (37) gives the relation
(38) 
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 nonsingular 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 nonincreasing function of time: