Steering droplets on substrates with plane-wave wettability patterns and deformations

Motivated by strategies for targeted microfluidic transport of droplets, we investigate how sessile droplets can be steered toward a preferred direction using travelling waves in substrate wettability or deformations of the substrate. To perform our numerical study, we implement the boundary-element method to solve the governing Stokes equations for the fluid flow field inside the moving droplet. In both cases we find two distinct modes of droplet motion. For small wave speed the droplet surfes with a constant velocity on the wave, while beyond a critical wave speed a periodic wobbling motion occurs, the period of which diverges at the transition. These observation can be rationalized by the nonuniform oscillator model and the transition described by a SNIPER bifurcation. For the travelling waves in wettability the mean droplet velocity in the wobbling state decays with the inverse wave speed. In contrast, for travelling-wave deformations of the substrate it is proportional to the wave speed at large speed values since the droplet always has to move up and down. To rationalize this behavior, the nonuniform oscillator model has to be extended. Since the critical wave speed of the bifurcation depends on the droplet radius, this dependence can be used to sort droplets by size.


Introduction
Liquid droplets on a solid substrate occur naturally during, e.g., rain and condensation [1,2] or artificially in printing of ink [3] or medical testing [4].Depending on their interaction with the substrate and gas phase, they either sit on a substrate, spread, contract, or even move laterally along the substrate [5].Lateral motion can be induced by a variety of mechanisms; through nonuniform electric fields (electrowetting) [6], gravity [7], Marangoni advection [8,9], gradients in wettability [4,10,11], or deformations of the substrate itself [12,13].Here, we focus on the latter two mechanisms.
In our previous work, we used the gradient in a wettability step to move a droplet forward using a feedback loop, which synchronizes step and droplet motion [14].However, due to surface imperfections at which droplets get pinned, such a motion can be difficult to realize in an experiment [5].To overcome the pinning, a large wettability gradient is needed [15], e.g., induced by structural changes of the substrate [16].
Deformations of the substrate itself also move a droplet [17].In Ref. [18] it has been demonstrated that capillary forces attract droplets towards regions with large inward curvature.An example of this are pipettes.It is possible to externally control such substrate deformations, e.g., using gels that reversibly swell and unswell in response to light with specific wavelengths [13,19,20].Already more than three decades ago wetting of curved rigid substrates was investigated [21].In more recent work, elastic substrates were studied that deform in contact with a liquid droplet [22][23][24].Internal tension builds up as a substrate deforms, which changes its surface tension what is known as Shuttleworth ef-fect [25,26].Here, we just consider curved rigid substrates, where this effect does not occur.
In this article we compare two mechanisms for inducing droplet motion by imposing the spatio-temporal pattern of a travelling wave onto the substrate; first in wettability and thereafter in the height profile of the substrate.This enables to steer the droplet into a prescribed direction.Because the travelling wave is spatially periodic, the droplet is continuously driven out of equilibrium and cannot settle into a sessile state.Instead, the droplet can be considered as a driven nonuniform oscillator [27], as we show, the driving force of which depends on its position relative to the travelling wave.In Fig. 1 we display two series of snapshots that illustrate the two possible modes of motion, which we term wobbling (left in Fig. 1) and surfing (right in Fig. 1), respectively.In the following we reformulate our formalism to apply the boundary element method (BEM) [28] to a wetting droplet, which we developed in Refs.[14,29], and adjust its boundary condition to account for either an imposed wettability pattern or a height profile of the substrate.Because in our previous work we only studied a plane substrate [14,29], the method needed to be extended to curved substrates.
First, in Section 2 we describe the theory of dynamic wetting for small droplets and our numerical approach based on the BEM.Next, in Section 3 we introduce the socalled nonuniform oscillator, which we use as an analytic model to interpret the observed droplet dynamics.We present and analyze the effect of travelling waves in wettability in Section 4 and travelling-wave deformations in Section 5. Finally, we comment on implications of our findings for the sorting of droplets by size in Section 6 and conclude with Section 7.

Boundary element method for dynamic wetting
We consider the motion of viscous droplets at small scales in the creeping flow regime, where viscous forces dominate, while inertia is negligible.Accordingly, fluid flow within droplets are described by the Stokes equations and the incompressibility condition, Here, we use fluid velocity field v, pressure p, and shear viscosity η.In addition, appropriate boundary conditions have to be chosen, most importantly, describing forces acting at the fluid boundary.Notably, the flow fields determined by the Stokes equations adapt instantaneously to the boundary conditions because the equations do not contain any time derivatives.As a result, when describing the dynamics of the droplet boundary by a surface velocity field R ̇, it is linked to a friction force field by the linear relation where G is a friction operator [30,31].Since the Stokes equations are linear in the fluid velocity field v, G is a linear operator that does not depend on v.After discretizing the droplet surface, G becomes a friction matrix and R and K are vectors.
Our goal is to determine K, R, and G explicitely for the mesh discretization of the droplet surface described below.This will give us a systematic approach for calculating the motion of the droplet.We start by introducing the discretized droplet surface, which specifies the vector R, then formulate the boundary condition at the interfaces to determine the generalized friction force K, and finally derive the friction matrix G.In addition, we nondimensionalize our quantities in Section 2.6.

