On the depinning of a driven drop on a heterogeneous substrate

We study the depinning of driven drops on a heterogeneous substrate using a long-wave evolution equation for the film thickness profile. The heterogeneity is incorporated into an additional pressure term describing the interaction of the liquid with the substrate or with an external field. A drop may be pinned by a hydrophilic defect at the back or a hydrophobic defect at the front. Two types of depinning occur depending on the strength of the driving. The first occurs via a saddle-node (sniper) bifurcation, while the second involves a Hopf bifurcation. The parameter dependence of the depinning process is studied using linear stability theory, and direct numerical simulations are used to explore the dynamical properties of the stick-slip motion of the drop after depinning.


Introduction
It is well known that liquid drops on an ideally smooth substrate move in response to external gradients. For example, a drop on an inclined substrate slides downslope in response to the gradient of potential energy [1,2,3]. Likewise a drop in a temperature gradient will move towards lower temperatures as a result of Marangoni forces caused by surface tension gradients [4,5]. A wettability gradient induced by chemical grading of the substrate also causes drop motion: in order to minimize its energy the drop moves towards the most wettable region [6,7,8,9].
Although on ideally smooth substrates drops will move even for arbitrarily small gradients, this is not the case for the 'real' substrates used in experiments. There the onset of drop motion is strongly influenced by spatial heterogeneities such as wettability defects [8,10,11,12,13,14,15,16,17,18], defects in the topography of the substrate or variations in the substrate temperature [19,20]. For drops in a condensor [21,22] spatial variations in the electrical field act in much the same way as wettability defects [23]. In all these cases a finite driving force is necessary to overcome the pinning influence of the heterogeneities [12,14,16,24,25,26,27]. In addition, substrate heterogeneities can suppress the coarsening of droplets in the late stages of dewetting [28,29,30]. Indeed, heterogeneities occurring on a micro-or mesoscale are known to affect the macroscopic movement of drops and are responsible, for instance, for the observed hysteresis between advancing and receding contact angles [31,32,33], and the roughening of contact lines [31,32,34,35,36]. Another signature of substrate defects is the so-called stick-slip motion of weakly driven contact lines encountered, for instance, in the late stages of capillary rise [26,37] or sessile drop experiments [38].
A common approach to drop pinning-depinning and the associated stick-slip motion is to regard the contact line motion as a transition between different metastable states as determined by the minimization of appropriate energies [16,17]. However, the resulting description is quasi-static and cannot account for the dynamics of the resulting stick-slip motion. In contrast, a dynamical model has to take into account not only the underlying energy landscape but also the dynamics of the relaxation towards the local minima. Other theories [27,39] predict contact angle hysteresis without substrate heterogeneities but cannot account for stick-slip motion.
In this paper we use a dynamical model to describe driven drops in a spatially heterogeneous system, and focus on understanding not only the processes that lead to pinning, but also the bifurcations that lead to depinning and the subsequent behavior of a moving drop. Pinning and depinning transitions have, of course, been studied extensively in a variety of fields. The simplest situation is described by the Adler equation [40] θ = µ − sin θ, where θ represents, for example, the drop position, and µ > 0 measures the strength of the driving force. When µ < 1 this equation has a pair of fixed points, one of which is stable and the other unstable. At µ = 1 these fixed points annihilate in a saddle-node bifurcation, but unlike the standard saddle-node bifurcation this bifurcation produces periodic motion for µ > 1. This result is simplest to understand if we write (1) asθ = −dV/dθ, V ≡ −µθ −cos θ. Evidently, (1) represents an overdamped particle in a cosinusoidal potential that is progressively tilted as µ increases. A 'particle' in a stable equilibrium at a local minimum of this potential 'spills out' once the tilt becomes large enough that its position no longer corresponds to a minimum. This occurs precisely at µ = 1. The periodic motion present for µ > 1 corresponds to the particle sliding down the resulting 'washboard' potential.
It is important to note that the period of this motion diverges as (µ−1) −1/2 [41]. The resulting bifurcation is sometimes called a Saddle-Node Infinite PERiod bifurcation, or 'sniper' for short.
Many depinning phenomena in physics may be understood using this simple picture. Usually this is so in systems with a continuous symmetry such as invariance under translations. In the absence of a heterogeneity spatially periodic structures may undergo a parity-breaking bifurcation that breaks the left-right symmetry of the pattern and produces a drift. The direction of the drift is determined by the associated tilt of the structure [42]. In this case the drift speed of the structure vanishes as the square root of the distance from the paritybreaking bifurcation. However, in the presence of spatial heterogeneities the situation changes dramatically because near the bifurcation even small amplitude heterogeneities suffice to pin the tilted structure. A detailed study of this regime [43] shows that while some depinning events are indeed analogous to the behavior described by the Adler equation, a quite different depinning mechanism is present as well. Here the tilted periodic state first undergoes a Hopf bifurcation that produces back-and-forth rocking motion of the tilted structure, but no net translation. As a parameter increases the amplitude of this oscillation increases leading to a global bifurcation involving an unstable fixed point and its translate by one period. This bifurcation generates oscillations with a nonzero mean drift, and this net drift increases with further increase in the parameter. In fact the production of oscillations with net drift is a common consequence of global (or gluing) bifurcations. Other examples of global bifurcations leading to drift can be found in the theory of the Josephson junction [44] and in spiral wave convection [45]. In contrast the pinning of fronts connecting homogeneous and structured states to the underlying structured state is described well by the sniper bifurcation [46,47]. This paper is devoted to the understanding of the depinning transition in a simple dynamical model for driven liquid drops or ridges in a spatially heterogeneous system. We choose a generic model describing driven drops on a solid substrate coexisting with a thin film or wetting layer [48], and introduce a heterogeneous interaction with the substrate with a well-defined spatial period. The heterogeneity might arise from chemical wettability defects, topographic defects or heterogeneous temperature or electrical fields. Random heterogeneities [34,35,36,49] are not considered. The lateral driving might be caused by gravitational or centrifugal forces. Other possibilities include gradients of wettability, temperature or electrical fields.
We use lubrication theory and model the effect of heterogeneities through an additional spatially periodic pressure term. This term may consist of a disjoining pressure [30,32,48] modelling wettability or contributions that account for the influence of electrical [21,23] or temperature [50,51,52] fields, or a combination thereof. This dynamical approach has the advantage that drop shape and contact angles are not a priori set, and permits us to study the evolution of the shape and stability of a drop trapped by a heterogeneity as a function of the driving. Three cases are considered in detail: a periodic sinusoidal modulation of the pressure term, and strongly localized 'hydrophilic' or 'hydrophobic' defects. As a result the drop may be pinned at the rear by a 'hydrophilic' region, or at its front by a 'hydrophobic' region. However, we do not seek an asymptotic description of macroscopic drops [9], focusing instead on the nature of the transitions involved in the process of depinning.
The paper is organized as follows. In the next section we summarize the thin film model we use for the study of drop depinning. In Section 3 we present the bulk of our results, and show that depinning occurs via both the sniper and Hopf bifurcations, and identify the conditions under which these mechanisms apply. These results are obtained by a combination of branch following methods, linear stability theory and direct numerical simulation of the model equation. Section 5 provides a brief outline of a theoretical understanding of the observed behavior using a second order equation valid near a codimension-three point in parameter space. We also include a brief summary of the corresponding results for thicker films for which the corresponding flat film is metastable, i.e., no primary linear surface instability exists (Section 6). All our calculations are two-dimensional. Our conclusions follow in Section 7. A two-dimensional liquid layer on an inhomogeneous solid substrate subject to a lateral driving force of strengthμ is modelled by an evolution equation for the film thickness profile h(x, t) derived from the Navier-Stokes equation using the long-wave approximation [48]:

