Emergence of Dewetting Patterns

Before going into the details of thin-film dewetting, we make some general remarks and simple observations. Assume that the shape of the fluid volume Ω(t)⊂R3 over a substrate at z=0 and time t is described by a function ht,x such thatΩ(t)={(x,z)∈R2×R+:0<z<h(t,x)}.[1]The wetted substrate area is ω(t)={x∈R2:h(t,x)>0}. Then, the surface energy driving the evolution is defined asE(h)=∫ω(γℓs−γgs)+γℓg1+|∇h|2 dx,with γ{ℓs,gs,ℓg} being the surface tensions of the liquid–solid, gas–solid, and liquid–gas interfaces, respectively. For incompressible fluids, the evolution of the film height conserves the total volume V(h)=∫ωh dx. For Newtonian fluids with velocity field u(t):Ω→R3, this energy decreases along solutions: that is,ddtE=−D=−∫Ω2μ|Du|2 dx dz−∫ωμb−1|u|z=02 dx≤0,decomposing the dissipation D into a bulk and substrate part with viscosity μ, Navier slip length b, and symmetric gradient Du=12(∇u+∇u⊤). Thereby, stationary states are energy minimizers, and a straight-forward minimization of E with constant V shows that minimizers are droplets (i.e., up to translational invariance ω is a disc, and h has constant curvature). Apparently, this assumes that contact lines can slide freely and that no pinning is present (2). However, it is also clear that arrays of droplets as shown in Fig. 1 are not global but local minimizers of E and thereby, also represent viable stationary states. While droplet patterns emerging from the spinodal dewetting feature seemingly random distributions of droplets with details depending on intermolecular interactions (30, 31), droplet distributions in other flow scenarios are reminiscent of the process that leads to their creation (e.g., heterogeneous nucleation, surface instabilities on patterned substrates, droplet production in confined environments, or even dynamically for sliding droplets) (32⇓⇓⇓⇓⇓–38). This raises the natural question: Which mechanism decides what pattern is generated?

This question is even more justified when realizing that the surface energies E leading to both patterns in Fig. 1 are qualitatively identical. Therefore, in the following, we show that physical systems with the same wetting energy but different magnitudes of dissipation (i.e., viscosity and Navier slip) will produce qualitatively different droplet patterns. While the exploration of the nonlinear pattern formation will mainly compare experiments and simulations, the general mechanism responsible for different patterns is the switching of different instabilities due to the dissipation (Movie S1).

Thin-Film Models and Instability

Problem Formulation.

After a film has ruptured by nucleation or by external forcing, forming a hole or in a planar-symmetric setting, a trench, the viscous fluid retracts to reduce the overall energy of the liquid–gas, solid–liquid, and solid–gas interfaces. The dewetting process is driven by the intermolecular potential ϕ between the film and the substrate. In the simplest case, it consists of a sum of attractive long-range van der Waals forces and short-range Born repulsion forces, the minimum of which yields the height h⋆ that remains in the dry region from which the film of uniform height h∞≫h⋆ has dewetted. Motivated by Lennard–Jones potentials, one often finds intermolecular potentials (31, 39) of the standard form ϕ(h)=ϕ¯h/h⋆, whereϕ¯(h)=Sn−mnh−m−mh−n,[2]with m=2,n=8, ϕ′(h⋆)=0, and ϕ(h⋆)=S. For partial wetting, the spreading coefficient S=γgs−γℓs−γℓg is negative.

Due to the slow dewetting rates of the polymer films with chain lengths below the entanglement length, the Navier–Stokes equations serve as the underlying model for the viscous fluid with a Navier slip boundary condition,t⋅(2μDu)n+μb−1t⋅u=0,[3]for tangential velocity implied by the previous energy dissipation balance. The scale separation of the characteristic height scale H=[z] and typical length scales L=[x] allows for a consistent thin-film approximation using the small parameterε=HL,[4]which leads to the reduction of the Stokes free boundary problem to a problem for the free boundary h(t,x) in closed form. In the following, all quantities are nondimensionalized with this scaling as explained along with details of derivation in ref. 22, which we here only summarize. In particular, we introduce the nondimensional slip length B=b/H.