Mesh discretization
As noted above, our aim is to reduce the droplet motion to the dynamics of its boundary consisting of the liquidsubstrate and liquid-gas interfaces as well as the threephase contact line, which requires special consideration.To treat the droplet dynamics numerically, we describe the entire droplet surface by a triangular mesh.Each triangle consists of three vertices and three edges, which it shares with its direct neighbour triangles.As long as the triangular mesh stays intact, the shape of the droplet is completely specified by the positions of the vertices r i , which define the configuration vector R = (r 1 , r 2 , r 3 , . ..).
In general, the part of the surface force field acting on the droplet at the liquid-solid and liquid-gas interactions derives from the functional derivative of a surface energy w.r.t. the configuration field R.However, since the droplet surface is discretized, the functional derivative becomes gradients w.r.t. the vertex positions.To continue, we construct an expression for the surface energy using surface integrals and discretize the integrals for our mesh in the following part.

Droplet free energy and discretization
The free energy of an incompressible liquid droplet at constant temperature can be written as where γ is the surface tension and the area A encompasses all the surfaces of the droplet.The bulk free energy is kept constant here.The surface integral can be split up into two parts: The liquid-gas (lg) interface and the solid-liquid interface (sl) with the respective surface tensions γ lg and γ sl .The solid-gas interface of the substrate not in contact with the droplet contributes to the free energy with surface tension γ sg .Whenever part of this interface is covered by the droplet, the surface energy changes by γ sl − γ sg .Thus we obtain where A lg and A lg are the areas of the respective interfaces.Young's law for the equilibrium contact angle θ eq is a force balance in the substrate plane that states that γ sl − γ sg = −γ lg cos θ eq . 1 Using it together with γ lg = const., we obtain where θ eq is determined by the local value of γ sg −γ sl .Because we want to study movable droplets on substrates with nonuniform wettability, the remaining integral cannot be further simplified here.We now turn our attention toward the constraints on the shape of the droplet.They are two-fold: First, the volume of the droplet is conserved and, second, the normal forces on the liquid-solid interface sum up to zero since any fluid forces are compensated by the forces exerted by the rigid substrate that resists any elastic deformation.Using the method of Langrangian multipliers, each constraint gives rise to an additional term in the energy functional: Here, the second term on the right-hand side constrains the droplet volume to the value V 0 and the pressure p 0 serves as the Langrange multiplier.The third term firmly connects the vertical position z of the liquid-solid interface to the height variations h(x) of the substrate written in Monge representation with x = (x, y).In the case of a flat substrate, h = 0, this term reduces to which guarantees that the solid-liquid interface A sl is defined by z = 0. Finally, note that the volume of the droplet can be expressed by the surface integral as one verifies immediately with Gauss's theorem.
To discretize all the surface integrals of the free energy F, we use the triangular mesh of the droplet surface and sum up the contributions from each triangle so that Here, we have to distinguish between functions f , the values of which are known for any position r and those functions, the values of which are only known at the vertices.In the first case, an example is cos θ eq , we use Gaussian quadrature with 4 sample values of the prescribed profile to calculate the integral on each triangle.The sample values are taken at positions (1/3, 1/3, 1/3), (3/5, 1/5, 1/5), (1/5, 3/5, 1/5), and (1/5, 1/5, 3/5) in barycentric coordinates of the triangle and weighted by -27/48, 25/48, 25/48, and 25/48, respectively [32].
In the second case, an example is the calculation of the volume V , we linearly interpolate between the values at the vertices.If f has a constant value f i on the area of the triangle, for example, in the case of γ lg , we simply use A i f i .

Method of Lagrange multipliers
After having discretized all contributions to F, we obtain terms as a function of the configuration vector R, including all the constraint terms, which we write as λ j g j (R) = 0.For convenience, we label λ 0 = p 0 and assign the remaining multipliers to the rigidity constraint.Thus, Eq. ( 6) becomes We discretize the surface integral of the rigidity constraint in Eq. ( 6), where we subsume the area A i surrounding each vertex into the λ i , respectively, and we sum over all N vertices on the solid-liquid interface.The force due to this constraint is guaranteed to be normal to the substrate, since ∇ i g i = n i , as shown in Appendix B, and therefore equivalent to the desired rigidity constraint. 2 This construction greatly simplifies the procedure below because ∇ R g i is easily accessible.The purpose of the Lagrange multipliers λ i is to constrain the force −∇ R F such that it does not have a component perpendicular to the manifold defined by the constraints g i (R) = 0, otherwise the constraints would be violated. 3Therefore, the condition must hold.Using Eq. (10) for F in this condition with we obtain a set of linear equations for the Lagrange multipliers λ i (including λ 0 ), which we solve numerically.With this, F is known up to a constant. 4