Model
Here γ is the surface tension, η is the dynamic viscosity, while P (h, x) is an additional pressure term accounting for wettability, heating, electrical fields etc. [23]. As a particular example we employ a disjoining pressure describing wetting properties [32,53,54,55], , where κ is a typical energy density scale and b > 0 denotes the ratio of two antagonistic interactions. Our choice corresponds to a combination of long-range power law and short-range exponential interactions [56,53,57], with the spatial heterogeneityξ(x), i.e., the influence of wettability defects, contained in the short-range part. Some examples of the dimensionless form of P (h) for a homogeneous substrate and the resulting stability diagram for a flat film are shown in figures 2 (a) and (b), respectively. Whenμ = 0 the resulting model allows for solutions representing static drops of arbitrarily large volume with a finite mesoscopic equilibrium contact angle sitting on a precursor film or wetting layer. We nondimensionalize Eq. (2) using the scales 3ηγ/κ 2 l for time and lγ/κ for the lateral coordinate. The film thickness is scaled by the interaction length l that provides a characteristic scale for the thickness of the wetting layer. The ratio of vertical and horizontal length scales lκ/γ, i.e., the smallness parameter of the underlying lubrication approximation, is closely related but not equivalent to the equilibrium contact angle for drops on a homogeneous substrate, cf. [9]. In addition we define the dimensionless driving standing, in general, for lateral forces that depend only weakly or not at all on film thickness like gravitational or centrifugal driving or a purely lateral temperature gradient. In the following the direction of driving and its opposite are called 'downstream' and 'upstream', respectively. We emphasize that any qualitatively similar pressure P yields like results, as shown for dewetting in [54,57,58,59,60,61], heated wetting films in [51], and for chemically driven running drops in [62,63,64].
To uncover the details of the depinning transition is suffices to restrict attention to drops that are not much higher than the wetting layer. This case includes not only nano-drops on ultrathin precursor films (l ≈ 1 . . . 10 nm), where P only incorporates wetting properties, but also micro-drops on macroscopic wetting layers (l ≈ 10 . . . 100 µm). In the latter case P may include wettability in the form of repulsive van der Waals interactions together with the influence of temperature and/or electrical fields [51,52,23]. The case of drops in a condenser is illustrated in Appendix A. The spatial modulation responsible for pinning can be included in any part of P . In the following we call any heterogeneity that attracts [ [19].
In the following we allow the short-range part of P to depend on the coordinate x, where ξ(x) is a periodic function of x, and take or where K(k) is the complete elliptic integral of the first kind and ∆ is such that the average of ξ(x) over a spatial period vanishes. The choice (5) allows us to vary the profile of the heterogeneity from sinusoidal, ξ(x) = sin(2π x/L het ) at k = 0, to a step profile as k → 1, with hydrophilic and hydrophobic stripes of equal width. In contrast, (6) yields profiles with unequal widths: as k → 1 (6) describes localized hydrophilic ( < 0) or hydrophobic ( > 0) defects. In the following it will be convenient to use the logarithmic measure s ≡ − log(1−k) to quantify the resulting heterogeneity profiles. The dimensionless film evolution equation then takes the form where P is the dimensionless pressure In the absence of lateral driving, i.e., µ = 0, one can assign an energy to every film profile h(x, t) using the Lyapunov functional [60,30] where L is the system size,h ≡ L 0 h(x)dx/L is the mean film thickness, f(x, h) = − P (x, h)dh is the local free energy, and C 1 is a Lagrange multiplier corresponding to mass conservation. The minima of this functional correspond to stationary solutions of Eq. (7), i.e., solutions of In writing this expression we have set the first integration constant C 0 to zero, as required by the symmetry of the problem with respect to reflection in x about a suitably chosen origin, determined by the form (5,6) of the heterogeneity. However, in a precursor film model with lateral driving, µ = 0, this reflection symmetry is absent, and in general C 0 = 0. In this case profiles h(x) that are stationary in the laboratory frame satisfy such profiles are accompanied by a nonzero flux of fluid through each spatial period. The resulting situation differs from slip models used to study drops of finite length with 'true' microscopic contact angles which permit a variational formulation. In these theories the onset of depinning is determined by comparing energies of metastable states, but the subsequent dynamics requires the specification by hand of the contact angles at the front and back (for examples involving contact line motion on homogeneous substrates, see [65,66]). This makes the slip Ansatz problematic since for moving drops both contact angles are determined by an interplay between heterogeneity and driving.

Results
In the following we focus on values of b for which the range of metastability extends to thick films, i.e., b 1.1 (cf. figure 2(b)). In this regime flattened (pancake) drops are unstable.

Unforced system
The solution structure of Eq. (7) with a sinusoidal wettability pattern and a disjoining pressure derived from diffuse interface theory [67] but no driving (µ = 0) is extensively studied in [29,30], focusing on the transition from long-time coarsening to pinning at the substrate heterogeneities that occurs with increasing heterogeneity strength or decreasing heterogeneity period. In contrast, the focus of the present paper is on the details of the depinning transition that occurs as the lateral driving increases. We recall first that thin films on a homogeneous substrate with ∂ hh f < 0 are unstable. This instability is a long wave instability and sets in when the system size or equivalently the spatial period L exceeds L = L c . A supercritical bifurcation then produces a branch of periodic drop-like solutions for L > L c . However, if the primary bifurcation at L = L c is subcritical such solutions extend to L s < L c , and connect to the flat film at L = L c via a branch of unstable nucleation solutions [2]. Solutions for different mean film thickness are shown in figure 3 in terms of their L 2 norm ||δh|| ≡ L 0 (h(x) −h) 2 dx/L and the relative energy per unit length E ≡ F [h]/L − f(h) as a function of the domain length L. The primary instability is supercritical when the film is thin (h = 1.5), but becomes subcritical in thicker films (h = 7). For even thicker films (h = 15) the flat film is linearly stable, and no primary bifurcation exists. However, the film is metastable and branches of drop and nucleation solutions exist for L > L s .
In general one expects that these solutions persist on a substrate with weak wettability variation provided their period is compatible with the period L het of the wettability pattern. These solutions can then be continued numerically towards stronger wettability variation. As a result we expect solutions not only with L = nL het , n = 1, 2, . . . (studied in [30]), but also some that result from solutions on the homogeneous substrate with L = L het /n for n = 2, 3, . . . and L > L c . Moreover, in the subcritical and metastable case L > L s , and additional solution branches result from the nucleation solutions [30]. For the wettability profiles (5) all of these branches come in pairs as splits the degeneracy between h(x) and h(x + L/2). Note that there is also one (nontrivial) solution branch corresponding to a flat film on a homogeneous substrate. Figure 4 shows the solution families for a film on a substrate with a sinusoidal wettability pattern of period L het and system size L = L het in the absence of driving. Here, the critical period for instability is L c = 15.5. When L = L het = 25 the heterogeneity selects two branches from the circle of drop-like solutions (related by translations modulo the period) present when = 0, one of which extends towards larger values of while the other annihilates with a branch arising from the uniform film solution (||δh|| = 0) on the homogeneous substrate. These branches all consist of a single drop per period. In contrast when L = L het = 50 branches of one, two and three drops per period are present, generated by the splitting of the corresponding solutions on the homogeneous substrate. As in the L = L het = 25 case these branches annihilate pairwise with increasing , leaving only the one-drop branch at larger heterogeneities. Examples of solutions on these branches are shown in figures 5(a,b). Figures 5(c,d) show the single drop profiles that survive to larger heterogeneities. We mention that on the smaller domain (L = 25) the mesoscopic contact angle increases with increasing amplitude of the heterogeneity, with the opposite trend for L = 50. Evidently, in the former case the drop is so large that it extends into the hydrophobic region and is therefore compressed, an effect that increases with . In contrast, on the larger domain the drop does not encroach on the hydrophobic area and so increasing lowers the mesoscopic contact angle. The picture is qualitatively similar for non-sinusoidal wettability patterns, with one difference: if the wettability pattern breaks the symmetry → − , x → x + L het /2 the solutions will depend on the sign of and the corresponding branches become asymmetric with respect to = 0. Figure 6 presents the overall characteristics of the families of solutions on heterogeneous substrates in the absence of driving. It gives results for different nonsinusoidal but symmetric wettability patterns with profiles given by Eq. (5), as well as for local wettability defects breaking that symmetry, given by Eq. (6).

Driven system
With lateral driving no steady drops exist on homogeneous substrates, i.e., substrates with no physical or chemical defect, and all drops slide with a velocity determined by the driving, liquid volume, wettability and the viscosity of the liquid. However, sliding drops that appear stationary in an appropriately moving frame, i.e., that slide without change of shape, are present [2,68]. The situation changes on heterogeneous substrates where drops can be pinned by inhomogeneities. Pinning arises when the downstream driving force is balanced by an upstream force along the substrate generated by the gradient of P (x). However, this balance requires the knowledge of both the drop shape and its position with respect to the heterogeneity. For very small driving it is always possible to find such drop-like solutions. For larger driving, however, no stationary drops may exist and instead sliding drops or a flow with a stationary profile modulated by the heterogeneity will be present. More complete information can be gathered by plotting families of stationary solutions as a function of the driving strength, as done in figure 7 for different wettability contrasts and heterogeneity profiles. Three main types of such families can be distinguished. (i) A monotonic transition with increasing µ between drop-like and film-like states, as occurs, for example, in figure 7 (b) for = 0.5. (ii) A non-monotonic transition between these states with increasing or decreasing µ, with multiple solutions present for µ sn2 < µ < µ sn1 , as in figure 7 (a) for = 0.5. (iii) A discontinuous case, consisting of three branches. Of these the branches of drop-like and nucleation solutions terminate at a finite µ sn1 in a saddle-node bifurcation, while the third branch consisting of film-like states is disconnected and exists for all µ. Consequently as µ increases one finds an abrupt transition from drop-like to film-like states, although the reverse transition is absent. An example is given by the = 0.2 curve in figure 7 (a) or the = ±0.4 curves in figure 7 (c,d).
The above discussion ignores the stability of these branches. One finds that for the supercritical values of L considered here the film-like state is always unstable, and hence the transition at the rightmost saddle-node necessarily leads to a time-dependent state. This is the depinning transition that is of particular interest in the present paper. of a prominent shoulder immediately prior to depinning. As this occurs the downstream (advancing) contact angle decreases, while the upstream (receding) angle also decreases, as the drop stretches markedly in extent. Panel (b) studies the influence of a hydrophobic defect that blocks the drop at its downstream end and leads to a strong increase of the downstream (advancing) contact angle. The upstream (receding) angle also increases as the drop becomes more and more compressed.

Linear stability of the steady states
To determine the nature of the transition leading to depinning we must examine the linear stability properties of the profiles h 0 (x) shown in figure 8. For this purpose we linearize the evolution equation (7) about h 0 (x) and look for solutions of the form h 1 (x, t) = h 1 (x) exp(βt). This procedure leads to the eigenvalue problem βh 1 (x) = L[h 0 (x)]h 1 (x) for the growth rates β and the associated instability modes h 1 (x). Note that L depends nonlinearly on h 0 and its derivatives (for a detailed example related to front instabilities see [69]). The solution of this problem yields the spectrum of eigenvalues as a function of the driving strength µ. Figure 9 shows the result for a sinusoidal heterogeneity corresponding to particular branches in figure 7 (a). When these are discontinuous only the results for the branch of droplike solutions are shown. For small values of the heterogeneity amplitude the upper drop-like state remains stable until the first saddle-node bifurcation where a real eigenvalue becomes positive ( figure 9(a,b)). For larger some of the real eigenvalues collide pairwise and move into the complex plane ( figure 9(b)). In this case only the real part of the complex pair is shown; such complex eigenvalues are indicated by dashed lines. For yet larger a pair of saddle-node bifurcations appears on the solution branch forcing the leading eigenvalues to undergo a 'drunken dance' as a function of µ. This is a consequence of the fact that at each saddle-node a real eigenvalue must pass through zero. Figure 9(c) shows that the collision of the two unstable real eigenvalues on the film-like portion of the solution branch leads to a pair of complex eigenvalues with positive real part, reflecting the oscillatory instability of the film-like portion of the branch at larger µ. For yet larger complex eigenvalues appear even on the upper drop-like portion of the branch (figure 9(d)), but this branch remains stable until the saddle-node bifurcation. As a result the leading complex eigenvalues must become real prior to the saddle-node bifurcation (figure 9(d)). Note that despite the presence of complex eigenvalues loss of stability (i.e., depinning) occurs via saddle-node bifurcations where a real eigenvalue passes through zero.
The situation becomes more interesting when the heterogeneity is non-symmetric as shown for a hydrophobic defect in figure 10. For moderate values of the wettability contrast (figure 10(a)) depinning continues to occur via a saddle-node bifurcation. However, with increasing these saddle-nodes annihilate in a so-called hysteresis bifurcation (figure 7), and the eigenvalues are no longer forced to become real. In this case a new phenomenon arises: the drop-like solutions on the upper branch lose stability at a Hopf bifurcation (figure 10(b)). Thus in this case the pinned drops start to 'rock' prior to the depinning transition. The nature of the depinning transition must therefore be quite different from the case where depinning occurs at a saddle-node bifurcation. Also of considerable interest is the behavior of the eigenvalues past the Hopf bifurcation (µ > µ H ): with increasing an interval of µ values appears with real eigenvalues (figure 10(c)), and these come closer and closer to the Hopf bifurcation as increases (figure 10(c-e)) and eventually result in a pair of saddle-node bifurcations (figure 10(f)). However, the Hopf bifurcation remains.
It is significant that one of the real eigenvalues involved in this process increases very abruptly with µ. This behavior reflects the sudden and abrupt loss of stability once the drop depins, and indicates that the drop accelerates rapidly from its equilibrium state prior to depinning, in accord with everyday observation.
The nature of the depinning process is captured in the eigenfunctions responsible for the loss of stability, i.e., the eigenfunctions at the onset of depinning, as well as the eigenfunctions associated with the large growth rates identified in figure 10. In figure 11(a) we show the eigenmode corresponding to a close to maximum (real) growth rate in figure 10(e). The mode peaks in the downstream direction of the drop, and has a prominent dip on the upstream side of the drop apex, clearly indicating that the loss of stability is associated with the transfer of liquid from the upstream to the downstream side, in other words, that the instability produces downstream motion. Indeed, the mode resembles in form the function h 0 (x) that corresponds to the translation mode on a homogeneous substrate. In contrast, figure 11(b) shows the oscillatory eigenmode at the first Hopf bifurcation in figure 10(d). Like the real eigenmode in figure 11(a) the modulus of the associated complex mode |h 1 (x)| also peaks on the downstream side of the drop. On the other hand there is a secondary (albeit smaller) maximum on the upstream side, suggesting the onset of asymmetric rocking motion of the drop. This motion is visible in the perturbation |h 1 (x)| cos(β i t + φ(x)) of the steady profile h 0 (x) shown in figure 12(b) in the form of a space-time plot. Figure 12(c) shows that farther from onset the rocking becomes more nonuniform and that its maximum amplitude moves downstream, while figure 12(d) shows the behavior associated with the second Hopf mode that sets in at µ ≈ 0.18 (see figure 10(d)). Finally, figure 12(a) shows the behavior at the Hopf bifurcation in figure 10(b). Figures 13 and 14 summarize the location of the saddle-node bifurcations µ sn1 , µ sn2 in the (µ, ) plane (heavy lines), together with the location of both primary and secondary Hopf bifurcations (dashed lines). The shaded regions correspond to three equilibria; only one equilibrium is present outside these regions. The regimes where stable pinned drops are present are indicated; these lose stability at a saddle-node bifurcation across the heavy line, and at a Hopf bifurcation across the dashed line. Note that contrary to appearance the line of Hopf bifurcations does not emerge from the cusp where the two saddle-node bifurcations annihilate. The figures also show that for a sinusoidal defect the saddle-nodes persist with increasing , and no Hopf bifurcations are present. Figures 15 shows the contact angles at the front and back of the pinned drop as a function of increasing driving strength µ. In the following we refer to these as the advancing and receding (equilibrium) angles, respectively. For a drop pinned at the back (left) the advancing [receding] angle decreases [increases] for small but increasing µ ( figure 15(a)). However, once the drop starts developing a shoulder at the back the receding angle decreases again until depinning occurs (at a saddle-node bifurcation). The situation is different for a drop pinned at the front (right). In this case both angles increase with µ almost until depinning but drop just prior to it ( figure 15(b)). The only decrease in the receding angle occurs at very small driving.
The advancing and receding angles at depinning shown in figure 16 provide a measure of (a) the contact angle hysteresis observed macroscopically. In the case of a hydrophobic defect at the front ( > 0) both angles increase nearly linearly, and continue to do so even for oscillatory depinning ( ¡ 0.6); the small hook visible in the figure near this transition is consistent with the fact that the Hopf bifurcation sets in prior to the disappearance of the saddle-node bifurcations. The behavior is more intricate when the pinning is by a hydrophilic defect at the back ( < 0). In this case the role of the two angles is reversed, and both decrease nearly linearly with slopes identical to those in the > 0 case. For −0.2, however, the receding angle reverses tendency and starts to increase again, while the advancing angle continues to decrease. This qualitative change in behavior is similar to that occurring with increasing driving in figure 8 (a): the force drags the main body of liquid downstream but the spot of higher wettability traps part of it upstream. For fixed µ the latter effect becomes more pronounced as | | increases.  6)).

