Self-organised dip-coating patterns of simple, partially wetting, nonvolatile liquids

TL;DR

This paper predicts a new self-organized pattern formation regime in dip-coating of simple, partially wetting, nonvolatile liquids, driven by a meniscus instability from competing dewetting and film deposition.

Contribution

It introduces a third coating regime characterized by regular liquid ridge patterns, supported by combined analytical and numerical analysis of the hydrodynamic thin-film model.

Findings

Existence of a pattern deposition regime in dip-coating.

Pattern formation driven by meniscus instability.

Mechanism involves competition between dewetting and Landau-Levich film deposition.

Abstract

When a solid substrate is withdrawn from a bath of simple, partially wetting, nonvolatile liquid, one typically distinguishes two regimes, namely, after withdrawal the substrate is macroscopically dry or homogeneously coated by a liquid film. In the latter case, the coating is called a Landau-Levich film. Its thickness depends on the angle and velocity of substrate withdrawal. We predict by means of a numerical and analytical investigation of a hydrodynamic thin-film model the existence of a third regime. It consists of the deposition of a regular pattern of liquid ridges oriented parallel to the meniscus. We establish that the mechanism of the underlying meniscus instability originates from competing film dewetting and Landau-Levich film deposition. Our analysis combines a marginal stability analysis, numerical time simulations and a numerical bifurcation study via path-continuation.

Equations

$${\partial }_{t}h(\mathbf{r},t)$$

$$\mathrm{where}\mathbf{\chi }(h)$$

$$f_{\mathrm{D}}(h)=-\frac{1}{2h^2}+\frac{1}{5h^5}$$

$$h=h_l,\nabla h=(-\alpha ,0{)}^T,\mathbf{n}\cdot \nabla \Delta h=0$$

$$\mathbf{n}\cdot \nabla h=0,\mathbf{n}\cdot \nabla \Delta h=0$$

$$h(x,y=0)=h(x,y=L_y).$$

$$f(h)=f_{\mathrm{D}}(h)+\frac{G}{2}{h}^2.$$

$$\omega (k)=-Q(h_0)k^2\left[k^2+f^{\prime \prime }(h_0)\right]+ikG\alpha Q^{\prime }(h_0).$$

$$h_{\mathrm{f}}\propto U^{2/3}.$$

$$h(x)=\frac{1}{2}(h_{l}-h_{\infty}+\alpha x)\big{(}1-\tanh(x+(h_{l}-h_{\infty})/\alpha)\big{)}+h_{\infty}.$$

$$\xi (h_0):=-f^{\prime \prime }(h_0)\mathrm{and}\eta (h_0):=\frac{\alpha G{Q}^{\prime }(h_0)}{Q(h_0)}$$

$$\tilde{\omega }(k):=\frac{\omega (k)}{Q(h_0)}=\xi k^2-k^4+\eta ik.$$

$$\frac{\mathrm{Re}\left\{\tilde{\omega }(k^{\ast })\right\}}{\mathrm{Im}\left\{k^{\ast }\right\}}=\tilde{c};\mathrm{Re}\left\{{\tilde{\omega }}^{\prime }(k^{\ast })\right\}=0;-\mathrm{Im}\left\{{\tilde{\omega }}^{\prime }(k^{\ast })\right\}=\tilde{c},$$

$${\tilde{c}}_{\pm }(h_0)=-\eta \pm \frac{2}{3}{\xi }^{3/2}\left[\left(2+\sqrt{7}\right)\sqrt{\frac{\sqrt{7}-1}{6}}\right].$$

$$k_s=\mathrm{Re}\left\{k^{\ast }\right\}+\frac{\mathrm{Im}\left\{\tilde{\omega }(k^{\ast })\right\}}{c}.$$

Full text

I Introduction

When a plate is withdrawn from a bath of simple, partially wetting, nonvolatile liquid, two distinct phenomena occur depending on the velocity of withdrawal $U$ Snoeijer and Andreotti (2013): At low velocities, no macroscopic liquid film is transferred onto the moving plate and a static meniscus forms at the bath surface Chan et al. (2012) (see sketch Fig. 1(a)). However, as studied in the seminal works by Landau and Levich Landau and Levich (1942) and by Derjaguin Derjaguin (1945), at higher velocities a liquid film of well-defined height is transferred onto the plate. This occurs above the so-called Landau–Levich transition (cf. Fig. 1(c)) and is the usual situation employed in coating applications Wilson (1982); Kheshgi et al. (1992); Weinstein and Ruschak (2004). The film height $h_{\mathrm{f}}$ of this so-called Landau–Levich film is determined by the velocity of withdrawal $U$, the viscosity of the liquid, and its surface tension. Recently, the behavior of simple, partially wetting, nonvolatile liquids on inclined plates was experimentally Snoeijer et al. (2008); Delon et al. (2008) and theoretically Snoeijer et al. (2006, 2007); Ziegler et al. (2009); Galvagno et al. (2014); Tseluiko et al. (2014); Gao et al. (2016) studied for velocities of withdrawal in the vicinity of the Landau–Levich transition. Further studies of Landau–Levich films for wetting liquids are found in Wilson (1982); Jin et al. (2005); Benilov and Zubkov (2008); Snoeijer et al. (2008); Benilov et al. (2010); Maleki et al. (2011). Landau–Levich films and the related transition are also relevant for other geometries and driving forces, e.g., they occur at menisci driven by surface-acoustic waves Morozov and Manor (2017), in rimming and coating flows, respectively, in and on rotating cylinders Wilson (1982); Melo (1993); Wilson and Williams (1997); Seiden and Thomas (2011) and for the Bretherton problem of a gas bubble moving through a liquid-filled tube Bretherton (1961); Teletzke et al. (1988).