Boundary condition and dynamic equation
The boundary condition of the Stokes equations, Eq. ( 1), is the balance of forces acting at each point on the interfaces of the droplet.There are three contributions to this force balance.First, the traction of the liquid on the interface, σn, where is the stress tensor, n the unit normal vector, and δ ij is the Kronecker symbol.For the following argumentation, we split the stress tensor into a part −p 0 n, where p 0 is the spatially uniform static pressure due to the volume constraint, and a spatially varying part σ ˜n = σn + p 0 n.Second, the traction of the substrate at the substrateliquid interface.It consists of the constraint force f rigid , with which the substrate resists deformations induced by the droplet, and f deform , which acts on the liquid when the substrate deforms according to the prescribed height profile h.And third, the traction due to the surface tension of both interfaces which includes the Young force acting on the contact line and which we denote here by −δF 0 /δr, i.e., the functional derivative of the surface free energy w.r.t. the surface parameterization.Note that we neglect the vapour pressure of the gas phase surrounding the droplet, here and in the following.In total, we have the force balance at each point of the droplet surface, where f rigid + f deform is zero at the gas-liquid interface.
Note that f deform is due to the boundary conditions at the liquid-substrate interface, which include continuity of the normal velocity component and a slip boundary condition of the tangential component, as we will specify in Section 2.5.2.The first two terms on the r.h.s. of Eq. ( 15) together with p 0 n can be rewritten as −δF/δr using the Lagrangian constraints introduced in Eq. ( 6).So we obtain To illustrate this force balance, we consider the equilibrium case of the plane surface, f deform = 0, and no spatially varying stress, σn = 0, such that δF/δr = 0.
On the liquid-gas interface, δF 0 /δr = −2κγ lg n [33] so that p 0 becomes the Laplace pressure where κ is the mean curvature of the droplet, which is here negative.We now rephrase boundary condition ( 16) for the discretized droplet surface.At each vertex i the traction σ ˜in i from the fluid is multiplied by area A i and collected in the friction vector Similarly, we collect the substrate traction f deform in a force vector deform A 3 , . ..) .
Our method to calculate F is described in detail in Section 2.5.As usual, the gradient of F describes a force −∇ R F, which takes the place of the functional derivative −δF/δr.Thus, in analogy to condition (16), we have Finally, with the ansatz K = −GR ̇from Eq. ( 2), which we will motivate further in the following section, the relaxational dynamics of the surface is given by Here, we observe that for a stationary substrate with F = 0, a local minimum of F corresponds to a stationary droplet since then Ṙ = 0.

Generalized friction
Now, we specify the friction matrix G such that all sources of friction related to the droplet are accounted for.Three types of friction are relevant in our setup.First, there is shear friction or viscous shear stresses in the droplet fluid, when neighboring fluid layers move relative to each other.Second, there is slip friction between the substrate and the directly adjacent fluid layer.And finally, there is friction of the moving contact line.The latter should naturally arise from the first and second mechanisms.However, so far this has proved difficult to show analytically [5] and also to implement numerically due to the flow singularity that arises at the contact line [34].Therefore, we use a hydrodynamically motivated friction law for the contact line.It is closely related to the Cox-Voinov law [35], which has successfully been matched to experiments, e.g., in Ref. [36].

Shear friction
We start with shear friction in the droplet fluid.The solution of the Stokes equations within a given domain is equivalent to a boundary-integral equation that relates flow velocity v(r) to surface forces σ(r)n and force/source dipoles at the domain surface [28], where α is the inward solid angle, O the Oseen tensor, and T the associated stress tensor.For the surface velocities relevant here, α = 2π at the fluid interface and α = 2θ at the contact line, while inside the droplet α = 4π.We note that σn in Eq. ( 22) can be replaced by σn since a constant pressure does not contribute to the integral as shown in Appendix A. To discretize the boundary-integral equation, we proceed as in our previous work [14].We assign to each surface vertex parts of the surrounding triangles such that the resulting polygonal cells with area A j do not overlap and cover the whole region of integration.The result is a set of linear equations for the vertex velocities v j , Recalling the definition of K in Eq. ( 18) and introducing block matrices, which we specify in Appendix C, or as where X −1 (C + Y ) is the shear contribution to the friction matrix.

Liquid-substrate friction
We now turn to friction between liquid and substrate due to a non-zero slip velocity.The conventional slip boundary condition is l s P t σn = ηP t v, where l s is the slip length and P t = 1 − n ⊗ n the projection operator on the tangential plane.Note that the surface normal n points out of the fluid according to Ref. [37] and that P t σn = P t σ ˜n.However, since we are interested in moving substrates, the fluid motion relativ to the local substrate velocity v subs and the deforming traction f deform become relevant and we have l s P t (σ ˜n − f deform ) = ηP t (v − v subs ).On the discretized solid/fluid interface of the substrate using Eqs.( 18) and ( 19), the boundary condition becomes Here, P || : R 3n → R 2m is the projection operator from all the n droplet vertices onto the tangential space of the substrate with m vertices.For consistency, we defined subs , v subs , . ..) as the vector of the substrate velocities at the vertices, which also includes the vertices on the liquid-gas interface with v is the slip contribution to the friction matrix, where δ ij is the Kronecker symbol and i, j refer to vertices on the substrate.For the parallel components of vector F in Eq. ( 21), we can now write P || F = QP || V subs so that K remains purely linear in Ṙ as we have stated in Eq. ( 2).So far, we have not addressed the continuity condition for the normal velocity component, v • n = v subs • n, at the liquid-substrate interface, which we already mentioned in Section 2.4.Similar to P || , it is convenient to introduce the projection operator onto the normal space of the substrate vertices P ⊥ : R 3n → R m .The boundary condition can then be restated for our discretized mesh as P ⊥ Ṙ = P ⊥ V subs .It directly determines the normal components of the vertex velocities at the substrate, which are contained in Ṙ.However, the condition also means that the substrate pushes or pulls on the liquid as it deforms its height profile h and the liquid resists with friction.We have already introduced the traction forces due to the substrate deformations.Their normal components are contained in the vector F ⊥ = P T ⊥ P ⊥ F , which has to be calculated numerically at the same time as those components of Ṙ which are still unknown.