Unsteady states
In this section we describe the results obtained from solving the time-dependent problem (7)-(8) for different choices of the wettability pattern ξ(x), focusing on the dynamics that result from the depinning transition. We write the partial differential equation for the film thickness profile as a system of coupled ordinary differential equations for the film thickness at equidistant points in space, and impose periodic boundary conditions in space. The system is integrated in time using an explicit variable order, variable time step Adams method (NAG library routine 'd02cjc' [70]). Simulations are either started from slightly perturbed flat films or from steady solutions at nearby parameter values. After an initial transient the profiles settle down on a steady state or approach a time-periodic state.

Sinusoidal heterogeneity
The pinned states for a sinusoidal wettability profile are studied in figure 7(a). Figure 9 shows that depinning occurs via a saddle-node bifurcation and that beyond depinning the growing eigenmode becomes oscillatory if the wettability amplitude and hence µ at depinning are large enough. In figure 17 we show space-time plots of time-dependent solutions of (7)-(8) corresponding to figure 9(d) for (a) µ = 0.042 and (b) µ = 0.08. Both motions are strictly periodic in time. The former is very close to the depinning transition, at a value of µ for which the resulting motion has a very large period. The figure shows that the drop remains near equilibrium for an extended period before showing a slight drift in the downstream direction. This drift then accelerates rapidly as the drop translates to the vicinity of the next equilibrium position one spatial period ahead. The overall dynamics resembles stick-slip motion. Figure 17(b) shows the corresponding behavior farther from the depinning transition, i.e., for a larger driving strength µ. We see that the drop leaves its equilibrium position more rapidly, and thins down as it flows to its next equilibrium position. The dynamics is more fluid, and no longer resembles stick-slip motion. Figure 18 extends the above results to different values of µ. We see that in terms of ||δh|| the time for the drop to leave its equilibrium position as µ becomes larger becomes comparable to the time it takes to return to such a position. As a result the depinning becomes less sudden. The same trend is revealed when we track the location of the maximum height of the drop ( figure 18(b)). The figure shows that in all cases the maximum moves downstream, rather abruptly for µ close to the depinning transition, more slowly for large µ. The vertical jumps in the figure indicate that the drop has moved one spatial period, and the motion continues with x max at x = 0. Figure 19 shows a bifurcation diagram presenting on the same plot the amplitude ||δh|| for the equilibria ( figure 7(a)) and the time average of this quantity for the time-dependent states after depinning. The figure shows that these time-dependent states do indeed bifurcate from the saddle-node in a sniper bifurcation, and confirms that the period of the resulting periodic motion does diverge as the inverse square root of the distance from depinning as in the Adler equation (1).

