Rupture of thin films with resonant substrate patterning

Journal of Colloid and Interface Science 303 (2006) 532–545 www.elsevier.com/locate/jcis

Rupture of thin films with resonant substrate patterning Justin C.-T. Kao, Alexander A. Golovin, Stephen H. Davis ∗ Department of Engineering Sciences and Applied Mathematics, Northwestern University, Evanston, IL 60208, USA Received 2 June 2006; accepted 4 August 2006 Available online 17 August 2006

Abstract We study the stability and rupture of thin liquid films on patterned substrates. It is shown that striped patterning on a length scale comparable to that of the spinodal instability leads to a resonance effect and an imperfect bifurcation of equilibrium film shapes. Weakly nonlinear analysis gives predictions for film shapes, stability, growth rates, and rupture times, which are confirmed by numerical solution of the thin-film equation. Film behavior is qualitatively different in the resonant patterning regime, but with sufficiently large domains rupture occurs on a spinodal length scale regardless of patterning. Instabilities transverse to the patterning are examined and shown to behave similarly as disturbances to films on uniform substrates. We explain some previously reported effects in terms of the imperfect bifurcation. © 2006 Elsevier Inc. All rights reserved. Keywords: Thin liquid films; Rupture; Dewetting; Patterning; Templating; Imperfect bifurcation; Nonlinear dynamics

1. Introduction The behavior of thin liquid films is of great importance to a variety of applications—for example, paint, lubrication, lithographic printing, foams, and most recently, semiconductor fabrication and microfluidic devices. For sufficiently thin films (100 nm), it is known that van der Waals forces may cause instability leading to film rupture and dewetting [1–5]. Given an appropriate functional form for the van der Waals potential, a critical wavelength and most dangerous mode can be found [6], and the evolution equation governing the nonlinear behavior derived [7]. Bifurcation analysis gives a Landau equation for perturbation amplitude and a nonlinear estimate for rupture time [8–11], while numerical simulations show the full evolution to rupture [12,13] and confirm self-similar solutions near rupture [14]. Effects such as thermocapillarity, evaporation, slip, and attractive–repulsive van der Waals potentials have also been considered [9,15]. Until recently however, theoretical investigations of dewetting have assumed homogeneous substrates, whereas in engineering applications it is likely that the solid surface is heterogeneous, with either deliberate patterning or natural imperfections. * Corresponding author. Fax: +1 847 491 2178.

E-mail address: [email protected] (S.H. Davis). 0021-9797/$ – see front matter © 2006 Elsevier Inc. All rights reserved. doi:10.1016/j.jcis.2006.08.015

Experimental studies of thin films on patterned surfaces demonstrate that dewetting is indeed affected by substrate inhomogeneities. In addition to film holes due to surface imperfections seen in experiments such as Stange and Evans [4], organized rupture from deliberately patterned substrates has also been shown. Meyer and Braun [16] and Zhang et al. [17–19] report “templating” of polystyrene films, while Higgins and Jones [20] and Sehgal et al. [21] observe that although substrate anisotropy or patterning may dictate the location and orientation of rupture, the morphological length scale remains close to that of spinodal dewetting on a homogeneous substrate— templating can only occur when the imposed pattern is sizecompatible with spinodal dewetting. The groups of both Sehgal and Zhang [18,19,21] further find that on striped substrates, the rupture location moves from stripe centers to stripe edges as the stripe width increases. A series of numerical simulations by Kargupta, Konnur, and Sharma [22–29] have clarified and in some cases inspired these experimental studies. Different intermolecular potentials and substrate patternings (stripes, dots, and checkerboards) are considered. As in experiment, it is found that film templating requires an appropriate patterning length scale and that rupture location changes as a function of length scale. It is demonstrated that rupture is accelerated on heterogeneous surfaces, while eventual film morphology follows spinodal dewetting length

J.C.-T. Kao et al. / Journal of Colloid and Interface Science 303 (2006) 532–545

scales when the corresponding homogeneous substrate would be destabilizing. For systems in which the film is stable on a homogeneous substrate, it is observed that a sufficiently large heterogeneity can engender rupture in the previously stable film. Morphological evolution and coarsening subsequent to rupture is also computed, with varying levels of self-organization found. There has also been interest in the behavior of flowing films [30–33] and droplets (fully ruptured films) [34–36] on patterned surfaces, as well as films on substrates with patterning of a length scale comparable to molecular radii [37,38]. These cases are not be considered further here; we restrict our attention to the simplest possible physical situation. Despite this array of experiments and simulations, there have been few systematic studies, with most authors focusing on exhibiting proper templating under specific sets of conditions. The most comprehensive theoretical examinations of film rupture are those of Thiele et al. [39,40], which investigate stationary states, stability, and post-rupture coarsening of the film, confirming previous experiments and simulations by means of bifurcation analysis. Though more informative than direct simulation, the results still rely on numerical solution of the equations. Furthermore, only a single heterogeneity wavelength is considered, far from the neutral stability threshold. To our knowledge there has been no analytic study of film stability and rupture on patterned substrates. This is our chief aim. In this paper, we examine the van der Waals instability of “(1 + 1)D” and “(2 + 1)D”1 thin films on stripe-patterned substrates and derive its dependence on the patterning wavenumber kp by analytic means. We find a resonance effect when the patterning wavenumber is close to the stability threshold kc for the homogeneous system, and show the stability and rupture of the film to be governed by an imperfect pitchfork bifurcation of the film thickness. Thus, we predict qualitatively different behavior for a film on a patterned substrate with kp ≈ kc , than would be expected for a film on homogeneous substrate, or for a film with kp far from kc . As kc is within an experimentally feasible range for many common physical systems, we expect that a consideration for the effects of resonance may be important in interpreting experiments and simulations. Our analytic results are confirmed by numerical calculations of stability and film evolution. We show also that real-world (2 + 1)D films on striped substrates, having disturbances transverse to the patterning, must be bounded in the transverse direction for the film to be stable. We begin by deriving the (1 + 1)D thin-film evolution equation for the 1D interface of a 2D liquid. Substrate heterogeneity is modeled by allowing the Hamaker constant to vary in space, A = A(x). When the Hamaker “constant” has small sinusoidal variation, bifurcation analysis enables derivation of a Landau equation for disturbance amplitude, giving rupture time estimates and the amplitude and stability of stationary states as a function of patterning wavelength. We find an imperfect pitchfork bifurcation, and a relationship between imperfection 1 One- and two-dimensional interfaces of two- and three-dimensional liquid films, respectively.

533

Fig. 1. Thin liquid film on solid substrate.

