{'fraction_non_alphanumeric': 0.0331696672122204, 'fraction_numerical': 0.009012547735951991, 'mean_word_length': 4.773718029482172, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Morphological transitions of sliding drops: Dynamics and bifurcations

We study fully three-dimensional droplets that slide down an incline by employing a thin-film equation that accounts for capillarity, wettability, and a lateral driving force in small-gradient (or long-wave) approximation. In particular, we focus on qualitative changes in the morphology and behavior of stationary sliding drops. We employ the inclination angle of the substrate as control parameter and use continuation techniques to analyze for several fixed droplet sizes the bifurcation diagram of stationary droplets, their linear stability, and relevant eigenmodes. The obtained predictions on existence ranges and instabilities are tested via direct numerical simulations that are also used to investigate a branch of time-periodic behavior (corresponding to repeated breakup-coalescence cycles, where the breakup is also denoted as pearling) which emerges at a global instability, the related hysteresis in behavior, and a period-doubling cascade. The nontrivial oscillatory behavior close to a Hopf bifurcation of drops with a finite-length tail is also studied. Finally, it is shown that the main features of the bifurcation diagram follow scaling laws over several decades of the droplet size.


I. INTRODUCTION
It is known from many experiments and everyday experience that drops on a smooth homogeneous solid substrate slide along the substrate when lateral driving forces are applied.
For heterogeneous substrates the driving force first needs to overcome pinning influences of the heterogeneities. Physical systems where individual or many sliding drops can be observed are, e.g., rain drops on a train window driven by the shear flow of the outside airstream, drops sliding or rolling down an incline due to gravitational force [1][2][3] or drops moving outwards on a spinning disk due to the centrifugal force. The motion of droplets may also be induced by chemical or thermal gradients along the substrate [4][5][6]. In all these settings, control parameters, such as, droplet volume, strength of the driving force, material properties of the liquid (viscosity, surface tension) and interaction properties of the combination liquid-substrate (equilibrium contact angle) determine the velocity and shape of the moving droplets.
Here, we mainly focus on sliding droplets of partially wetting liquids on an incline, i.e., droplets that form at equilibrium finite contact angles with the substrate and slide down under the influence of gravity. Both, experimental [3,[7][8][9] and theoretical [10][11][12][13][14][15][16][17][18] studies, analyzed the dynamic behavior and revealed interesting morphological changes with increasing driving force and/or droplet volume. Reference [3] employs small-contact-angle silicon oil droplets on an inclined fluoro-polymer coated substrate and describes a sequence of transitions that droplets of fixed volume undergo with increasing substrate inclination (see Fig. 1 of Ref. [3]): At small inclinations the droplets slide with constant speed and a constant oval-like shape of their base with a small front-back asymmetry (while at zero inclination they are shallow spherical caps). With increasing inclination angle the asymmetry between front and back increases, as does the droplet velocity. The front still resembles a circular arc but the back becomes increasingly pointed and develops a cusp at a critical capillary number (or non-dimensional droplet speed). Beyond the corresponding substrate inclination an instability occurs, where the sliding drop emits smaller satellite droplets at its back. This effect is also denoted as pearling [3], a notion we also adopt here. In parallel, the shape of the unstable droplet develops an elongated protrusion at its back. The emitted droplets can be of identical size or follow a rich variety of dynamical periodic patterns. In this context, Ref. [3] mentions a cascade of bifurcations that involve frequency divisions.
It is found that the dependence of the sliding speed on inclination is nearly linear with a change of the corresponding slope at the onset of pearling. Power law dependencies on volume are also studied. Similarly, Ref. [7] derives and experimentally verifies for small droplet speeds (assumption that the droplet nearly remains a spherical cap) a power law that also includes dependencies on other experimental parameters. The limit of nearly non-wetting drops rolling down an incline is investigated in Refs. [2,19].
It is important to mention that the experiments in Ref. [3] are actually not performed with individual sliding drops but with a continuous sequence of identical drops that are placed at a defined constant time interval at the top of the inclined plate. This does not affect the behavior below the threshold for the pearling instability. However, it implies that beyond the threshold one has pearling events where satellite drops split off the main drop and also coalescence events where main drops absorb satellite droplets left behind by the preceding main drop. This allows one to consider the pearling process together with coalescence events as a periodic process either in the fixed laboratory frame or in a frame moving with the main drops. This observation will below be used where we employ a long-wave model to analyze the morphological changes of the droplet shape and also the pearling-coalescence cycle. Also note that a similar change in droplet morphology (up to the appearance of a cusp) had earlier been observed for Marangoni-contracted water drops sliding down inclined hydrophilic substrates [1].
To model and simulate the sliding droplets, different approaches might be followed: One may employ the fully nonlinear Navier-Stokes equations describing the evolution of the full three-dimensional velocity field [20]. In most cases, describing slowly sliding drops, the linearized small-Reynolds-number form -the Stokes equation -may be used [21]. Both approaches have to be complemented with appropriate boundary conditions at the solidliquid and the liquid-air interface. Alternatively, one may combine the bulk equations with an additional phase-field dynamics, thereby modeling the liquid-air interface as diffuse [22,23].
Usually, these approaches involve a rather high computational cost, in particular, if threedimensional drops are considered. Other methods include dissipative particle dynamics [24] and lattice-Boltzmann methods [25].
Alternatively, one may employ an asymptotic model derived from the Navier-Stokes and continuity equations and boundary conditions for free surface flows for which all typical length-scales parallel to the substrate are large as compared to the film height, i.e., for droplets with small contact angles and surface slopes. This small-gradient or long-wave approximation results for three-dimensional drops in a highly nonlinear evolution equation for the film thickness profile h(r, t) as a function of the two coordinates r = (x, y) along the planar substrate [26,27]. In the present case the resulting equation accounts for capillarity through a Laplace pressure and for wettability through a Derjaguin (or disjoining) pressure [28,29], that combines two antagonistic power laws [30,31]. When adding the lateral driving force due to gravity, they serve to study droplets sliding down an inclined homogeneous substrate. Such long-wave models with various forms of Derjaguin pressure were already used to thoroughly analyze morphological transitions of stationary sliding two-dimensional drops (2d), i.e., with one-dimensional profiles h(x, t) [11,12]; and for a few time simulations of three-dimensional droplets (3d) that show the pearling instability [12,16,32] or states close to pearling [33]. Such models are also employed to study the combined spreading and sliding of perfectly wetting drops [26,34].
Note that we neglect any roughness of the substrate beneath the sliding droplets. Therefore, any contact angle hysteresis known from experiments is excluded from the dynamics that we observe, and there is no finite onset inclination angle at which the droplets start to move. The related depinning behavior is observed with similar long-wave models where modulations of the Derjaguin pressure represent regular arrays of wettability defects in Refs. [35] (2d) and [36,37] (3d). Long-wave equations are also employed to study droplet spreading over random topographical substrates [38]. The influence of regular and random substrate modulations on the onset of droplet sliding is also studied with lattice-Boltzmann methods in Refs. [39] and [25], respectively. Here, we always assume an ideally smooth substrate as can, e.g., be realized through atomically flat single crystal surfaces or through liquid-infused substrates.
The present work employs such an asymptotic long-wave model to investigate the morphological changes of fully three-dimensional sliding drops as observed in Ref. [3]. This is done employing pde2path, a continuation and bifurcation package for elliptic systems of partial differential equations (PDEs) [40,41]. In particular, we study the behavior of sliding droplets at several fixed drop volumes as the inclination angle of the substrate is changed that represents the strength of the driving force. Our analysis shows that a family of stationary sliding droplet solutions exists that undergoes steady and oscillatory bifurcations.
These bifurcations are related to morphology changes as, e.g., the development of protru-sions at the droplet back (related to a pair of saddle-node bifurcations) and the emergence of satellite droplets in a pearling instability that technically is a global (homoclinic) bifurcation. Beside the bifurcation diagram of stationary droplets we obtain their linear stability and the relevant eigenmodes. We obtain predictions on existence ranges and instabilities that are tested via direct numerical simulations that are also used to investigate the branch of time-periodic behavior (corresponding to pearling-coalescence cycles) which emerges at the global instability, the related hysteresis in behavior, and a period-doubling cascade that occurs when further increasing the inclination. The non trivial oscillatory behavior close to a Hopf bifurcation of drops with finite-length tails is also studied. Finally, it is shown that the main features of the bifurcation diagram are universal as they follow scaling laws over several decades of the droplet size. The quantitative changes are condensed in the form of scaling laws. Finally, we consider the total dissipation of the different droplet types and analyze where dissipation is localized within the sliding droplets.
The paper is structured as follows. In Sec. II we present our model and discuss the numerical approaches of time-stepping and path-continuation. The sub-sections of Sec. III present our main results, first, in the form of the overall bifurcation diagram and, then, detailing the behavior on the various found sub branches of stationary sliding drops. Sec. IV analyses the dependence of main features on the droplet volume while Sec. V discusses the velocity field and dissipation within the drops. The work concludes with a discussion in Sec. VI.

A. Governing equation
We describe the dynamics of sliding drops using an asymptotic long-wave model, i.e., an evolution equation for the film thickness profile. The model accounts for capillarity through a Laplace pressure, wettability through a Derjaguin (or disjoining) pressure and the down-hill component of gravitation as lateral driving force [12,16,26]. It is derived from the Navier-Stokes equations using a long-wave approximation [26,27] employing a no-slip boundary condition at the smooth solid substrate and stress-free boundary conditions at the liquid-air interface, i.e., the free surface of the film. This approach is valid in the case of thin liquid layers that have small surface slopes, i.e., for liquid films and drops that are described by a height profile h(x, y, t) which only exhibits small gradients in the lateral directions. In other words, the surface deflections are small (or the entire film is thin) relative to the lateral dimensions of height modulations.
For a general lateral driving force χ(h), the non-dimensional evolution equation for the film thickness profile reads: where Q(h) = h 3 is the mobility function and Π(h) = −df (h)/dh is the Derjaguin pressure modeling wettability, i.e., the effective interaction between the free surface of the film and the substrate [28,42]. The term ∆h represents the Laplace pressure due to curvature. We employ the wetting potential f (h) = 1 5h 5 − 1 2h 2 that corresponds to a precursor film model, i.e., a macroscopic drop coexists at equilibrium with an thin adsorption layer of thickness h p = 1 determined by Π(h p ) = 0 [30,31]. Note that there are many other possible choices for the wetting potential and resulting Derjaguin pressure. Another common ansatz employs, e.g., exponential functions [12,32] or a combination of a power law and an exponential [43,44], however the resulting overall shape of Π(h) is rather similar. Instead of using a Derjaguin pressure that ensures the existence of a wetting layer, other approaches are possible, such as a slip-model that allows for a finite slip velocity of the liquid at the liquid-solid interface [45].
Although there exist quantitative differences to the model used here, time simulations show qualitatively similar results for individual sliding drops (see, e.g., Ref. [33]). Furthermore one can show that asymptotically these models are equivalent in describing moving contact lines [46]. However, for the current work, a slip model is less suitable because additional assumptions are needed to treat topological changes of drop solutions with compact support.
Such changes frequently occur during the studied pearling-coalescence cycles.
In the non dimensional form of the thin film equation in Eq. (1), all material parameters, such as viscosity, surface tension, and Hamaker constant, are already absorbed into the scaling of height, the coordinates, and time, resulting in the remaining dimensionless parameters (as described in note [47]). For details of the used scaling, we refer to Refs. [12,31,48,49].
Note that the scaling is chosen in such a way that each of the most important control parameters only enters the equation through one dimensionless parameter but is not contained in the scales. Here, these is the inclination of the substrate and the volume of the drop.
In the case of an inclined substrate, the influence of gravity is incorporated into the driving force χ(h) through the scaled long-wave inclination angleα: Here, G 0 = 10 −3 presents a dimensionless gravity parameter and the plane is tilted in the x direction; i.e., there is no driving force in the y direction. For simplicity, we always use α = G 0α as an effective inclination.
We note, that Eq. (1) with driving (2) may be expressed in the general variational form of a gradient dynamics on a free energy functional F [h]. Namely, The final term in F [h] corresponds to the lateral part of the potential energy. Note, that (4) has no lower bound, as the global minimum is only reached if all the liquid volume approaches x → ∞. Therefore, we have the interesting case of a system that remains out of equilibrium forever, but is still described by a gradient dynamics that evolves towards an energy minimum (it can be shown that F [h] is a Lyapunov functional). This dramatically distinguishes the case treated here from the relaxational dynamics discussed in Refs. [31,50] where the system moves towards equilibrium and approaches it.
To investigate droplets that slide with a constant shape and velocity U, i.e., stationary that determines the stationary states h = h 0 (x, y). To determine their stability, we add a small perturbation h 1 (x, y), i.e., use the ansatz h(x, y, t) = h 0 (x, y) + h 1 (x, y) exp(λt).
Introducing the ansatz into Eq. (1) and linearizing in h 1 gives (in the comoving frame) the linear eigenvalue problem where λ and h 1 represent the eigenvalue and eigenfunction, respectively. The linear operator The presented model has a variety of different solution types. A homogeneous film of any height is always a (trivial) solution. However, for the employed Derjaguin pressure flat films are unstable in most cases (linearly unstable for h > 2 1/3 ). Also, effectively one-dimensional solutions are possible, such as transversally invariant ridges that slide down the incline.
Their cross sections in the x direction are one-dimensional sliding drops, i.e., h = h(x, t).
Another option are rivulets that represent free surface channels guiding the liquid downslope.
They have a cross section in the ydirection independent of x, i.e., h = h(y, t). Here, we mainly analyze truly three-dimensional droplets, i.e., h = h(x, y, t), using both, continuation techniques and direct numerical simulations, which are briefly described next. scheme [54]. The resulting nonlinear problems are solved with a Newton method, using for the linear problems a biconjugate-gradient-stabilized method (BiCGStab) with a symmetric successive overrelaxation (SSOR) as preconditioner. To implement Eq. (1) using this ansatz, we split the model into two equations second order in space: w =∆h .
In a weak formulation for general test functions φ h and φ w the equations read The boundary integrals occurring due to the performed partial integrations vanish due to the used periodic boundary conditions. As initial conditions we use a parabolic cap of prescribed volume V above the precursor layer. When the time-dependent behavior of unstable solutions is studied, solutions obtained by continuation (Sec. II B 2) are used as initial conditions.
To obtain the velocity of the sliding droplets, we track the position of maximal height over time and differentiate numerically using fourth order finite differences. To better represent the results in the forthcoming plots of height profiles, the x coordinate of the solutions is shifted such that the maximum is always at the same position. This corresponds to a transformation into the frame moving with the numerically determined droplet velocity.

Path-continuation
To analyze the model (1) with continuation methods [55,56], i.e., to directly track solutions in parameter space, we make use of the numerical pseudo-arclength continuation package pde2path [40,57] which is based on the finite-element methods of matlab's pdetoolbox.
Basically, path continuation determines steady states of an ordinary differential equation PDEs including boundary conditions and side conditions as, e.g., volume conservation can always be approximated by an ODE system of large dimension. The primary continuation parameter is in our case, e.g., the inclination α of the substrate in Eq. (1). Naturally the continuation of sliding drops requires the use of the system in the comoving frame as given by Eq. (5), implying that beside the primary continuation parameter also the velocity of the frame has to be determined. In this way one may start at an analytically or numerically given solution, continue it in parameter space and obtain a broad range of solution families including their bifurcations and accompanying changes in morphology [55,56]. However, it depends on the continuation package whether only time-independent solutions can be continued or whether this is also possible for time-periodic states. As the latter is often not the case, direct numerical simulations are needed to complete the bifurcation diagram of a system.
In time simulations of (1), the volume is automatically conserved because the equation respectively. Here h old denotes the solution obtained in the previous continuation step.
Taking a single static drop with α = U = 0 as starting solution, we conduct an pseudoarclength continuation for increasing inclination angles α. The corresponding velocity U of the comoving frame is used as an additional free parameter, as the drops slide with varying velocities that depend on the inclination and are automatically adapted to the corresponding α during continuation.  From the main branch of stationary states, a side branch of pearling solutions bifurcates, that represents the emission of small satellite drops from the larger sliding drops. Due to the periodic boundary conditions in the direction of droplet motion, these solutions are timeperiodic in the comoving frame, which corresponds to the experimental setup in Ref. [3] where a continuous sequence of identical drops slides down an incline. The distance between the drops corresponds to the periodicity of the domain considered here.
Starting with a droplet at rest for a vanishing inclination (α = 0), the first part of the curve U(α) (denoted in the following as sub-branch I ) exhibits a linear increase of velocity with increasing inclination angle. The overall shape of the drop profile does not change significantly (cf. Fig. 1, snapshot I; for details see Sec. III A). For a specific value of α = α SN1 , which depends on the drop volume, a saddle-node bifurcation (or fold) occurs where the stable sub-branch I connects with a linearly unstable sub-branch II. The latter continues towards smaller α; i.e., it coexists with sub-branch I for a range of inclination angles and always exhibits smaller velocities. The shape of the droplets changes significantly along this sub-branch, as a protrusion is formed at the rear of the drop (cf. Fig. 1, snapshot II) that is the mesoscopic equivalent of cusps observed for macroscopic drops. This is further discussed in Sec. III C. Near the saddle-node bifurcation, on sub-branch II also a global bifurcation occurs, that is responsible for the pearling instability and the emission of small satellite droplets. The resulting non-stationary solutions are time-periodic and show a repeated breakup of the main drop into a large main droplet and a small satellite droplet and the subsequent coalescence of the two. This type of behavior is determined by direct numerical simulations and is denoted as side branch A. It is further analyzed in Sec. III B.
Following sub-branch II with decreasing α, at a lower inclination threshold α = α SN2 , a second saddle-node bifurcation occurs. There, a second real eigenvalue becomes positive.
Therefore beyond the fold, "sub-branch III" consists of doubly unstable drops and continues towards larger α. The two real unstable eigenvalues approach each other, become identical already slightly after the saddle-node bifurcation. Then they transform into a complex conjugate pair. Beyond this point, time simulations that start with such unstable steady drops on sub-branch III show breathing-like oscillatory destabilization and then approach a drop on the stable upper sub-branch I. Increasing the inclination angle further, the pair of complex conjugate eigenvalues crosses the imaginary axis, i.e., a Hopf bifurcation occurs. At this bifurcation, we expect a branch of time-periodic solutions to emerge, e.g., periodically oscillating droplets. However, the employed continuation package is not yet able to follow such a branch of time-periodic solutions. As our direct time simulations in this parameter range did not show such solutions, we conjecture that the branch of time-periodic solutions is a short (small α range) unstable branch that starts subcritically at the found Hopf bifurcation and ends in a homoclinic bifurcation (as sometimes seen in related systems; see e.g. Ref. [58]).
Note, that this hypothetical short branch should not be confused with the aforementioned side-branch A of periodic pearling and coalescence that emerges at a different homoclinic bifurcation. Beyond the Hopf bifurcation, we are on sub-branch IV that is formed by linearly stable droplets. The protrusion at the back of the droplets becomes even more significant in this region as an elongated tail is formed whose length increases with increasing inclination angle (cf. Fig. 1 snapshots III/IV). Note that depending on α this non-trivial droplet solutions coexists with either the stable drops of sub-branch I or with the stable time-periodic behavior on side-branch A. Increasing α even further, on branch IV the droplet elongates until it forms a rivulet solution in the x direction (see Sec. III D for details). In the next sections we further analyze the phenomena on the individual sub-branches.
A. Stable stationary sliding drops: Sub-branch I First we give some more details on sub-branch I, i.e., the part of the bifurcation curve in Fig. 1, that emerges from the origin. Due to the absence of substrate inhomogeneities or other pinning mechanisms, a drop at rest begins to slide for arbitrarily small finite values of α (cf. Refs. [48,59] and the present case in comparison to Refs. [35,37]). Inclining the substrate further, the drop velocity increases linearly while the shape of the drop does not vary significantly. The radially symmetric spherical cap shape at α = 0 becomes slightly oval as the droplet prolongs a bit in the direction of the incline. Further, one can observe that the x → −x symmetry is broken as the trailing back of the drop becomes a bit more pointed than the advancing front (cf. Fig. 1, snapshot I). The value of the slope of the linear dependence U(α) depends on the volume of the drop, and can be brought into the form of a simple power law. This is quantitatively analyzed below in Sec. IV. The final part of sub-branch I at α 1.5 × 10 −3 deviates from the linear behavior towards lower velocities until it becomes vertical where it ends in the first saddle node at α SN1 ≈ 1.85 × 10 −3 . In this final part of the sub-branch the shape asymmetry of the droplets gets stronger. The dependencies of the critical angle α SN1 and the velocity U SN1 at the bifurcation on droplet volume can also be brought into the form of power laws as discussed below in section IV. As initial condition for the DNS, we use the stable solutions on sub-branch I below but close to α SN1 , and increase the inclination angle to an α > α SN1 above the critical one. i.e., the behavior corresponds to a limit cycle solution. At first view, this periodic behavior might seem artificial as it is induced by the periodic boundary conditions. However, these conditions are appropriate also for the real experimental system [3] where one periodically To be able to include the branch of these time-periodic solutions into the same bifurcation diagram in (Fig. 1) that shows the stationary sliding drops, we compute the velocity of the larger drop averaged over ten periods T per , for which we track the position of maximal drop height over time and differentiate it numerically.  which is typical for an onset of a periodic cycle via a homoclinic bifurcation [60]. This hypothesis is also supported by the fact that there exists some hysteresis, i.e., in a small interval of inclination angles α bif < α < α SN1 both, the time-periodic solutions and the stationary drops of sub-branch I are stable (see also Fig. 7 below). This zoom also clarifies that the time-periodic branch ends on sub-branch II of unstable droplets as expected for a homoclinic bifurcation.
For larger inclination angles α 2.1 × 10 −3 , the time-period T per increases again. Also the standard deviation starts to increase (Fig. 4). This indicates that the period of the subsequent pearling-coalescence cycles is not identical anymore, but begins to vary. In fact, from α ≈ 2.15 × 10 −3 period doublings occur: First, the period begins to alternate between two different values but for larger α more periods occur.
This period doubling can best be illustrated by studying the volumes of all drops after pearling, i.e., in the split-up state. To determine the volume of the drops, we identify connected areas in the simulation domain, that are elevated above a certain threshold height h threshold = 1.05 and integrate over this subdomain: The threshold height is subtracted in the integration to ensure that drops of identical real volume, i.e., volume above the adsorption layer, but different base area give the same value    For inclination angles above α ≈ 2.9 × 10 −3 , the drops become long enough to span the whole domain; i.e., they form rivulets and then no pearling or coalescence occurs. But still, these rivulets can exhibit a time-periodic behavior, e.g., its local height varies along its length and these surface-wave-like modulations travel downward. Further details of the behavior of the rivulets are not in our present scope. Note that the inclination at which the domainspanning rivulets form depends, of course, on domain size. Another option that we do not pursue here further, is to consider not one drop on a domain of length L x , but two identical drops on a domain of size 2L x . They may show a different pearling behavior, because the volume conservation now holds for the two drops and their satellite drops together and not individually for each drop and its own satellite drop. Therefore also other time-periodic behaviors are possible, like two drops that share one satellite drop.   Fig. 1 focusing on the monotonously linearly unstable sub-branch II and the oscillatory linearly unstable sub-branch III, parts of the neighboring sub-branches I and IV and the beginning of side-branch A. The arrows indicate the dynamics that occurs when unstable drops on sub-branches II and III are disturbed. Normally, they all relax to the stable sub-branch I (cf. Fig. 9).

C. Unstable stationary solutions: Sub-branches II and III
Next we give some more details on sub-branch II, i.e., the part of the bifurcation curve in Fig. 1 that connects the first two saddle-node bifurcations. We start at the first saddle-node bifurcation at α SN1 and follow sub-branch II with decreasing inclination angle; cf. Fig. 7 for a zoom of the bifurcation diagram in this region. One clearly sees that side branch A (discussed in Sec. III B) emerges from sub-branch II at an α that is slightly smaller than α SN1 . This point corresponds to the homoclinic bifurcation discussed in Sec. III B. All solutions on sub-branch II are linearly unstable drop solutions that along the branch undergo a morphological change. Namely they develop a protrusion at their back which increases in length with decreasing α. This is shown in the snapshots III and IV in Fig. 1. When performing time simulations using these drops as initial conditions, their instability always leads to a dynamical change of the drop shape and velocity. Normally, the protrusion at the end of the drop shrinks until the drop has the shape of the oval drop on the stable subbranch I at identical α. Alternatively, the drop may undergo a single pearling-coalescence cycle before approaching sub-branch I. At the inclination α SN2 ≈ 1.266 × 10 −3 , a second saddle-node bifurcation occurs in which one could expect the unstable eigenvalue to become stable again. This is, however, not the case illustrated in Fig. 8  the imaginary part is shown as the gray dashed line). By increasing α further, the unstable pair crosses the imaginary axis, i.e., it stabilizes in a Hopf bifurcation at α hopf = 1.392×10 −3 .
There the stable sub-branch IV starts. i.e., its front and back interact [cf. Fig. 11(c)]. This results in a number of saddle-node bifurcations as shown in the bifurcation diagram in the left panel of Fig. 11. However, here only the first of these at α ≈ 3.5 × 10 −3 is of interest as it marks the end of sub-branch IV.
The part beyond (sub-branch V, dashed line in Fig. 11) is unstable and non generic as its particular location depends of course on the chosen domain size. The interaction and subsequent fusion of drop front and back results in a rivulet solution that is first modulated but becomes translationary invariant in the x direction [cf. Fig. 11(d)] in the final bifurcation (filled circle in Fig. 11) that corresponds to a Hopf bifurcation. Note that for a perfect rivulet no velocity can be defined, therefore the dotted line in Fig. 11 represents the velocity of small traveling modulations on the rivulet that are picked up by the numerical continuation.
For larger domain sizes, one needs to go to a larger inclination angle to sufficiently lengthen a drop such that it spans the whole domain and becomes a rivulet. Therefore, solutions as in Fig. 11(c) would occur for larger inclination angles. This would extend the bifurcation diagram towards larger inclination angles.

IV. DEPENDENCE ON DROP VOLUME: SCALING LAWS
This section focuses on the power laws that have been mentioned in the previous sections.
Quantifying these power laws allows us to generalize our observations (that are up to here mainly obtained for a single drop volume) to a broad range of drop volumes. Figure 12 shows bifurcation diagrams of stationary drops for several different drop volumes V . It shows that its overall appearance does not change over the considered range of drop volumes.
The first part, sub-branch I (as discussed in Sec. III A), shows an increasing initial slope as the volume increases. This corresponds to the fact that larger spherical or oval drops slide with higher velocities because the gravitational pulling increases more strongly than the viscous friction (also see Ref. [18]). Fitting the linear parts for all analyzed drop volumes V reveals how the slope U/α depends on V . These slopes are shown in Fig. 13(a) for a volume range of about three orders of magnitude. Using a power law ansatz for this dependency we i.e., the exponent β 0 is larger than 1/2 and smaller than 2/3. This can be understood from a simple consideration: The velocity U of the sliding drop results from a force balance between gravitational pulling and the viscous friction force. The gravitational force F G is proportional to the inclination angle and the volume of the drop, F G ∼ αV . For the friction force F F one may discuss two limiting cases: It might either be proportional to the size of the footprint of the drop or to the length of the contact line. In the former case, one expects Assuming the force balance F G = F F holds, one will obtain either U α lin ∼ V 1−2/3 = V 1/3 or U α lin ∼ V 1−1/3 = V 2/3 . The numerically obtained exponent β 0 ≈ 0.569 is closer to 2/3 than to 1/3. This indicates that the viscous friction is mostly localized close to the contact line and is less related to the footprint of the drop. This is further discussed in Sec. V below, where we analyze the spatial distribution of the viscous dissipation within the drops.
The next power law gives the dependence of the angle of the first saddle node bifurcation α SN1 on drop volume V [cf. Fig. 13(b)]: and a 1 = (2.9 ± 0.1) × 10 3 .
This locus nearly coincides with the angle at which the pearling sets in (see discussion of details in Sec. III B). Studying the dependence of the corresponding velocity U SN1 at this bifurcation on V reveals a slightly more complicated power law shown in Fig. 13(b) that includes an offset: On the one hand, the relations Eqs. (14) and (15)

V. ENERGY DISSIPATION IN A SLIDING DROP
Finally, we investigate how the dissipation mechanisms differ for the various described stationary sliding drops. We determine the local and total dissipation for stationary sliding drops based on the dissipation per volume that is in long-wave approximation given by   We now turn to a discussion of the dissipation inside the drop, where the local persubstrate-area dissipation in long-wave approximation is obtained by integration across the film height. In dimensionless form it is which is equivalent to Note that the factor 1/3 in the dissipation corresponds to a factor of 1/3, which we absorbed into the scaling of the mobility. The total dissipation is then given bȳ We use the x-translation invariant rivulet h r (y) as reference state. It exists for any inclination and is approached by the studied drop solutions at large inclination angles (cf. Sec. III D, Fig. 11). Writing the total dissipation as a functionalD[h] we have for the relative total withD where we have used that for such rivulets there are no lateral mean flows, i.e., ∂ yyy h + ∂ y Π(h) = 0. Note that the reference state is mainly subtracted to ensure readability of the graphical representation. The relative dissipation for droplets on the main branch of Fig. 1 is given as solid and dashed lines in Fig. 15. On sub-branch I the relative total dissipation monotonously increases with increasing α. Beyond the first saddle-node bifurcation at α SN1 , the dissipationD rel decreases along the unstable sub-branch II (with decreasing α). Beyond the second saddle-node bifurcation at α SN2 the dissipation first increases again with increasing α along sub-branches III and IV. It reaches a local maximum at about α ≈ 2.2 × 10 −3 before decreasing towards the reference rivulet state.
The spatial distribution of the dissipation within the oval sliding droplets on sub-branch I is significantly stronger close to the advancing and receding parts of the moving contact line as compared to the central part of the footprint or, indeed, the lateral parts of the contact line. This contrast gets even stronger with increasing droplet volume; cf. i.e., ridges with imposed cylinder symmetry are studied in Ref. [61]. They also find that dissipation at the contact lines dominates.
Interestingly, the dissipation pattern is quite different in the second type of stable sliding drops, i.e., the droplets with elongated tails on sub-branch IV. They show some dissipation at the leading part of the moving contact line, however, most of the dissipation is localized at the center of the tail along its entire length. The absolute value of the total relative dissipation is at identical α always larger for the droplets on sub-branch I than those on subbranch IV. The time-averaged dissipation for the time-periodic pearling-coalescence cycles on side branch A lies always between the values for sub-branches I and IV.
In this work we have presented a theoretical investigation of the bifurcation structure of three-dimensional sliding droplets on a smooth homogeneous inclined substrate. Focusing on small inclination angles and small contact angles and free surface slopes, we have employed an asymptotic long-wave evolution equation for the film thickness profile h(x, y, t) that accounts for capillarity and wettability through Laplace and Derjaguin pressure terms. The latter consists of antagonistic power laws, ensuring that at equilibrium drops of finite equilibrium contact angle coexist with a thin adsorption layer (sometimes called a precursor film). The driving by gravity has been incorporated through an additional lateral force term. We have analyzed the model with both a numerical path-continuation technique and direct numerical simulations. The former has allowed us to determine stable and unstable stationary sliding droplets of various kinds as summarized in a comprehensive bifurcation diagram, while the latter has allowed us to determine time-periodic pearling-coalescence cycles of droplets and, in general, to determine the dynamic evolution of the unstable solution types.
The bifurcation diagram has been presented in the form of a dependence of the velocity of sliding drops on inclination angle. It shows which droplet types exist at which parameter values, how the various sub-branches are connected, and where stabilities change at the various bifurcations. Overlaps of linearly stable sub-branches indicate multistability, e.g., between simple stable oval sliding drops and stable sliding drops with a prolonged tail.
The different sub-branches are connected by saddle-node and Hopf bifurcations that are responsible for the change in stability along the branch.
Interestingly, the side branch of time-periodic pearling-coalescence cycles does not emerge at one of these local bifurcations but emerges from an intermediate point on an unstable sub-branch in a global homoclinic bifurcation, i.e., without change of linear stability. The behavior on this sub-branch has been found to correspond to interesting dynamics: first, a repeated cycle of splits into separate drops and their subsequent coalescence. For increasing α, these cycles undergo a cascade of period doublings resulting in irregular behavior at larger angles. This closely resembles the behavior found for dripping faucets in dependence of the flow rate [62]. One may even argue that our transition between individual oval drops to the stable drops with prolonged tails bears a passing resemblance to the dripping-jetting transition observed for a dripping faucet [63].
We have generalized our results by studying the influence of the droplet volume on the encountered transitions in behavior. We have found that the general appearance of the bifurcation diagram is universal over several orders of magnitude in volume (also cf. note [64]). Main features, such as, the volume-dependence of the slope of the linear relation of sliding velocity and inclination angle for simple oval droplets and the dependence of the locus of the first two saddle-node bifurcations can be brought into the forms of power laws.
In fact, these power laws are similar to experiment results in Refs. [3,[7][8][9]. Here, we found the same functional dependency on the drop volume with an exponent of about 0.57 whereas the results in the references conclude an exponent of 2/3.
The time evolution of the droplets in the pearling-coalescence regime can be well compared to earlier numerical approaches to similar models. Reference [16] presented a number of time simulations also using a precursor film model based on a disjoining pressure combining two power laws. They described the transition to pearling when increasing the driving force and mentioned that a large droplet may decay into an almost chaotic pattern. Reference [65] explored the influence of the precursor height, mesh resolution, and numerical accuracy on the morphologies. Shapes of sliding droplets are also investigated in Ref. [33] employing a slip model and studying droplets of finite support; this allows these authors to observe a spectrum of droplet shapes from ovals to drops with monotonic or non monotonic tails; however, the pearling itself can not be observed. Our results significantly add to these partial pictures by providing the entire bifurcation structure, including stable and unstable droplets and their transitions. Numerical and analytical approaches focusing on the cusp formation [9,14] also obtain a universal dependence of the sliding velocity on the inclination angle when staying below the pearling instability. However, asymptotic approaches to the cusp properties do not allow for the identification and further analysis of the pearling instability as a bifurcation, which is what we have presented here.
Another important aspect is our analysis of local and global dissipation in the sliding drops. Calculating the dissipation by inner velocity gradients of the drop we have found that for the simple oval drops the dissipation is localized in the vicinity of the advancing and receding parts of the moving contact line of the drop. The first power law with a rather simple friction ansatz for a sliding spherical cap suggests the the friction is localized at the edge of the cap, in total agreement with the spatial distribution of the dissipation as it results from internal friction. Interestingly, the dissipation pattern is quite different in the droplets with elongated tails, where it is localized at the center of the tail along its entire length.
A possible future extension of our present analysis could be a consideration of inhomogeneities, e.g., substrates with regular or random modulations of their topographical or chemical properties. Such systems have already attracted much experimental [39,66] and theoretical [25,35,37] interest that highlighted interesting phenomena, especially concerning the pinning of sliding drops on heterogeneities. As the multiplicity of steady and traveling states is much larger on heterogeneous substrates, further effort is needed to understand the various competing influences. Advancing towards real-world applications, such as the movement of water with pesticides on plant leaves [67,68], such considerations are obviously necessary.
Another important future step in the theoretical analysis of sliding drops is the analysis of large scale dynamics, i.e., the time evolution of large ensembles of droplets where splitting and coalescence events continuously occur. Interestingly, the here-discussed power law for the onset of pearling not only provides the critical inclination angle at fixed drop volume but also a critical volume at a fixed angle. In this way, it should be possible to employ the information from our single-drop bifurcation study to predict the behavior of ensembles of interacting drops of various volumes. Such a connection of small-scale and large-scale analysis will be pursued elsewhere.
Finally, we emphasize that the approach pursued here of combining path-continuation and direct time-simulation to establish the bifurcation structure for sliding drops on a twodimensional substrate that is described by a thin-film equation is equally well suited to closely related equations. Examples are (i) the convective Cahn-Hilliard equation describing phase decomposition under the additional influence of an external lateral driving force or faceting dynamics of growing crystals [69] and (ii) the Kuramoto-Sivashinsky equation whose variants describe surface waves and deposition patterns [70,71]. We are not aware of any path-continuation study that would elucidate the bifurcation structure of their fully twodimensional solutions.