Hydrophilic defect
As shown in figure 13(b) depinning from a hydrophilic defect also occurs via a sniper bifurcation. It is not surprising therefore that the qualitative results for this case resemble those for the sinusoidal wettability profile. Figure 20 shows a typical bifurcation diagram: the large amplitude drop-like states remain stable until a saddle-node bifurcation at which the drop depins and periodic stick-slip motion ensues. The period of this motion diverges as Figure 21(a) shows the dynamics of the profile for µ beyond depinning. The figure confirms that for µ near critical the translation of the drop occupies a much smaller fraction of the period than for larger driving. Figure 21

Hydrophobic defect
As already mentioned the hydrophobic heterogeneity permits a distinct mechanism for depinning. In this case linear theory shows that the static drop may lose stability at a Hopf bifurcation as the driving strength µ increases. Figure    frequency at onset remains finite, and figure 23(a) shows that near µ H the drop spends a long time near its equilibrium configuration before an abrupt translation to its next equilibrium position takes place. As expected the time spent near the equilibrium position decreases with increasing µ, so much so that figure 23(d) shows a drop travelling at almost constant speed, only slightly modulated by the heterogeneity. Figure 24 shows these results in a different but more quantitative form. The figure shows that the position of the global maximum changes discontinuously at two points indicating that some liquid remains trapped behind the defect. that the Hopf bifurcation is in fact slightly subcritical. As a result we cannot follow the branch of time-dependent solutions all the way to the Hopf bifurcation. It should be mentioned that with both increasing and decreasing the branch of static drops develops a pair of saddle-node bifurcations via a hysteresis bifurcation (Figures 13, 14). However, in both cases the Hopf bifurcation occurs above the upper saddle-node bifurcation. With decreasing wettability contrast the Hopf bifurcation moves towards this bifurcation, and disappears on colliding with it. Thus the oscillatory instability disappears in a Bogdanov bifurcation. This Bogdanov bifurcation is non-standard, however, since it involves the collision of a Hopf bifurcation not with a saddle-node but a sniper. This is because the spatial periodicity of the problem requires that the problem be posed in a phase space that is a cylinder (see Section 5). For smaller wettability contrasts depinning occurs via the sniper bifurcation, as shown in figure 25 (a) for = 0.4.