size and patterning amplitude. The bifurcation is subcritical for van der Waals potentials like h−3 , while for layered (attractive/repulsive) potentials, a supercritical bifurcation is possible. Steady states are found numerically for sinusoidal and square-wave patterning. Using Floquet theory, their stability and growth rates are investigated. Rupture time is calculated via numerical simulation of the full evolution equation, and it is shown that flat-film initial conditions are far enough from equilibrium for patterned substrates as to give significantly different rupturing behavior than seen for near-stationary initial conditions on both patterned and uniform substrates. Finally, we examine the stability of the (2 + 1)D film on striped substrates and observe that all stationary states found for the (1 + 1)D film are subject to long-wave transverse instability, as known for homogeneous substrates. 2. Thin-film evolution equation in (1 + 1)D Consider a thin liquid film on a flat solid substrate (z = 0) with the surface of the film at z = h(x, t), as shown in Fig. 1. For an incompressible liquid and conditions of no-slip at the liquid–solid interface, zero tangential stress at the liquid–air interface, and constant surface tension, the governing Navier– Stokes system in the fluid (0 < z < h) is ∇ · u = 0,

(1)

ρ(ut + u · ∇u) = −∇p − ∇φ + μ∇ u. 2

(2)

At the substrate (z = 0), u = 0. At the surface of the film (z = h),  −1  −p + 2μ ux h2x − (uz + wx )hx 1 + h2x  −3/2 = σ hxx 1 + h2x , (3) 2(ux − wz )hx − uz − wx + (wx + uz )h2x = 0,

(4)

w = ht + uhx ,

(5)