Contact-line friction
Finally, we address the friction of the contact line.In 1964 Moffatt solved the Stokes flow at the contact line of a moving wedge filled with a viscous liquid and found an analytic expression for the traction σn along the liquid-gas interface [38], Here, H is the distance of the interface from the substrate, θ the contact angle, and v cl the tip velocity of the wedge along the substrate.Again, if the substrate moves, we have v cl = v−v subs and its deforming traction f deform , which modify Eq. ( 30).The velocity lies in the tangential plane of the substrate and is normal to the contact line.To proceed, we integrate the traction over height H along the liquid/gas interface from a microscopic length l s , i.e., the slip length, up to a mesoscopic length ζ In the discretized representation of the droplet surface, we choose for the integral the surface force σ i n i assigned to the vertex i on the contact line, which is contained in the friction vector K, i.e., For the contact line, we therefore have the additional contribution where P cl : R 3n → R k is the projection operator on the direction normal to the contact line and tangential to the substrate, where k is the number of vertices on the contact line.Obviously, these vertices contribute with to the friction matrix.Here we use the Kronecker symbol δ ij , which means that M is diagonal.As in Section 2.5.2, we identify P cl F = M P cl V subs as contribution to F in Eq. (21).Therefore, in total we find where T indicates matrix transposition.Finally, we are able to formulate the full set of equations.Collecting the contributions from Eqs. ( 27), (28), and (33) linear in Ṙ, the total friction matrix G is thus given by By inverting G, we arrive at the dynamical equation for R, with K given by the force balance, Eq. ( 20), which we further specify here as using Eq.(35).We solve Eq. ( 37) numerically to find the unknown components of Ṙ and F ⊥ .They include all components of Ṙ corresponding to vertices on the liquid-gas interface and to the tangential motion of the substrate vertices, as well as all non-zero components of F ⊥ .They are connected to the normal motion of the substrate vertices.

Nondimensionalization and parameters
We choose characteristic units, which derive from the motion of the droplet surfaces.First, as a unit of length l, we choose the initial radius R 0 of the liquid-solid interface of the droplet.The dominant time scale for the deformation of the droplet is determined by the motion of the contact line.A good estimate for the speed of the contact line was derived by Voinov in Ref. [35] using Eqs.(30) and (15).He ultimately arrives the Cox-Voinov law From the prefactor of the Cox-Voinov law and 0 , we obtain an inherent time scale for the motion of the contact line.Furthermore, a characteristic force follows from combining the liquid-gas surface tension γ lg and R 0 as f c = γ lg R 0 .Now, writing Eq. ( 1) in nondimensional form using these characteristic units, gives rise to a dimensionless viscosity η ˜= [9 ln(ζ/l s )] −1 .Thus, our results become independent from specific values of γ lg and R 0 .They only depend on the ratios η ˜, ζ/R 0 and l s /R 0 , for which we choose appropriate values.As in Ref. [14] the parameters of our system are matched to droplets from a 90%-glycerol/10%-water mixture for which de Ruijter et al. [36] measured in experiments ln(ζ/l s ) = 44, η = 209 mPa s, γ lg = 65.3 mN m −1 , l s = 1 nm, and mass density ρ = 1.24 g ml −1 .To determine a value for ζ needed in the friction matrix of Eq. ( 34), we fitted the experimental data for the droplet relaxation in Ref. [36] in our simulations and obtained ζ = 0.17 R 0 .All simulations are performed using a surface mesh of 1199 vertices which, in their initial configuration, have an average distance of 0.094 R 0 from the nearest neighbor.In the following, we use R 0 = 100 µm.
From the characteristic parameters, we determine the Reynolds number as which matches the assumption of negligible inertia inherent in Eq. ( 1).