Dependence on loading
Finally, we also explore the dependence of the depinning transition on the loading of the defect. We use the spatial period L as a measure of the volume of liquid in the system, and hence of the amount of liquid that has to be held in the pinned state. Figure 26 shows the loci of the sniper bifurcation for a slightly hydrophobic defect in the (µ, L) plane. The figure shows that as the loading L increases the driving strength µ at which depinning occurs decreases (µ ∼ L −1 , at large L), and that µ is smaller for smaller wettability contrast. The presence of the cusps is notable, and indicates very high sensitivity of the driving strength at depinning to the loading in this regime. When L is too small (for the present parameter values L 20, a value that depends only weakly on the wettability contrast) no drops form, and instead the solution takes the form of a thin but stationary film whose thickness is modulated by the underlying heterogeneity.

Depinning mechanism
The depinning transitions via the sniper and Hopf bifurcations interact in parameter space near the location where the two saddle-nodes annihilate. In this region the Hopf bifurcation interacts with one of the saddle-node bifurcations producing the codimension-two Bogdanov bifurcation; in addition the second saddle-node is annihilated by the first in a hysteresis bifurcation, indicating that a third parameter must be tuned as well. Thus it is natural to consider a codimension three point in parameter space where the two saddle-nodes annihilate and at the same time interact with a Hopf bifurcation with zero frequency. The vicinity of this codimension-three point provides the key to the various transitions observed in our computations. These are captured by the system where A > 0 is a constant, and the parameters ν, µ and λ are unfolding parameters and hence small. Here the variable θ is a periodic variable representing for example the behavior of the maximum of the film profile modulo the spatial period L of the system. Since µ and λ are small the equilibria of interest will have small θ. Consequently we may locally replace equation (12) bÿ In this equation it is helpful to think of θ as O( ), and take the time derivative to be also O( ), . Of these we think of ν, λ as fixed, and take µ as the distinguished bifurcation parameter; this parameter plays the role of the driving strength, and increases through zero from left to right with µ = 0 at the codimension-three point.
In the following we denote the equilibria by θ 0 , and note that there may be as many as three. These annihilate pairwise in saddle-node bifurcations at λ + 3θ 2 0 = 0, and the saddlenode bifurcations themselves annihilate at µ = λ = 0. Thus three equilibria are present for λ < 0 and one for λ > 0. Of these the leftmost equilibrium is stable. When λ < 0 a real eigenvalue passes through zero at the rightmost saddle-node bifurcation, followed by a further loss of stability at the leftmost saddle-node. Thus the rightmost equilibrium is always doubly unstable. When λ > 0 the stability assignments for large |µ| are preserved but since no saddle-node bifurcations are present it follows that the change in stability must occur via a Hopf bifurcation. To understand how the latter fits into the picture we examine the stability of the equilibria in more detail.
The linear stability of these is given by the eigenvalue equation Steady state bifurcations (s = 0) correspond to the saddle-node bifurcations already mentioned; a Hopf bifurcation (s = ±iω) occurs when ν + Aθ 0 = 0 and requires λ + 3θ 2 0 > 0. The latter condition is satisfied in the region λ + 3ν 2 /A 2 > 0, i.e., for λ < 0. The Hopf bifurcation has zero frequency when λ + 3θ 2 0 = 0. Thus the curve λ + 3ν 2 /A 2 = 0 represents the location of the Bogdanov bifurcations. The corresponding values of the bifurcation parameter µ follow from µ+λθ 0 +θ 3 0 = 0. Note that the Bogdanov bifurcation coincides with the annihilation of the saddle-nodes only when ν = 0; when ν > 0 the Bogdanov bifurcation occurs on the lower equilibrium branch and vice versa. The calculation presented in figure 13 indicates that the value of ν in our problem is quite small; the degenerate Hopf bifurcation (ω = 0) almost coincides with the cusp at which the saddle-node bifurcations annihilate. Equation (12) can be used to follow the solutions of equation (13) away from µ = 0, λ = 0. One finds that as µ increases for λ < 0 the rightmost saddle-node bifurcation is indeed a sniper, provided that ν > 0. This bifurcation produces a limit cycle that is a rotation, i.e., θ advances by 2π each period. In contrast, the Hopf bifurcation creates oscillations with no net change in θ in one period. Oscillations of this type correspond to libration. Since sliding motion necessarily corresponds to rotations, it follows that a global bifurcation must take place that changes librations into rotations; such a bifurcation must involve an unstable fixed point. In the case with only one fixed point the most likely candidate for such a bifurcation involves the collision of the zero-mean oscillation with the now unstable fixed point. Such a global bifurcation requires a three-dimensional system, however, and hence cannot occur within (12). It does occur in other systems such as the Takens-Bogdanov bifurcation with O(2) symmetry [71], case II−.

