# Signatures of slip in dewetting polymer films

See allHide authors and affiliations

Edited by David A. Weitz, Harvard University, Cambridge, MA, and approved March 18, 2019 (received for review December 3, 2018)

## Significance

Dewetting is the hydrodynamic process where a uniform layer of liquid destabilizes and decays into distinct patterns of stationary droplets by virtue of interfacial and intermolecular energies. These patterns can now be predicted theoretically, and their evolution can be followed numerically in striking similarity to experimental results. The droplet arrays observed during the quasistatic evolution are generated by two distinct types of instabilities: droplet shedding of moving liquid rims and satellite droplet pinch-off during the decay of stationary ridges. The mechanism causing the emergence of different droplet patterns is the dissipation at the liquid–solid interface modeled by an effective Navier slip.

## Abstract

Thin polymer films on hydrophobic substrates are susceptible to rupture and hole formation. This, in turn, initiates a complex dewetting process, which ultimately leads to characteristic droplet patterns. Experimental and theoretical studies suggest that the type of droplet pattern depends on the specific interfacial condition between the polymer and the substrate. Predicting the morphological evolution over long timescales and on the different length scales involved is a major computational challenge. In this study, a highly adaptive numerical scheme is presented, which allows for following the dewetting process deep into the nonlinear regime of the model equations and captures the complex dynamics, including the shedding of droplets. In addition, our numerical results predict the previously unknown shedding of satellite droplets during the destabilization of liquid ridges that form during the late stages of the dewetting process. While the formation of satellite droplets is well known in the context of elongating fluid filaments and jets, we show here that, for dewetting liquid ridges, this property can be dramatically altered by the interfacial condition between polymer and substrate, namely slip. This work shows how dissipative processes can be used to systematically tune the formation of patterns.

The no-slip condition is widely accepted as an appropriate boundary condition for flows of Newtonian liquids sheared along a solid surface. A notable exception arises in the presence of a moving contact line between a viscous liquid and a rigid solid substrate, where the use of the no-slip condition leads to a nonintegrable singularity in the stress field (1, 2). In the past decades, however, it has been shown that thin films of polymer melts can exhibit significant slip when sheared along a substrate, where slip lengths much larger than the film thickness have been observed (3⇓⇓⇓⇓⇓–9).

For retracting rims as they emerge after a hole or trench has opened, the magnitude of slip has a direct impact on the dewetting dynamics. When the slip length is very small or zero, the retraction rate is independent of the size of the growing rim and hence, approximately constant, except for logarithmic corrections (10). For slip that is large compared with the film thickness, viscous dissipation increases with the rim size and gives rise to a

The evolution of the polymer melts during these dewetting regimes from early stages after hole formation to the late stages of rupture of the ridges, where the hole boundaries meet, is most appropriately modeled by thin-film models. A systematic asymptotic derivation from full Navier–Stokes equations (22, 23) exploiting the separation between the lateral and normal length scales revealed that the resulting dimension-reduced thin-film model depends on the order of magnitude of slip, leading to two asymptotic distinguished limits: a weak-slip regime and a strong-slip regime.

While thin-film models have been shown to be of great advantage for the analysis of free boundary problems, predicting the evolution over long time and large spatial scales deep into the nonlinear regimes is still a major computational challenge. The primary problem is to resolve the length scales associated with nanoscopic residual layers that remain after the film has dewetted, typically about

In this article, we present a numerical algorithm that is able to answer this need featuring a strategy for local adaptivity and an optimized treatment of the intermolecular potential. We will show that the difference in slip lengths indeed leads to the instability patterns seen in experiments. Our numerical solutions also confirm the Rayleigh–Plateau-type instability of the residual ridges during the late phases of the dewetting for both cases, the no-slip and the intermediate-slip cases, with similar wavelengths at the onset, which had been predicted previously based on a linear stability analysis (17, 20). Similar studies were concerned with cases of infinite ridges with (24) and without (25) gravity followed by a broad range of investigations in the literature for this situation using different contact line models and approximations. Numerically, the work by Diez et al. (26) focuses on finite-length ridges but also, includes a review and elucidating comparison of results for the infinite case. Here, most unstable modes are due to varicose perturbations, and the preferred wavelength of the instability is set by the balance of the destabilizing capillary forces, contact line conditions, and for sufficiently viscous liquids, the viscous dissipation. Stability analysis of the effect of slip vs. no-slip on the instability of a stationary ridge shows that both are linearly unstable and have similar wavelengths (17, 20).

