Patterns formed in a thin film with spatially homogeneous and non-homogeneous Derjaguin disjoining pressure
aa r X i v : . [ n li n . PS ] J a n Pattern formation in a thin-film equation withspatially homogeneous and non-homogeneousDerjaguin disjoining pressure
A. S. ALSHAIKHI , M. GRINFELD , and S. K. WILSON Department of Mathematics and Statistics, University of Strathclyde,Livingstone Tower, 26 Richmond Street, Glasgow G1 1XH, UK email : [email protected] Abstract
We consider pattern formation in a two-dimensional thin film on a planar substratewith a Derjaguin disjoining pressure and periodic wettability stripes in one of the direc-tions. We rigorously clarify some of the results obtained numerically by Honisch et al. [ Langmuir
31: 10618–10631, 2015] and embed them in the general theory of thin-filmequations. For the case of constant wettability, we elucidate the change in the globalstructure of branches of stationary solutions as the average film thickness and the sur-face tension are varied. Specifically we find, by using methods of local bifurcation theoryand the continuation software package AUTO, both nucleation and metastable regimes.We discuss admissible forms of spatially non-homogeneous disjoining pressure, arguingfor a form that differs from the one used by Honisch et al. , and study the dependenceof the stationary solutions on the wettability contrast in that case.
Keywords: thin films, disjoining pressure, non-homogeneous substrates, pattern formation
Thin liquid films on solid substrates occur in many natural situations. For example, theyappear in tear films in the eye which protect the cornea [5] or in streams of lava from avolcanic eruption [15]. Moreover, thin liquid films occur in many technological applications,such as coatings [18] and lubricants (e.g. oil films which lubricate the piston in a car engine[36], drying paint layers [14], and in the manufacture of microelectronic devices [17]). Forextensive reviews of thin-film flow see, for example, Oron et al. [25] and Craster and Matar[8].As these liquid films are thin, the Navier–Stokes equation governing their flow can be re-duced to a single degenerate fourth-order quasi-linear parabolic partial differential equationusually known as a thin-film equation. In many applications a choice of a disjoining pres-sure, which we denote by Π, must be made. Such a term describes the action of surfaceforces on the film [32]. In different situations, different forms of disjoining pressure areappropriate; these may incorporate long-range van der Waals forces and/or various types1f short-range interaction terms such as Born repulsion; inclusion of a particular type ofinteraction can have significant effects on the wettability of the surface and the evolutionof the film film, sometimes leading to dewetting phenomena, i.e. the rupture of the filmand the appearance of dry spots. (Here and subsequently by “wettability” of the surfacewe mean the ability of a solid surface to reduce the surface tension of a liquid on contactwith it such that it spreads over the surface and wets it.)Witelski and Bernoff [38] were one of the first authors to analyse mathematically the ruptureof three-dimensional thin films. In particular, considering a disjoining pressure of the formΠ = − / (3 h ) (we use the sign convention adopted in Honisch et al. [13]), they analysedplanar and axisymmetric equilibrium solutions on a finite domain. They showed that adisjoining pressure of this form leads to finite-time rupture singularities, that is, the filmthickness approaches zero in finite time at a point (or a line or a ring) in the domain. Ina related more recent paper, Ji and Witelski [16] considered a different choice of disjoiningpressure and investigated the finite-time rupture solutions in a model of thin film of liquidwith evaporative effects. They observed different types of finite-time singularities due tothe non-conservative terms in the model. In particular, they showed that the inclusion ofnon-conservative term can prevent the disjoining pressure from causing finite-time rupture.A pioneering theoretical study of a thin-film equation with a disjoining pressure term givenby a combination of negative powers of the height of the thin film is that by Bertozzi et al. [3], who studied the formation of dewetting patterns and the rupture of thin liquid filmsdue to long-range attractive and short-range Born repulsive forces, and determined thestructure of the bifurcation diagram for stationary solutions, both with and without therepulsive term.Aiming to quantify the temporal coarsening in a thin film, Glasner and Witelski [11] exam-ined two coarsening mechanisms that arise in dewetting films: mass exchange that influencesthe breakdown of individual droplets and spatial motion that results in droplet rupture aswell as merging events. They provided a simple model with a disjoining pressure whichcombines the effects of both short- and long-range forces acting on the film. Kitavtsev etal. [19] analysed the long-time dynamics of dewetting in a thin-film equation by using adisjoining pressure similar to that used by Bertozzi et al. [3]. They applied centre manifoldtheory to derive and analysed an ordinary differentual equation model for the dynamics ofdewetting.The recent article by Witelski [37] presents a review of the various stages of dewettingfor a film of liquid spreading on a hydrophobic substrate. Different types of behaviour ofthe film are observed depending on the form of the disjoining pressure: finite-time sin-gularities, self-similar solutions and coarsening. In particular, he divides the evolution ofdewetting processes into three phases: an initial linear instability that leads to finite-timerupture (short time dynamics), which is followed by the propagation of dewetting rims andinstabilities of liquid ridges (intermediate time dynamics), and the eventual formation ofquasi-steady droplets (long time dynamics).Most of the previous studies of thin liquid films focussed on films on homogeneous substrates.However, thin liquid films on non-homogeneous chemically patterned substrates are also ofinterest. These have many practical applications, such as in the construction of microfluidicdevices and creating soft materials with a particular pattern [27]. Chemically patternedsubstrates are an efficient way to obtain microstructures of different shapes by using different2ypes of substrate patterning [30]. Chemical modification of substrates can also be used toavoid spontaneous breakup of thin films, which is often highly undesirable, as, for example,in printing technology [4, 24].Due to their many applications briefly described above, films on non-homogeneous sub-strates have been the object of a number of previous theoretical studies which motivatethe present work. To mention some of them, Konnur et al. [20] found that in the caseof an isolated circular patch with wetting properties different from that of the rest of thesubstrate that chemical non-homogeneity of the substrate can greatly accelerate the growthof surface instabilities. Sharma et al. [31] studied instabilities of a liquid film on a substratecontaining a single heterogeneous patch and a substrate with stripes of alternating less andmore wettable regions. The main concern of that paper was to investigate how substratepatterns are reproduced in the liquid film, and to determine conditions for best templating.Thiele et al. [34] performed a bifurcation analysis using the continuation software packageAUTO [9] to study dewetting on a chemically patterned substrate by solving a thin-filmequation with a disjoining pressure, using the wettability contrast as a control parameter.The wettability contrast measures the degree of heterogeneity of the substrate; it is intro-duced and defined rigorously in (5.1) in Section 5. Honisch et al. [13] modelled an experimentin which organic molecules were deposited on a chemically non-homogeneous silicon oxidesubstrates with gold stripes and discuss the redistribution of the liquid following deposition.In a recent preprint, Liu and Witelski [23] studied thin films on chemically heterogeneoussubstrates. They claim that in some applications such as digital microfluidics, substrateswith alternate hydrophilic and hydrophobic rectangular areas are better described by apiecewise constant function than by a sinusoidal one. Therefore, in contrast to other studies,including the present one, they study substrates with wettability characteristics describedby such a function. Based on the structure of the bifurcation diagram, they divide thesteady-state solutions into six distinct but connected branches and show that the onlyunstable branch corresponds to confined droplets, while the rest of the branches are stable.In the present work, we build on the work of Thiele et al. [34] and Honisch et al. [13]. Partof our motivation is to clarify and explain rigorously some of the numerical results reportedin these papers. In the sinusoidally striped non-homogeneous substrate case, we offer ajustification for using a form of the disjoining pressure that differs from the one used inthese two papers. A detailed plan of the paper is given in the last paragraph of Section 2. Denoting the height of the thin liquid film by z = h ( x, y, t ), where ( x, y, z ) are the usualCartesian coordinates and t is time, Honisch et al. [13] consider the thin-film equation h t = ∇ · { Q ( h ) ∇ P ( h, x, y ) } , t > , ( x, y ) ∈ R , (2.1)where Q ( h ) = h / (3 η ) is the mobility coefficient with η being the dynamic viscosity. Thegeneralized pressure P ( h, x, y ) is given by P ( h, x, y ) = − γ ∆ h − Π( h, x, y ) , γ is the surface tension and we follow [13] in taking the Derjaguin disjoining pressureΠ( h, x, y ) in the spatially homogeneous case to be of the formΠ( h, x, y ) = − Ah + Bh (2.2)suggested, for example, by Pismen [26]. Here A and B are positive parameters that measurethe relative contributions of the short-range forces (the h − term) and the long-range ones(the h − term). However, we will see that both of these constants can be scaled out of themathematical problem.In what follows, we study thin films on both homogeneous and non-homogeneous substrates.In the non-homogeneous case, we will modify (2.2) by assuming that the Derjaguin pressureterm Π changes periodically in the x -direction with period L . The appropriate forms of Πin the non-homogeneous case are discussed in Section 5.Hence, in order better to understand solutions of (2.1), we start by characterising y -independent stationary solutions of that equation for 0 < x < L subject to periodic bound-ary conditions at x = 0 and x = L . In other words, we seek solutions of (2.1) in the form h ( x, y, t ) = h ( x ), satisfying in dimensional variables the non-local boundary value problem γh xx + Bh − Ah − L Z L (cid:20) Bh − Ah (cid:21) d x = 0 , < x < L, (2.3)subject to the constraint 1 L Z L h ( x ) d x = h ∗ , (2.4)where the constant h ∗ ( >
0) denotes the (scaled) average film thickness, and the periodicboundary conditions h (0) = h ( L ) , h x (0) = h x ( L ) . (2.5)Now we non-dimensionalise. Setting H = (cid:18) BA (cid:19) / , h = H ˜ h, and x = L ˜ x, in (2.3) and removing the tildes, we obtain ǫ h xx + f ( h ) − Z f ( h ) d x = 0 , < x < , (2.6)where f ( h ) = 1 h − h , (2.7)and ǫ = γB / L A / , (2.8)subject to the periodic boundary conditions h (0) = h (1) , h x (0) = h x (1) , (2.9)4nd the volume constraint Z h ( x ) d x = ¯ h := h ∗ A / B / . (2.10)Note that the problem (2.6)–(2.10) is very similar to the corresponding stationary problemfor the Cahn–Hilliard equation considered as a bifurcation problem in the parameters ¯ h and ǫ by Eilbeck et al. [10]. The boundary conditions considered in that work were thephysically natural double Neumann conditions. The periodic boundary conditions (2.9) inthe present problem slightly change the analysis, but our general approach in characterisingdifferent bifurcation regimes still follows that of Eilbeck et al. [10], though the correctinterpretation of the limit as ǫ → + is that now we let the surface tension γ go to zero.In particular, we perform a Liapunov–Schmidt reduction to determine the local behaviourclose to bifurcation points and then use AUTO (in the present work we use the AUTO-07pversion [9]) to explore the global structure of branches of stationary solutions both for thespatially homogeneous case and for the spatially non-homogeneous case in the case of an x -periodically patterned substrate.We first investigate the homogeneous case and, having elucidated the structure of the bifur-cations of non-trivial solutions from the constant solution h = ¯ h in that case in Sections 3and 4, we study forced rotational ( O (2)) symmetry breaking in the non-homogeneous casein Section 5. In Appendix A, we present a general result about O (2) symmetry breaking inthe spatially non-homogeneous case. It shows that in the spatially non-homogeneous case,only two stationary solutions remain from the orbit of solutions of (2.6)–(2.10) induced byits O (2) invariance. We concentrate on the simplest stationary solutions of (2.6)–(2.10), asby a result of Laugesen and Pugh [21, Theorem 1] only such solutions, that is, constantsolutions and those having only one extremum point, are linearly stable in the homogeneouscase.Below we use k · k to denote L ([0 , We start by performing an analysis of the dependence of the global structure of branches ofstationary solutions of the problem in the spatially homogeneous case given by (2.6)–(2.10)on the parameters ¯ h and ǫ .In what follows, we do not indicate explicitly the dependence of the operators on the pa-rameters ¯ h and ǫ , and all of the calculations are performed for a fixed value of ¯ h and closeto a bifurcation point ǫ = ǫ k for k = 1 , , , . . . defined below.We set v = h − ¯ h , so that v = v ( x ) has zero mean, and rewrite (2.6) as G ( v ) = 0 , where G ( v ) = ǫ v xx + f ( v + ¯ h ) − Z f ( v ( x ) + ¯ h ) d x.
5f we set H = (cid:26) w ∈ C (0 ,
1) : Z w ( x ) d x = 0 (cid:27) , where G is an operator from D ( G ) ⊂ H → H , then D ( G ) is given by D ( G ) = (cid:26) v ∈ C (0 ,
1) : v (0) = v (1) , v x (0) = v x (1) , Z v ( x ) d x = 0 (cid:27) . The linearisation of G at v applied to w is defined by dG ( v ) w = lim τ → G ( v + τ w ) − G ( v ) τ . We denote dG (0) by S , so that S applied to w is given by Sw = ǫ w xx + f ′ (¯ h ) w. (3.1)To locate bifurcation points, we have to find the nontrivial solutions of the equation Sw = 0subject to w (0) = w (1) , w x (0) = w x (1) . (3.2)The kernel of S is non-empty and two dimensional when ǫ = ǫ k = p f ′ (¯ h )2 kπ for k = 1 , , , . . . , (3.3)and is spanned by cos(2 kπx ) and sin(2 kπx ). That these values of ǫ correspond to bifurcationpoints follows from two theorems of Vanderbauwhede [35, Theorems 2 and 3, pp. 361–363].In a neighbourhood of a bifurcation point ( ǫ k ,
0) in ( ǫ, v ) space, solutions of G ( v ) = 0 on H are in one-to-one correspondence with solutions of the reduced system of equations on R , g ( x, y, ǫ ) = 0 , g ( x, y, ǫ ) = 0 , (3.4)for some functions g and g to be obtained through the Liapunov–Schmidt reduction [12].To set up the Liapunov–Schmidt reduction, we decompose D ( G ) and H as follows: D ( G ) = ker S ⊕ M and H = N ⊕ range S. Since S is self-adjoint with respect to the L -inner product denoted by h· , ·i , we can choose M = N = span { cos(2 kx ) , sin(2 kx ) } , and denote the above basis for M by { w , w } and for N by { w ∗ , w ∗ } . We also denote theprojection of H onto range S by E .Since the present problem is invariant with respect to the group O (2), the functions g and g must have the form g ( x, y, ǫ ) = xp ( x + y , ǫ ) , g ( x, y, ǫ ) = yp ( x + y , ǫ ) , (3.5)6or some function p ( · , · ) [7], which means that in order to determine the bifurcation structure,the only terms that need to be computed are g ,xǫ and g ,xxx , as these immediately give g ,yǫ and g ,yyy and all of the other second and third partial derivatives of g and g areidentically zero.Following Golubitsky and Schaeffer [12], we have g ,xǫ = h w ∗ , dG ǫ ( w ) − d G ( w , S − EG ǫ (0)) i , (3.6) g ,xxx = h w ∗ , d G ( w , w , w ) − d G ( w , S − E [ d G ( w , w )]) i , (3.7)where d r G ( z , . . . , z r ) = ∂ r ∂t . . . ∂t r G ( t z + . . . + t r z r ) (cid:12)(cid:12)(cid:12)(cid:12) t = ... = t r =0 for r = 1 , , , . . . , (3.8)and we choose w = w ∗ = cos(2 kπx ) , where w ∈ ker S and w ∗ ∈ (range S ) ⊥ . In particular, from (3.8) we have d G ( z , z ) = ∂ ∂t ∂t G ( t z + t z ) (cid:12)(cid:12)(cid:12)(cid:12) t = t =0 = ∂ ∂t ∂t h ǫ k ( t z ,xx + t z ,xx ) + f ( t z + t z + ¯ h ) − Z f ( t z + t z + ¯ h ) d x i(cid:12)(cid:12)(cid:12)(cid:12) t = t =0 = f ′′ (¯ h ) z z − Z f ′′ (¯ h ) z z d x, and so d G (cos(2 kπx ) , cos(2 kπx )) = f ′′ (¯ h ) cos (2 kπx ) − Z f ′′ (¯ h ) cos (2 kπx ) d x = f ′′ (¯ h ) cos (2 kπx ) − f ′′ (¯ h ) . To obtain S − E [ d G ( w , w )], which we denote by R ( x ), so that SR = E [ d G ( w , w )],we use the definition of ǫ k given in (3.3) and solve the second order ordinary differentialequation satisfied by R ( x ), R xx + 4 k π R = 4 k π f ′′ (¯ h ) f ′ (¯ h ) cos (2 kπx ) − k π f ′′ (¯ h ) f ′ (¯ h ) , subject to R (0) = R (1) and R x (0) = R x (1)in ker S . We obtain R ( x ) = S − E [ d G ( w , w )] = − f ′′ (¯ h ) f ′ (¯ h ) cos(4 kπx ) , d G ( w , S − E [ d G ( w , w )]) = d G (cid:18) cos(2 kπx ) , − f ′′ (¯ h ) f ′ (¯ h ) cos(4 kπx ) (cid:19) = f ′′ (¯ h ) cos(2 kπx ) (cid:18) − f ′′ (¯ h ) f ′ (¯ h ) cos(4 kπx ) (cid:19) − Z f ′′ (¯ h ) cos(2 kπx ) (cid:18) − f ′′ (¯ h ) f ′ (¯ h ) cos(4 kπx ) (cid:19) d x = −
16 [ f ′′ (¯ h )] f ′ (¯ h ) cos(2 kπx ) cos(4 kπx ) . (3.9)In addition, from (3.8) we have d G ( z , z , z ) = ∂ ∂t ∂t ∂t G ( t z + t z + t z ) (cid:12)(cid:12)(cid:12)(cid:12) t = t = t =0 = f ′′′ (¯ h ) z z z − Z f ′′′ (¯ h ) z z z d x, and therefore d G (cos(2 kπx ) , cos(2 kπx ) , cos(2 kπx )) = f ′′′ (¯ h ) cos (2 kπx ) − Z f ′′′ (¯ h ) cos (2 kπx ) d x = f ′′′ (¯ h ) cos (2 kπx ) . (3.10)Putting all of the information in (3.9) and (3.10) into (3.7) we obtain g ,xxx = h w ∗ , d G ( w , w , w ) − d G ( w , S − E [ d G ( w , w )]) i = Z cos(2 kπx ) (cid:20) f ′′′ (¯ h ) cos (2 kπx ) − (cid:18) −
16 [ f ′′ (¯ h )] f ′ (¯ h ) cos(2 kπx ) cos(4 kπx ) (cid:19)(cid:21) d x = 38 f ′′′ (¯ h ) + 18 [ f ′′ (¯ h )] f ′ (¯ h ) . (3.11)In addition, G ǫ ( v ) = v xx , so that G ǫ (0) = 0 at v = 0, and hence we have d G ( w k , S − EG ǫ (0)) = 0 . Furthermore, since dG ǫ ( w ) = w xx , from (3.6) we obtain g ,xǫ = h w ∗ , dG ǫ ( w ) − d G ( w , S − EG ǫ (0)) i = Z cos(2 kπx ) (cid:0) − π k cos(2 kπx ) (cid:1) d x = − k π . (3.12)Referring to (3.5) and the argument following that equation, the above analysis showsthat as long as f ′ (¯ h ) > ǫ = ǫ k a circle of equilibria bifurcates from the constantsolution h ≡ ¯ h . The direction of bifurcation is locally determined by the sign of g ,xxx given by (3.11). Hence, using 1 /ǫ as the bifurcation parameter, the bifurcation of nontrivialequilibria is supercritical if g ,xxx is negative and subcritical if it is positive.By finding the values of ¯ h where g ,xxx given by (3.11) with f ( h ) given by (2.7) is zero, wefinally obtain the following proposition: 8 roposition 1. Bifurcations of nontrivial solutions from the constant solution h = ¯ h ofthe problem (2 . – (2 . are supercritical if . < ¯ h < . and subcritical if . < ¯ h < . or if ¯ h > . . Proof
The constant solution h ≡ ¯ h will lose stability as ǫ is decreased only if f ′ (¯ h ) > − / ¯ h + 3 / ¯ h >
0, for ¯ h > / ≈ . g ,xxx = 57¯ h − h + 6512¯ h (¯ h − , so that g ,xxx < h ∈ (1 . , h / there are no bifurcations from the constant solution h = ¯ h . Furthermore, wehave the following proposition: Proposition 2.
The problem (2 . – (2 . has no nontrivial solutions when ¯ h ≤ . Proof
Assume that such a nontrivial solution exists. Then, since ¯ h ≤
1, its global mini-mum, achieved at some point x ∈ (0 , x to be an interior point by translation invariance.) But then ǫ h xx ( x ) = Z f ( h ) d x − f ( h ( x )) < , so the point x cannot be a minimum. To describe the change in the global structure of branches of stationary solutions of theproblem (2.6)–(2.10) as ¯ h and ǫ are varied, we use AUTO [9] and our results are summarisedin Figure 1.As Figure 1 shows, a curve of saddle-node (SN) bifurcations which originates from ¯ h ≈ . /ǫ ≈ .
432 satisfies ¯ h → + as 1 /ǫ → ∞ , while a curve of SN bifurcations whichoriginates from ¯ h ≈ . /ǫ ≈ .
998 satisfies ¯ h → ∞ as 1 /ǫ → ∞ .Figure 1 identifies three different bifurcation regimes, denoted by I, II and III, with differingbifurcation behaviour occurring in the different regimes, namely (using the terminology of[10] in the context of the Cahn-Hilliard equation): • a “nucleation” regime for 1 < ¯ h < / ≈ .
259 (Regime I), • a “metastable” regime for 2 / < ¯ h < .
289 and ¯ h > .
747 (Regime II), and • an “unstable” regime for 1 . < ¯ h < .
747 (Regime III).9 h . . . /ǫ . . . Figure 1: The global structure of branches of stationary solutions with a single maximum,including both saddle-node (SN) (shown with dash-dotted curves) and pitchfork (PF) bi-furcation branches (shown with solid curves). The nucleation regime 1 < ¯ h < / ≈ . / < ¯ h < .
289 and ¯ h > .
747 (Regime II), and theunstable regime 1 . < ¯ h < .
747 (Regime III) are also indicated.In regime I, the constant solution h ( x ) ≡ ¯ h is linearly stable which follows from analysingthe spectrum of the operator S for f ′ (¯ h ) < ǫ is decreased the constant solution h ( x ) ≡ ¯ h loses stability through asubcritical bifurcation.In the regime III, as ǫ is decreased the constant solution h ( x ) ≡ ¯ h loses stability through asupercritical bifurcation. Honisch et al. [13] chose the Derjaguin disjoining pressure Π( h, x, y ) to be of the formΠ( h, x, y ) = (cid:18) h − h (cid:19) (1 + ρ G ( x, y )) , (5.1)where the function G ( x, y ) models the non-homogeneity of the substrate and the parameter ρ , which can be either positive or negative, is called the “wettability contrast”. FollowingHonisch et al. [13], in the remainder of the present work we consider the specific case G ( x, y ) = sin (2 πx ) := G ( x ) , (5.2)corresponding to an x -periodically patterned (i.e. striped) substrate.10here are, however, some difficulties in accepting (5.1) as a physically realistic form ofthe disjoining pressure for a non-homogeneous substrate. The problems arise because thetwo terms in (5.1) represent rather different physical effects. Specifically, since the 1 /h term models the short-range interaction amongst the molecules of the liquid and the 1 /h term models the long-range interaction, assuming that both terms reflect the patterningin the substrate in exactly the same way through their dependence on the same function G ( x, y ) does not seem very plausible. Moreover, there are other studies which assumethat the wettability of the substrate is incorporated in either the short-range interactionterm (see, for example, Thiele and Knobloch [33] and Bertozzi et al. [3]) or the long-rangeinteraction term (see, for example, Ajaev et al. [1] and Sharma et al. [31]), but not bothsimultaneously. Hence in what follows we will consider the two cases Π( h, x ) = Π LR ( h, x )and Π( h, x ) = Π SR ( h, x ), where LR stands for “long range” and LR stands for “short range”where Π LR ( h, x ) = 1 h − h (1 + ρG ( x )) (5.3)and Π SR ( h, x ) = 1 h (1 + ρG ( x )) − h , (5.4)in both of which G ( x ) is given by (5.2) and ρ is the wettability contrast.For small wettability contrast, | ρ | ≪
1, we do not expect there to be significant differencesbetween the influence of Π LR or Π SR of the bifurcation diagrams, as these results dependonly on the nature of the bifurcation in the homogeneous case ρ = 0 and on the symmetrygroups under which the equations are invariant. To see this, consider the spatially non-homogeneous version of (2.6), that is, the boundary value problem ǫ h xx + f ( h, x ) − Z f ( h, x ) d x = 0 , < x < , (5.5)where now f ( h, x ) = Π LR ( h, x ) or f ( h, x ) = Π SR ( h, x ) . (5.6)subject to the periodic boundary conditions and the volume constraint, h (0) = h (1) , h x (0) = h x (1) , and Z h ( x ) d x = ¯ h. (5.7)Seeking an asymptotic solution to (5.5)–(5.7) in the form h ( x ) = ¯ h + ρh ( x ) + O ( ρ )in the limit ρ →
0, by substituting this anzatz into (5.5) we find that in the case of Π LR ( h, x )we have h ( x ) = − ¯ h sin (2 πx )4 π ¯ h ǫ − h + 6 , (5.8)while in the case of Π SR ( h, x ) the equivalent result is h ( x ) = ¯ h sin (2 πx )4 π ¯ h ǫ − h + 6 . (5.9)For non-zero values of ρ , in both the Π LR and Π SR cases, the changes in the bifurcationdiagrams obtained in the homogeneous case ( ρ = 0) are an example of forced symmetry11reaking (see, for example, Chillingworth [6]), which we discuss further in Appendix A.More precisely, we show there that when ρ = 0, out of the entire O (2) orbit only twoequilibria are left after symmetry breaking.Figure 2 shows how for small wettability contrast | ρ | ≪
1, the resulting spatial non-homogeneity introduces imperfections [12] in the bifurcation diagrams of the homogeneouscase ρ = 0 discussed in Section 4. It presents bifurcation diagrams in regimes I, II and IIIwhen the disjoining pressure Π LR is given by (5.3) for a range of small values of ρ togetherwith the corresponding diagrams in the case ρ = 0. The corresponding bifurcation diagramswhen the disjoining pressure Π SR is given by (5.4) are very similar and hence are not shownhere. /ǫ k h − ¯ h k (a) PSfrag replacements1.5 /ǫ k h − ¯ h k (b) PSfrag replacements1.5 /ǫ k h − ¯ h k (c) Figure 2: Bifurcation diagrams of solutions with a single maximum, showing k h − ¯ h k plotted as a function of 1 /ǫ when the disjoining pressure is Π LR for ρ = 0 (dashed curves), ρ = 0 .
005 (doted curves) and ρ = 0 .
05 (solid curves) for (a) ¯ h = 1 .
24, (b) ¯ h = 1 .
3, and (c)¯ h = 2, corresponding to regimes I, II, and III, respectively.12or large wettability contrast, specifically for | ρ | ≥
1, significant differences between the twoforms of the disjoining pressure are to be expected. When using Π LR , one expects globalexistence of positive solutions for all values of | ρ | ; see for example Wu and Zheng [39]. Onthe other hand, when using Π SR , there is the possibility of rupture of the liquid film for | ρ | ≥
1; see for example [3, 39], which means in this case we do not expect positive solutionsfor sufficiently large values of | ρ | . ρ k h − ¯ h k Figure 3: Bifurcation diagram showing k h − ¯ h k plotted as a function of ρ when the disjoiningpressure is Π LR , ¯ h = 3 and 1 /ǫ = 50. The leading-order dependence of k h − ¯ h k on ρ as ρ →
0, computed using (5.8) is shown with dashed lines.In Figure 3 we plot the branches of the positive solutions of (5.5)–(5.7) with a single maxi-mum when the disjoining pressure is Π LR for ¯ h = 3 and 1 /ǫ = 50, so that when ρ = 0, weare in Regime II above the curve PF of pitchfork bifurcations from the constant solution(see Figure 1). The work of Bertozzi [3] and of Wu and Zheng [39], shows that strictlypositive solutions exist for all values of | ρ | , beyond the range [ − ,
2] for which we have weperformed the continuation.Figure 4 shows that the situation is different when the disjoining pressure is Π SR (with same¯ h and ǫ ). At | ρ | = 1, one of the branches of smooth solutions disappears due to ruptureof the film, so that at some point x ∈ [0 ,
1] we have h = 0 and a strictly positive solutionno longer exists, while the other two branches disappear at a saddle node bifurcation at | ρ | ≈ . ρ = 0 is in fact a whole O (2)-symmetric orbit of solutions predicted by the analysis leading to Figure 1.13 k h − ¯ h k Figure 4: Bifurcation diagram with k h − ¯ h k plotted as a function of ρ when the disjoiningpressure is Π SR , for ¯ h = 3 and 1 /ǫ = 50. The leading order dependence of k h − ¯ h | on ρ as ρ →
0, computed using (5.9), is shown with dashed lines.Note that when the disjoining pressure is Π SR , given by (5.4), we were unable to use AUTOto continue branches of solutions beyond the rupture of the film. xh Figure 5: Solutions h ( x ) when the disjoining pressure is Π SR for ¯ h = 2 and 1 /ǫ = 30 for ρ = 0, 0.97, 0.98, 0.99 and 1, denoted by “1”, “2”, “3”, “4” and “5”, respectively, showingconvergence of strictly positive solutions to a non-strictly positive one as ρ → − .14igure 5 shows the film approaching rupture as ρ → − at the point where the coefficientof the short-range interaction term disappears when ρ = 1, i.e. when 1 + sin(2 πx ) = 0 andhence at x = 3 /
4. These numerical results are consistent with the arguments of Bertozzi etal. [3].Investigation of the possible leading-order balance in (2.6) when the disjoining pressure isΠ SR and ρ = 1 in the limit x → / h ( x ) = O ( x − / / ; continuing theanalysis to higher order yields h = (2 π ) (cid:18) x − (cid:19) − ǫ (cid:0) π (cid:1) / (cid:18) x − (cid:19) + O (cid:18) x − (cid:19) ! . (5.10)Figure 6 shows the detail of the excellent agreement between the solution h ( x ) when ρ = 1and the two-term asymptotic solution given by (5.10) in a neighbourhood of x = 3 / xh h ( x ) AsymptoticSolution Figure 6: Details near x = 3 / h ( x ) shown with solid line when thedisjoining pressure is Π SR and ρ = 1, and the two-term asymptotic solution given by (5.10)shown with dashed lines for ¯ h = 2 and ǫ = 1 / LR and Π SR , respectively, in the (1 /ǫ, ρ )-plane in the three regimesshown in Figure 1.In Figure 8 the horizontal dashed line at ρ = 1 indicates rupture of the film and loss of asmooth strictly positive solution, implying that there are regions in parameter space whereno such solutions exist. 15 /ǫρ (a) PSfrag replacements1.5 /ǫρ (b) PSfrag replacements1.5 /ǫρ (c) Figure 7: Multiplicity of strictly positive steady state solutions with a single maximumin the (1 /ǫ, ρ )-plane when the disjoining pressure is Π LR for (a) ¯ h = 1 .
24 (Regime I), (b)¯ h = 1 . h = 2 (Regime III).16 /ǫρ (a) PSfrag replacements1.5 /ǫρ (b) PSfrag replacements1.5 /ǫρ (c) Figure 8: Multiplicity of strictly positive steady state solutions with a single maximumin the (1 /ǫ, ρ )-plane when the disjoining pressure is Π SR for (a) ¯ h = 1 .
24 (Regime I), (b)¯ h = 1 . h = 2 (Regime III). In the present work we have investigated the steady state solutions of the thin-film evolu-tion equation (2.1) both in the spatially homogeneous case (2.6)–(2.9) and in the spatially17on-homogeneous case for two choices of spatially non-homogeneous Derjaguin disjoiningpressure given by (5.3) and (5.4). We have given a physical motivation for our choicesof the disjoining pressure. For reasons explained in the last paragraph of Section 2, weconcentrated on branches of solutions with a single maximum.In the spatially homogeneous case (2.6)–(2.10), we used the Liapunov-Schmidt reduction ofan equation invariant under the action of the O (2) symmetry group to obtain local bifurca-tion results and to determine the dependence of the direction and nature of bifurcation inbifurcation parameter 1 /ǫ = l on the average film thickness ¯ h ; our results on the existenceof three different bifurcation regimes, (namely nucleation, metastable, and unstable) aresummarised in Propositions 1 and 2 and in the phase diagram shown in Figure 1 obtainedusing AUTO.In the spatially non-homogeneous case (5.5)–(5.7), we clarified the O (2) symmetry breakingphenomenon (see Appendix A) and presented imperfect bifurcation diagrams in Figure 2and global bifurcation diagrams using the wettability contrast ρ as a bifurcation parameterfor fixed ǫ and ¯ h in Figures 3 and 4.Our results can be described using the language of global compact attractors [39]. Let usdiscuss Figure 3 in more detail; it is reproduced in Figure 9 with labelling added to thedifferent branches of strictly positive steady state solutions with a single maximum. Belowwe explain what these different labels mean. ρ X O BAC k h − ¯ h k Figure 9: Figure 3 with the different solution branches labelled.18
Figure 10: Sketch of the global attractor for ρ = 0. The thicker circle represents the O (2)orbit of steady states and O represents the constant solution h ( x ) = ¯ h . BAC
Figure 11: Sketch of the global attractor for small non-zero values of | ρ | . The points A , B , C correspond to the steady state solutions labelled in Figure 9.When ρ = 0, for 1 /ǫ = 50, the attractor is two-dimensional; the constant state h ≡ ¯ h denoted by O has a two-dimensional unstable manifold and X corresponds to a whole O (2)orbit of steady state solutions. A sketch of the attractor in this case is shown in Figure 10.When | ρ | takes small positive values, only two steady state solutions, denoted by A and C remain from the entire O (2) orbit, as discussed in Appendix A, while the constant steadystate continues to B without change of stability. The resulting two-dimensional attractoris sketched in Figure 11.Increasing | ρ | causes equilibria B and C to collide in a saddle node bifurcation and disappear,so that the attractor degenerates to a single asymptotically stable equilibrium. It would beinteresting to understand why this collision of equilibrium branches has to occur.19e have also explored the geometry of film rupture which occurs as ρ → − when thedisjoining pressure is given by Π SR ; this phenomenon is shown in detail in Figures 5 and 6.Finally, in Figures 7 and 8, we showed the results of a two-parameter continuation study inthe (1 /ǫ, ρ ) plane, showing how the multiplicity of positive steady state solutions changes asparameters are varied, and in particular, indicating in the case of Derjaguin pressure Π SR shown in Figure 8 regions in parameter space where no such solutions exist. We conjecturethat in these regions the solution of the unsteady problem with any positive initial conditionconverges to a weak solution of the thin-film equation with regions in which h ( x ) = 0, i.e.solutions with film rupture (dry spots). For a discussion of such (weak) solutions of thin-film equations in the homogeneous case the reader is referred to the work of Laugesen andPugh [22].In the case of disjoining pressure (5.4), we could not use the AUTO-07p version of AUTOto continue branches of solutions beyond touchdown. It would be an interesting project todevelop such a capability for this powerful and versatile piece of software.Figures 8(b) and 8(c) provide numerical evidence for the existence of a curve of saddle-node bifurcations converging to the point (0 ,
1) in the (1 /ǫ, ρ ) plane; an explanation for thisfeature of the global bifurcation diagrams requires further study.To summarise: our study was primarily motivated by [13]. While we have clarified themathematical properties of (2.6)–(2.9) and (5.5)–(5.7), so that the structure of bifurcationsin Figure 3(a) of that paper for non-zero values of ρ is now understood, many of their othernumerical findings are still to be explored mathematically, for example, the stability of ridgesolutions shown in their Figure 5.A final remark that might be of interest to the reader is that the solutions of equations(5.5) – (5.7) are the steady state solutions of the one-dimensional spatially inhomogeneousversion of the evolution equation (2.1), which is a degenerate quasi-linear fourth orderequation of the form h t = ( Q ( h ) P ( x, h ) x ) x , < x < L. (6.1)These solutions can also be thought of as the steady state solutions of a much simpler(Rubinstein-Sternberg type [28]), second order semi-linear non-local equation, h t = γh xx + Π( h, x ) − L Z L Π( h, x ) d x, < x < L. (6.2)It would be interesting to compare the dynamics of (6.1) and (6.2), for example using thespectral comparison principles of Bates and Fife [2]. Acknowledgement
We are grateful to Prof. U. Thiele (University of M¨unster) for clarifications concerning thework of Honisch et al. [13] and for sharing with us the AUTO codes used in that work whichformed the basis of our continuation analysis.20 O (2) Symmetry Breaking by Spatial Non-homogeneity
In this Appendix, we present an argument that shows that when the wettability contrast ispresent, i.e. when ρ = 0, the breaking of the O (2) symmetry which equation (5.5) with theperiodic boundary conditions (5.7) has for ρ = 0 , leaves only two stationary solutions.This is, in principle, a known result (see, for example, Chillingworth [6]), but, since weare not aware of an easily accessible reference, we give the details here. As before, we set G ( x ) = sin(2 πx ). We provide the proof for Π SR given by (5.4), the proof for Π LR given by(5.3) is similar.For the case of Π SR , let us rewrite the boundary value problem (5.5) in the form ǫ h xx + f ( h ) + ρf ( h ) G ( x ) − Z [ f ( h ) + ρf ( h ) G ( x )] d x = 0 , < x < , (A.1)where f ( h ) = 1 h − h , and f ( h ) = 1 h , i.e. we separate the spatially homogeneous and spatially non-homogeneous components ofthe disjoining pressure. Equation (A.1) is subject to the periodic boundary conditions (5.7).Suppose that when ρ = 0 there is an orbit of stationary solutions, i.e. a continuous closedcurve of solutions h ,s ( x ), parameterised by s ∈ R / [0 , h ,s ( x ) := h ( x + s ),for some function h ( x ), i.e. all these solutions are related by translation. We aim tounderstand what remains of this orbit for small non-zero ρ .Fix s ∈ R / [0 , h ( x ) = h ,s ( x ) + ρh ( x ) + O ( ρ ) . Substituting this expansion into (A.1) and collecting the O ( ρ ) terms, we have ǫ h ,xx + ( f ′ ( h ,s ) + f ′ ( h ,s )) h − Z (cid:2) f ′ ( h ,s ) + f ′ ( h ,s ) (cid:3) h d x = − f ( h ,s ) G + Z f ( h ,s ) G d x = 0 , (A.2)where, just like h ,s ( x ), h ( x ) also satisfies the periodic boundary conditions (5.7).Now set Ku := ǫ u ,xx + ( f ′ ( h ,s ) + f ′ ( h ,s )) u − Z [ f ′ ( h ,s ) + f ′ ( h ,s )] u d x, and let D ( K ), the domain of K , be D ( K ) = (cid:8) f ∈ C ([0 , | f (0) = f (1) , f ′ (0) = f ′ (1) (cid:9) . K is self-adjoint with respect to the L ([0 , L ([0 , Ku = 0.Next, we show that u := h ′ ,s solves Ku = 0. Indeed, by differentiating (A.1) with ρ = 0with respect to x , we find that u solves the equation ǫ u xx + ( f ′ ( h ,s ) + f ′ ( h ,s )) u = 0 . Integrating this equation over the interval [0 , Z ( f ′ ( h ,s ) + f ′ ( h ,s )) u d x = 0 . Hence 0 = ǫ u xx + ( f ′ ( h ,s ) + f ′ ( h ,s )) u = ǫ u xx + ( f ′ ( h ,s ) + f ′ ( h ,s )) u + Z ( f ′ ( h ,s ) + f ′ ( h ,s )) u d x = Ku.
Also note that as h ,s ( x ) satisfies periodic boundary conditions, Z h ′ ,s ( x )d x = 0 . (A.3)Hence the solvability condition for (A.2) is Z h ′ ,s ( r ) (cid:20) − f ( h ,s ) G + Z f ( h ,s ) G d x (cid:21) d r = 0 . (A.4)By (A.3), this condition reduces to Z f ( h ,s ) h ′ ,s G d x = 0 . (A.5)Now recall that h ,s ( x ) = h ( x + s ), so if we write F ( x + s ) = f ( h ( x + s )) h ′ ( x + s ), thefunction F ( · ) is 1-periodic in x with zero mean. Hence F ( z ) = ∞ X k =1 α k sin(2 kπz ) + β k cos(2 kπz ) . Therefore for G ( x ) = sin(2 πx ), the solvability condition for (A.2) becomes α sin(2 kπs ) − β cos(2 πs ) = 0 , (A.6)which has two solutions s ∈ R / [0 , h ( x ) onlyfor two choices of s ∈ R / [0 , O (2) orbit that exists for ρ = 0. 22 eferences [1] Ajaev, V. S., Gatapova, E. Y. & Kabov, O. A. (2011) Rupture of thin liquid films onstructured surfaces. Physical Review E , 041606.[2] Bates, P. W. & Fife, P. C. (1990) Spectral comparison principles for the Cahn-Hilliardand phase-field equations, and time scales for coarsening. Physica D: Nonlinear Phe-nomena , 335–48.[3] Bertozzi, A. L., Gr¨un, G. & Witelski, T. P. (2001) Dewetting films: bifurcations andconcentrations. Nonlinearity , 1569–92.[4] Brasjen, B. J. & Darhuber, A. A. (2011) Dry-spot nucleation in thin liquid films onchemically patterned surfaces. Microfluidics and Nanofluidics , 703–16.[5] Braun, R. J., King-Smith, P., Begley, C., Li, L. & Gewecke, N. (2015) Dynamics andfunction of the tear film in relation to the blink cycle. Progress in Retinal and EyeResearch , 132–64.[6] Chillingworth, D. (1985) Bifurcation from an orbit of symmetry. In: Pnevmatikos, S. N.(ed.) Singularities and Dynamical Systems . Amsterdam. Elsevier, pp. 285–94.[7] Chossat, P. & Lauterbach, R. 2000
Methods in Equivariant Bifurcations and DynamicalSystems . Singapore, World Scientific.[8] Craster, R. V. & Matar, O. K. (2009) Dynamics and stability of thin liquid films.
Reviews of Modern Physics , 1131–98.[9] Doedel, E. J. & Oldeman, B. E. 2009 Auto07p: Continuation and Bifurcation forOrdinary Differential Equations . Montreal. Concordia University.[10] Eilbeck, J. C., Furter, J. E. & Grinfeld, M. (1989) On a stationary state characterizationof transition from spinodal decomposition to nucleation behaviour in the Cahn–Hilliardmodel of phase separation.
Physics Letters A , 272–75.[11] Glasner, K. B. & Witelski, T. P. Collision versus collapse of droplets in coarsening ofdewetting thin films.
Physica D: Nonlinear Phenomena , 80–104.[12] Golubitsky, M. & Schaeffer, D. G. 1985
Singularities and Groups in Bifurcation Theory,Vol. I . New York. Springer.[13] Honisch, C., Lin, T.-S., Heuer, A., Thiele, U. & Gurevich, S. V. (2015) Instabilities oflayers of deposited molecules on chemically stripe patterned substrates: ridges versusdrops.
Langmuir , 10618–31.[14] Howison, S. D., Moriarty, J. A., Ockendon, J. R., Terrill, E. L. & Wilson, S. K. (1997)A mathematical model for drying paint layers. Journal of Engineering Mathematics , 377–94.[15] Huppert, H. E. (1982) The propagation of two-dimensional and axisymmetric viscousgravity currents over a rigid horizontal surface. Journal of Fluid Mechanics , 43–58.2316] Ji, H. & Witelski, T. P. (2017) Finite-time thin film rupture driven by modified evap-orative loss.
Physica D: Nonlinear Phenomena , 1–15.[17] Karnaushenko, D., Kang, T., Bandari, V. K., Zhu, F. & Schmidt, O. G. (2020) 3 D Self-assembled microelectronic devices: concepts, materials, applications.
Advances inMaterials , 1902994.[18] Kistler, S. F. & Schweizer, P. M. (eds). 1997 Liquid Film Coating: Scientific Principlesand Their Technological Implications . New York. Chapman & Hall.[19] Kitavtsev, G., Lutz, R. & Wagner, B. (2011) Centre manifold reduction approach forthe lubrication equation.
Nonlinearity , 2347–69.[20] Konnur, R., Kargupta, K. & Sharma, A. (2000) Instability and morphology of thinliquid films on chemically heterogeneous substrates. Physical Review Letters , 931–34.[21] Laugesen, R. S. & Pugh, M. C. (2000a) Linear stability of steady states for thin filmand Cahn-Hilliard type equations. Archive for Rational Mechanics and Analysis ,3–51.[22] Laugesen, R. S. & Pugh, M. C. (2000b) Properties of steady states for thin film equa-tions.
European Journal of Applied Mathematics , 293–351.[23] Liu, W. & Witelski, T. P. (2020) Steady-states of thin film droplets on chemicallyheterogeneous substrates. arXiv preprint arXiv:2002.11286 .[24] Ondar¸cuhu, T. & Aim´e, J. P. 2013 Nanoscale Liquid Interfaces: Wetting, Patterningand Force Microscopy at the Molecular Scale . Stanford. Pan.[25] Oron, A., Davis, S. H. & Bankoff, S. G. (1997) Long-scale evolution of thin liquid films.
Reviews of Modern Physics , 931–80.[26] Pismen, L. M. (2002) Mesoscopic hydrodynamics of contact line motion. Colloids andSurfaces A , 11–30.[27] Quake, S. R. & Scherer, A. (2000) From micro-to nanofabrication with soft materials.
Science , 1536–40.[28] Rubinstein, J. & Sternberg, P. (1992) Nonlocal reaction–diffusion equations and nucle-ation.
IMA Journal of Applied Mathematics , 249–64.[29] Rynne, B. P. & Youngson, M. A. 2008 Linear Functional Analysis , 2nd edition. Berlin.Springer.[30] Sehgal, A., Ferreiro, V., Douglas, J. F., Amis, E. J. & Karim, A. (2002) Pattern-directeddewetting of ultrathin polymer films.
Langmuir , 7041–48.[31] Sharma, A., Konnur, R., & Kargupta, K. (2003) Thin liquid films on chemically het-erogeneous substrates: self-organization, dynamics and patterns in systems displayinga secondary minimum. Physica A: Statistical Mechanics and its Applications , 262–78. 2432] Starov, V. M. (2013) Disjoining pressure. In: Li, D. (ed).
Encyclopedia of Microfluidicsand Nanofluidics . Berlin. Springer.[33] Thiele, U. & Knobloch, E. (2006) On the depinning of a driven drop on a heterogeneoussubstrate.
New Journal of Physics , 313.[34] Thiele, U., Brusch, L., Bestehorn, M. & B¨ar, M. (2003) Modelling thin-film dewettingon structured substrates and templates: bifurcation analysis and numerical simulations. European Physical Journal E Ordinary and Partial Differential Equations ,Berlin. Springer, pp. 356–65.[36] Wang, D., Song, Y., Tian, J., Shigu, E. & Haidak, G. (2018) Research on the fluid filmlubrication between the piston-cylinder interface.
AIP Advances , 105330.[37] Witelski, T. P. (2020) Nonlinear dynamics of dewetting thin films. AIMS Mathematics Physica D: Nonlinear Phenomena , 155–76.[39] Wu, H. & Zheng, S. (2007) Global attractor for the 1-D thin film equation.
AsymptoticAnalysis , 101–11.[40] Xue, L., Zhang, J. & Han, Y. (2012) Phase separation induced ordered patterns in thinpolymer blend films. Progress in Polymer Science37