In the present work, we focus on a scenario which has not yet been described, namely the self-organised deposition of a periodic array of liquid ridges from a meniscus onto the moving plate (see Fig. 1(b)). We show that this type of behavior can be found for liquids which exhibit spinodal dewetting Mitlin (1993); Bestehorn and Neuffer (2001) at film heights in the range of the thicknesses $h_{\mathrm{f}}$ of the dynamically produced Landau–Levich films. To investigate the effect, we employ a well-studied mesoscopic hydrodynamic thin-film model for the withdrawal of a plate from a simple, partially wetting, nonvolatile liquid Galvagno et al. (2014). At equilibrium, it allows for the coexistence of a macroscopic drop and a microscopic adsorption layer (precursor film) and is therefore well suited to describe the dynamics behaviour in the vicinity of the Landau–Levich transition as there the deposited coating changes from a microscopic adsorption layer that is macroscopically “dry” to a Landau–Levich coating film. The implied continuous character of the film surface profile across the contact-line region allows for detailed bifurcation studies and the identification of a variety of transition scenarios Galvagno et al. (2014); Tseluiko et al. (2014); Wilczek et al. (2016).

Conceptually, such a deposition of periodic arrays of stripes or ridges belongs to a class of pattern formation processes triggered by a moving front as in the frame of the moving plate the meniscus can be considered as a front that moves with velocity $-U$. This externally imposed front velocity is a control parameter of the system. Examples for investigations of such triggered pattern formation range from the experimentally and theoretically investigated structure formation in Langmuir–Blodgett films Köpf et al. (2010, 2012a, 2012b); Köpf and Thiele (2014); Wilczek and Gurevich (2014); Wilczek et al. (2015), over the study of Cahn–Hilliard-type model equations in one (1D) and two (2D) dimensions for externally quenched phase separation (e.g., by a moving temperature jump for films of polymer blends or binary mixtures) Foard and Wagner (2009, 2012); Krekhov (2009) to the rigorous mathematical analysis of trigger fronts in a complex Ginzburg–Landau equation as well as in a Cahn–Hilliard and an Allen–Cahn equation Goh and Scheel (2014, 2015); Monteiro and Scheel (2017). In the aforementioned systems, a switch from a linearly stable to an unstable state takes place at a certain position within the considered domain. This position is moved with a constant velocity which is (directly or indirectly) controlled, leading to intriguing pattern formation phenomena.

The outline of the present work is as follows: First, in Section II we briefly present the thin-film model including the pertinent boundary conditions and outline the employed numerical methods. Next, in Section III we present the relevant bifurcation diagrams for the Landau–Levich transition and discuss the various branches of different types of linearly stable and unstable surface profiles that exist near the transition. Then, Section IV presents numerical time simulations and shows that there exists a regime where periodic structures form at the meniscus, i.e., a periodic array of liquid ridges parallel to the meniscus is deposited onto the moving plate. A phase diagram is obtained for these states that shows their existence in an extended region in parameter space. Section V employs a marginal stability analysis to determine the upper velocity threshold for the formation of periodic structures as the linear propagation velocity of a dewetting front into spinodally unstable Landau–Levich films. The subsequent Section VI employs an advanced path continuation technique to determine the full bifurcation structure including steady and time-periodic states. These are related to the results of the marginal stability analysis and numerical time simulations for 2D substrates. For the latter, Section VII compares the structures formed by meniscus-guided dewetting to dewetting structures found on horizontal substrates. Section VIII concludes with an outlook.

II Mathematical Model

We model the evolution of the local height $h(\mathbf{r},t)$, $\mathbf{r}=(x,y)$ of the film of simple, nonvolatile liquid in the laboratory frame for the dip-coating geometry shown in Fig. 1 with the well-established long-wave, thin-film or lubrication equation Oron et al. (1997); Craster and Matar (2009) in its nondimensional form, including a Derjaguin (or disjoining) pressure accounting for partial wettability Starov and Velarde (2009); Thiele (2010), gravitational contributions and an advection term accounting for the withdrawal of the substrate Galvagno et al. (2014); Wilczek et al. (2016):

[TABLE]

Here, $\alpha$ is the inclination angle of the substrate in long-wave scaling, $G$ a dimensionless gravity parameter and $Q(h)$ is the mobility resulting for a no-slip boundary condition at the substrate Oron et al. (1997). For the specific scaling used in (1) and (2) see Galvagno et al. (2014). We use the specific wetting potential Pismen (2001); Thiele (2010)

[TABLE]

resulting in the Derjaguin pressure $-f_{\mathrm{D}}^{\prime}(h)=-1/h^{3}+1/h^{6}$.

Equation (1) is in general defined on a two-dimensional domain $\left[0,L_{x}\right]\times\left[0,L_{y}\right]$, where we choose the first dimension to correspond to the direction of withdrawal (cf. Fig. 1). The boundary conditions, which are also employed in Wilczek et al. (2016), read

[TABLE]

The downstream boundary conditions (4) at $x=0$ model the smooth transition to the liquid bath, while the upstream boundary condition at $x=L_{x}$ should permit the liquid to flow out of the considered domain. The value of $h_{l}$ in (4) is chosen large enough that one can safely neglect the next terms of the asymptotic series derived via central manifold projection in Tseluiko et al. (2014).