Driven droplets as nonuniform oscillators
Our droplets are driven by travelling waves either in wettability or in deformations of the substrate.For small wave speeds they are strongly connected to and surf on the wave, while for large speeds they can no longer follow and perform a wobbling motion.Before we describe the phenomenology of both substrate types, we introduce and study a simplified model for the motion of a droplet under a spatially periodic external stimulus.Neglecting inertia, we start with a simplified force balance for the center of mass y of the droplet, −Γy ̇+F (y, t) = 0, where Γ is a friction constant and F is the effective force exerted by the substrate, on which the droplet is moving.
The travelling wave pattern is modeled by the effective force with wavelength λ and wave speed v wave .It models that there are specific positions in the travelling wave, y − v wave t = nπ with n = 0, 1, . .., where no force is excerted on the droplet.For our specific wave patterns, these are the minima and maxima.In the co-moving reference frame of the substrate pattern, y com = y − v wave t, the droplet's velocity y ̇com = y ̇− v wave generally oscillates around −v wave , i.e., with amplitude v c = F 0 /Γ.Depending on its position on the pattern, the droplet's speed is increased or reduced.In particular, for y com = nλ/2 the force on the droplet is zero, which means y ̇= 0 or y ̇com = −v wave .Although y ̇com oscillates around −v wave , that is not its time average.Instead −v wave is the spatial average of y ̇com .We can map Eq. ( 43) onto the so-called ideal nonuniform oscillator for phase variable φ with constants a and c [27] This results into the expression φ ̇= a − c sin φ (45) also called the Adler equation [39].Unlike the uniform (harmonic) oscillator, the phase of which has a constant rate of change, φ determines its own rate of change.In particular, the oscillation stops for |c| ≥ |a| when the two terms on the r.h.s.balance each other [40].
Strogatz derives the period T of the nonuniform oscillator in Ref. [27] and following his derivation we find with our parameters To predict the long-time behavior of the droplet, for example, its position, the time averaged droplet speed v īs useful, which we derive now.Generally, in steady state the droplet position oscillates relative to the substrate pattern with a period where the difference between v wave and v ¯gives the droplet's speed relative to the substrate pattern with wavelength λ.Differently speaking, T is the period of the forcing experienced by the droplet through the travelling wave.Setting Eq. ( 46) equal to Eq. ( 47), we obtain which relates v ¯to v wave with only a single free parameter, v c .Consequently, in this model v c contains all material properties of the driving mechanism and dropletsubstrate interaction.The limiting case for large v wave is v ¯≈ v 2 c /(2v wave ) and therefore which is useful to find v in the simulations but also in experiments.
However, v c is not just a fitting parameter because it also marks a change in the dynamics of the droplet.For v wave ≤ v c the oscillation stops, as predicted by the non-uniform oscillator, which means y ̇com = 0 or, equivalently, y ̇= v wave = v ¯.In our following analysis, we will refer to this stationary dynamics as 'surfing' and to the oscillatory dynamics as 'wobbling'.

Travelling waves in wettability
We first study how travelling waves in wettability can move droplets on a completely flat substrate, similar to our earlier study of moving steps in wettability in Ref. [14].However, there are clear differences, which we explain in the following.

Substrate dynamic
Here, we impose a travelling wave in wettability on a flat substrate such that the prescribed equilibrium contact angle varies according to with wavelength λ = 2 R 0 , v wave the speed of the wave in y-direction, the Heaviside or Θ-function, the amplitude ∆θ = 30 • , and θ 0 = 90 • . 5he droplet moves toward regions of high wettability, which means the valleys of the travelling wave in θ eq .Because these valleys are steadily moving forward in y-direction, the droplet is biased to also move in that direction.However, as noted in Section 3, the droplet can perform different kinds of motion under the influence of this driving mechanism, namely surfing and wobbling motion (see Fig. 1), which we will analyse in detail now.

Surfing and Wobbling
We first review the basic phenomenology.A surfing droplet moves at the same speed as the travelling wave.Consequently, its shape remains constant and a stationary flow field exists inside the droplet (see Movie M1 in the ESI).This stationary state is identical to the surfing dynamics we observed in Ref. [14].
A wobbling droplet moves at a slower speed than the travelling wave pattern.It oscillates in speed and shape as it repeatedly passes over the peaks and valleys of the wave (see Movie M2 in the ESI).The wobbling dynamics was not present in our previous study of moving steps in wettability [14] because it requires a periodic wettability pattern.
According to the nonuniform oscillator model, we outlined in Section 3, wobbling occurs above a critical speed of the travelling wave pattern.As the speed of the travelling wave pattern decreases, the period at which the droplet shape oscillates grows and it diverges exactly at the transition.So, the resulting surfing state can be considered as a wobbling which is frozen in time.This bifurcation is analyzed in detail in Section 4.4.