Metastable flat film
The above results were obtained for liquid volumes, as measured by the spatial period L, for which the flat film is unstable. For these values of L no steady profile exists other than pinned drops. As a result beyond depinning the thickness profile must vary in space and time. The situation may be different, however, for situations where the mean thickness corresponds to a metastable flat film, i.e., the locally stable flat film coexists with stable pinned drops, and is separated from them by unstable nucleation solutions. This is the situation we focus on in the present section.
In the absence of driving the typical situation on a homogeneous substrate is as shown in figure 3 (dashed line). For L > L s two branches are present, the upper corresponding to stable drops, and the lower to unstable nucleation solutions. In addition there is a branch of stable flat states given by ||δh|| = 0 with zero energy. Small heterogeneity with L het = L splits the former solutions into two each, and perturbs the flat state into a corrugated state. Thus, for = 0 there will now be five distinct states, instead of the three states when the primary instability is supercritical. The situation is illustrated in figure 27 which shows the resulting five branches as a function of for both sinusoidal (figure 27a,b) and localized (figure 27c,d) wettability defects. In a driven system these solutions annihilate pairwise as the driving strength µ increases. Figure 28 shows that for a sinusoidal wettability pattern with a small the nucleation branches are the first to disappear, followed at larger µ by the drop-like states. As a result for large µ only the nearly flat state remains. However, when is larger the larger amplitude nucleation state annihilates the lower drop-like state, reducing the number of solutions at given µ first to three, and then with increasing µ to one (see figure 28(a) for = 0.6). Indeed when localized defect (see figure 28(b)) five solutions are present at given > 0 not too large (see figure 27(c)), of which the middle two annihilate first as µ increases; the next two can persist to quite large driving strength. In contrast, when > 0 is large ( ¡ 2.1) only three solution branches are present, and the top two annihilate with increasing µ leaving the nearly flat film (not shown).
The most important difference from the case studied in the preceding sections lies therefore in the presence of a stable, almost flat film at larger driving strength. As a result the depinning at the rightmost saddle-node does not automatically lead to a periodic orbit; instead the solution may jump to the stable nearly flat state. In this case the drop does not slide downstream but is adsorbed into a flowing film with a steady surface profile or one that is modulated by surface waves. States of this type are, however, difficult to compute numerically since they are connected to very large values of and/or µ. For instance, the modulated flat film state and the 'last' nucleation state at µ = 0 annihilate, for the parameters of figure 28 (b), at = 2211! Such large values lie outside the range of validity of evolution equations derived using lubrication theory.
In the present case, however, the above possibility does not arise and one again finds stick-slip motion beyond the saddle-node bifurcation. Figure 29 shows the corresponding bifurcation diagram as a function of µ, while figure 30 shows the position of the drop maximum, the L 2 norm and the profiles as a function of time over slightly more than one period in time. The sharp first peak in the L 2 norm (figure 30(a)) coincides with the accelerating motion due to the passage of the receding contact line over the hydrophobic defect; the main peak arises when the advancing contact line becomes trapped behind the next defect. The small oscillations superposed on the main peak are a numerical artifact related to 'stick-slip' motion on the numerical grid.