In the presence of gravity, the hydrostatic term in $\boldsymbol{\chi}(h)$ can be combined with the wetting potential to give the effective local energy (cf. Fig. 2):

[TABLE]

It exhibits a minimum at a small film height $h_{a}$ and thus for a horizontal substrate a macroscopic drop can coexist with a stable thin adsorption layer of height $h_{a}$. The corresponding region outside the drop can be considered as macroscopically “dry” (see, e.g., the case of pure liquid in Thiele et al. (2018)). This corresponds to the “moist” case in De Gennes (1985). This implies for the present dragged meniscus setting that there such an adsorption layer always exists upstream of the meniscus, even below the Landau–Levich transition Galvagno et al. (2014). This is in contrast to the slip model in Snoeijer et al. (2007); Ziegler et al. (2009). Here, the equilibrium contact angle $\theta_{\mathrm{eq}}$ in the long-wave scaling is related to the energy $f(h)$ through $\theta_{\mathrm{eq}}=\sqrt{2\left|f(h_{a})\right|}$ and is for the employed $G\ll 1$ practically identical to the case without gravity.

Next, we discuss the linear stability of a flat film of thickness $h_{0}$ on a resting inclined homogeneous substrate. This is done by inserting the ansatz $h(x,t)=h_{0}+\varepsilon~{}\mathrm{e}^{ikx+\omega(k)t}$ into (1) for $U=0$, expanding in $\varepsilon\ll 1$. To order $\mathcal{O}(\varepsilon)$, we obtain the dispersion relation

[TABLE]

that determines the development of harmonic perturbation modes of small amplitude. For $f^{\prime\prime}(h_{0})<0$, the film is unstable and develops a surface instability that results in spinodal dewetting Mitlin (1993); Thiele et al. (2001a). For the particular local energy (7), liquid films are linearly stable below a critical film height slightly larger than the adsorption layer height $h_{a}$. Then, short-range stabilising contributions dominate over the long-range destabilising terms. Mesoscopic film heights above this threshold are unstable (cf. Fig. 2). Only at larger film heights $h\approx(3/G)^{1/4}$, flat films get eventually stabilised by the hydrostatic contribution.

Our numerical investigations of the equations (1)-(5) are conducted by two complementary approaches. Stable and unstable steady profiles as well as time-periodic profiles together with their bifurcation diagrams are calculated by numerical pseudo-arclength continuation Dijkstra et al. (2014); Engelnkemper et al. (2019) using the package auto07p Doedel et al. (1991a, b). Formulating the 1D version of (1) with $\partial_{t}h=0$ as a small system of three first order ordinary differential equations (ODE) in space (after one integration in space), the parameter continuation allows us to obtain bifurcation diagrams of the steady film profiles which are briefly discussed in Section III. For a recent review of such techniques cf. Engelnkemper et al. (2019), tutorials can be found in Thiele et al. (2014). However, if we employ a spatial finite difference approximation of (1)-(5), we are able to formulate the full time-dependent problem as a rather large system of first order ODE in time. This allows us to also determine stable and unstable time-periodic states by continuation although for a reduced domain size (cf. Section VI). For further implementation details see Köpf and Thiele (2014); Lin et al. (2016) where similar techniques are employed for thin-film equations and closely related driven Cahn–Hilliard equations.

To investigate the time-periodic states, we furthermore conduct numerical time simulations using a finite element method. The approach uses quadratic elements with linear (Q1) ansatz and test functions as well as an implicit second-order Runge–Kutta time stepping scheme. They are implemented using the open-source framework DUNE-PDELab Bastian et al. (2008a, b, 2010). For further implementation details, see also Wilczek et al. (2016), where a similar numerical approach is employed. One-dimensional simulations are performed on a domain $\Omega_{\mathrm{1D}}=[0,400]$ discretised on 400 elements, while two dimensional simulations are performed on $\Omega_{\mathrm{2D}}=[0,600]\times[0,300]$ discretised on $300\times 150$ elements.

III Steady States: The Landau–Levich Transition

First, we discuss steady solutions of equations (1)-(5) in 1D as obtained by numerical continuation. In Fig. 3, the emergence of Landau–Levich film solutions at a critical velocity of plate withdrawal is illustrated. In particular, we show the height of the transferred film $h_{\mathrm{f}}$ in dependence of the velocity of plate withdrawal for three different inclinations $\alpha$. Above the Landau–Levich transition that occurs at about $U=0.04$, $h_{\mathrm{f}}$ monotonically increases following the well-known scaling law derived in Landau and Levich (1942):

[TABLE]

