Impact of energy dissipation on interface shapes and on rates for dewetting from liquid substrates

We revisit the fundamental problem of liquid-liquid dewetting and perform a detailed comparison of theoretical predictions based on thin-film models with experimental measurements obtained by atomic force microscopy. Specifically, we consider the dewetting of a liquid polystyrene layer from a liquid polymethyl methacrylate layer, where the thicknesses and the viscosities of both layers are similar. Using experimentally determined system parameters like viscosity and surface tension, an excellent agreement of experimentally and theoretically obtained rim profile shapes are obtained including the liquid-liquid interface and even dewetting rates. Our new energetic approach additionally allows to assess the physical importance of different contributions to the energy-dissipation mechanism, for which we analyze the local flow fields and the local dissipation rates. Using this approach, we explain why dewetting rates for liquid-liquid systems follow no universal power law, despite the fact that experimental velocities are almost constant. This is in contrast to dewetting scenarios on solid substrates and in contrast to previous results for liquid-liquid substrates using heuristic approaches.

The evolution of many physical systems is governed by thermodynamical or mechanical energetic principles [1][2][3][4] . Such principles are versatile instruments that allow the derivation of underlying physical equations 5 . For flows of incompressible liquids, energy-dissipation principles are known for a long time [6][7][8] . In particular for thin-film flows, the great success in the quantitative understanding of viscous flows with contact-line motion has supplied a universal tool that enables the nano-and microstructuring and functionalization of surfaces, but moreover allows to relate flow patterns with liquid properties and substrate chemistry 9 . Typical phenomena governed by such principles are the dewetting of liquids from solid substrates and from liquid substrates, or general wetting and spreading phenomena [10][11][12][13][14] , where the balance of the decline of energy and the dissipation = − ≤  0 E D can be used to derive power-law rates for the velocity of moving contact lines, see Fig. 1.
When such a power-law rate x c ~ t β exists, its exponent β reveals the dominant physical effect, e.g., gravity, surface tension, viscous dissipation in bulk and on interfaces, and the geometry of the problem 15 . One basic assumption behind such rate estimates is that there exists a simple relationship between the rate of change of the energy and the shape of the time-dependent domain  Ω ⊂ t ( ) 3 occupied by the liquid layer, often also including assertions about the self-similarity of the evolution. For instance, for a large class of free boundary problems where a liquid  ( ) dewetts from a substrate (s) with a straight contact line, the change of surface energy can be approximated by is the spreading coefficient of the system constructed from the corresponding surface-and interface-tensions γ α of interfaces Γ α with surface area |Γ α | and  x c is the contact line velocity. However, it is challenging or even impossible to find a general and similarly simple closed-form approximation for the energy-dissipation rate since the flow field  Ω → u: 3 and the corresponding shear stress τ(u) have a complicated local structure that depends on fine details of the shape of Ω. Therefore, one requires a deeper understanding of the specific dissipation mechanisms and of the domain evolution in order to understand the dynamics of the corresponding processes.
In the pioneering works for liquid-liquid dewetting by Joanny 13 and Brochard-Wyart et al. 14 dewetting rates for small equilibrium contact angles and for limiting regimes of the liquid-liquid viscosity ratios were predicted. While both works combine valid hydrodynamical and dissipation arguments to derive expressions for contact line velocities, the impact of non-trivial interface shapes on the flow and dissipation remained unclear. Since then, many theoretical studies are concerned with the derivation of appropriate thin-film models to study the long-time morphological evolution of the liquid layers. Apart from investigations into stationary states and how they are approached 16,17 , a number of studies focussed on modes of instability in liquid-liquid dewetting using stability analysis and numerical simulations of the thin-film models [18][19][20][21][22] even with additional surfactants 23 .
On the experimental side dewetting rates and morphologies for liquid-liquid model systems such as polystyrene on polymethyl methacrylate have been investigated systematically by Krausch et al. 24,25 by varying the heights and viscosities of the liquid layers. Similar experimental studies were performed by Pan et al. 26 for further layer viscosities and heights. However, the shapes observed by Krausch et al. 25 differ considerably from the empirical predictions used to derive dewetting rates 13,14 , which were found to be constant.
To the best of our knowledge, fundamental dynamic properties like dewetting rates have not been settled up to now. The main reason for this is certainly the absence of theoretical confirmations for the observed shapes of dewetting rims, which then might help to understand the mechanisms behind certain dissipation balances and dewetting rates. Additionally, a quantitative study also requires the key parameters of the experimental system, i.e., surface tensions, viscosities, and layer thicknesses, to be determined sufficiently precise. The focus of this study is thus to supply a quantitative understanding of the dewetting mechanics by detailed comparisons of experimentally obtained rim shapes, their evolution, and their dewetting dynamics with those computed from thin-film equations. Additionally, we examine the underlying mechanisms by discussing flow fields and energy dissipation mechanisms.