Conclusions
In this paper we have explored the process of pinning and depinning of driven drops on heterogeneous substrates. For this purpose we have formulated the depinning process as a bifurcation problem, and focused on a generic problem of this type. This approach permitted us to restrict attention to drops not much higher than the wetting layer. For simplicity we considered the case of a heterogeneous interaction between the free interface and substrate with a well-defined spatial period. Such an interaction might arise from spatially varying wetting properties, or heterogeneous temperature or electrical fields, while the driving might be caused by gravitational or centrifugal forces, or gradients of wettability, temperature or electrical fields. Appendix A lists selected results for a drop pinned in a heterogeneous electric field in a condenser that highlights the close resemblance between this system and the generic problem studied in the main body of the paper. We have examined two types of pinning: pinning by a hydrophilic defect at the back of the drop, and pinning by a hydrophobic defect in front of it, and identified two mechanisms whereby depinning takes place. In the case of a hydrophilic defect the drop stretches markedly just prior to depinning as the driving strength increases; the pinned drop loses stability at a saddle-node bifurcation, resulting in periodic motion as the drop slides over the periodic array of hydrophilic defects. We have referred to this type of bifurcation as the 'sniper'. The periodic motion that results is slow when the drop is stretching, and fast once the drop breaks away from a defect and spills onto the next one. The situation is richer for hydrophobic defects that pin the drop by blocking it. In this case in addition to the steady state sniper bifurcation a new depinning mechanism was observed: the drop loses stability to an oscillatory mode leading to back and forth oscillation of the drop. We believe (but have been unable to demonstrate) that this bifurcation must be supercritical in order to form a homoclinic connection with an unstable drop, and allow the drop to spill over the blocking defect. However, we have been unable to find the (stable) zero-mean oscillations that must be created at such a Hopf bifurcation, possibly because the anticipated global bifurcation occurs very close to the primary Hopf bifurcation. Note that such a Hopf bifurcation is only present because of the flow through the precursor. For larger drops this flow becomes negligible and consequently the dynamics associated with the Hopf bifurcation must all occur in a small range of µ near this bifurcation. Thus the motion that results from oscillatory depinning is in fact very similar to that produced via the sniper bifurcation; the main difference between these two depinning scenarios is confined to behavior very close to the depinning transition. In the sniper scenario the average speed of the drop vanishes like (µ − µ c ) 1/2 , while in the Hopf scenario it drops to zero logarithmically at the global bifurcation separating the rocking motion from sliding. Of course, the average speed at the saddle-node bifurcation that occurs as µ decreases ( figure 25(b)) remains finite.
In the paper we have examined the dependence of the above scenarios on the parameters of the system specifying the profile of the heterogeneity, as well as on the drop volume as measured by the mean thicknessh, and the heterogeneity period L. Not surprisingly we found that for larger volumes depinning occurs for smaller heterogeneity amplitude | | (figure 26). This is simply an effect of the weight of the liquid in the drop.
The results of this paper clarify several issues in the theory of liquid films within the continuum lubrication description. On a homogeneous substrate the use of a disjoining pressure defines the mesoscopic contact angle; this is the contact angle that is measured in a laboratory setting. This angle is not specified a priori in the theory, and can only be determined by solving the liquid film equation. As a consequence the downstream and upstream contact angles in driven systems are in general different. Since the problem posed on a homogeneous substrate is invariant under translations, generic spontaneous symmetry-breaking (i.e., dropforming) instabilities are necessarily Hopf bifurcations and lead to traveling (i.e., sliding) drops. Thus on a homogeneous substrate all nonuniform states slide with respect to the substrate, and one concludes that the precursor film formulation of the problem does not permit pinned drops under lateral driving [2,52]. In other approaches the precursor is absent and the downstream contact line location and angle is specified a priori. Such theories are singular in the sense that the film thickness vanishes outside of the drop, and it is this singularity that mathematically speaking violates the theorem and permits stationary drops. In our formulation such singularities are absent and we must rely on heterogeneities to pin driven drops. In this paper we have examined in detail the processes that lead to such pinning, as well as the processes that lead to depinning as the strength of the driving is progressively increased. It is unclear to us how depinning would be described within theories that prescribe the contact line.
Although our results are applicable, in principle, to a wide spectrum of systems, in some of these the required parameter values may be hard to realize experimentally. The spectrum of possible applications ranges from (i) nano-drops on ultrathin precursor films of thicknesses l ≈ 1 . . . 10 nm where the film-substrate interaction is dominated by wettability (i.e., a disjoining pressure) to (ii) micro-drops on macroscopic wetting layers with l ≈ 10 . . . 100 µm. In the latter case the film-substrate interaction may combine repulsive van der Waals interactions and the influence of temperature and/or electrical fields [23,51,52]. In principle, the spatial modulation responsible for pinning can be present in any part of this interaction. We estimate the nondimensional driving force µ using Eq. (3) based on the wetting layer thicknesses l and the typical energy density κ = A/l 3 for a film of thickness l, where A is the Hamaker constant characterizing the van der Waals interaction (note that at the precursor film thickness the typical energy densities of the two antagonistic terms are roughly equal). For a droplet on a substrate inclined at a (small) angle α to horizontal µ = l 5 γ 1/2 A −3/2μ [see Eq. (3)], whereμ = ρgα ≈ 10 4 α×N/m 3 , A ≈ 10 −20 Nm and γ ≈ 10 −1 N/m, we find that µ = l 5 10 33 α × m −5 . Thus for l = 1nm, 10nm, 100nm and 1µm the quantity µ/α is 10 −12 , 10 −7 , 10 −2 and 10 3 , respectively. This result should be compared with the typical value µ ≈ 10 −2 at the transition between the two types of depinning from the heterogeneities considered in this paper. It follows that gravity-driven depinning via a Hopf bifurcation can only be observed for wetting layer thickness of about 100 nm or more, i.e., in case (ii). A specific example is given in Appendix A. In case (i) gravity-driven depinning of droplets will always take place via a sniper bifurcation.