Fig. 4 further illustrates that due to the interplay of wetting potential, plate inclination and dragging by the moving plate the details of the transition to the Landau–Levich film can be rather complex and strongly depend on the inclination angle (cf. Galvagno et al. (2014) for a discussion in the context of unbinding transitions). The lower panels of Fig. 4 give a selection of different stable and unstable steady profile types occurring for $\alpha=2.0$ in the vicinity of the Landau–Levich transition. The upper panel gives the corresponding bifurcation diagram where as solution measure we use the mean amount of liquid per area. In contrast to the cases shown in Galvagno et al. (2014), the branch of Landau–Levich films (cf. profiles I and II) is not connected to the meniscus states but to a snaking branch of detached-foot solutions where an extended foot is connected to the meniscus by a finite length of a thin nearly flat film (e.g., profiles III and IV in Fig. 4). In experimental studies Snoeijer et al. (2008), similar profiles were found to detach from the meniscus as the speed $U$ is increased above the Landau–Levich threshold. Furthermore, another snaking branch of standard attached-foot solutions exists without deep depression between foot and meniscus. Note, however that the film height does not need to be monotonic (see profiles V-VIII). This branch is not connected to the Landau–Levich film branch but to the branch of meniscus solutions. This branch and the corresponding attached-foot solutions are extensively investigated with slip and precursor models in Refs. Galvagno et al. (2014); Tseluiko et al. (2014); Snoeijer et al. (2008); Ziegler et al. (2009). For a detailed study of solution branches and their reconnections with variation of the inclination see Galvagno (2015). Here, the study of the steady states sets the stage for our main focus that lies on time-periodic states. The latter have to our knowledge not yet been described in the literature. In particular, the steady Landau–Levich films discussed at Fig. 3 play a crucial role in the theoretical analysis of the periodic deposition structures and are furthermore used as initial conditions for the numerical time simulations.

IV Periodic Depositions

Next, we perform time simulations of equations (1)-(5) in 1D. As initial condition, we either use a steady meniscus profile or a steady Landau–Levich film obtained in Section III by numerical continuation. At given fixed inclination angle $\alpha$, the Landau–Levich film is always stable at large velocities $U$. Similarly, for sufficiently small velocities $U$, a stable meniscus solution is observed. However, in a certain intermediate velocity range at about and not too far above the Landau–Levich transition, we find that the system converges towards stable time-periodic states. They correspond to a time-periodic detachment of liquid ridges from the meniscus, i.e., at constant frequency. Then, the ridges move as a spatially periodic stripe-like pattern up with the moving substrate (but not at the same velocity). Fig. 5 shows space-time plots illustrating the process for $U=0.07$ and $U=0.08$. Note that the respective arrow faithfully indicates the motion of a point fixed to the moving substrate. This clearly indicates that the ridges travel slower than the plate once they are detached from the meniscus. This is because the ridges slowly slide down the plate due to gravity.

An overview of the behaviour in the $(\alpha,U)$ parameter plane is presented in the form of a morphological phase or state diagram in Fig. 6. The red squares indicate the patterning region situated between regions of steady Landau–Levich film deposition (blue filled circles) and of steady meniscus states (green filled triangles). Each symbol represents a time simulation initialised by a steady solutions at the chosen pair $(\alpha,U)$ as obtained by continuation. However, the borders of the patterning region depend on initial conditions indicating that the system exhibits bistability and hysteresis between steady and time-periodic states. For instance, the grey shaded area in Fig. 6 indicates the patterning regime obtained for a simple initial profile that interpolates between a straight free surfaces for both, the dragged film on the moving substrate and the liquid bath, i.e., the initial condition is

[TABLE]

This initial condition is designed to match the boundary conditions at the left boundary, where $h(x=0)=h_{l}$ and $\left.\nabla h\right|_{x=0}=(\alpha,0)$ is demanded, as well as at the right boundary, where a flat film with height $h_{\infty}=1$ is imposed.

Close inspection of the dynamics in the numerical time simulations leads us to the following interpretation of the mechanism of the observed pattern formation: Above but close to the Landau–Levich transition, linearly unstable flat films are transferred onto the moving plate. These linearly unstable films exhibit spinodal dewetting if placed on a horizontal substrate Mitlin (1993); Thiele et al. (2001a, b); Lam et al. (2018). Here, the same surface instability is active if its timescale is comparable to the one determined by the velocity of withdrawal $U$ and the size of the meniscus region. Thus, the periodic deposition occurs as a result of the instantaneous dewetting of the deposited liquid layer near the meniscus region. Next, we more closely analyse the upper bound of the patterning region in the case where a well-defined steady Landau–Levich film is employed as initial condition of the time simulations.

V Marginal Stability Analysis

A well-established technique for the analytical calculation of the propagation speed of instability fronts that invade linearly unstable states is the marginal stability analysis van Saarloos (1988); Huerre and Monkewitz (1990); van Saarloos (2003). In the context of thin-film flows similar analyses have been employed, e.g., in Limat et al. (1992); Lin and Kondic (2010); Lam et al. (2018).

Here, we employ it to calculate the upper bound of the patterning region as it is given by the plate velocity at which the surface instability starts to invade the steady Landau–Levich film from the upstream end. First, we consider the propagation of a surface instability resulting in dewetting along flat unstable films $h_{0}(x)=h_{0}$ as described by the one-dimensional version of equation (1) omitting the advection ($U=0$). The basic assumption of the marginal stability analysis is that the propagation of an instability front into an unstable state is controlled by the dynamics linearised about this state. With other words, the instability front is pulled by its leading edge that is well described by a linearization of the governing equations.

We introduce the constants

[TABLE]

and the rescaled dispersion relation (8)

[TABLE]

Note that for $\alpha\neq 0$ one has a non-zero imaginary part of the dispersion relation that indicates the occurrence of an oscillatory or travelling-wave instability. Here, it simply refers to a drift of the harmonic spinodal modes. The equations determining the scaled marginally stable linear propagation speed $\tilde{c}$ of a perturbation read van Saarloos (2003):

[TABLE]