We seek the scalar function h:[0,T]×R2→[0,∞) as in Eq. 1, where the film height h(t,x) depends on time t∈[0,T] and on space x=(x,y)∈R2 with given initial data h0(x)=h(t = 0,x) (Fig. 2). It has been shown that, depending on the magnitude of slip length B, there exist two asymptotic distinguished limits (22) for ε≪1. The first distinguished limit assumes B=O(1) and leads to the weak-slip thin-film model given by the fourth-order parabolic partial differential equation (PDE)∂th−∇⋅(m(h)∇π)=0,[5a]where the mobility ism(h)=13h3+Bh2,[5b]and the generalized pressure π is defined as the functional derivative of the energy functionalE(h)=∫R212|∇h|2+ϕ(h) dx[5c]with respect to h: that is,π=δEδh=−∇2h+Π(h),[5d]and Π(h)=∂hϕ(h) is the derivative of the nondimensional intermolecular potential. Effectively, by introducing Π, one enforces positivity of h:[0,T]×ω¯→[0,∞) and replaces the complement of the set ω(t)⊂ω¯⊂R2 by the region where h≈h⋆ as it is sketched in Fig. 2.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

Sketch of solution h(t,x) in the x−z plane.

The second distinguished limit assumes that B=ε−2β and leads to the strong-slip thin-film model given by the system of PDEs for the film height h:[0,T]×ω¯→[0,∞) and the lateral velocity u¯:[0,T]×ω¯→R2:∂th+∇⋅hu¯ = 0,[6a]Re∂tu¯+(u¯⋅∇)u¯ = 1h∇⋅σ−∇π−β−1u¯h,[6b]where Re is the rescaled Reynolds number and the effective shear stress isσ=h∇u¯+(∇u¯)⊤+2(∇⋅u¯)I.[6c]Two important limiting cases are the no-slip and the intermediate-slip thin-film models. The no-slip model is obtained for B→0 so that the degenerate mobility becomes m(h)=13h3. The latter intermediate-slip model is obtained in the limit B→∞ but ε−2B→0, which on rescaling of time, produces the mobility m(h)=h2. In this sense, no-slip m=13h3 mobility and intermediate-slip mobility m=h2 are two extremal limits in the slip length B.

It has been shown experimentally and near the linear regime it was also demonstrated theoretically that details of the dynamic and morphological evolution strongly depend on the magnitude of slip at the solid–liquid interface (40). For no slip or weak slip, the liquid accumulates in a growing rim in front of the contact line and destabilizes. The case of strong slip is characterized by very asymmetric retracting rims with a monotone spatial decay toward the unperturbed film for particularly large slip, which undergoes a transition to oscillatory decay as the slip length decreases below a critical value (22, 41). In fact, this observation has been used to identify the strong-slip regime in experiments and also to determine the slip length quantitatively (42⇓–44). The molecular origin of slip in polymer melts was investigated in ref. 45. More details about this topic can be found in the work by Bäumchen and Jacobs (46).

Experimental Setup and Methods

For the experiments, atactic polystyrene (PS; purchased from PSS; molecular weights are as listed in the experiments) is used as a model viscous liquid. The films were produced by spin casting a toluene solution (Selectipur or LiChrosolv; Merck) of PS on freshly cleaved mica sheets. The glassy thin films were then floated onto an ultrapure water (organic impurities of <6 ppb, resistance at 25○C: <18.2 M Ω cm) surface and were then picked up with hydrophobic Si wafers.

Hydrophobic Si wafers were achieved by two different preparation methods: (i) on the cleaned Si surface, a self-assembled monolayer of silane molecules [dodecyltrichlorosilane (DTS); Sigma Aldrich/Merck] was prepared (52), or (ii) the cleaned Si wafer was dipped into a solution of a fluoropolymer layer (AF1600; Sigma Aldrich/Merck).