Experimental and Theoretical Methods
fields in the x-z-plane. The layered polymer system was prepared using standard thin film preparation techniques 17 as shown in Fig. 2. In a first step, the underlying liquid PMMA substrate is directly spun from a toluene solution onto a silicon substrate that was previously cleaned with Piranha etch. Simultaneously, the upper PS layer is spun from a toluene solution onto freshly cleaved muscovite mica. In a second step, the PS layer is transferred onto the surface of ultra clean water and picked up from there with the PMMA coated silicon substrate. During the transfer of the thin PS layer onto the water, the PS layer ruptures into smaller pieces, which are subsequently transferred onto the spin coated PMMA layer. Straight boundaries of these patches that are sufficiently remote not to be disturbed from neighboring patches are selected to observe the dewetting process. The cross section of these patches are almost ideal rectangular steps and thereby correspond to the start configurations at t = 0 introduced above. The typical film-thicknesses  h h , s used in our dewetting experiments range from 45 nm to 250 nm and we performed experiments for various ratios  h h : s . To remove potential nanoscopic air bubbles that might have been been trapped between the PMMA and the PS layer during the transfer process, the samples were allowed to set after preparation for at least 24 h. The dewetting process is then started by heating the materials above the glass transition temperature and monitored by in situ atomic force microscopy (AFM). The dewetting experiments were conducted at a temperature of T = 140 °C. The shape of the PS-air and PMMA-air interface can be determined in situ using AFM in soft tapping mode. Quenching the sample to room temperature the shape of the buried PS-PMMA interface can be additionally determined by AFM in tapping mode after stripping the upper PS layer with a selective solvent (cyclohexane, Sigma Aldrich) 17 . The full shape of all polymer interfaces is obtained by composing PS-air, PMMA-air and PS-PMMA surfaces, a procedure which generates shapes as shown in Fig. 3. This composition of the 3-D image requires rotation, shift and tilt of upper and lower AFM scan as postprocessing for a perfect match. The contact line is also aligned parallel to the y-axis, so that cross sections can be averaged over a few scan lines in the y-direction. Both polymers were purchased from Polymer Standard Service Mainz (PSS-Mainz, Germany) with polydispersity of M w /M n = 1.05 and molecular weights of M w = 64 kg/mol and M w = 9.9 kg/mol for PS and PMMA, respectively. The glass transition temperatures are T g = 100 ± 5 °C for PS(64 k) and T g = 115 ± 5 °C for PMMA(9.9 k). The viscosities of these polymers at T = 140 °C were determined for PS as μ ≈  700kPas and for  PMMA as μ s ≈ 700 kPas applying the method of self-similarity profiles of stepped polymer films 27,28 . Using stationary droplet profiles 17 we experimentally determined the involved surface tensions to γ = . ± . The fluid flow at time t is described by the continuous velocity field  Ω → u: 3 and satisfies a no-slip condition at the solid substrate, i.e., u(x, y, z = 0) = 0. In particular, this flow is incompressible ∇ ⋅ = u 0 and according to (2) it dissipates the energy s s s s where we introduced the shear stress for a two-phase system consisting of Newtonian liquids as , so that viscoelastic effects can be safely neglected 32 and the treatment of the PS-PMMA system as Newtonian liquids is justified.
The energy from (1)   where we introduced the wetted area as the time-dependent support set of The apparent problem of an infinite energy  in (7) for unbounded ω in (8) can be circumvented by formally restricting the whole domain to a sufficiently big finite box. For given dissipation  and energy  it is basically known since the works of Helmholtz 6 and Rayleigh 7 , that the evolution of the system is given by a minimal dissipation principle v or equivalently upon differentiation by the weak formulation a(u, v) = f(v) for all divergence-free v, where using (5) we have the bilinear form The computation of f requires the formal variation of surface measures |Γ i | with respect to perturbations of the surface by a flow v. This variation involves the Laplace-Beltrami operator 33 and for our energy  can be expressed as where κν is the mean curvature vector on Γ  and ∇ the tangential gradient. The Young force ν γ = ∑ f i i i cl appears at the contact line ∂Γ and is generated by performing integration-by-parts using the Laplace-Beltrami on surfaces. In the context of finite elements discretization for non-parametrized surfaces this method was introduced by Dziuk 34 . This finishes the construction of the variational structure behind the energy-dissipation equality E D = −  that holds by construction. When one derives the partial differential equation that formally corresponds to (9), then one usually introduces an additional pressure  Ω → p: variable that acts as a Lagrange multiplier to enforce the incompressibility condition. The resulting model for the flow of highly viscous Newtonian fluids is the coupled Stokes system (6)  , and corresponding jump conditions on the interface Γ  s , generated by (10). This implies a condition for the pressure-jump at the interface Γ  s, . At the contact line the Young force f cl = 0 imposes further conditions on the triple junction using the Neumann triangle. When the velocity field is computed by solving (11), the domain Ω(t) is evolved according to the velocity field u, so that in particular the velocity u Γ of points on the free boundary satisfy the kinematic condition on the free interfaces and surfaces Γ i . The fact that the domain shape Ω(t) is part of the unknowns makes the problem a free boundary problem. In the following we understand the evolution of the domain shapes parametrized using non-negative functions  (3) and (4) and shown in Fig. 3. For simplicitly the dependence of  h h , s on y will be dropped, since the solutions are assumed translation invariant in y-direction due to the particular experimental setup.
The system of Eq. (11) is now non-dimensionalized using   . On the complement  ω \ only h s is unknown and solves the standard thin-film equation  In order to be able to predict interface shapes for this model, we developed a novel finite-element based numerical scheme 35 , which uses advanced energetic arguments to discretize the contact line motion with the natural and essential interface conditions at x c mentioned above.