where the wavenumber $k^{*}$ is defined in the complex plane and the three unknowns are $\tilde{c},\ \mathrm{Im}\left\{k^{*}\right\}$ and $\mathrm{Re}\left\{k^{*}\right\}$. Since parity of the system ($x\to-x$ symmetry) is broken due to the inclination-related parameter $\eta$, we obtain two distinct marginal velocities $\tilde{c}_{+}$ and $\tilde{c}_{-}$. This is in contrast to systems where parity holds as, e.g., in the case of the decomposition of a binary mixture described by the Cahn–Hilliard equation Scheel (2017), for a film underneath a horizontal substrate destabilised by gravity described by a thin-film equation Limat et al. (1992) or dewetting described by a thin-film equation Lam et al. (2018). The velocities are given by

[TABLE]

Here, $c_{-}=Q(h_{0})\tilde{c}_{-}$ corresponds to a propagation towards the bath while $c_{+}=Q(h_{0})\tilde{c}_{+}$ represents a propagation in the direction away from the bath. Returning to the full system (1) at fixed inclination $\alpha$, where the plate is withdrawn with velocity $U>0$, we can now identify a threshold velocity $U_{\mathrm{lim}}$. It is the velocity of withdrawal which exactly compensates for the speed $c_{-}$ of the instability moving towards the bath along the flat deposited film of height $h_{\mathrm{f}}$. However, here the height $h_{\mathrm{f}}$ is not an independent parameter but given as the thickness of the transferred Landau–Levich film at $U=U_{\mathrm{lim}}$.

In practice, to obtain the threshold velocity $U_{\mathrm{lim}}$ at fixed $\alpha$, we first calculate $h_{\mathrm{f}}$ as a function of $U$ by numerical continuation (see Section III). Then, we employ eq. (14) to calculate $c_{-}(U)=c_{-}\left(h_{\mathrm{f}}(U)\right)$ and finally obtain the zero crossing of $g(U):=|c_{-}(U)|-U$ that gives $U_{\mathrm{lim}}$ as illustrated in Fig. 7 (left).

Fig. 7 (right) shows that these zero crossings at $U_{\mathrm{lim}}$ when calculated as a function of $\alpha$ (shown as a solid line) indeed very accurately predict the upper limiting velocity of the regime of periodic deposition. Note that an analogeous approach for the calculation of such limiting velocities has been employed in Köpf et al. (2012a) for a reduced model system for pattern formation during Langmuir–Blodgett transfer.

In addition, one can calculate the wavenumber $k_{s}$ of the pattern behind the front which is selected by the leading edge of the patterning front van Saarloos (1988), for the case without phase slips at the front. The calculation of the selected wavenumber $k_{s}$ behind the front is based on the assumption that the flux of nodes present in the oscillatory leading edge of the instability front with the wavenumber $k^{*}$ obtained from (11) is conserved across the front Dee and Langer (1983). From this consideration, one obtains

[TABLE]

Here, in the dip-coating geometry, the instability front moves into the system from the domain boundary on the upstram side of the system, i.e., where the Landau–Levich film is dragged out of the domain. The analytically determined wavenumber can only be expected to describe the periodic structures formed close to the onset of the instability well away from the meniscus region. This is only the case at the upper limit of the patterning regime, which is well described by the marginal stability analysis. We have extracted the emerging wavelength $\lambda$ from our numerical time simulations and present them in Fig. 8 (left) in the $(\alpha,U)$-plane where the black solid line again indicates the analytically obtained instability threshold. The right panel of Fig. 8 compares the wavelength $\lambda$ obtained from for two fixed values of $\alpha$ as a function of $U$ with the full curve of results $\lambda_{s}(U)=2\pi/k_{s}(U)$ of the marginal stability analysis. Indeed, the two wavelengths well agree at onset, i.e., on the boundary of the patterning region, but deviate at smaller $U<U_{\mathrm{lim}}$. Note, that the curve corresponding to time simuations for $\alpha=2.6$ hints at a strong increase in $\lambda$ for smaller $U$. This could indicate a divergence of $\lambda$ as seen at a global bifurcation. To further discuss this hypothesis and to better understand the hysteresis between steady and time-periodic states, we next perform a numerical bifurcation analysis that includes time-periodic states.

VI Bifurcation Structure of Periodic Solutions

Up to here, we have beside numerical time simulations and analytical linear considerations employed numerical pseudo-arclength continuation to determine bifurcation diagrams of steady states as presented, e.g., in Fig. 3. To do so, we have reformulated the 1D steady version of (1) as a system of three first order ODEs in space. This has allowed us to use the auto07p-mode for a boundary value problem (BVP) that comes with a fully adaptive spatial discretisation and, therefore, facilitates consideration of relatively large domains. However, to obtain a more complete picture of the behaviour close to the Landau–Levich transition, in particular, close to the onset of periodic deposition where bistability and hysteresis can occur (cf. Fig. 6), we now pursue the numerical continuation of time-periodic solutions. As the toolbox auto07p supports continuation of time-periodic solutions only for systems of ODEs, we discretise Eq. (1) “by hand” with a finite difference discretisation employing a fixed equidistant grid of $N$ nodes in space resulting in $N$ coupled ODEs in time. The full auto07p-intrinsic adaptive grid is then employed for the discretisation in time for the time-periodic solutions only. This approach has also been employed in Köpf and Thiele (2014) for the numerical continuation of periodic solutions of a Cahn–Hilliard model for Langmuir–Blodgett transfer experiments. For the continuation of periodic solutions of a similar problem with periodic boundary conditions, a spatial discretisation by pseudo-spectral methods is performed in Lin et al. (2016). For a recent review see Engelnkemper et al. (2019). Since the continuation of periodic solutions of $N$ ODEs is for large $N$ computationally very expensive, we adapt our system parameters accordingly to $h_{l}=40,~{}L=200$, i.e., move to a smaller system than the one studied before. The system is still large enough that we expect to find all features observed before. We choose a grid with $N=300$ points and tested that a larger value for $N$ does not change the location of the bifurcations. The steady solutions obtained with the BVP approach are reproduced in a good approximation.