Dewetting is initiated by heating the glassy polymer film above its glass transition temperature. The dewetting of the retracting straight fronts was monitored in situ by optical microscopy on a heating plate (Linkam) or by atomic force microscopy (AFM; Dimension ICON; Bruker).

The dewetted distance was typically obtained from optical micrographs. In AFM experiments, the dewetted distance can also be calculated from 3D scans of the rim on the basis of volume preservation. The values resulting from both approaches were checked for consistency.

Slip lengths have been calculated using the rim profile analysis method (53, 54); structural details, surface roughness values, and wetting properties of the coatings are given in the supplementary material of ref. 40. Polymer film thicknesses have been determined by ellipsometry on the glassy film or by AFM on the edge of a film. The equilibrium contact angles are 𝜗=(88±2)○ and 𝜗=(66±2)○ for AF1600- and DTS-coated substrates, respectively. Comparing typical values for rim height and width, we find typical values in the range ε=0.15−0.3 for the dewetting systems.

In Fig. 3, we show retracting rims with four different chain lengths of the PS film, each thereby corresponding to drastically different slip length b. For the smallest chain length PS(65k), the slip length b is smaller than the rim height H and hence, should correspond to a no-slip situation. As chain length and correspondingly, slip length b increase, we observe a more pronounced instability, as we would expect for a transition from a no-slip scenario to an intermediate-slip scenario. However, if b is increased even further by two orders of magnitudes for PS(390k), the instability is suppressed again. Interestingly, this was predicted for the transition to a strong-slip regime.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 3.

Comparisons of the dewetting behavior of h∞=110-nm-thick PS films for different chain lengths dewetting from AF1600-coated substrate after the rim traveled the distance given above the images. Films are annealed to different temperatures to speed up dewetting in case of films with larger viscosities: PS(65k) at T=140○C, PS(101k) at T=140○C, PS(186k) at T=150○C, and PS(390k) at T=150○C. In the chosen systems, the slip length increases with increasing molecular weight, leading to the values indicated next to the corresponding images obtained by fitting the rim profile as reported by Bäumchen et al. (45) or McGraw et al. (55).

For the remaining part of the paper, we will only consider experiments that can be captured with a no-slip (AF1600 coating) or intermediate-slip (DTS coating) regime with PS(13.7k) or PS(10.3k) and use corresponding theoretical models with cubic m=13 h3 or quadratic m=h2 mobility.

Numerical Methods

We now explain the intricacies of numerically resolving the multiple length scales of the problem. We focus on the two cases of no slip and intermediate slip, as most of the experimental results regarding the various instabilities fall into either of these two regimes.

Our solution of the thin-film model in Eq. 5 is based on a P2 finite element method (FEM), where the fourth-order equation is split into a system of two second-order equations. The FEM uses piecewise quadratic elements and local mesh refinement. We use a semiimplicit time discretization, where only the highest-order derivative is treated implicitly. Therefore, the thin-film equation (Eq. 5) is multiplied with a test function v and integrated by parts to obtain∫ω¯∂thv+m(h)∇π⋅∇v dx=0,[10a]where boundary terms vanish due to the no-flux boundary condition n⋅∇π=0 imposed on ∂ω¯. We also rewrite the pressure π in the weak form as∫ω¯πv dx=∫ω¯∇h⋅∇v+Π(h)v dx,[10b]where again, we used integration by parts and n⋅∇h=0 on the boundary ∂ω¯. In particular, this statement of the PDE implies conservation of volumeddt∫ω¯h(t,x) dx=0,which can be seen when selecting v=1 in the continuous or discrete weak formulation (Eq. 10a). Evaluating the solution at discrete times hn(x)=h(nτ,x), we use the time discretization ∂th=τ−1(hn−hn−1). This allows us to rewrite the weak form of the thin-film model in Eq. 10 so that we seek (hn,π)∈W such that∫ω¯hnv+τ m* ∇π⋅∇v dx=∫ω¯hn−1v dx,[11a]∫ω¯πw−∇hn⋅∇w dx=∫ω¯Π*w dx,[11b]which needs to hold for all (v,w) from a suitable function space W and with initial data h0. When we define the mobility m*=m(hn−1) and Π*=Π(hn−1), this becomes a semiimplicit time discretization.