Discussion of Shapes and Rates
For a fixed thickness ratio (13) shows that, due to the absence of other intrinsic time and length scales for Newtonian liquids, the influence of the absolute height is to scale time proportionally without changing the rescaled rim profiles. To check this prediction two experimental liquid-liquid systems with thickness ratio =  Fig. 4 for corresponding dewetting distances. The good reproducibility of the characteristic rim profiles within experimental errors confirms the previously made assumption that the fluids can be considered Newtonian and allows us to focus our study on different aspect ratios. Comparing experimentally measured and theoretically computed interface profiles, we find an excellent agreement of both the characteristic shapes of the liquid-air/substrate-air interfaces measured by in-situ AFM in . Also some material of the liquid substrate is dragged along generating a depletion near x < x c and an accumulation of substrate material near x > x c . The contact line itself is elevated by the flow, a dynamic feature quite common for soft substrates 36 but not observed in stationary droplets for sufficiently thick substrates 17 . Right next to the contact line, the liquid-liquid interface extends deeply into the substrate and generates a trench which generates additional resistance against the dewetting motion. The size of this trench depends only weakly on the size of the dewetting rim, i.e., the dewetting distance.
Compared to the ratio =  h h : 1:1 s in Fig. 6(a,b), thickness ratios of 2:1 or 1:2 do not lead to qualitatively new features. For smaller aspect ratio = For small dewetting distances, the dewetting rates in Fig. 7 suggest a linear behavior x c ~ t for all thickness ratios in agreement with the results by Lambooy and coworkers 24 . For fixed substrate film thickness h s , the dewetting rates are larger for liquid layers thinner than the substrate, <  h h s , and smaller for thicker liquid layers, >  h h s . But, a close inspection of the seemingly constant dewetting rates in Fig. 7 (left) indicates that the dewetting velocity slowly decreases over time. This fact is most apparent for aspect ratio 2:1, while for aspect ratio 1:2 the velocity even appears to increase. However, the experimental accuracy is not sufficient to fully clarify this claim.
To clarify the dependence of the dewetting rates, results from simulation are plotted in Fig. 7 (right) for physical dewetting times of several month, which are not accessible experimentally together with further results for other film thickness ratios. Note the small variation in the velocities during dewetting, which explains why dewetting rates appear almost constant. However, the intricate transient behavior of the velocity  x c featuring inflection points in the simulations coincides with the before mentioned experimental observation. For instance, for an aspect ratio of 2:1 and the experimentally accessible (normalized) times t = 10 −1 … 10 0 h nm −1 , cf. Fig. 7 (left), the dewetting rate decreases, while for an aspect ratio of 1:2 the rate slightly increases within the observed dewetting interval. Note the striking agreement of all experimental and theoretical rescaled rates x c (t) also suggests the validity of the introduced parameters, i.e., viscosities and surface tensions. Furthermore, for all simulated parameters we find that for large times the velocity slowly decays to zero. Next, we are going to discuss the physical dissipation mechanism behind the observed transient dewetting rates.