The resulting bifurcation diagram for $\alpha=2.0$ is shown in Fig. 9. It well illustrates the intricate pattern formed by the branch of steady states and a number of branches of time-periodic states. Starting at large velocity $U$, we decrease $U$ and follow the branch of steady Landau–Levich films. Before reaching the saddle-node bifurcation at $U\approx 0.058$, we find a series of 17 Hopf bifurcations. The 17 bifurcations form the starting points of 9 branches of time-periodic states. Each of the 8 branches (red dashed lines in Fig. 9) connects two of the Hopf bifurcations and consists of unstable time-periodic states. In contrast, the branch emerging at the first Hopf bifurcation (at $U\approx 0.11$, blue dashed line in Fig. 10) connects nearly vertically with the unstable branch of steady states connecting the two saddle-node bifurcations that corresponds to the detached foot solutions discussed in Sec. III and shown in Fig. 4. The time period diverges upon approach of the steady state indicating that this branch of time-periodic states ends in a homoclinic global bifurcation Strogatz (2014). A partly similar general picture is also found for the aforementioned Cahn–Hilliard model of Langmuir–Blodgett transfer Köpf and Thiele (2014), where the authors termed it a harp-like bifurcation structure.

The branch of time-periodic states bifurcating at the highest value of $U$ emerges subcritically towards the stable Landau–Levich films at higher $U$, i.e., it consists of unstable states. It then passes a first saddle-node bifurcation where it turns back towards smaller $U$ and becomes linearly stable. A second saddle-node bifurcation at slightly lower $U$ destabilises the branch again that now continues again towards larger $U$. This brief interlude becomes slightly more pronounced at larger $\alpha$ (see Fig. 10 (top)). At a third saddle-node bifurcation, the branch again turns back towards smaller $U$ and becomes linearly stable. It remains linearly stable until ending in the global bifurcation. The stable part of the branch corresponds to the line deposition commonly found in the discussed time simulations. The extended subcritical region between the primary Hopf bifurcation at $U=0.1088$ and the third saddle-node at $U=0.1121$ explains the bistability and hysteresis found in the time simulations.

In Fig. 10, we show focus on the region where hysteresis occurs between steady and time-periodic states and show the bifurcation diagram for $\alpha=2.5$ together with examples of space-time plots of unstable and stable time-periodic solutions. Close to the primary Hopf bifurcation, the time-periodic states correspond to small-amplitude oscillations close to the right hand end of the domain (set I). By eye, only two spatial periods can be distinguished. This behaviour is consistent with the interpretation that at the Hopf bifurcation point, the leading edge of an instability starts to propagate from the right hand boundary towards the meniscus counteracted by the advection term. With other words, the Hopf bifurcation obtained through numerical continuation should nearly coincide with the stability border obtained via the marginal criterion outlined in Sec. V evaluated for the film height of the deposited Landau–Levich film. Fig. 11 shows that this is indeed a very good approximation across a range of inclination angles $\alpha$.

However, the time-periodic solutions on the first subbranch, i.e., the first subcritical part of the branch emerging at the Hopf bifurcation, are linearly unstable. This might be explained by considering a marginal propagating front into an inhomogeneous steady state (corresponding to the film height away from the meniscus that decays towards the Landau–Levich film height). Similar theoretical concepts are used in hydrodynamics in the context of so called global modes that occur, in particular, in open flow systems with broken Galilean invariance Couairon and Chomaz (1997, 1999), a property shared by the present dragged film system. In time simulations initialised with individual snapshots from an unstable time-periodic state such as state I in Fig. 10, the instability front typically propagates towards the meniscus and the system subsequently relaxes to a time-periodic state corresponding to the uppermost subbranch shown in Fig. 10, i.e., where states III and IV are located.