The FEM constitutes a method for the construction of a finite-dimensional subspace Wh⊂W, where we have an admissible decomposition ω¯=⋃k=1NTk of the domain into triangles Tk, on which we define Wh as the space of continuous functions that are piecewise quadratic on each triangle (i.e., P2 finite elements). Then, we seek a discrete solution (hn,π)∈Wh of Eq. 11 valid for all (v,w)∈Wh. The integrals appearing in Eq. 11 are solved exactly or using a seven-point Gauss quadrature. A typical droplet pattern emerging from the simulation with m=h2 is shown in Figs. 4 and 5.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 4.

Simulated pattern formation with intermediate-slip m=h2: 3D view with light areas showing dry regions, whereas elevated dark areas show the liquid–air interface.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 5.

(Left) Dewetting rim during pinch-off of a single droplet highlighted with dashed lines and (Right) magnification of this droplet and corresponding locally refined triangulation. The mesh consists of 92, 272 vertices, which require us to solve for 2×368, 423=736, 846 unknowns (hn,π) for the P2 FEM discretization at each time step.

The precursor thickness h⋆ in the potential of Eq. 2 needs to be chosen much smaller than any droplet size we intend to resolve. Then, the solution h will feature large regions, where the solution is either almost constant h≈h⋆ or the solution is smooth and h≫h⋆ as shown in Fig. 2. However, where those regions meet, the solution features a kink, which should be resolved in the triangulation (compare with Fig. 5). We perform a heuristic local mesh refinement in these connecting regions where the contact line is situated. When constructing a new mesh, we start from a coarse base mesh and determine which of its elements are crossed by the contact line of the previous solution. Those elements are refined by inserting additional vertices to an extent that the kink is resolved again. Furthermore, we make sure that neighboring elements are also refined, so that the contact line remains in the refined region of the mesh for a number of time steps. Based on the set of newly created vertices, we perform a new triangulation and interpolate the old solution onto the new triangulation.

To allow the control of the minimum and the derivatives of ϕ separately via h⋆ and ε⋆, it is advantageous to work with an alternative potential representation of the form ϕ(h)=ϕ^(h−h⋆)/ε⋆ withϕ^(s)=12γ1+s−(1+γ)exp(−s2),[12]which has similar properties as the rescaled version of Eq. 2. In particular, with this potential, we still have the same Γ convergence property as h⋆,ϵ⋆→0. This abstract statement ensures that equilibrium contact angles are maintained and that the energetic contributions of ϕ from Eq. 12 and Eq. 2 coincide in this limit. We note that, for 0<γ≪1, the minimum slightly shifts away from h⋆; however, one gains a slightly stabilizing potential with ϕ″ > 0 for h≫h⋆. In addition, we point out that we monitor that the minimum of the solution minxh(t,x) never violates the nonnegativity requirement for the given choice of ϕ(h) and parameters h⋆,ε⋆ (Fig. 6), in particular, by keeping the time step size τ sufficiently small. We also refer to previous mathematical studies (56, 57) that incorporate the nonnegativity property into their numerical scheme.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 6.

Lennard–Jones potential Eq. 2 compared with the much shorter-ranged exponential potential Eq. 12, both with h⋆=1/20 and ε⋆=1/16. Note that, in the numerics we rescaled so that at the contact line the slope is |∇h|=1, leading to a nondimensional potential with ϕ(h⋆)=−1/2.

Initial Data for Rims