Quantitative analysis
One characteristic quantity of the surfing and wobbling dynamics is the droplet speed and we ask how it relates to the speed of the wave pattern.The droplet speed needs to be defined carefully, especially for wobbling, because it varies periodically over time.It is appropriate to average the droplet speed over one period of oscillation T .The period of oscillation refers to the shortest duration between two identical shapes of the droplet, albeit not identical positions in the lab reference frame.The average droplet speed is then formally given by where the starting time t is set after the decay of any transients and the number of periods n is chosen according to the available simulation data.Note that, for the case of a surfing droplet, we have y(t) = v wave t + y 0 with position y 0 at t = 0. Here, the time intervall T is arbitrary and the definition above reduces to v ¯= v wave .
In Fig. 2(a) we plot the speed ratio v ¯/v wave versus v wave as determined from our BEM simlations.The speed ratio is 1 in the surfing state and it decays as v −2 wave in the wobbling state (dashed line), as predicted from the nonlinear oscillator model and Eq. ( 48) in the limit v wave ≫ v c .This means v ¯is inversely proportional to v wave for wobbling and equal to v wave for surfing.The time-averaged speed of the droplet decays to zero which also matches our observation from an earlier article that droplets are less susceptible to rapid changes in wettability [29].For a full quantitative comparison with our model, we need to determine the velocity v c .Following Eq. (49), we plot √ 2v wave v ¯versus v wave in Fig. 2(b) and read off v c = 2.65 R 0 τ −1 at large wave speeds (dashdotted red line).The dashed line in Fig. 2(a) is the prediction from Eq. (48) using the determined value for v c .
The transition from surfing to wobbling appears continuous, meaning there is no jump in droplet speed, as predicted by the non-linear oscillator model.However, the transient dynamics decays slowly close to the tran- sition so that we lack data close to v wave = v c .If the transition is indeed continuous, the droplet moves fastest exactly at the critical wave speed v c .The period T of wobbling oscillations, displayed in Figure 2(c), grows inversely proportional to v wave − v c as v wave approaches v c from above.These observations for the wobbling state again quantitatively match the prediction of the nonlinear oscillator model.

Comparison with travelling steps in wettability
In Ref. [14] we observed a droplet surfing along a substrate on a wettability step that moves with velocity v step and a sigmoidal profile σ(y).Similar to Eq. ( 43) we can describe the droplet velocity in the co-moving frame of the step with where for the sigmoidal function we chose σ(y) = (1 + e −y ) −1 , and δ is the step width.Most importantly, for v step < v c two surfing states exist.They sit outside the center of the step and are stable and unstable fixed points, which eventually meet in a saddle-node bifurcation at v step = v c .Furthermore, there is a sessile steady state where the droplet sits outside of the step.It occurs for any value of v step in the limit |y com /δ| ≫ 0 so that the contribution from σ ′ vanishes.The saddle node bifurcation with the stable and unstable fixed points as well as the stable fixed point of the sessile droplet are schematically illustrated in Fig. 3(a) using an appropriate shape variable.
Comparing the moving wettability step treated in Ref. [14] to the travelling wettability wave considered here and modeled by Eq. ( 43), we note that a sessile state is impossible due to the periodicity of the wave pattern.Instead, the wobbling dynamic takes the place of the sessile state for v wave > v c .In Fig. 3(b) the resulting bifurcation diagram is sketched.The saddle-node bifurcation of the moving wettability step is still present but for travelling waves it collides with the limit cycle associated with wobbling.This collision comprises a global bifurcation which is called Saddle-Node-Infinite-PERiod (SNIPER) bifurcation [27].
A SNIPER bifurcation was previously observed for droplets exposed to an oscillating external force on a static heterogeneous substrate [39] and for a droplet subject to a constant external force on a rotating cylinder [41,42] with the force acting perpendicular to the cylinder's axis of rotation.In both cases, the bifurcation marks the transition from a pinned to a sliding droplet.

Travelling-wave deformations
We now turn to the second driving mechanism introduced in Section 2. The substrate deforms such that its height profile follows a travelling wave.The wettability is kept uniform and we solely rely on the travelling height profile to move the droplet.The influence of gravity is negligible because the capillary length l c = √︁ γ lg ρ −1 g −1 ≈ 23 R 0 is much larger than the droplet radius.This raises the question of the driving mechanism that moves the droplet forward.
For an equilibrium contact angle of θ eq = 90 • , which according to Young's law occurs at surface tensions γ sl = γ sg , only the liquid-gas interface is relevant and the droplet deforms to minimize the area of this interface.Hence, on a curved surface the droplet moves away from peaks and toward valleys because there the liquid-gas interface occupies less area compared to a droplet with the same volume sitting on a flat part or a peak of the substrate.For θ eq < 90 • the droplet is even more attracted to valleys in the substrate, as it is energetically advantageous to cover more of the solid with liquid, while for θ eq > 90 • the situation is less clear.In Ref. [18] Lv et al. demonstrate how a droplet moves on the inside and outside of a conical substrate without any additional external stimulus.On the outside of the cone the liquid-substrate curvature is negative and on the inside of the cone it is positive.Notably, for all cases studied with θ eq ranging from 30 • to 120 • , the droplets move in the direction of increasing curvature.The same behavior was recently confirmed for droplets on cylinders of varying size [43].Indeed, when we place a droplet onto a substrate peak with θ eq = 120 • in our BEM simulations and then perturb the droplet, the position is unstable.Lv et al. term this mechanism curvi-propulsion and in this section we investigate if it generates the same wobbling and surfing dynamics we found for the wettability waves in Section 4 and if there are any differences between the two cases.