Following the branch of time-periodic states from state I through the three saddle-node bifurcations and then towards lower velocities of withdrawal (passing states II, III, and IV of Fig. 10, the height modulations become stronger transforming into a periodic array of ridges. In parallel, the spatial onset of the modulations moves continuously closer to the meniscus. Finally, at sufficiently small velocities $U$, the ridges are formed directly at the meniscus (see state IV in Fig. 10). In this limit one can say that an oscillating meniscus emits ridges that are transported away by the moving plate. This is another reason why a marginal stability analysis of a homogeneous state is not able to give the lower instability onset.

VII Two-Dimensional Flows

Up to here, all our considerations have concerned 1D film flows, i.e., we have assumed that all structures are translationally invariant in the direction orthogonal to the direction of dragging ($x$-coordinate). In the final section we consider 2D dragged films.

It is known that on horizontal substrates in the absence of lateral driving terms, fronts of spinodal dewetting are often transversally unstable Reiter and Sharma (2001); Neto et al. (2003); Beltrame and Thiele (2010). Our numerical time simulations suggest that this is not the case in the considered dip-coating geometry as the dragging of the plate stabilises the transversal instability modes. In Fig. 12, we show snapshots of time simulations for (bottom panels) 2D substrates for a dewetting front on a horizontal substrate and (top panels) for a dragged film. For the dip-coating geometry, i.e., for non-zero inclination angles and dragging velocities, the deposition of regular straight stripes still persists after a long time (Fig. 12 (top)) indicating that transversal ridge instabilities Thiele and Knobloch (2003) upon formation of the stripes are either completely stabilised or so weak that we do not observe them in the finite computation time.

In contrast, a standard dewetting front on a horizontal substrate is clearly transversally unstable, see Fig. 12 (bottom) where such a front is initiated at the left hand end of the domain and moves towards the right. Note that both simulations in the bottom row start with an initially flat film of a height corresponding approximately to the height $h_{\mathrm{f}}$ obtained for the respective set of parameters in the dip-coating geometry, i.e., the two dragging velocities in the top row of Fig. 12 result in two specific $h_{\mathrm{f}}$ that are also employed in the bottom row. For both initial film heights, first the front is nearly transversally invariant but soon develops a strong transversal instability. This occurs even without the addition of explicit noise to the initial conditions. Even without conducting a transient transversal stability analysis of stripe-forming fronts, the time simulations clearly suggest that the stripe deposition process in dip-coating is significantly less unstable due to the presence of strongly symmetry-breaking directional forcing. Note, however, that each deposited stripe is, in general, unstable with respect to a Plateau–Rayleigh-type instability (cf. Thiele et al. (2003); Diez et al. (2009)) that acts, however, on a larger time scale.

VIII Conclusion and outlook

In the present work, we have employed a hydrodynamic thin-film model to numerically and analytically investigate the withdrawal of a solid substrate from a bath of nonvolatile partially wetting simple liquid. We have predicted that in this classical dip-coating geometry, there exists a third withdrawal regime beside the two well-known ones Teletzke et al. (1988) where the substrate is either (i) macroscopically dry (covered by a microscopic adsorption layer) or (ii) homogeneously coated by a macroscopic liquid Landau–Levich film whose thickness depends on the angle and velocity of substrate withdrawal. The third regime consists in (iii) the deposition of a regular pattern of liquid ridges that are oriented parallel to the meniscus. Although the possibility of an instability of the deposited Landau–Levich film is mentioned in the conclusion of Ref. Teletzke et al. (1988), to our knowledge the regime was not analysed before as part of the solution and bifurcation structure of the long-wave description of dip-coating processes.

In contrast to thin-film slip models of dip-coating Snoeijer et al. (2008); Ziegler et al. (2009) that incorporate partial wettability directly through an imposed finite contact angle at the contact line (in slip models the point where the film height strictly goes to zero), here we have employed a precursor-film model Galvagno et al. (2014). It includes a Derjaguin (or disjoining) pressure that guarantees that even a macroscopically dry region is covered by a microscopic adsorption layer and also accounts for the correct macroscopic equilibrium contact angle. Employing such a model has allowed us to investigate the interplay of the dynamic deposition of the film and a surface instability of the deposited Landau–Levich film that is related to spinodal dewetting Mitlin (1993); Thiele et al. (2001a); Lam et al. (2018). This interplay provides the mechanism of the meniscus instability that results in the periodic shedding of ridges that form a regular deposition pattern. One may say that the pattern results due to instantaneous dewetting of the coating film as it emerges from the meniscus.

Our detailed analysis has combined a marginal stability analysis of a dewetting front advancing into an unstable flat film, time simulations in the fully nonlinear regime and a bifurcation study combining two types of numerical path-continuation. Thereby, we have shown that the marginal stability analysis determines accurately the upper limiting velocity observed in numerical time simulations. This analysis determines the threshold substrate velocity below which a linear dewetting instability front that moves from the upstream boundary of the domain into the deposited Landau–Levich film (i.e., towards the meniscus) is not anymore advected away by the moving plate. Furthermore, the fully nonlinear simulations have shown that regular patterns of deposited ridges can also be found in a certain region above this upper threshold velocity. This implies that there is a range of velocities where the system shows bistability between homogeneous and periodic deposition, i.e., hysteresis effects.

The detailed numerical bifurcation analysis revealed a number of Hopf bifurcations on the solution branch describing a homogeneous transfer of which the one occurring at highest velocity corresponds to the upper limiting velocity determined by the marginal stability analysis. Together they result in the emergence of an intricate structure formed by branches of time-periodic states. This solution structure exhibiting sub-critical Hopf bifurcations well explains the existence of bistability and hysteresis observed in the numerical time simulations. Note that somewhat similar structures are found in other related systems where a spatial heterogeneity (here, the meniscus where the bath ends) interacts with a lateral driving force (here, the withdrawal of the plate) in the vicinity of a phase transition (here, the wetting transition occurring when the equilibrium contact angle goes to zero and the dewetting instability ceases to exist). For instance, a more complicated bifurcation structure is found in a study of Langmuir–Blodgett transfer of a layer of surfactant employing a model based on a convective heterogeneous Cahn–Hilliard equation Köpf et al. (2012b); Köpf and Thiele (2014). There, many branches of time-periodic concentration profiles emerge from an intricately snaking branch of steady profiles. The branches of time-periodic states represent the periodic transfer of stripes of different surfactant phases and normally connect a Hopf bifurcation and a global (sniper or homoclinic) bifurcation. Overall, they form a harp-like structure that is topologically similar to our Fig. 9. Note, that in both cases the number of Hopf bifurcation points will depend on system size as more modes fit into larger domains.

Note, that the understanding of the bifurcation scenario of deposition processes is not only important for the dip-coating process at hand and the Langmuir–Blodgett transfer Köpf et al. (2012b); Köpf and Thiele (2014). It is also relevant for line deposition from solutions and suspensions with volatile solvent where other instability mechanisms dominate, see, e.g., the case of evaporative dewetting Frastia et al. (2012); Doumenc and Guerrier (2013). It was pointed out in Frastia et al. (2012) and more extensively discussed in Köpf and Thiele (2014) that the onset of periodic deposition at the lower limiting velocity may be considered a depinning transition as beyond this critical velocity part of the steady meniscus profile depins and is dragged away from the bath as liquid ridge. This relates our result to other depinning transitions, e.g, of driven drops on substrates with defects Beltrame et al. (2011); Varagnolo et al. (2013), of drops on rotating cylinders Lin et al. (2016), and of interface undulations of an air finger in a liquid-filled channel from its advancing tip Thompson et al. (2014). In most cases, depinning is triggered at global (sniper or homoclinic) bifurcations and the emerging branches of time-periodic states normally end in Hopf bifurcations at large driving velocities. Also note that the combination of an imposed plate withdrawal velocity and a heterogeneity (meniscus fixed in the lab system) in our dip-coating system shows similarities with a quenching front moving with an imposed velocity. The latter may produce a variety of patterns (for a Cahn–Hilliard-type system, see Ref. Foard and Wagner (2009, 2012)) as mathematically investigated in detail in Ref. Goh and Scheel (2015).

Finally, we focus on implications of our results for dip-coating, other coating processes and related hydrodynamic systems. The underlying mechanism of the ridge deposition is universal and will apply to all dragged-film and other coating systems when the deposited coating film is linearly unstable, independently of the instability mechanism and geometry. This includes, e.g., films of dielectric liquids destabilised by an electrical field (see appendix of Ref. Thiele and Knobloch (2006)), liquid films deposited on a heated substrate that are unstable w.r.t. a long-wave Marangoni instability (as discussed in Refs. Oron et al. (1997); Thiele and Knobloch (2004)), partially wetting liquids driven by surface accoustic waves (systems as in Refs. Mhatre et al. (2016); Morozov and Manor (2017) but using partially wetting liquids) and coating processes for a fiber Quéré (1999).

Further note that our results are also relevant for the original Bretherton problem of a gas bubble that moves through a liquid filled tube Bretherton (1961); Teletzke et al. (1988). There, recent experiments Kreutzer et al. (2018) investigate long bubbles that move in rectangular channels filled with partially wetting liquid (also cf. Dreyfus et al. (2003); Jose and Cubaud (2014)). They describe three distinct regimes: (i) a fully dewetted regime without a visible film, (ii) a fully wetted regime where the deposited film is stable on the convective timescale of the bubble and (iii) an intermediate regime where the deposited film ruptures and a sequence of droplets emerges. The three regimes closely correspond to the three regimes described here, i.e., our study also provides an outline of the basic bifurcation structure for the Bretherton problem for partially wetting liquids. A somewhat similar transition related to the transition between our regimes (ii) and (iii) was recently described in connection with the dynamics of relatively thick viscous liquid films flowing down a cylindrical fiber Ji et al. (2019). Of interest in the context of the present study is their transition between a Rayleigh–Plateau regime, where a regular sequence of relatively small drops connected by short stable liquid films slides down the fiber (stable travelling wave), and an isolated droplet regime where widely spaced larger drops slide down the fiber connected by an unstable coating film that develops into a number of small slower droplets. In the frame moving with the large drops small droplets are shed from their trailing meniscus. The transition occurs with increasing drop size and distance as larger drops slide faster and leave a thicker coating layer behind. This layer becomes linearly unstable above a critical thickness and the small droplets develop in analogy to the time-periodic shedding of ridges in our case. This implies that the here developed methodology based on the combination of marginal stability analysis and numerical bifurcation study could prove fruitful for the mentioned systems.

The main part of our study has focused on basic mechanisms and the bifurcation structure for one-dimensional substrates, i.e., for an imposed translational invariance in the direction orthogonal to the direction of withdrawal. However, we have also employed numerical time simulations to test the validity of our main results in the full two-dimensional geometry. We have indeed found that the symmetry-breaking forcing introduced by the plate withdrawal in the dip-coating geometry significantly stabilises the formation of stripes as compared to the case of an initially straight dewetting front that recedes without additional forcing on an isotropic homogenous substrate. In the future, it would be desirable to deeper investigate the two-dimensional geometry, i.e., to perform transient transversal linear stability analyses of the one-dimensional solutions as well as numerical bifurcation analyses via continuation techniques. The linear analysis should establish whether there exists a parameter region where the formed ridge patterns are indeed stable with respect to transversal perturbations on any time scale. An investigation of the bifurcation structure should be embedded into a more general analysis of the bifurcation structures of spatially one- and two-dimensional time-periodic states for several of the above described systems resulting in pattern deposition and depinning.

Acknowledgements.

SG and WT acknowledge support by Deutsche Forschungsgemeinschaft (DFG) within PAK 943 (Project No. 332704749). UT thanks the DFG for funding under Grant No. TH781/8-1.