where ρ, μ, σ are constant density, viscosity and surface tension, respectively, u(x, z, t) = (u, w) is fluid velocity, and A ¯ φ(x, h) = 6π φ(h)f (x) is the van der Waals potential. Following the derivation in Oron et al. [15] with a modified lubrication scaling for the long-wave small parameter η  1,   √ ηx z h (X, Z) = AC , H= , (6) , h0 h0 h0   ημ ρh0 ρh0 (U, W ) = u, w , T = A2 C 2 t, (7) μ ημ ρh0 ¯ ηρh2 φ(h) f (x), ϕ(X, H ) = P = 2 0 p, (8) ¯ 0) μ φ(h

534

J.C.-T. Kao et al. / Journal of Colloid and Interface Science 303 (2006) 532–545

Table 1 Cutoff and fastest-growing wavelengths (λc , λ∗ ) for various film thicknesses (h0 ), based on σPS = 0.033 N/m, σPVC = 0.039 N/m, APS = 5.0 × 10−22 J, and APVC = 6.7 × 10−21 J (estimated from tables in [41]) h0 (nm) 2 5 10 20 50 Fig. 2. Dispersion relation for ϕ = H −3 .

we expand in η and obtain the scaled evolution equation for H (X, T ),  3 H HT = (9) (−HXX + ϕ)X , 3 X where C = μ2 /η3 ρh0 σ is the capillary number, A = ηA ρh20 × ¯ 0 )/6πμ2 is the nondimensional Hamaker parameter, ϕ is φ(h the normalized van der Waals potential (ϕ = H −3 F (X) for instance), and h0 is taken to be the mean film thickness. 3. Homogeneous substrate We first summarize previous results for thin films on homogeneous substrates, where the Hamaker parameter is spatially uniform. 3.1. Stability Consider Eq. (9) with a purely attractive and unpatterned van der Waals potential as used in Ruckenstein and Jain [6] and Williams and Davis [7], A , φ= 6πh3

(10)

which corresponds to ϕ = H −3 . Linear stability analysis (δ  1) with H = 1 + δeikX+ωT yields the dispersion relation  1 (11) 3 − k2 k2, 3 as shown in Fig. 2. The √flat film state is stable to short-wave perperturbaturbations (k > kc ≡ 3 ) and unstable to long-wave √ tions (k < kc ), with largest growth rate at k ∗ ≡ 3/2. This corresponds to a dimensional stability threshold λc = (2π)3/2 h20 ×



σ/A and fastest-growing wavelength λ∗ = 4π 3/2 h20 σ/A . To stabilize a film on a homogeneous substrate, the domain must be smaller than λc , while for larger domains the film will rupture and display a characteristic length scale λ∗ . Table 1 gives estimates for λc and λ∗ for several liquids on a quartz substrate, calculated from data in Israelachvili [41]. In [9] and [11], a layered van der Waals potential (including a repulsive term) is used, ω=

φ=

B A − , 6πh3 hn

(12)

Polystyrene

Polyvinyl chloride

λc (µm)

λ∗ (µm)

λc (µm)

λ∗ (µm)

0.36 2.3 9.0 36 230

0.11 0.67 2.7 11 67

0.51 3.2 13 51 32

0.15 0.95 3.8 15 95

which, after rescaling, gives ϕ = H −3 − BH −n , where B =  6πB  / hn−3 0 A . Linear stability analysis gives a similar dispersion relation ω=

 1 2 k − k2 k2, 3 c

(13)

where now the cutoff wavenumber is kc2 = 3 − nB. The instability can occur when kc2 > 0, i.e., for B < 3/n; otherwise, the repulsive term dominates and all modes are linearly stable. The fastest growing mode corresponds to k ∗2 = (3 − nB)/2. Many other forms of the intermolecular potential are possible [42], and are possibly more realistic, but these are sufficient for our purposes—examining stability and initial rupture. 3.2. Bifurcation analysis Using bifurcation analysis [8–11], one obtains a Landau equation for the amplitude Q1 and phase ϑ1 of a film H (X) = 1 + Q1 cos(kX + ϑ1 ) close to the stability threshold, k = (1 +  2 k2 )kc : Q˙ 1 = αQ1 + κQ31 , ϑ˙1 = 0.

(14) (15)

Note that the parameters , k2 , and Q1 are determined by two equations so we may take k2 = 1 or pick an arbitrary small ; this will not be the case when the substrate has patterning. The constants α and κ are 2 α = − ϕ  (1)2 k2 , (16) 3  1   2 ϕ (1) + 3ϕ  (1)ϕ  (1) , κ= (17) 72 consistent with [9–11]. In steady state, (14) describes a pitchfork bifurcation for the amplitude of the perturbed solution. For the purely attractive (ϕ = H −3 ) van der Waals potential, κ = 19 2 and the subcritical pitchfork opens to the right as √ shown in Fig. 3a. Stationary solutions are Q1 = 0 or Q1 = ± 12k2 /19. On the other hand, with the attractive/repulsive van der Waals potential ϕ = H −3 − BH −n , it is possible for κ to change sign, giving a supercritical bifurcation [9,11]. In particular, there is a range of B for which κ < 0, at B slightly less than B for which the flat film is stable (Fig. 4). In this region the instability saturates, giving a stationary corrugated film as shown in [9] and [11].

J.C.-T. Kao et al. / Journal of Colloid and Interface Science 303 (2006) 532–545

(a)

535

(b)

Fig. 3. Steady state perturbation amplitude Q1 vs wavenumber deviation k2 from the stability threshold. (a) Subcritical bifurcation for ϕ = H −3 , (b) supercritical bifurcation for ϕ = H −3 − BH −n with B = 0.6, n = 4. Solid lines correspond to stable branches and dashed lines to unstable ones.

(a)

(b)

Fig. 4. (a) κ vs B for n = 4, 5, . . . , 9 (shortest to longest dashes). The bifurcation is supercritical for κ negative. Curves are truncated at B above which the flat film is linearly stable. (b) Parameter space for attractive/repulsive van der Waals. Supercritical bifurcation in the black region, subcritical in the white regions, and no instability exists in the gray region.

See Fig. 3b. The stationary amplitudes are Q1 = 0 or Q21 = 48kc4 k2 /[n2 B 2 (4n2 + 11n + 7) − 3nB(3n2 + 17n + 74) + 684]. Equation (14) can be solved in implicit form for τ (Q1 ) to find a nonlinear estimate of rupture time [8–11]:  1 α ln 1 + τr = lim τ = (18) . Q1 →∞ 2α 2κQ1 (0)2 In the limit α → 0, τr → [2κQ1 (0)2 ]−1 . Leshansky and Rubinstein [11] found a good agreement between this estimate near kc and numerical simulations of the evolution equation (9) far from kc . 4. Heterogeneous substrate: sinusoidal patterning in (1 + 1)D Consider a van der Waals potential with small sinusoidal variations,   ϕ(H, X) = ϕ0 (H ) 1 +  3 cos kp X , (19) where   1. This is a simplified model of a real substrate with imperfections of a length scale 2π/kp and strength  3 .

4.1. Bifurcation analysis When kp is far from kc , stationary solutions are found by regular perturbation, expanding H (X) = 1 +  3 H1 (X) + · · · to give H1 (X) − ϕ0 (1)H1 (X) = −ϕ0 (1)kp sin kp X. Since

kc2

≡ −ϕ0 (1)

(we assume

ϕ0 (1) < 0

H1 (X) = c1 cos kc x + c2 sin kc x −

so that

(20) kc2

ϕ0 (1) cos kp x. kp2 − kc2

> 0), (21)

This regular perturbation solution is also derived in Thiele et al. [39,40]. It describes well the effect of a modulation in A when the patterning wavenumber kp is not near kc , but becomes invalid for kp → kc , indicating resonance between the spatial forcing of the Hamaker parameter and the neutrally stable stationary solution. In order to examine the system near the instability threshold kc , we apply the method of strained coordinates with ξ = kp X, τ =  2 T and H (X, T ) = 1 + H1 (ξ, τ ) + · · · , where the patterning wavenumber is expanded as kp = (1 +  2 k2 )kc (see,

536

J.C.-T. Kao et al. / Journal of Colloid and Interface Science 303 (2006) 532–545

(a)

(b)

Fig. 5. An imperfect bifurcation occurs for spatially periodic Hamaker patterning. (a) ϕ = H −3 [1 +  3 cos ξ ]. The turning point of the lower curve is at √ √ 3 3  2 kp = kc ≡ (1 + 19 /4)kc , Q1 = −1/ 19. (b) ϕ0 = H −3 − BH −n . Here n = 4, B = 0.6. Solid lines correspond to stable branches and dashed lines to unstable ones.

e.g., [43]). At O(), the solution has wavenumber kp ,   H1 (ξ, τ ) = Q1 (τ ) cos ξ + ϑ1 (τ ) ,

(22)

with Q1 (τ ) and ϑ1 (τ ) as yet unknown. At O( 2 ),   H2 (ξ, τ ) = Q2 (τ ) cos ξ + ϑ2 (τ )    Q1 (τ )2 ϕ0 (1) + cos 2 ξ + ϑ1 (τ ) .  12ϕ0 (1)

kp2 =

(23)

At O( 3 ), the elimination of secular terms gives an equation for the amplitude Q1 , Q˙ 1 (τ ) = γ + αQ1 + κQ31 ≡ g(Q1 ),

1 ∓ d, and HXX ≈ ±dkp2 (as H ≈ 1 + d sin kp X), we find d is a root of the equation

(24)

and we find that ϑ1 = 0. Here α and κ remain the same as in Eq. (16), and γ = 13 ϕ0 (1)ϕ0 (1) < 0, creating an imperfection of the bifurcation [44]. 4.2. Stationary states In steady state, Eq. (24) describes an imperfect pitchfork bifurcation for the amplitude of the deviation from a flat film. Since α and κ are the same as in the homogeneous Hamaker parameter case, the main effect of the heterogeneity on the stationary solutions is to smooth the transition across kc . For the purely attractive van der Waals potential ϕ0 (H ) = H −3 , the equation for stationary state amplitude becomes 19Q31 = 12k2 Q1 + 2 (shown in Fig. 5a), and for ϕ0 (H ) = H −3 − BH −n , we have Q31 = [48kc4 k2 Q1 − 24kc2 (B − 1)] × [n2 B 2 (4n2 + 11n + 7) − 3nB(3n2 + 17n + 74) + 684]−1 (shown in Fig. 5b). The new stability threshold for the attractive van √ der Waals potential is at kc ≡ (1 + 3 19 2 /4)kc , and we use the purely attractive potential exclusively henceforth. As an alternative to bifurcation analysis, the amplitude of the stationary film corrugations can be estimated as follows. The evolution equation in steady state is integrated twice to give the static pressure balance at the interface in the long-wave approximation, −HXX + ϕ = P¯ , (25) where P¯ is constant across the entire domain. Then, assuming that the Hamaker parameter has extrema F0 ± δF , Hmin / max ≈

F0 (3 + d 2 ) − δF (d −1 + 3d) . (1 − d 2 )3

(26)

This gives a fully nonlinear estimate of the stationary film shape, shown as the dashed line in Fig. 6. (For the comparison, F0 = 1 and δF =  3 .) We verify these estimates by numerical solution of the evolution equation (9) in steady state. The full evolution equation with ϕ(H, X) = F (X)H −3 is spatially discretized in conservative form:  (H 3 )m+1/2 dHm D+ Fm = D− − D+ D+ D− Hm + dT 3 3  −1  − Fm+1/2 H (27) D H . m+1/2 + m D+ and D− are forward and backward spatial finite differences, respectively. The stationary states are found by a simple collocation requiring that the right-hand side of Eq. (27) be zero at all grid points in a periodic domain of size 2π/kp . A comparison is shown in Fig. 6, with steady states found via weakly nonlinear analysis (Eq. (24), solid line), nonlinear estimate (Eq. (26), dashed line), and numerically (Eq. (27), dots and squares). We find that Eq. (26) gives a rather good approximation to the numerical solution, even far from kc . Also, observe that for kp ≈ kc , the stationary film shape on patterned substrates differs substantially from the stationary film shape on a uniform substrate. The purely attractive potential used in Fig. 6 has  3 = 18 . It is worth noting that Fig. 5a of Thiele et al. [40] is in a plane perpendicular to the imperfect bifurcation shown here, in the three-dimensional bifurcation space of wavenumber (kp ), patterning strength () and equilibrium film amplitude. We have chosen kp for the bifurcation parameter while [40] uses . The diagrams are consistent—the result of Thiele can be recovered by varying  in the equations (valid for arbitrary 0 <   1). Also note that the mean film thickness in Eq. (9) has been normalized to unity. Increasing or decreasing the dimensional film thickness h0 for a fixed substrate patterning—as often oc-

J.C.-T. Kao et al. / Journal of Colloid and Interface Science 303 (2006) 532–545

Fig. 6. Stationary state film amplitude vs wavenumber (kp ) for perturbation solution (solid line), nonlinear estimate (dashed line), numerical solution for sinusoidal Hamaker patterning (dotted) and numerical solution for square-wave Hamaker patterning (open squares). The bifurcation diagram for the homogeneous substrate is shown in gray. F0 = 1, δF =  3 = 18 .

curs in experiment—has the effect of altering the nondimensionalized patterning wavenumber kp by a factor 1/ h20 , moving the system right or left respectively on the bifurcation diagram (Fig. 6), and sampling qualitatively distinct film behaviors [21]. 4.3. Stability of the branches Near the bifurcation point (Q ≈ 0 and kp ≈ kc ) the film is well described by the bifurcation analysis presented above and stability is easily found from the Landau equation (24). Farther from the threshold, numerical solution of Eq. (27) is required to find an accurate stationary state. In this section, we explore the stability of the numerically obtained bifurcation branches by means of the Floquet theory. Let H (X, T ) = H¯ (X) + δeωt eiqX h(X), where H¯ is the stationary solution, h is 2π/kp periodic, q is restricted to −kp /2  q  kp /2, and δ  1. Then linearizing in δ, we find that h and ω must satisfy the equation   H¯ 3 4 F Dq h + H¯ 2 H¯  Dq3 h + Dq2 h 3 H¯    F H¯  F + 2 − 5 2 Dq h H¯ H¯   F F H¯  F H¯  2 F  H¯  + (28) + ω h = 0, −5 2 −4 2 +8 H¯ H¯ H¯ H¯ 3 ∂ ∂ + iq and (·) = ∂X (·). We solve this numeriwhere Dq = ∂X cally as the eigenvalue problem for the finite difference matrix representing (28), and obtain the growth rates ω and eigenmodes h(X)eiqX of each stationary state. Fig. 7a shows the growth rates of disturbances to the three branches of stationary states, which are displayed in Fig. 7b, plotted against the wavenumber kd of the disturbances.2 We consider first spatially synchronous disturbances (kd = kp ) to the stable small-amplitude middle branch. These always decay,

537

indicating stability. Disturbances with longer wavelength than the patterning (kd = kp /n, n  2) follow roughly the dispersion relation for homogeneous substrates and as kp /kd = n → ∞ the system approaches the homogenization limit. (In Fig. 7a, the linear dispersion and kd = kp /10 growth rates are visually indistinguishable except at very small wavenumber.) Hence, to avoid rupture of the middle branch of √ equilibrium √ films, one must limit the system size to λc = 2π/ 3(1 + 3 19 2 /4); for larger systems, rupture will occur with a characteristic √ length √ scale near that of the fastest-growing mode, λ∗ ≡ 2 2π/ 3. This is in agreement with previous experiments by Higgins and Jones [20] (see Fig. 2), Sehgal et al. [21] (see Figs. 2 and 3), and simulations by Kargupta, Konnur, and Sharma [22–29]. We observe that rupture of the small-amplitude state may occur at every alternate stripe, or at every nth stripe, depending on the relative sizes of kp , kc , and allowable kd , which also confirms previous work [18,19,21–29]. The large-amplitude upper and lower branches are always unstable, and the film will rupture at every stripe, or at alternate stripes, again depending on kp , kc , and kd . Note that for very long wavelength disturbances to the unstable branches, the effective wavenumber may be at a multiple of the nominal wavenumber kd , and hence have a greater growth rate. 4.4. Rupture We now examine the rupturing of the film over a single wavelength of patterning, characterizing the rupture time and touchdown location(s) as functions of patterning wavelength and initial conditions. As with Eq. (14), the amplitude equation (24) can be solved in implicit form for τ (Q1 ), Q1 (τ )

τ=

dx g(x)

(29)

Q1 (0) Q1 (τ )

= Q1 (0)

=

3 i=1

c2 c3 c1 + + dx x − r1 x − r2 x − r3

  ci 3 

Q1 − ri Q1 − ri ci log , = log Q1 (0) − ri Q1 (0) − ri

(30)

(31)

i=1

where ri are the roots of g(r) = 0 such that Q1 (τ ) = ri is a stationary solution, ci = 1/(ri − rj )(ri − rk ) = 1/g  (ri ), and c1 + c2 + c3 = 0. This gives estimates for the nonlinear rupture time τr , in the case when rupture occurs synchronously with patterning, kd = kp . In contrast to the homogeneous Hamaker parameter case, here τr depends on the direction of rupture—whether Q1 → ∞ or Q1 → −∞: ci 3 

Q1 − ri τr+ = lim τ = lim log Q1 →∞ Q1 →∞ Q1 (0) − ri i=1

2 Alternatively, k can be thought of as the wavenumber of the domain, a d

lower bound on the wavenumber of allowable disturbances.

= log

3



i=1

Q1 (0) − ri

−1/g  (ri )

,

(32)

538

J.C.-T. Kao et al. / Journal of Colloid and Interface Science 303 (2006) 532–545

(a)

(b)

Fig. 7. (a) Growth rates for periodic disturbances to (1 + 1)D stationary films versus the disturbance wavenumber kd . Curves are shown for patterning at various multiples of the disturbance wavenumber: harmonic (kp = kd ⇔ q = 0, in black), subharmonic (kp = 2kd , or q = kp /2, in dark gray), and also kp = 10kd (q = kp /10, light gray). Line styles in (a) correspond to stationary states in (b): small-amplitude middle branch (thick solid lines), upper and lower branches (dashed), and dispersion relation for the flat film on homogeneous substrate (thin solid line).

(a) k2 = −1

(b) k2 ≈ 0

(c) k2 = 14

Fig. 8. Rupture time vs initial film amplitude for k2 = {−1, 0, 14 }. The heavy line is for periodic Hamaker patterning and the light line for constant Hamaker parameter, with τ =  2 T .

τr− = log

3



−1/g  (ri ) ri − Q1 (0) .

(33)

i=1

Plots of τr are shown in Fig. 8. In Fig. 8a, at kd below kc , the rupture time and stationary state are merely shifted slightly compared to those found for the film on a uniform substrate. Once kd  kp however (Figs. 8b, 8c), the stationary state differs significantly from the flat film. Rupture time on the heterogeneous surface is much reduced compared to rupture time on a homogeneous substrate for near-flat initial conditions, similar to Kargupta et al.’s simulations using piecewise-constant patterning [22–29]. However, we find this accelerated rupture to be inherent to spinodal dewetting on heterogeneous surfaces, rather than a separate instability. The full evolution equation is solved by applying the method of lines to the finite difference equation (27) on a periodic domain of size 2π/kp . Time-marching is handled by MATLAB’s ode15s, an implicit variable-order adaptive time stepping method based on numerical differentiation formulas [45]. Sufficient convergence is obtained with 200 spatial grid points, the default error tolerances,3 and a maximum time step size of 0.1 units. We verify that numerically calculated rupture times are consistent with the analytic predictions shown in Fig. 8. 3 0.1% relative and 10−6 absolute error.

Some typical time-evolution results for various initial conditions (near-flat films, and films near the unstable upper and lower bifurcation branches), are shown in Fig. 9. For longwavelength patterning kp < kc , the fastest-growing mode (λ = λ∗ ) dominates (kp = 0.25 plots in Figs. 9b and 9c), whereas for kp ≈ kc there is an interaction between the instability and the patterning causing the film to remain spatially synchronous with the patterning (kp = 1, 2, and 2.1 plots in Fig. 9). With kp > kc , the upper and lower branches may either rupture or evolve to the stable small-amplitude middle branch, depending on initial conditions (kp = 2.1, 2.5, 3 in Figs. 9a and 9b). The time and location of rupture depend strongly on initial conditions, as shown in Fig. 10. When kp < k ∗ and the film evolution begins near an equilibrium corresponding to the upper branch in Fig. 6, the location of the first touchdown point (black and gray circles in Fig. 10a) is not strongly correlated with the patterning, while for initial conditions of approximately flat films (black and gray squares in Fig. 10a), rupture tends toward the point of largest Hamaker parameter, Amax . On the other hand, if kp > k ∗ , the upper branch ruptures at either Amax or Amin , the lower branch ruptures only at Amax (black and gray triangles in Fig. 10a), and the flat film continues to rupture at Amax , provided kp < kc . The change in rupture location preferences for near-equilibrium films can be attributed to the fact

J.C.-T. Kao et al. / Journal of Colloid and Interface Science 303 (2006) 532–545

539

(a)

(b)

(c) Fig. 9. Typical time evolution (gray lines at T = 0.1) of initial conditions (thin black line), with sinusoidal patterning (dashed line). (a) Initial conditions near lower branch, (b) initial conditions near upper branch, (c) flat initial conditions. Final film shape at touchdown or equilibrium is shown (thick black line).

that when kp approaches k ∗ , the fastest-growing mode becomes synchronous with the patterning. Fig. 10b displays the corresponding rupture times. Below kc , the flat film (squares) ruptures considerably faster than the near-equilibrium states (circles). Furthermore, the effect of perturbations (black vs. gray points) in the initial conditions is significant for near-equilibrium conditions, whereas the nearflat film is so far from equilibrium that disturbances do not affect the rupture time much except at very small kp . These rupture times are consistent with the analytic predictions above and suggest that the accelerated rupture seen in experiments and previous simulations may be due to starting far from equilibrium in the first place, as seen from Fig. 6 in the region kp ≈ kc . Finally, note that above kc , rupture from the unstable equilibria occurs extremely quickly.

5. Heterogeneous substrate: sinusoidal patterning in (2 + 1)D Real-world films are three-dimensional with a two-dimensional interface. In this case the evolution equation becomes 

H3 2 ∇(−∇ H + ϕ) . HT = ∇ · 3

(34)

Here we consider only one-dimensional patterning, i.e., ϕ = ϕ(H, X) and the stationary states found for the (1 + 1)D film still apply. As before, Floquet theory gives the stability of each branch. Let H (X, Y, T ) = H¯ (X) + δeωt eilY eiqX h(X), where the new term eilY represents a sinusoidal disturbance in the transverse

540

J.C.-T. Kao et al. / Journal of Colloid and Interface Science 303 (2006) 532–545

(a)

Fig. 11. Neutral stability curves for spatially synchronous (black) and subharmonic (gray) disturbances to equilibrium films. Line styles indicate branch of the equilibrium state, as in Fig. 7b.

(b) Fig. 10. (a) Rupture location and (b) rupture time versus patterning wavenumber. Dot shape indicates initial conditions: flat-film (square), upper-branch (round), or lower-branch (triangular). Dot color indicates the amount of noise in the initial conditions: 10% (black) or 0.1% (gray).

direction. Then h and ω are solutions of the equation   F 2 ¯3 2 H¯ 3 4 2 ¯ 3 ¯ D h + H H Dq h + − H l Dq2 h 3 q 3 H¯     ¯ F FH + 2 − 5 2 − H¯ 2 H¯  l 2 Dq h ¯ H H¯   F F H¯  F H¯  2 F  H¯  + −4 +8 −5 2 2 H¯ H¯ H¯ H¯ 3   1 ¯3 2 F 2 + H l − l + ω h = 0, 3 H¯

(35)

∂ ∂ + iq and (·) = ∂X (·). We solve this numerically Dq = ∂X as the eigensystem of the finite difference matrix representing (35), and obtain the growth rates ω and eigenmodes h(X)eiqX+ilY of each stationary state. It is found that for all values of disturbance wavenumber kd , there is a transverse instability, spatially synchronous with the patterning. However, the spatially subharmonic disturbance is stable for sufficiently short-wavelength patterning, k > kc . This is shown in Fig. 11, which displays the neutral stability curves in kd and l for disturbances to the three stationary branches. For a given kd , transverse instability occurs for wavenumbers l below the neutral-stability curves at lc . The lower bound on lc is √ 3, as for kc on a uniform substrate. Although not shown in the

figure, disturbances with kd = kp /n, n > 2, are also decaying at short wavelengths. The corresponding growth rates are shown in Fig. 12. For fixed transverse disturbance wavenumber l = 1 (Fig. 12a), all three branches are unstable even at large kd , while at l = 2 (Fig. 12b), the small-amplitude middle branch is always stable and the upper branch is stable for kd < kc . For fixed longitudinal disturbance wavenumber kd = 1.8, just below the nose kc of the imperfect bifurcation (Fig. 12c), the upper branch is unstable to synchronous disturbances with sufficiently long transverse wavelength, while the subharmonic disturbance of the middle branch is just passing into stability. Subharmonic disturbances of the other two branches are extremely unstable and do not fit on the graph. Finally, at kd = 2.1 (Fig. 12d), synchronous disturbances to all three branches display long-wavelength instability in the transverse direction, while subharmonic disturbances of the middle branch are stable. Again, subharmonic disturbances of the other two branches are extremely unstable. 6. Heterogeneous substrate: square-wave patterning in (1 + 1)D For an engineered patterned substrate, a piecewise-constant Hamaker parameter would be a better approximation than a sinusoidal Hamaker parameter. Consider ϕ = ϕ0 (H )F (ξ ), with ξ = kp X, and ⎧ π 3 ⎨1 +  , 0  ξ  2 , π F (ξ ) = 1 −  3 , 2  ξ  3π 2 , ⎩  ξ  2π, 1 +  3 , 3π 2

(36)

(37)

where the mean Hamaker parameter is normalized to 1 and  3 is the amplitude of the patterning.

J.C.-T. Kao et al. / Journal of Colloid and Interface Science 303 (2006) 532–545

(a)

(b)

(c)

(d)

541

Fig. 12. Growth rates of transverse perturbations to the steady state. (a) Transverse wavenumber is fixed at l = 1, (b) l = 2, (c) longitudinal wavenumber is fixed at kd = 1.8, (d) kd = 2.1. Line color indicates synchronous (kd = kp , black) and subharmonic (kd = 2kp , gray), while line style indicates the branch of the stationary state as in Fig. 7b.

The patterning amplitude  3 may be found from experimental observations. For two patterned-substrate patches whose fastest growing wavelengths in the homogeneous case would be λ∗+ and λ∗− , assuming the form (10) for the potential gives 3 =

∗2 λ∗2 + − λ− ∗2 λ∗2 + + λ−

(38)

,

and hence λc = (1 +

√ 3 19 2 /4)−1 λc is known.

6.1. Stationary states For films on uniform substrates, trivial bifurcation branches exist at integer fractions of the cutoff wavenumber, kc /n, corresponding to solutions with periodicity n in a single nominal wavelength. These remain trivial in the case of purely sinusoidal patterning of the substrate. However, for more general patterning, e.g., square-wave, the higher harmonics of the patterning shape lead to imperfections of the previously trivial branches. For a patterning of the form F (X) =

∞

cm cos(mkp X),

(39)

Fig. 13. Stationary state film amplitude vs wavenumber (kp ) for film on square-wave patterned substrate with  3 = 18 . Dots represent numerical solution of the thin-film equation (9) in steady state. Grey lines represent film amplitudes from Eq. (40), around the five largest bifurcation points.

the sinusoidally patterned substrate, as shown in Fig. 6. For simplicity, we will neglect these small-kp branches for the remainder of this paper and consider only kp  0.882, above the n = 2 branch.

m=1

an analysis near the nth bifurcation point, kp = (1 +  2 k2 )kc /n, gives H1 = Q1 cos(nkp X) with the amplitude Q1 satisfying 19Q31

− 12Q1 k2 − 2cn = 0

(40)

in steady state. A comparison of this estimate and numerical results is shown in Fig. 13. Note, though, that these branches occur in long-wavelength domains, and so are extremely unstable and not physically significant. At more relevant wavenumbers, the bifurcation diagram and stationary states resemble those of

6.2. Stability of branches As before, stability is examined via Eq. (28). Fig. 14a shows the growth rates of disturbances to the film, versus the wavenumber kd of the disturbance, with line styles corresponding to stationary branches in Fig. 14b. The results here are essentially the same as in the case of sinusoidal pa√ √ Hamaker rameter. Perturbations longer than λc = 2π/ 3(1 + 3 19 2 /4) will grow, causing rupture at every nth stripe, with n depending on kp , kc , and kd .

542

J.C.-T. Kao et al. / Journal of Colloid and Interface Science 303 (2006) 532–545

(a)

(b) Fig. 14. Square-wave Hamaker patterning,  3 = 18 . (a) Growth rates for periodic disturbances to (1 + 1)D stationary films versus the disturbance wavenumber kd . Curves are shown for patterning at various multiples of the disturbance wavenumber: harmonic (kp = kd ⇔ q = 0, in black), subharmonic (kp = 2kd , or q = kp /2, in dark gray), and also kp = 10kd (q = kp /10, light gray). Line styles in (a) indicate stationary states in (b): small-amplitude middle branch (thick solid lines), upper and lower branches (dashed), and dispersion relation for the flat film on homogeneous substrate (thin solid line).

6.3. Rupture

edges of a stripe to rupture at the center of a stripe, consistent with Fig. 4 of Kargupta et al. [23] (abscissa is inverted). For larger kp (shorter wavelength patterning), the flat film ruptures at the middle of the region of larger Hamaker parameter, as in the sinusoidal patterning case. The stationary state shows only slight rupture location preference until kp ≈ kc . Rupture times show that the flat film is far from equilibrium —its rupture time is roughly independent of the amplitude of the initial disturbance, whereas the stationary state initial conditions show a clear rupture time dependence on disturbance size. For long patterning wavelengths, the flat-film rupture time does not rise as in the sinusoidal patterning and near-equilibrium cases. This is because the flat film remains far from equilibrium at the edges of the square wave, whereas in the sinusoidal case, the flat film is increasingly close to stationary for long wavelengths. However, at spinodal wavelengths, rupture times with flat initial conditions on square-wave patterning are comparable to rupture times on sinusoidal patterning, further supporting the idea that despite the preferred rupture location, the accelerated rupture observed here and in simulations by Kargupta et al. [22–29] may be attributed to far-from-equilibrium initial conditions rather than the gradient in the potential as previously suggested. 7. Heterogeneous substrate: square-wave patterning in (2 + 1)D The stability of (2 + 1)D films on a substrate with onedimensional square-wave patterning is examined via Eq. (35), as in the sinusoidal patterning case. Fig. 17 shows the neutral stability curves in kd and l, with transverse instability of the stationary states occurring below each curve. Growth rates are shown in Fig. 18. As with sinusoidal patterning, all branches are subject to some transverse instability though the smallamplitude middle branch is stable with respect to spatially subharmonic disturbances. 8. Conclusion and discussion

For the case of square-wave patterning, we compute rupture time numerically only. Some sample time evolutions are shown in Fig. 15, and rupture time and locations plotted in Fig. 16. Unlike the sinusoidal patterning case, here there are welldefined preferred rupture locations for given initial conditions of both flat films and near-equilibrium films. Examination of film evolution suggests that for unstable films, the square-wave shape of the Hamaker parameter leads to greatest disequilibrium (and hence rupture alignment) near the pattern’s edges, provided the wavelength of patterning is long enough. This is in contrast to sinusoidal patterning, where the strength of the Hamaker parameter always favors rupture on the centerline. Furthermore, when kp  kc /2, two wavelengths of the instability fit to a single wavelength of the patterning, leading to near-simultaneous rupture at the edges of the pattern. This geometric argument may explain the different rupture locations and resulting isolated droplets seen in previous experiments [18,19, 21] and simulations [23,24,26–29]. The bump in rupture time near kp ≈ 0.8 corresponds to the switch from rupture at the

We have studied systematically the effects of striped (1D) substrate patterning on the stability and rupture of thin liquid films subject to van der Waals forces, with a focus on resonant patterning. We have shown both analytically and numerically that solutions of the evolution equation (9) are governed by an imperfect pitchfork bifurcation in the patterning wavelength (Fig. 6). Thus an inhomogeneous Hamaker parameter with O( 3 ) variation results in an O( 2 ) change in the stability threshold from kc to kc and an O() change in the stationary states, destabilizing the film as compared to films on equivalent homogeneous substrates. The bifurcation is insensitive to the shape of the patterning (sinusoidal versus square-wave variation of the Hamaker parameter), though higher harmonics in the patterning can induce imperfections in the small-kp bifurcation branches that were trivial for sinusoidal patterning and for homogeneous substrates (Fig. 13). We have investigated the growth rates of small disturbances to the stationary states (Figs. 7 and 14). We have confirmed

J.C.-T. Kao et al. / Journal of Colloid and Interface Science 303 (2006) 532–545

543

(a)

(b)

(c) Fig. 15. Typical time evolution (gray lines at T = 0.1) of initial conditions (thin black line), with square-wave patterning (dashed line). (a) Initial conditions near lower branch, (b) initial conditions near upper branch, (c) flat initial conditions. Film shape at touchdown or equilibrium is shown (thick black line).

the stability of stationary films on the bifurcation branches as determined by weakly nonlinear analysis, and further found that regardless of patterning wavelength, rupture tends to occur on a length scale comparable to λ∗ , corresponding to the fastest-growing mode on a homogeneous substrate. This is in accordance with Higgins and Jones [20], Sehgal et al. [21], and Kargupta et al. [24,26–29]. Thus, while one may control the location and alignment of rupture via substrate patterning alone, control of the length scale of rupture requires manipulation of the overall liquid or substrate properties. We have found that on striped substrates, transverse (3D) disturbances to the stationary state obey the linear dispersion relation for films on homogeneous substrates and cause rupture on sufficiently wide domains. Therefore, real-world film templating requires either fully two-dimensional patterning or transverse bounding of the film in order to avoid uncontrolled rupture. We have also performed numerical simulations of the evolution equation (9) and recorded rupture times and locations for a range of initial conditions (Figs. 10 and 16). Although bifurcation curves and rupture times are only slightly affected by the shape of the patterning, we observed a strong rupturelocation dependence on the patterning shape and size, particularly for near-flat initial conditions. Our results on rupture location within a stripe confirm findings of Sehgal et al. [21]

and Zhang et al. [18,19]. In addition we have found both analytically and numerically that accelerated rupture on patterned substrates may be due to initial conditions that are very far from the equilibrium states, rather than a new instability mechanism as proposed by Kargupta, Konnur, and Sharma [22,23]. Many of these results reproduce previously known experimental and numerical findings, including: appropriate patterning length scales for templating [21,24,26–29], rupture location transition [18,19,21,23,24,26–29], accelerated rupture [22–29], and bifurcation [39,40]. However, we have explained some of these phenomena differently than in previous work. The significance of this work is that the results here are derived by more transparent methods and placed in the unified theoretical setting of imperfect bifurcation. Thus we can describe the effects of patterning in terms of three regimes: (1) kp < kc , (2) kc < kp < kc , and (3) kc < kp , where kp is the wavenumber of both the patterning and the bounded or periodic film, kc is the neutral stability threshold for homogeneous substrates, and kc is the modified threshold due to substrate patterning. In region (1), films on both uniform and patterned substrates are unstable. However, near-flat films on patterned substrates display accelerated rupture due to a nonflat equilibrium state. In (2), we find that patterning causes previously stable films to become unstable—the region of instability expands because of imperfect bifurcation. The equilibrium state in this region

544

J.C.-T. Kao et al. / Journal of Colloid and Interface Science 303 (2006) 532–545

(a)

Fig. 17. Neutral stability curves for spatially synchronous (black) and subharmonic (gray) disturbances to equilibrium films. Line styles indicate branch of the equilibrium state, as in Fig. 7b.

(b) Fig. 16. (a) Rupture location and (b) rupture time versus patterning wavenumber. Dot shape indicates initial conditions: flat-film (square), upper-branch (round), or lower-branch (triangular). Dot color indicates the amount of noise in the initial conditions: 10% (black) or 0.1% (gray).

is also strongly affected by patterning, resulting in quick rupture for previously stable films. Finally, films in (3) are stable regardless, but patterned substrates induce a modulation in the equilibrium state. For substrate patternings where Amin and Amax both generate rupture and (Amax − Amin ) is small, the values of kc and kc are known analytically. For other cases, outside the validity of the perturbation analysis, we nonetheless expect the qualitative behavior to remain. In particular, it may be that while a

(a)

(c)

(b)

(d)

Fig. 18. Growth rates of transverse perturbations to the steady state. (a) Transverse wavenumber is fixed at l = 1, (b) l = 2, (c) longitudinal wavenumber is fixed at kd = 1.8, (d) kd = 2.1. Line color indicates synchronous (kd = kp , black) and subharmonic (kd = 2kp , gray), while line style indicates the branch of the stationary state as in Fig. 14a.

J.C.-T. Kao et al. / Journal of Colloid and Interface Science 303 (2006) 532–545

homogeneous system is spinodally stable at all wavenumbers (kc2 < 0), the corresponding heterogeneous kc exists and is positive, leading to instability at sufficiently long wavelengths. We believe that this is the case in the simulations by Kargupta et al. which exhibit rupture of “spinodally stable” films on substrates with a strong inhomogeneity [22,23,25,28]. Similarly, for films not constrained to have wavenumber kp , as on larger domains where rupture occurs on a spinodal length scale regardless of patterning, rupture time continues to be governed by distance from an equilibrium which can be strongly affected by patterning. We believe this accounts for the frequently observed accelerated rupture of films on heterogeneous substrates. Finally, we explain the rupture locations and preferred length scales seen in previous work [18–24,26–29] in terms of the growth rates for periodic disturbances as shown in Figs. 7 and 14. The varying location of rupture within a stripe is a geometric consequence of being able to fit multiple wavelengths of the instability in a single period of the patterning. That films on patterned substrates continue to rupture at spinodal wavelengths whenever possible is due to the fact that although patterning changes the stability threshold kc to kc , the wavenumber k ∗ of the fastest-growing mode remains the same. These findings should prove useful in the interpretation of existing and future thin-film experiments and simulations. Other factors that remain to be considered include the effects of large patterning amplitude, the ratio of stripe width to patterning period, more general functional forms of the intermolecular potential, nonperiodic films on longer domains, and postrupture film evolution and coarsening. Acknowledgments J.C.T.K. gratefully acknowledges support of the NSF-IGERT fellowship program, #DGE-9987577. A.A.G. acknowledges the support of the US Department of Energy, Grant #DE-FG0203ER46069. We thank the reviewers of this paper for many helpful comments and suggestions. References [1] [2] [3] [4] [5]

A. Vrij, Discuss. Faraday Soc. 42 (1966) 23–33. A. Sheludko, Adv. Colloid Interface Sci. 1 (1967) 391–463. G. Reiter, Phys. Rev. Lett. 68 (1) (1992) 75–78. T.G. Stange, D.F. Evans, Langmuir 13 (16) (1997) 4459–4465. R. Xie, A. Karim, J.F. Douglas, C.C. Han, R.A. Weiss, Phys. Rev. Lett. 81 (6) (1998) 1251–1254.

545

[6] E. Ruckenstein, R.K. Jain, J. Chem. Soc. Faraday Trans. 70 (1974) 132– 147. [7] M.B. Williams, S.H. Davis, J. Colloid Interface Sci. 90 (1) (1982) 220– 228. [8] T. Erneux, S.H. Davis, Phys. Fluids A 5 (5) (1993) 1117–1122. [9] A. Oron, S.G. Bankoff, J. Colloid Interface Sci. 218 (1999) 152–166. [10] B.Y. Rubinstein, A.M. Leshansky, Langmuir 16 (2000) 2049–2051. [11] A.M. Leshansky, B.Y. Rubinstein, Phys. Rev. E 71 (2005) 040601. [12] M.P. Ida, M.J. Miksis, Appl. Math. Lett. 9 (3) (1996) 35–40. [13] J. Becker, G. Grün, R. Seemann, H. Mantz, K. Jacobs, K.R. Mecke, R. Blossey, Nat. Mater. 2 (2002) 59–63. [14] W.W. Zhang, J.R. Lister, Phys. Fluids 11 (9) (1999) 2454–2462. [15] A. Oron, S.H. Davis, S.G. Bankoff, Rev. Mod. Phys. 69 (3) (1997) 931– 980. [16] E. Meyer, H.-G. Braun, Macromol. Mater. Eng. 276 (2000) 44–50. [17] Z. Zhang, Z. Wang, R. Xing, Y. Han, Surf. Sci. 539 (2003) 129–136. [18] Z. Zhang, Z. Wang, R. Xing, Y. Han, Polymer 44 (2003) 3737–3743. [19] C. Luo, R. Xing, Z. Zhang, J. Fu, Y. Han, J. Colloid Interface Sci. 269 (2004) 158–163. [20] A.M. Higgins, R.A.L. Jones, Nature 404 (2000) 476–478. [21] A. Sehgal, V. Ferreiro, J.F. Douglas, E.J. Amis, A. Karim, Langmuir 18 (2002) 7041–7048. [22] R. Konnur, K. Kargupta, A. Sharma, Phys. Rev. Lett. 84 (5) (2000) 931– 934. [23] K. Kargupta, R. Konnur, A. Sharma, Langmuir 16 (2000) 10243–10253. [24] K. Kargupta, A. Sharma, Phys. Rev. Lett. 86 (20) (2001) 4536–4539. [25] K. Kargupta, A. Sharma, J. Colloid Interface Sci. 245 (2002) 99–115. [26] K. Kargupta, A. Sharma, J. Chem. Phys. 116 (7) (2002) 3042–3051. [27] K. Kargupta, A. Sharma, Langmuir 18 (2002) 1893–1903. [28] A. Sharma, R. Konnur, K. Kargupta, Physica A 318 (2003) 262–278. [29] K. Kargupta, A. Sharma, Langmuir 19 (2003) 5153–5163. [30] D.E. Kataoka, S.M. Troian, Nature 402 (1999) 794–797. [31] L. Kondic, J. Diez, Phys. Rev. E 65 (4) (2002) 045301. [32] L. Kondic, J. Diez, Colloids Surf. 214 (2003) 1–11. [33] L. Kondic, J. Diez, Phys. Fluids 16 (9) (2004) 3341–3360. [34] H. Gau, S. Herminghaus, P. Lenz, R. Lipowsky, Science 283 (1999) 46– 49. [35] A.A. Darhuber, S.M. Troian, S.M. Miller, J. Appl. Phys. 87 (11) (2000) 7768–7775. [36] A.A. Darhuber, S.M. Troian, Annu. Rev. Fluid Mech. 37 (2005) 425–455. [37] L. Rockford, Y. Liu, P. Mansky, T.P. Russell, Phys. Rev. Lett. 82 (12) (1999) 2602–2605. [38] N. Rehse, C. Wang, M. Hund, M. Geoghegan, R. Magerle, G. Krausch, Eur. Phys. J. E 4 (2001) 69–76. [39] L. Brusch, H. Kühne, U. Thiele, M. Bär, Phys. Rev. E 66 (2002) 011602. [40] U. Thiele, L. Brusch, M. Bestehorn, M. Bär, Eur. Phys. J. E 11 (2003) 255–271. [41] J.N. Israelachvili, Intermolecular and Surface Forces, Academic Press, 1992. [42] A. Sharma, Langmuir 9 (3) (1993) 861–869. [43] D.W. Jordan, P. Smith, Nonlinear Ordinary Differential Equations, Oxford Univ. Press, 1999. [44] B.J. Matkowsky, E.L. Reiss, SIAM J. Appl. Math. 33 (2) (1977) 230–255. [45] L.F. Shampine, M.W. Reichelt, SIAM J. Sci. Comput. 18 (1997) 1–22.