We still need to specify the initial data h0(x) for the numerical simulation to describe the various stages of the dewetting process. In experiments, the dewetting is initiated from a uniform flat layer bounded by a nearly straight edge. When the sample is heated, the layer liquefies, and the edge becomes a moving contact line of the dewetting process. In simulations, we choose the supporting domain sufficiently large ω¯=[0,500]×[0,600]; using h⋆=1/20, we represent the uniform layer with the nearly straight contact line at x0 with initial heighth0(x)=h⋆+12(h∞−h⋆)1+tanhx−x0(y)Δ,[13]and the smooth initial contact line position is represented by x0(y)=20+∑n=150an⁡cos(nπy/600). With Δ=1/2, one can interpret h0 to be an approximation of a step-like profile, where slight corrugations of x0 along the y direction are introduced using the normally distributed coefficients an with zero mean and SD 1/10. This choice is then maintained throughout all dewetting simulations for both mobilities. Due to the introduced scaling, we have h∞=1.

To study the systematic effect of the slip boundary condition on dewetting patterns, we study the contact line instability of moving rims and stationary ridges for the mobility functions m(h)=13h3 (no slip) and m(h)=h2 (intermediate slip). In Fig. 7, the evolution of the rims for no slip and intermediate slip is shown at times where they have accumulated a similar volume, starting with the same initial film thickness h∞ as given by the experiments. This corresponds to the experimental results in Fig. 8, where the evolution of the rims is shown at the same distance measured from the initial coordinates. The distinctive dewetting patterns for the intermediate-slip case are observed, while for the no-slip case, only shallow oscillations will occur in excellent comparison with the experimental results in Fig. 8.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 7.

Numerical solutions on ω¯=[0,500]×[0,600] (h⋆=1/20, γ=1/10) for (Upper) m=h2 (intermediate-slip B→∞) and (Lower) m=13h3 (no-slip B=0) for initially slightly perturbed straight front and time progressing from left to right shown for similar rim progression (Movie S1).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 8.

Series of micrographs from experiments with thin films (Upper) dewetting from a hydrophobized Si wafer covered with a silane monolayer DTS with b=1 μm (intermediate slip) and (Lower) dewetting from an AF1600-covered Si wafer with b=40 nm (no slip). In both series, a h∞=100 nm thin PS(13.7k) film dewets at T=120○C. In both series, undulations are formed along the rim; however, only in the intermediate-slip case above, budding is observed that later leads to fingers and a pinch-off of droplets. Comparing rims that have traveled a similar distance ensures that only rims of similar volumes are evaluated. Reprinted with permission from ref. 19. Copyright 2015 by the American Physical Society.

Breakup of Liquid Ridges

Toward the final stages of dewetting process, for long times and domain size far beyond ω¯, the rims approach each other and merge, yielding a polygonal network of ridges. As is known from previous studies, stationary ridges are susceptible to a Rayleigh–Plateau instability for both no-slip and intermediate-slip boundary conditions at the substrate. Both the no-slip case and the intermediate-slip case have been extensively investigated in the literature (20, 24, 26).

To study the instability of a capillary ridge, we first compute a 1D stationary solution hstat(x) (i.e., a time-independent 2D solution of Eq. 5 with ∂th≡0, ∂yh≡0, and constant thickness h→h⋆ as x→±∞ as shown in Fig. 9, with a symmetric and approximately parabolic shape around the maximum). Then, we study the linear stability with respect to perturbationsh(t,x)=hstat(x)+δh1(x;k)exp(iky+σt)[14]with fixed wavenumber k, which returns an eigenproblem for the perturbation h1 of the base state (20, 51). We select the most unstable mode kmax and define as initial datah0(x)=hstat(x)+δh1(x;k)cos(ky)[15]with a sufficiently small δ so that the nonnegativity of h0 is also not violated. In SI Appendix , SI Text and Figs. S1–S3, we present the details of the linear stability analysis and show that solutions of the nonlinear problem follow the predicted exponential amplification while the linear regime is valid. Furthermore, we checked that the alternative choice of potential Eq. 12 compared with the standard choice Eq. 2 has no visible impact on the linear stability as long as the minimal value ϕ(h⋆) and the minimal height h⋆ are the same.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 9.

The 1D steady-state hstat with ε⋆=1/20 and h⋆=1/80.

However, we are interested in using the numerical simulation as a tool to observe features in the nonlinear regime that are distinctive features of the mobility law. Therefore, to resolve finer features of the ridge shape (e.g., secondary or tertiary droplets), the general form of the intermolecular potential is the same, but we use a smaller h⋆=1/80.