Surfing and wobbling
Instead of a wettability profile, we now impose a travelling height profile of the substrate with amplitude h 0 , h(y, t) = h 0 q(t) sin 2 y − ξ − y wave (t) 2λ , where ξ determines the initial phase shift or position of the wave relative to the droplet at y = 0. Especially, in the surfing dynamic we use a value of ξ = 0.3 R 0 to place the droplet close to its surfing position on the wave profile, while ξ is irrelevant for the wobbling dynamics.We introduce the factor and The purpose of the factor q(t) is to allow a continuous evolution of the surface deformation from a flat to the sinusoidal shape to which the droplet can adjust gradually starting at a specific position relative to the wave.Then, the deformation wave starts to travel at t = t 0 , where we choose t 0 in the interval 0.1 τ ≤ t 0 ≤ 0.2 τ .To be specific, we set λ = 4 R 0 for the wavelength, h 0 = 0.1 R 0 for the amplitude, and an equilibrium contact angle of θ eq = 90 • . 6In Section 6, we will explore travelling deformation waves for a range of wavelengths.With time the droplet again settles into either surfing or wobbling motion depending on the wave speed v wave .Because it is attracted to the valleys in the substrate, it tries to follow them as the travelling deformation wave progresses.If the travelling wave moves sufficiently slowly, the droplet can surf on a position behind the valley (see Movie M3 in the ESI).If the travelling wave moves faster than a critical speed, the droplet is always overtaken by the peak behind it, relaxes into the closest valley, and falls back again giving rise to the wobbling motion (see Movie M4 in the ESI).

Quantitative analysis
In analogy to the travelling wettability pattern, we analyze the time-averaged speed of the droplet as defined in Eq. (51).In Fig. 4(a) we again plot v ¯/v wave versus v wave .In contrast to the wettability wave, the speed of the wobbling droplet now grows linearly with v wave at large v wave .Due to this clear difference, we realize that our model needs to be adjusted.
To start, we re-examine the nonlinear oscillator model of Eq. (43), where a spatially periodic force drives the droplet with a constant amplitude v c .For the deforming substrate we introduced a traction force F ⊥ in Eq. ( 38), which imposes the deformation of the substrate onto the droplet shape.However, this force depends on v wave : It is zero for v wave ≤ v c , where the droplet's shape is constant and beyond v c it grows with v wave as the droplet is deformed more rapidly.As before, the force is also spatially periodic since it originates from the wave pattern of the substrate's height profile.Therefore, to include the dependence of the amplitude of the driving force on v wave , we extend Eq. ( 43) and write (56) The dimensionless coefficient ε ≪ 1 quantifies the relative strength of the wave-speed depending part of the amplitude.In analogy to Eq. ( 46), we find for the time period, which together with Eq. ( 47) gives the mean droplet speed, Here, we have used ).The small parameter ε means that v ¯increases proportional to v wave for large v wave .Note that ε does not affect the scaling of the wobbling period T with the inverse wave speed well above v c .Now, our procedure to determine ε and v c , the central parameters of the model, is as follows.First, from Eq. (58) we determine the limiting value of v ¯/v wave = ε for large v wave , which we can directly read off from Fig. 4(a) as ε = 3.565 • 10 −3 .Then, we subtract εv wave from both sides of Eq. (58), and expand the r.h.s. for large v wave to obtain Solving for v c yields which for ε = 0 indeed simplifies to Eq. (49).Accordingly, in Fig. 4(b) we plot v * as defined in Eq. (61).Because ε appears as coefficient of v 2 wave , v * is sensitive to small errors in ε at large v wave .However, around v wave /R 0 τ −1 = 2 it takes on a constant value, from which we determine v c = 0.185 R 0 τ −1 .The dashed black line labelled model in Fig. 4(a) shows that with these parameters our revised model quantitatively matches the BEM data (blue dots).Similarly, it also reproduces the wobbling period displayed in Fig. 4(c).Because around v wave = v c , the original and the revised model behave the same, the latter also predicts a SNIPER bifurcation at v wave = v c for droplets driven by travelling substrate deformations.
Our revised model fits the BEM data well and thus provides an interpretation of the asymptotic behavior of v ¯for large v wave .To capture the asymptote, we introduced an amplitude of the driving force, which depends on v wave .It accounts for the fact that a travelling substrate deformation always needs to displace the droplet, which requires a larger driving force with inreasing v wave .This is the clear difference to travelling wettability patterns treated in Section 4, where the pattern can just pass beneath the droplet without the need to lift it up.

A note on particle sorting
In Sections 4 and 5 we described that droplets move in response to spatial gradients in wettability or curvature of the substrate, respectively.Both, gradient and curvature, change in magnitude when we adjust the wavelength λ of the travelling wave pattern.Therefore, we can expect that the speed of the droplet v ¯also depends on λ.
In Fig. 5 we display for both driving mechanisms the mean droplet velocity v ¯as a function of λ for v wave = 50 R 0 τ −1 , which places the droplet far into the wobbling dynamic.For travelling deformation waves (green crosses in Fig. 5) no clear trend is discernible.However, for travelling waves in wettability (blue dots in Fig. 5) we observe v ¯decays as (λ/R 0 ) −2 .Recall from Eq. (49) that the asymptotic behavior at large v wave is v ¯= v 2 c /2v wave .Since v wave is independent of λ, we can infer that v c ∼ (λ/R 0 ) −1 .Consequentially, the ratio λ/R 0 can be used to tune at which speed v c the transition from surfing to wobbling occurs.Now consider a poly-disperse collection of droplets, such that R 0 varies among them.This means, for each droplet size there is a specific v c and some droplets might be in the surfing regime and others in the wobbling regime at the same value of v wave .Since droplets generally move more slowly in the wobbling regime than in the surfing regime, this naturally leads to a sorting mechanism for poly-disperse collections of droplets.