Our numerical method allows us to follow the evolution until rupture and predicts significant morphological differences. In particular, it reveals that, for the no-slip case, the breakup is accompanied by the formation of a cascade of satellite droplets that has never been studied before in this context. Moreover, for the intermediate-slip regime, no satellite droplets appear. Interestingly, a closer look at experimental results confirms these predictions. In other contexts, such as liquid jets or fluid filaments, formation of satellite droplets during rupture is well known, and their destabilization was observed to have a rich structure of intermediate asymptotic regimes [e.g., the works by Tjahjadi et al. (27), Eggers and Villermaux (28), and Castrejón-Pita et al. (29)].

We start by introducing the dissipative free surface flow and the relation to pattern formation. Then, we introduce the nonlinear model for dewetting flow and revisit results for dewetting rates and their impact on instabilities in the early stages of the process. We introduce a highly adaptive numerical method that bridges the multiple length scales and timescales and is able to resolve the dewetting with droplet pinch-off. The solutions are discussed, and we compare with experiments.

## 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

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

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,

We seek the scalar function **1**, where the film height

The second distinguished limit assumes that

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).

### Dewetting Rates.

We first discuss the dewetting dynamics of a straight rim, where we assume translational invariance

For intermediate slip, where

Rims with capillary humps, like the ones that appear here, are known to be subject to Rayleigh–Plateau-like instabilities (20, 24, 26), where the higher capillary pressure in thinner parts squeezes even more liquid into the thicker parts, hence promoting the growth of undulations along the rim. The linear stability of dewetting rims is complicated by the fact that the base state itself grows in time, giving rise to a linearized PDE that cannot be solved exactly using separation of variables. Instead, the linearized PDE can be solved numerically, and the amplification of an initial perturbation can be tracked over time (48). Interestingly, the perturbation evolves into a universal long-time shape that is not sensitive to the initial perturbation. Comparing these shapes reveals an important difference between the no-slip and intermediate-slip cases. The former is much more symmetric and closer to the classical varicose mode observed in the Rayleigh–Plateau instability than the latter (48). Moreover, the maximum amplification is significantly higher for the intermediate-slip case.

These results were analyzed further using an asymptotic sharp interface approach for large rims and a Wentzel–Kramers–Brillouin analysis (also known as WKB analysis, see ref. 49), which established that the long-time dominant wavenumber is given by an equal area rule and is shorter than the prediction by a frozen-mode analysis (50, 51). Furthermore, theory indicates that perturbations of the rim remain small in the no-slip case and for moderate values of B, while in the intermediate-slip case

For slip lengths that are much larger than the film height, the dynamics of the evolution change yet again. In the strong-slip regime, the fluid flow is a plug flow, and the evolution is described by a system of PDEs for the film height and the lateral velocity (6), where the contribution from elongational stresses enters to the same order as the effects from the friction due to slip. Dewetting rim solutions of this model were explored numerically and asymptotically (22, 41). In this regime, the shape of the profile becomes highly asymmetric, with a steep side facing the dewetted area and a much flatter decay to the unperturbed film

## 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

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

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.

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

## 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 **5**) is multiplied with a test function v and integrated by parts to obtain**10a**). Evaluating the solution at discrete times **10** so that we seek

The FEM constitutes a method for the construction of a finite-dimensional subspace **11** valid for all **11** are solved exactly or using a seven-point Gauss quadrature. A typical droplet pattern emerging from the simulation with