We now consider nonlinear simulations using the stationary ridge with very small monochromatic perturbations as initial data. Specifically, we consider an unstable perturbation with kmax and carry out simulations with the mobilities h2, 13h3.

As long as the perturbations are small enough, the evolution of the corrugations in the nonlinear model closely follows the predictions from the linear stability analysis. After the perturbations become comparable with the size of the rim, nonlinear terms become relevant, and the growth rate changes.

Regarding the morphological evolution of the ridges for the no-slip and intermediate-slip cases, our numerical simulations show qualitatively and quantitatively very similar behavior during the linear regime. However, deep into the nonlinear regime, the breakup into droplets follows different scenarios. In the no-slip case, a cascade of satellite droplets emerges, while for the intermediate-slip case, they disappear. Starting with the same initial condition, the evolution for both cases is shown in Fig. 10. Interestingly, this different behavior is also observed in our experimental results as seen in Fig. 11. The instability mechanisms causing these different patterns are encoded in the behavior near the self-similar pinch-off (e.g., refs. 58 and 59 are related theoretical studies of thread breakup).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 10.

Numerical solutions showing ridges with (Upper) no slip m(h)=13h3 and (Lower) intermediate slip m(h)=h2. The initial data are h0(x)=hstat(x)+δh1(kmax;x)cos(ky), with time increasing from left to right. The wavenumber kmax is the one with the largest amplification from Eq. 14, and the domain ω¯=[−12,0]×[0,π/kmax] is then extended to [−12,12]×[−π/kmax,3π/kmax] using the implied symmetry of the solution (Movie S1).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 11.

Experimental results. Close-up view of the late-stage ridge of PS(10.3k) on AF1600 with time progressing from left to right showing decay into satellite and subsatellite droplets. The initial film height is h∞=115 nm, and b=40 nm. With B=b/H and H corresponding to rim height, this gives rise to the no-slip regime (21).

Conclusion and Discussion

We have introduced a highly adaptive finite element-based numerical approach that correctly captures the complex dewetting process described by a class of thin-film models with degenerate mobilities. We showed that, for the no-slip condition, the droplet pinch-off is absent during the retraction of the rim, while for the intermediate-slip case, self-repeating droplet pinch-off occurs in excellent agreement with experimental observations. The ability to resolve the different length scales for long timescales also enables the prediction of phenomena, such as the formation of satellite droplets, as a function of the mobility. The emergence of satellite droplets is well known during the break up of liquid jets and the related problem of liquid filaments. For the latter problem, destabilization is due to the difference in the axial contribution to the capillary pressure between thicker and thinner parts. In this system, the pressure is higher in the thinner parts and squeezes the liquid into the bulges, thus increasing the perturbation until the filament breaks up. Apart from the huge literature on experimental studies, the problem has sparked numerous numerical and analytical investigations (28, 60, 61). In particular, the work by Tjahjadi et al. (27), where the emergence of satellite droplets was captured numerically, and the highly the accurate numerical schemes developed by Kim et al. (62) improved the understanding the underlying physical processes considerably.

For the situation of a ridge that destabilizes on a solid substrate, the additional influence of the substrate enters the linear stability analysis, in particular through the contact angle dynamics. It will be interesting to investigate analytically the rupture behavior for this problem to help understand the influence of the boundary condition at the substrate. Similarly, it remains an interesting open question how the dissipation due to a dynamic contact angle would affect the whole pattern formation process. Such a study would certainly require even more sophisticated finite element techniques (63, 64).

The decision about the pathway leading to specific droplet patterns is then often decided by the specific influence that the dissipation has on flow instabilities. In this study, these instabilities are the shedding of droplets from moving rims and the symmetry of the Rayleigh–Plateau instability leading to satellite droplets, and they are influenced by the magnitude of the interface dissipation encoded in the Navier slip condition.

We conclude that, by controlling dissipative effects in this dewetting films system, we can steer pattern formation without changing the driving forces.