Conclusions
We have studied liquid-droplet motion on substrates using two different mechanisms: travelling waves in wettability and travelling-wave deformations of the substrate.
To that end, we first implemented the boundary element method (BEM) for the fluid velocity field inside the droplet starting from the free energy of the droplet, the balance of forces at the droplet interface, and the relevant friction contributions.
To interpret our numerical findings from applying the BEM, we first introduced an analytic model for a particle driven by a travelling-wave force, which maps onto the so-called nonuniform oscillator.For increasing wave velocity, the analytic model exhibits a bifurcation between two states, which we call surfing and wobbling and which directly connect to the phenomenology observed in our BEM simulations.
For travelling waves in wettability the results from BEM simulations for wobbling period and mean droplet speed are consistent with the analytic model and indeed confirm the bifurcation in the full dynamical system.This further implies that the motion of the droplet's center of mass can be characterized solely by the critical wave speed for the bifurcation, which is determined by the material properties of the droplet and substrate.For a specific set of parameters we determined the critical wave speed by varying the wave speed in the BEM simulations.
We also investigated travelling-wave deformations of the substrate, which give rise to what is called curvipropulsion in Ref. [18].Again, we could identify a transition from surfing to wobbling with increasing wave speed but, interestingly, the simulation results reveal that the droplet speed does not approach zero for large wave speeds.By including a speed-dependent amplitude of the driving force, the modified analytic model could account for this new behavior, as we demonstrated for a specific set of parameters in the BEM simulations.
Finally, we also performed simulations where we varied the wavelength λ of the travelling wave for both driving mechanisms in the wobbling regime.For waves in wettability the droplet speed decays as λ −2 , which implies that the critical wave speed is proportional to droplet radius over wavelength.This suggests a potential sorting mechanism for a collection of droplets by droplet size.For traveling-wave deformations of the substrate, we could not observe a clear trend in droplet speed w.r.t.wavelength.
Our findings compare and contrast two different mechanisms for steering droplets toward a preferred direction, one based on switchable wettability and one based on externally-induced deformations of the substrate.They demonstrate that the relevant quantities such as droplet speed and wobbling period can be reproduced and interpreted with an analytic model using only one or two fitting parameters, respectively.Moving droplets is relevant for lab-on-a-chip applications.Our work suggest two methods for realizing this.Another attractive possibility to explore are substrates with switchable softness [44] and the occurence of durotaxis.

Figure 1 :
Figure 1: Time series (top to bottom) of snapshots for wettability waves travelling at different speeds: (left) wobbling motion at large speeds and (right) surfing motion of the droplet at small speeds.The greyscale shading indicates wettability and the red lines indicate the contact line of the droplet.

Figure 2 :
Figure 2: Diagrams of (a) the mean droplet speed v ¯, (b) the helper quantity √ 2v wave v ¯, and (c) the wobbling oscillation period T plotted versus v wave for a travelling wave in wettability with wavelength λ = 2 R 0 .Blue dots indicate our results from BEM simulation, the dashed black lines indicate predictions from the nonuniform oscillator model and the dash-dotted red line indicates the critical wave speed v c , which is the sole fitting parameter of the model.

Figure 3 :
Figure 3: Schematic bifurcation diagrams for an appropriate droplet-shape variable plotted versus pattern velocity.(a) The droplet dynamics in the step wettability pattern (c.f.Ref. [14]) corresponds to a saddle-node bifurcation at the critical velocity v c .(b) The travelling wave pattern generates an additional limit cycle that collides with the saddle and node, giving rise to a saddle-node-infinite-period (SNIPER) bifurcation.

Figure 4 :
Figure 4: Diagrams of (a) the mean droplet speed v ¯, (b) the helper quantity v * , and (c) the wobbling oscillation period T as functions of v wave for a travelling deformation wave with wavelength λ = 4 R 0 .Blue dots indicate our results from BEM simulation, the dashed black line indicates the revised nonuniform oscillator model, Eq. (58), the dash-dotted red line indicates the critical wave speed v c , and the dotted red line indicates the parameter ε.The latter two are the fitting parameters of the model.

Figure 5 :
Figure 5: Diagram of mean droplet speed v ¯as a function of λ for travelling waves with speed v wave = 50 R 0 τ −1 .Blue dots indicate results from BEM simulation with waves in wettability (left axis), and green crosses indicate results from BEM simulation with deformation waves (right axis).Connecting dotted and dashed lines are added to guide the eye.The solid black line indicates an inverse-square-law scaling.