Discussing the Role of Dissipation
Since the thin-film model accurately predicts shapes and speeds of the liquid-liquid dewetting, we extend our approach and discuss local flow features that are experimentally inaccessible in order to explain the observed dewetting dynamics. For instance, the rescaled dissipation balances with the driving surface tension in a 2-D cross section according to z z . While the driving force is straightforward to understand, the dissipation depends on local details of the flow field. As a representative example to discuss the qualitative behavior we use numerical solutions with thickness ratio 1:1 and show rim profiles at different times overlapped with the dominant horizontal component of the velocity reconstructed using (17) in the left panel of Fig. 8 and the corresponding dissipation ∂ the dissipation vanishes at the liquid/air and substrate/air interfaces, whereas the flow field is zero at the solid/substrate interface z = 0. The latter results in a large shear rate and a large energy dissipation at the solid interface, cf. right panels of Fig. 8. Close to the backflow and close to the contact line the maximal dissipation density is reached. However, due to the small size of these regions the integrated dissipation near the contact line and in the remaining rim are of the same order, at least for the transient times and moderately large rims considered here. To visualize this fact, we additionally show the cumulative dissipation inside the (liquid) substrate and inside the (dewetting) liquid in Fig. 8 for different times. Since the shear rate is large at the solid interface where z = 0, clearly the dissipation for an aspect ratio 1:1 is large in the substrate for the short and intermediate times considered experimentally. Nevertheless, with the volume of the liquid rim increasing in time, ultimately the dissipation in the liquid layer will dominate for large times or for higher aspect ratios. A slightly more detailed visual description of this dynamics and the corresponding experiments is provided in the attached supplemental video.
Accordingly, one can identify two different zones where the energy is dissipated in the liquid and the substrate. A significant amount of the dissipation is produced in a small region near the contact line. This can be seen in the steep increase of the cumulative dissipation in the right panels of Fig. 8. The remaining contribution to the dissipation is more or less evenly distributed over the rim width resulting in a moderate and constant increase of the cumulative dissipation over the width of the rim. For large times this bulk contribution will dominate the dissipation and forces the velocity to decay to zero. This can also be seen in the temporal evolution of the dissipation profiles, which is decreasing due to the quadratic dependence on the velocity scale in (18). This variable contribution to the energy dissipation  ∼ directly impacts the observed dewetting rates. This can be best explained on the basis of the known dewetting behavior on solid substrates, which qualitatively also applies for dewetting from liquid substrates. When the contact line position is x c (t) and starts at x c (t = 0) = 0, then asymptotically the rim cross-sectional area grows proportionally to ×  x h c by volume uptake from the unperturbed liquid layer. Consequently, assuming self-similar growth, typical rim geometry such as rim width or like rim height grow proportionally to x c . Additionally, we assume self-similar growth of the dissipation according to , where longer simulation times (right) reveal that rates  x c decrease depending on details such as the aspect ratio.  where it remains to specify α using the dominating dissipation mechanism and its dependence on rim geometry.
In the intermediate slip model 11 , the dominant contribution comes from a substrate dissipation, so that the total dissipation is proportional to the rim width and thereby α = 1/2 in (19). Then, the energy-dissipation balance produces c leading to an asymptotic x c (t) ~ t β dewetting law with β = 2/3. Another example 11 is the no-slip model, where the dissipation is predominantly localized near the contact line and in the bulk domain. The contact line area does not scale with the volume, whereas the gradients in the dissipation in (18) cancel the growth of the cross-sectional area, thereby both leading to α = 0 in (19). This produces a linear dewetting law β = 1, however with logarithmic corrections 37 . Similarly, the power-law dewetting rates on liquid substrates predicted by Brochard et al. 14 rely on the assumptions that the dissipation is generated in only one such localized zone together with a nearly self-similar growth of rim shapes. However, these assumptions fail in the considered situation of liquid-liquid dewetting since dissipation is clearly not generated in one single zone but accumulates in the substrate, in the liquid and near contact lines on a similar order of magnitude. This explains why in our setting the liquid-liquid dewetting process is not in a regime dominated by a specific dissipation mechanism that would admit a simplification to a power-law rate, and thereby challenges the applicability of the theoretical results by Joanny 13 and Bochard-Wyart et al. 14 to experimental systems considered in this paper in early stages. However, the weak scaling of the dissipation with increasing x c theoretically explains the nearly linear dewetting rate, that was observed experimentally by Krausch 24 and in the present work.
The consideration of the liquid-liquid dewetting using thin film models with explicit contact line dynamics, conducted here, allows to describe the variable energy dissipation in a liquid-liquid system and to quantitatively derive rim shapes and dewetting rates. Nevertheless, the predicted slowdown of the dewetting velocity is expected when the dissipated energy is not soley confined to the contact line. The exact slopes in the log-log plot of Fig. 7 (right) depend on details such as thickness ratio and viscosities, and thereby do not support a universal power-law behavior. This observation confirms previous speculations by Krausch et al. 24 that were based on experimental findings about the transient nature of the experimentally measured dewetting dynamics.