The precursor thickness **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

To allow the control of the minimum and the derivatives of ϕ separately via **2**. In particular, with this potential, we still have the same Γ convergence property as **12** and Eq. **2** coincide in this limit. We note that, for *ϕ*″ > 0 for

## Initial Data for Rims

We still need to specify the initial data

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

## Breakup of Liquid Ridges

Toward the final stages of dewetting process, for long times and domain size far beyond

To study the instability of a capillary ridge, we first compute a 1D stationary solution **5** with *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

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

We now consider nonlinear simulations using the stationary ridge with very small monochromatic perturbations as initial data. Specifically, we consider an unstable perturbation with

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).

## 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.

## Acknowledgments

We thank Kirstin Kochems for technical assistance. D.P. acknowledges the financial support of the Einstein Foundation Berlin in the framework of the research center Matheon.

## Footnotes

- ↵
^{1}To whom correspondence should be addressed. Email: dirk.peschka{at}wias-berlin.de.

Author contributions: D.P., K.J., A.M., and B.W. designed research; D.P., S.H., and L.M. performed research; and D.P., K.J., A.M., and B.W. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1820487116/-/DCSupplemental.

- Copyright © 2019 the Author(s). Published by PNAS.

This open access article is distributed under Creative Commons Attribution License 4.0 (CC BY).

## References

- ↵
- ↵
- ↵
- ↵
- Brochard-Wyart F,
- de Gennes PG

- ↵
- ↵
- ↵
- ↵
- Lauga E,
- Brenner MP,
- Stone HA

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Masson JL,
- Olufokunbi O,
- Green PF

- ↵
- Brochard-Wyart F,
- Redon C

- ↵
- Reiter G,
- Sharma A

- ↵
- ↵
- ↵
- Marquant L

- ↵
- Münch A,
- Wagner B,
- Witelski TP

- ↵
- ↵
- ↵
- Sekimoto K,
- Oguma R,
- Kawasaki K

- ↵
- Diez JA,
- González AG,
- Kondic L

- ↵
- Tjahjadi M,
- Stone HA,
- Ottino JM

- ↵
- ↵
- Castrejón-Pita JR, et al.

- ↵
- ↵
- Seemann R,
- Herminghaus S,
- Jacobs K

- ↵
- ↵
- ↵
- Seemann R,
- Brinkmann M,
- Kramer EJ,
- Lange FF,
- Lipowsky R

- ↵
- Dangla R,
- Kayi SC,
- Baroud CN

- ↵
- Eggersdorfer ML,
- Seybold H,
- Ofner A,
- Weitz DA,
- Studart AR

- ↵
- Cuellar I, et al.

- ↵
- Wilczek M,
- Tewes W,
- Engelnkemper S,
- Svetlana V

- ↵
- ↵
- ↵
- Evans PL,
- King JR,
- Münch A

- ↵
- ↵
- ↵
- Fetzer R,
- Rauscher M,
- Münch A,
- Wagner BA,
- Jacobs K

- ↵
- ↵
- ↵
- Flitton JC,
- King JR

- ↵
- ↵
- Holmes MH

- ↵
- Dziwnik M,
- Korzec M,
- Münch A,
- Wagner B

- ↵
- King J,
- Münch A,
- Wagner B

- ↵
- Lessel O, et al.

- ↵
- ↵
- Bäumchen O,
- Fetzer R,
- Münch A,
- Wagner B,
- Jacobs K

- ↵
- McGraw JD, et al.

- ↵
- Zhornitskaya L,
- Bertozzi AL

- ↵
- Grün G,
- Rumpf M

- ↵
- ↵
- Dallaston MC,
- Fontelos MA,
- Tseluiko D,
- Kalliadasis S

- ↵
- ↵
- Bertozzi AL

- ↵
- Kim J,
- Kang K,
- Lowengrub J

- ↵
- Peschka D

- ↵
- Peschka D

## Citation Manager Formats

## Article Classifications

- Physical Sciences
- Physics