Appendix A. Dielectric liquid in a condenser
In this Appendix we present selected results for a drop of dielectric liquid in a condenser of gap width d. The drop is pinned by a spatial heterogeneity. The evolution equation for the thickness profile h(x, t) of a film of dielectric liquid of relative permittivity ε in air (relative permittivity equal to one) can be derived from the general two-layer result in Ref. [23]. In the nondimensional variables used in the body of the paper the equation for h(x, t) is (A.1) The pressure P (h, x) contains van der Waals interactions with both condenser walls and the electrostatic pressure. The heterogeneity is assumed to be in the electrical field only, i.e., P (h, x) = [p up vdW + p low vdW + b (1 + ξ(x)) p el ]. (A. 2) The disjoining pressure models wetting interactions that prevent rupture of the liquid film at the lower plate and rupture of the air film at the upper plate, respectively.
For simplicity we assume these to be of equal strength. As in the main body of the paper we base the energy density scale κ on van der Waals interactions, i.e., κ = A/l 3 for wetting layer thickness l. For a vertical electrical field and a layer of dielectric liquid of permittivity ε under air the dimensional vertical field components in the liquid and the air are and , (A. 6) respectively. In the lubrication approximation the horizontal components can be neglected. The electrostatic pressurep el ≡ ε 0 (ε − 1)E liq E air [21,23], where ε 0 denotes the permittivity of vacuum, provides the dimensionless pressure and the dimensionless relative interaction strength where ε > 1. The direction of the applied voltage has no influence on the result. With ε = 2 (a typical value for silicon oil [23]) and the estimate for κ from section 7 the relative interaction strength is b = 10 9 U 2 l×V −2 m −1 , i.e., for wetting layer thickness between l = 100nm and 1µm and U = 1V, b is in the range 100 to 1000. Figure A1 shows a family of steady solutions for b = 200. The transition between depinning via sniper and depinning via Hopf bifurcation occurs at a driving strength µ ≈ 1, well in the range of gravitational driving. This is also the case for two layers of immiscible liquids of different permittivities and viscosities.