Conclusion
Motivated by the long-standing puzzle between theoretically predicted and experimentally observed rates for liquid-liquid dewetting, we performed a combined theoretical and experimental investigation of the transient interface shapes and dewetting rates. Conducting a full simulation of the sharp interface thin-film model for Newtonian liquids without any a priori assumptions on rim shape development or energy dissipation we obtained a full agreement with experimentally determined interface shapes and dewetting dynamics using the relevant experimental parameters like viscosities, aspect ratios, and surface energies. As the main tool to assess the transient nature of the flow, we reconstructed local flow and dissipation fields. Such a detailed analysis of a local energy balance provides deep insights into underlying mechanisms driving such a process.
By analyzing the local energy dissipation, we have found that the liquid-liquid dewetting system is in a transient state with no self-similar behavior and the dissipation is not distributed exclusively at the contact line, in the substrate, or in the bulk. While the dewetting rate in the observed experimental regime is almost constant and thereby of powerlaw type in the strict sense, the absence of self-similarity and localization of dissipation underlines the absence of a dominant mechanism behind this rate. A similar energetic argument provided the explanation why, for very large times beyond experimental reach, the dewetting velocity slowly decreases to zero. Such predictions would be impossible using heuristic approaches, since the transient internal flow is rather complex and results from a complex interaction of substrate and liquid. Without such a theoretical toolbox, the observed dewetting rates might otherwise be misinterpreted as a regime with potential dominant physical effects. The demonstrated ability to use energetic arguments to quantitatively describe liquid-liquid systems thus set grounds for a similarly complete understanding as already obtained for liquid-solid dewetting systems. In particular, the analysis of the local dissipation distribution provides a powerful tool to identify dominant physical regimes or their absence. It might be in particular possible to extend this approach also to fluids with complex rheological behavior.