Steering droplets on substrates using moving steps in wettability

Droplets move on substrates with a spatio-temporal wettability pattern as generated, for example, on light-switchable surfaces. To study such cases, we implement the boundary-element method to solve the governing Stokes equations for the fluid flow field inside and on the surface of a droplet and supplement it by the Cox-Voinov law for the dynamics of the contact line. Our approach reproduces the relaxation of an axisymmetric droplet in experiments, which we initiate by instantaneously switching the uniform wettability of a substrate quantified by the equilibrium contact angle. In a step profile of wettability the droplet moves towards higher wettability. Using a feedback loop to keep the distance or offset between step and droplet center constant, induces a constant velocity with which the droplet surfes on the wettability step. We analyze the velocity in terms of droplet offset and step width for typical wetting parameters. Moving instead the wettability step with constant speed, we determine the maximally possible droplet velocities under various conditions. The observed droplet speeds agree with the values from the feedback study for the same positive droplet offset.


Introduction
If a droplet can go uphill and move in response to light, where are the limits of its motility? In their seminal experiment, Chaudhury and Whitesides [1] demonstrated that liquid droplets can be driven up a tilted plane against gravity by chemically treating the plane so it gradually becomes more wettable with height. Later, Ichimura et al. [2] showed that sessile droplets start moving in response to gradients in wettability, which are produced by a photo-chemical reaction. These experiments essentially demonstrated an early use of structured light to create motion on the micron scale, an approach which has recently become the focus of intense research [3,4,5]. Light-driven fluid motion [6,7,8,9,10] is an especially favorable control mechanism because of its high precision and controllability through the established experimental methods in optics [11]. Notably, it can be used in combination with or as an alternative to electrowetting techniques [12]. Precise control of fluid motion is foundational for advanced labon-a-chip devices [13,14,15], as well as self-cleaning surfaces [16,17], and printing with sub-droplet precision [18,4].
Recently, significant advances toward lightswitchable substrates with large wettability gradients have been made [19]. First, Zhu et al. [20] demonstrated light-induced droplet motion on a substrate modified by azobenzene-calix [4]arene where the equilibrium contact angle, as a measure for wettability, ranged from 33 • to 110 • . Second, several groups [21,22] have designed arrays of micropillars made from light-switchable polymer material.
Switching the pillars between an upright and buckled shape offers the possibility to reversibly switch the substrate into a superhydrophic state where wettability is vanishingly small.
In this article we implement the boundary-element method [23,24,25,26,27,28] to solve the full threedimensional equations of Stokes flow in order to study droplet motion induced by non-uniform and dynamic wettability patterns. In particular, we apply our method to moving step profiles in wettability and explore the limits of droplet motion induced by such wettability gradients, i.e., we determine the maximally possible droplet speeds under various conditions. Literature provides two research lines relevant to us. First, in two articles McGraw et al. implemented the boundary element method to gain insights into the internal flow fields of axisymmetric droplets [29,30]. They combined theory and experiment to study how droplets relax towards their new equilibrium shape on a substrate after an instantaneous change in the spatially uniform wettability. Second, Glasner [31] implemented the boundary element method to study the dynamics of droplets placed on substrates with static but non-uniform wettability profiles. The study employed a quasi-static approximation where the gas-liquid interface of the droplet is always assumed to be equilibrated w.r.t. to the shape of the contact line. In a very recent example of the same research line, Savva et al. [32] extended Glasner's method by describing the gas-liquid interface with the so-called thin-film approximation for Stokes flow for highly wettable substrates.
Our approach goes beyond both research lines because it combines dynamic and spatially non-uniorm wettability patterns and thereby offers the possibility to study continuous droplet motion. It explores more widely the research avenue recently opened in an experiment by Gao et al. [33] who continuously drove a droplet up a tilted plane using a thermally induced wettability gradient that moved along with the droplet.
We apply our method to investigate a droplet surfing on a moving step in wettablity in two ways: First, by positioning the droplet at a constant distance or offset from the center of the wettability step using a feedback loop, we induce a self-regulated motion of the droplet. Its velocity depends subtly on the properties of the wettability step and its offset from the droplet. Second, by letting the wettability step approach the droplet with a constant velocity, we recover a subset of the steady states observed in the previous study with a feedback loop. This helped us to understand the role of convexity in designing wettability profiles to induce steady droplet motion.
The article is structured as follows. We present the theory of the boundary element method and its implementation in Sec. 2. Our numerical approach is validated in Sec. 3 by studying a sessile droplet on a substrate with homogeneous wettability. We study the droplet surfing on a wettability step using a feedback loop in Sec. 4 and by moving the wettability step with a constant speed in Sec. 5. Finally, in Sec. 6 we present our conclusions.

Boundary element method for dynamic wetting: Theory and implementation
We consider a liquid droplet sitting on a solid substrate with switchable wettability (Fig. 1). The droplet is bounded by two surfaces: Its interface with the substrate and its free interface with the gaseous atmosphere. The intersection where the surfaces meet is the contact

Boundary integral equation for Stokes flow
The boundary element method (BEM) approximates solutions to linear partial differential equations (PDEs) based on their Green's functions [23,25]. For droplets moving on substrates we need to study Stokes flow in incompressible fluids, which obeys two linear PDEs for velocity v and pressure p, [26] where µ is the shear viscosity and r the position vector. Green's function for Stokes flow in 3 is the Oseen tensor O(r) = 1 8πµ|r| where I is the unit matrix. Its associated stress field is a tensor of rank three [26], We are interested in the flow field v inside and, in particular, on the surface of the closed volume of a droplet. On any compact region D ⊂ 3 with boundary ∂ D and outward normal described by unit vector n, the flow field fulfills the following integral equation [26,23,25]: Here, σ is the stress tensor of the fluid, T is the stress field of the Oseen tensor introduced above, and f (r) is a dimensionless coefficient, which assumes different values inside D and at the boundary ∂ D [25]: where ∂ D has a corner with inward solid angle a α . (5) The first and the second term on the right-hand side of Eq. (4) describe the flow initiated by surface forces σn and a layer of force/source dipoles at the surface, respectively. Therefore, they are also called the singlelayer and double-layer potentials. To close Eq. (4), one either needs to eliminate the surface stress force σn or the velocity v by appropriate boundary conditions. We introduce the boundary conditions in the next paragraph. Furthermore, in Sec. 2.3 we describe our discretization of the integral equations so that approximate values for v or σn can be calculated on each piece of ∂ D. Thereafter, we can calculate the velocity field v(r) anywhere in D [24].
The Stokes flow problem is only fully determined with boundary conditions. The droplet boundary consists of two parts, one interface with a solid substrate and a second with a gas phase. This means each part has its own boundary condition. At the interface with the substrate there is a Robin boundary condition (for fluids called the Navier condition), µv + λσn = 0, with a slip length λ [34], which connects tangential velocity and tangential stress, while the normal component of v is zero. A typical value for the slip length is λ ≈ 1 nm [35] much smaller than the droplet dimension. At the interface towards the gas phase we assume zero viscosity of the a The solid angle is calculated locally by considering a small sphere centered at the position r on the droplet surface. The surface cuts out a volume from the sphere directed towards the droplet, which encloses the solid angle α. Thus, α = 2π on smooth parts of the droplet surface and α = 2θ at the contact line, where θ is the contact angle. gas phase so that the tangential stress is zero. There remains a Neumann boundary condition (in hydrostatics the Laplace pressure), n · σn = −2γκ, which locally relates the normal stress component n · σn to the mean curvature κ, where γ is the local surface tension.

Cox-Voinov law for the contact line
Droplets in equilibrium minimize the sum of their surface energies from the interfaces towards the gas phase and substrate, respectively, while preserving their volume [35]. If wetting is energetically more favorable than a dry substrate, the liquid phase covers the substrate with a thin film. Otherwise, a droplet forms with an energetically optimal ratio of the areas between liquid-gas and liquid-substrate interfaces. On homogeneous and planar substrates, the optimal shape obviously is a spherical cap. The balance of all forces acting on the three-phase contact line determines the contact angle θ eq in equilibrium It results in Young's equation, which relates θ eq to the surface tensions γ ij between liquid (l), gas (g), and substrate (s), respectively: Higher wettability means smaller contact angle and according to Eq. (6) that the difference γ sg − γ sl of both surface tensions at the substrate increases, since then it is energetically more favorable to cover the substrate with liquid.
Here, when we talk about dynamic wettability patterns, we will vary ∆γ s = γ sg − γ sl and thereby θ eq in time. On substrates with heterogeneous wettability the optimal droplet shape can be much more complicated including complex shapes of the contact line and varying curvature of the liquid-gas interface [31].
The dynamics of the contact line and the contact angle are not directly determined within the approach outlined so far. Because the droplet surface is not smooth at the contact line, it is not clear which boundary condition to use here. For the Neumann boundary condition towards the gas phase the mean curvature κ diverges, while the Navier condition at the substrate does not contain any surface tension necessary to obtain Young's equation (6) in equilibrium. Therefore, the contact line dynamics has to be determined by an additional relation, which translates the fluid motion in the microscopic contact line region to the macroscopic scale. Several such relations for dynamic wetting exist [35,36,37,38,39,40]. The most well-known among them is the Cox-Voinov law [39,40] derived from hydrody-namic considerations: It relates the difference of the cubes of dynamic and equilibrium contact angles, θ dyn and θ eq , to the velocity of the contact line v contact . The dimensionless coefficient ln(h/λ) is given by the ratio of slip length λ to a macroscopic length scale h. Voinov defines h as the height above the substrate at which θ dyn is measured. However, experiments have demonstrated that ln(h/λ) should be treated as a free parameter characterizing the mobility of the contact line [41]. The velocity v contact determines the component of the fluid velocity perpendicular to the contact line and in the plane of the substrate. We evaluate θ dyn using the surface normal n and the substrate normal e z at the contact line: The resulting velocity v contact of Eq. (7) is included in the discretized BEM problem as a boundary condition by introducing a corresponding stress at the contact line as an additional unknown. Furthermore, in the direction of the substrate normal we prescribe a vanishing velocity and along the tangent of the contact line we set the stress to zero. Together, these three conditions fully determine fluid velocity and stress along the contact line. b On a light-switchable substrate the light patterns affect the surface tension difference ∆γ s = γ sg − γ sl in Eq. (6). Therefore, as mentioned already, we will treat ∆γ s as a heterogeneous and dynamic quantity, which locally determines θ eq and, thereby, deforms the droplet.

Discretization
We discretize Eq. (4) using a triangular mesh to cover all droplet surfaces (solid-liquid and vapour-liquid). From the triangular mesh we define a dual mesh in which each vertex at position r (i) is surrounded by a compact polygonal region C i = C(r (i) ) ( Fig. 2). Together, all the regions C i cover the entire surface. For a sufficiently fine and well-adapted mesh, field variables are approximately constant on each C i and integrals can thus be evaluated piece-wisely. As a result, Eq. (4) becomes a system of linear equations of the form: b A similar approach was previously applied to bubbles on a submerged substrate [42] although the dynamics of submerged contact lines are necessarily distinct from our case. Figure 2: Schematic of triangular mesh (black) in the neighborhood of vertex r (i) and its associated polygonal region C i (red). The region C i is defined by and subdivided into triangles, one of which is indicated in green. The three corners of the triangles are the central vertex r (i) , the point halfway toward a neighboring vertex r ( j) , and the centroid (barycenter) of one of the black mesh triangles.
where v ( j) and t ( j) = σ ( j) n ( j) are the respective values of the flow field and the stress or surface force at vertex j. The tensors A (i j) and B (i j) are piece-wise integrals, i.e., with Kronecker symbol δ i j and identity matrix I, and The linear algebra problem in Eq. (9) is fully solvable by eliminating either v ( j) or t ( j) with the help of boundary conditions at each vertex. It will be convenient to collect all matrices A (i j) and B (i j) into block matrices A and B and all vertex velocities and stress vectors into combined vectors v and t, so that Eq. (9) can be written simply as Av = Bt .
To perform the integrations in Eqs. (10) and (11), we subdivide the polygonal region C i into triangles (Fig. 2) and numerically integrate over each triangle using a Gaussian quadrature stencil with 9 points for nonsingular integrands and with 400 points for singular integrands (i.e. when i = j in Eqs. (10) and (11)) [23,28]. Note that we never evaluate the integrand at the vertex r (i) , where the tensors O or T have integrable singularities, but always choose points inside the triangles for performing the Gaussian quadrature, which is sufficient.
To solve Eq. (12), we first apply the boundary conditions introduced above to remove one of the variables at the vertex, either velocity or surface force. Specifically, at the interface towards the gas phase the tangential component of the surface force is zero while we replace the normal component by the Laplace pressure −2γκ.
To calculate the principal curvatures at each vertex, we use the discrete Laplace-Beltrami operator described in Method C of the review of Guckenberger et al. [43] and the normal vector n (i) follows by taking the sum of the normal vectors of all neighboring triangles and normalizing it. At the interface with the solid substrate we set the normal velocity component to zero and eliminate the tangential component of the surface force by the Robin boundary condition. Finally, at the contact line the velocity component normal to the substrate is zero. The in-plane velocity component normal to the contact line obtains the value prescribed by the Cox-Voinov law, and the tangential component of the surface force is zero. We then collect all the remaining unknown variables on the left-hand side and all the known quantites on the right-hand side and solve the resulting system of linear equations by inversion of the coefficient matrix. Finally, the missing velocity and surface-force variables are calculated from the boundary conditions. In total, solving the linear problem of Eq. (12) yields approximations for v and σn anywhere on the droplet surface. Once these quantities are known, v(r) can be calculated anywhere within the droplet by numerical integration from Eq. (4).

Time evolution
In the discretized version the dynamics of the free droplet surface is specified by the velocity of each vertex with position vector r i . Concretely, the vertex moves with the local normal component of the fluid velocity and an artificial tangential component w (i) t , which is due to mesh optimization and introduced rigorously in Sec. 2.6, Motion tangential to the droplet's free surface is physically irrelevant because it does not affect the shape of the droplet. However, tangential motion can be used to improve mesh fitness, i.e., adaptation of the mesh to the droplet shape.
Note that time appears only through Eq. (13) in our problem because Stokes flow adapts instantaneously to the droplet's shape. Numerically, the dynamics of each r (i) is calculated using a state-of-the-art Runge-Kutta algorithm [44] once the velocity vectors v (i) are calculated with the BEM as described above.
We refine the numerical setting described so far with two additional techniques described in Secs. 2.5 and 2.6: a volume constraint and a mesh optimization method.

Volume constraint
We follow a method first suggested by Alinovi et al. [45] to preserve the volume of the droplet during motion through a side condition for the velocity vectors of the droplet surface. Since volume conservation can be formulated by the chain rule, we obtain We call c (i) = ∇ r (i) V the local constraint vector and for further use we concatenate all c (i) into a vector c, in the same order as v and t in Sec. 2.3). Alinovi et al. estimate c (i) as the product of the local normal n (i) vector and the surface element, which requires a smooth surface around each vertex. c In order to also treat the kink in the droplet surface at the contact line, we instead calculate the gradients ∇ r (i) V directly by numerical differentiation, which gives well-defined constraints for the contact line. Now, with the introduction of a Lagrange multiplier Λ for the constraint of Eq. (14), Eq. (12) is extended to Any solution for the flow field on the droplet surface is now guaranteed to conserve volume for an infinitesimal step in direction of v.

Mesh optimization
During the simulation we optimize the mesh representing the droplet meaning we avoid acute-angled triangles and keep the triangle areas uniform across the whole droplet surface. This guarantees that discretization errors remain minimal. According to Zinchenko et al. [27], a mesh fitness is introduced that is maximized for equilateral triangles. It is a function of all the edge lengths of the triangles and their areas. Thus, improving the mesh fitness w.r.t. the triangle shape minimizes discretization errors. Following Zinchenko et al. [27], we apply two measures at different phases of time integration. First, when calculating the r.h.s. of Eq. (13), we choose the tangential velocity of each vertex on the free droplet surface including vertices on the contact line such that the rate of change of the overall mesh fitness is minimal. d We then move all the vertices on the free droplet surface by one Runge-Kutta time step. Note that the normal velocity, which affects droplet shape and motion, remains unchanged. Second, as a next step we check the fitness of the mesh facing the substrate and if the maximum edge lengths exceeds 150 % of the minimal edge length, the positions of all substrate vertices (without the ones on the contact line) are adjusted by tangential displacements. e Since this technique is completely derived from literature and by its very nature must not affect the dynamics of the droplet, we do not reproduce the algorithm in detail here. In Appendix A, Fig. 11 we display mesh fitness over time for a benchmark case to show the effectiveness of the implemented method.

Nondimensionalization and material parameters
We are able to nondimensionalize our continuum equations by introducing a characteristic length scale R 0 , a time scale τ, and a force scale f c . Using these three parameters, we will then rewrite all quantities a i as dimensionless quantitiesã i . For R 0 we choose the radius of the initial circular base area of the droplet. Furthermore, we remove dynamic and kinematic viscosities from our equations, which then implies the characteristic time and force scales Note that f c is the intrinsic force scale of a Newtonian fluid and the limit of small Reynolds number means that all acting forces are smaller than f c . All relevant quantities can now be written in dimensionless form: curvatureκ = R 0 κ, slip lengthλ = λ/R 0 , stress tensor componentsσ i j = σ i j R 2 0 / f c , and surface tensionsγ ij = γ ij R 0 / f c . Since the reduced time and coordinates aret = t/τ and x = x/R 0 , the respective reduced derivatives becomẽ In Ref. [41] Ruijter et al. considered the relaxation of a droplet on a uniform substrate towards its equilibrium shape assuming that the droplet is always a spherical cap (spherical cap model). Then, the constant volume V or its reduced valueṼ is always strictly related to the momentary contact angle. In particular, with the initial d [41] fitted θ dyn (t) as a solution of Eq. (19) to their experimental data and observed the value 44 for ln(h/λ). They interpret this large value as the result of "an additional source of energy dissipation within the three-phase zone" [41].
contact angle θ 0 we havẽ For the half-sphere droplet with θ 0 = π/2, this givesṼ = 2π/3, as it should. Assuming that the dynamics of the relaxing droplet is completely determined by the motion of the contact line, the authors identify a characteristic relaxation time τ c , which in units of τ reads In Table 1 we collect the material parameters and their dimensionless counterparts for a 90 % glycerine/10 % water mixture for various length values R 0 and use them in the following.

Validation: Switching a homogeneous substrate
We apply our method to a simple case to validate its performance. The simplest applicable test case is the relaxation of a sessile droplet on a homogeneous substrate that switches from one wettability to another (see movie M01 in the ESI † ). Thus the equilibrium contact angle is switched at a specific time and we observe how the droplet relaxes towards the new equilibrium contact angle. Note that our method includes a valid model for the motion of the contact line because we impose the Cox-Voinov law explicitly. However, the motion of the contact line is coupled to the motion of the free surface and the fluid inside the droplet. In the following we illustrate the overall dynamics of the droplet at a specific system. We rely on experimental data published by de Ruijter et al. [41] who observed the spreading of ca. 5-10 mm 3 of a 90 % glycerol/10 % water mixture on PET. All material properties and parameters are collected in Table 1   describe the droplet dynamics by the dynamic contact angle θ dyn . It obeys the following differential equation: The characteristic time τ c introduced in Eq. (18) follows by linearizing Eq. (19) in small changes of θ dyn and including only the coefficients independent of θ dyn or θ eq .
In Fig. 3 we present results for the relaxing glycerin droplet using the material parameters of row 2 in Table 1. We compare our results for the contact angle from the BEM to numerical solutions of Eq. (19) in Fig. 3(a). We find qualitative agreement between both approaches. However, the observable deviations of the two curves, the BEM shows a slower relaxation in the contact angle, reveals that the fluid flow included in our method slows the relaxation further down. Our simulated dynamics also matches the experimental results in both the shape of the relaxation curve of the contact angle and relaxation time, as the experimental points in Fig. 3(a) show. We achieve this despite several uncertainties: We know neither the exact initial shape of the droplet in the experiments of de Ruijter et al. nor the exact volume of liquid they used. We also note during our simulations fluctuations in the droplet volume do not exceed 0.15 % of the initial volume while the droplet base radius and height, displayed in Fig. 3(b), change by over 20 % of their initial values. This cearly demonstrates the robustness of our implementation. Finally, in Fig. 4 we display the stream lines of the fluid at t = 7τ, along which material is transported from the top to the sides of the droplet as it spreads.
We also studied a number of wettability switches on homogeneous substrates to test the technical limitations of our mesh stabilization. Specifically, we were able to realize changes in contact angle of ±15 • when starting from 60 • and ±30 • when starting from 90 • . Notably, the latter case means that our method is able to cross into regions where droplets overhang their base area without any special modifications to treat this case. For example, in Fig. (5) we demonstrate the relaxation of the contact angle from θ dyn = 60 • to θ eq = 120 • . (A corresponding animation is included as movie M01 in the ESI. † ) A clear deviation from the spherical-cap model is visible.

Surfing on a wettability step profile with feedback 4.1 Setup
A droplet placed on a substrate with a gradient in wettability will move in the direction of increasing wettability meaning smaller contact angle. In the following we consider a smooth step profile in wettability (cf. Fig. 6) and move it with the droplet so that the distance between wettability step and the droplet's center of mass stays constant. Thus, the droplet quasi "surfes" on the wettability step (see movie M02 in ESI † ). Using such a feedback between droplet motion and wettability step, we will determine the maximum speed achievable by droplets moving in a gradient of fixed steepness. In our approach we quantify the wettability of the substrate by the local equilibrium contact angle. To realize a step profile in wettability, θ eq ( y, t) (displayed in Fig. 6), we use the logistic step function and write Thus, the wettability step is characterized by its maxmimum and mimimum contact angles θ max eq and θ min eq and by the width ∆ y as a measure for its maximum steepness. Via feedback control we keep the position of the step center at y m (t) = y c (t) − s, where y c (t) is the posi- tion of the droplet center of mass and s the fixed offset between the step and droplet centers. Here, s < 0 means that the droplet lags behind the wettability step. In the following, we study how the four parameters θ max eq , θ min eq , ∆ y, and s influence the droplet speed. In addition, we analyse the deformation of the droplet to a nonequilibrium but steady shape while in motion and the steadystate flow field inside the droplet.
Here, all droplets are simulated using the material properties set out in row 1 of Table 1 which describes a 90 % glycerol/10 % water mixture with R 0 = 100 µm.

Flow field
Regardless of the specific values for each parameter for the wettability step, the droplet eventually settles into a steady state with constant shape and constant flow field. The moving droplet never adjusts to the local equilibrium angles exactly. Rather a difference remains between θ dyn and θ eq , which drives the droplet and determines its steady-state speed Notably, the droplets never enter into a limit cycle with oscillating speed but rather approach a speed value constant in time.
As a representative example, in Fig. 7 we display the internal flow field of a droplet surfing on a wettability step characterized by θ min eq = 90 • , θ max eq = 120 • , and ∆ y = 1/3 at an offset s = 0, i.e., directly on the steepest location. The flow field forms one single vortex filling the droplet in the plane parallel to its direction of motion and two vortices in the plane perpendicular to this direction. The fluid is transported predominantly along the free surface and so, in effect, it rolls on the substrate as the droplet moves forward.

Speed
The steady-state droplet speed v 0 depends on the various parameters of the wettability profile. In Fig. 8(a) we plot the droplet speed versus the offset s for several parameter sets structured by (i) the limiting equilibrium contact angles on both sides of the wettability step and (ii) the step width ∆ y.
First of all, we observe that v 0 increases overall with decreasing ∆ y which is consistent with the expectation that steepness of the wettability step drives the droplet forward. Furthermore, as ∆ y decreases for the same combinations of contact angles, an optimal offset s becomes prominent close to s = −0.2 R 0 (circle and triangle symbols in Fig. 8(a)), meaning the droplet lags behind the step and is pushed forward. Since the wettability gradient is largest at s = 0, we conclude that the steepness of the wettability step does not solely determine droplet speed. Rather, there is a competing influence from the fact that the equilibrium contact angles θ max eq and θ min eq contribute nonlinearly to droplet speed via the contact line. According to the Cox-Voinov law the speed of the contact line depends on the difference of the cubes of θ dyn and θ eq . Thus, v 0 increases with θ eq even when the difference θ dyn − θ eq remains the same. As a result, the droplet speed increases for larger θ eq , which the droplet encounters if it lags behind at s < 0. This is nicely visible in Fig. 8(a), when we compare the 60 • -90 • combination in equilibrium contact angles (red symbols/dash-dotted lines) to the 90 • -120 • combination (purple symbols/full line). Even though the difference between the limiting angles is 30 • in both cases, the droplet speed is larger for the larger equilibrium contact angles.

Deformation
Since the droplets are placed on a non-uniform wettability pattern, their shapes adapt and become asymmetric. There is a competition between the Laplace pressure depending on the local curvature of the free surface, which drives the droplet toward a uniform curvature, and the forces acting on the contact line, which drive the local contact angle toward its equilibrium value. Therefore, the free surface will deviate from a spherical cap. The result of these competing forces can be observed in the variation of the dynamic contact angle θ dyn along the contact line in steady state. To quantify this variation, we plot in Fig. 8(b) the difference of the largest and smallest contact angle, ∆θ dyn = θ max dyn − θ min dyn , versus the offset s for the different parameter sets. Simply put, a larger difference ∆θ dyn indicates a more asymmetric droplet shape. Note that θ min dyn and θ max dyn are realized at the front and the back of the droplet w.r.t. to its direction of motion as indicated in the inset of Fig. 8(b), where the equilibrium contact angles are extremal. At these locations the contact line velocity v contact must equal the droplet speed v 0 in steady state. So the difference between the extrema of θ dyn and θ eq along the contact line determines v 0 .
We observe that ∆θ dyn is primarily determined by large variations in θ eq , which makes sense. For example, the wettability profiles with θ min eq = 45 • and θ max eq = 90 • [green symbols/dashed lines in Fig. 8(b)] cause larger variations ∆θ dyn than step profiles with θ min eq = 90 • and θ max eq = 120 • [purple symbols/solid lines in Fig. 8(b)] for otherwise same parameters.
In addition, for ∆ y/R 0 = 1/3 and 1/2, ∆θ dyn shows a strong dependence on the offset s. Interestingly, for these cases the droplet asymmetry decreases notably as s approaches the speed optimimum s = −0.2 R 0 . Apparently, a more symmetric droplet has larger differences between θ dyn and θ eq along its contact line and therefore achieves higher speeds v 0 .
In summary, we observe droplets settle into a rolling motion when they are a placed on a wettability step profile with a constant offset to the location of steepest wettability gradient. Their limiting speed increases with the steepness of the gradient. For large steepness there is an optimal offset so that the droplets are pushed forward by regions of lower overall wettability. Their shapes become more asymmetric in response to steeper gradients, while the optimal offset for speed is correlated with a more symmetric shape.

Setup
As an alternative to the feedback method described above, we study the simpler setup where the wettability step of Eq. (20) moves with constant speed. Here, two scenarios are possible: If the droplet can keep up with the moving wettability profile, it will surf on the pattern, otherwise it will fall behind, stop moving, and eventually adapt its shape to the lower wettability of the substrate.
To gain insight into these scenarios, we re-examine one particular pattern from the previous part: a wettability step with θ min eq = 90 • and θ max eq = 120 • , and ∆ y = 1/3. We place a glycerol droplet with parameters according to row 1 of Table 1 well before the step and, unlike before, let the step approach the droplet with a constant speed v s independent of the position of the droplet.

Long-time dynamics
We first position the droplet at t = 0 with an initial lead s 0 = 2 R 0 before the steepest point of the step, which means the droplet is on the more wettable side of the step. This gives the droplet some time to adapt to a surfing shape under the influence of the approaching step. We then observe how the droplet lead,  if the droplet can follow the wettability step, the droplet lead approaches a constant value.
In Fig. 9 we display trajectories leading to either droplets surfing on the wettability step or droplets that cannot follow. Above a step speed of roughly v s = 1.1 · 10 −3 R 0 /τ droplets cannot keep up with the step meaning their droplet lead tends to minus infinity. Below, droplets match the step speed and surf, meaning their droplet lead approaches a constant value. The transition between both scenarios occurs close to the maximum droplet speed of 1.1·10 −3 R 0 /τ, which we observed using the feedback mechanism in Sec. 4.3, Fig. 8(a). Furthermore, we note that the droplet lead of 0.24 R 0 observed for v s = 10 −3 R 0 /τ is close to the lead 0.2 R 0 , which occured for a similar speed of v 0 = 1.01·10 −3 R 0 /τ in Sec. 4.3 [purple triangles in Fig. 8(a)].
So far, we can draw two conclusions from our investigations in comparison to the feedback studies of the previous section. First, for a given step speed a droplet sitting ahead of the approaching wettability step assumes the offset s given in the upper curve of Fig. 8(a) by the branch to the right of the maximum. Second, droplets cannot follow wettability steps with larger speeds.

Preferred surfing state
How does the final surfing state of the droplet depend on the initial droplet lead s 0 ? To answer this question, we placed a droplet with spherical cap shape and contact angle θ dyn = 90 • at different positions relative to the wettability step, which moved with different veloc-  Fig. 8(a)]. Same parameters as in Fig. 9 are used.
ities. The result is presented in the state diagram s 0 vs. v s of Fig. 10. The black solid-dashed line indicates the result from the feedback analysis in Fig. 8(a). The yellow triangles mean that the droplet lead decreases. For v s ≤ 1.1 R 0 /τ, when the lead reaches the upper branch of the black curve, the droplet is taken along by the step and moves with constant speed. For v s > 1.1 R 0 /τ the droplet cannot follow the wettability step and is left behind (see movie M03 in ESI † ). The blue triangles inside the black curve indicate droplets that initially move faster than the wettability step. The droplet lead s increases until it hits again the solid black line and the droplet again moves with constant speed v s (see movie M04 in ESI † ).
So the droplet speed determined in our feedback analysis as a function of the offset or droplet lead s gives exactly the outcome when we use a constant step speed. However, the black curve in Fig. 10 has a stable branch (solid line) indicating that the droplet will mostly move in front of the wettability step, while the positions given by the dashed line are unstable. This behavior is reminiscent of a saddle-node bifurcation near v s ≈ 1.1 R 0 /τ. Our results also show that for droplets on wettability steps feedback has a stabilizing effect on an otherwise unstable state (see, e.g., Ref. [46]).
The (in)stability of either branch is intuitively explained by the steepness and local curvature of the sigmoid curve, which we use for the wettability profile (see Fig. 6): By definition, a leading droplet is on the convex part of the sigmoid which means, as it falls behind it is exposed to an increasing gradient which speeds it up again. Conversely, a lagging droplet is on the concave part of the sigmoid and as it falls behind it is exposed to a decreasing gradient which slows it down further. Only close to the turning point also lagging positions are stable. We explain this by the fact that the droplet is an extended object. The unstable and stable positions meet at the bifurcation point, which is near the steepest part of the sigmoid. There, the droplet is exposed to the largest possible gradient and therefore achieves its maximum speed.

Conclusions
In this article we have applied the boundary element method and developed a numerical scheme to simulate droplets on substrates with switchable and non-uniform wettability realized, for example, by light-switchable surfaces. A strong emphasis was put on developing a stable grid on the droplet surface. We are able to simulate droplets on the whole range of possible contact angles and, in particular, investigated speed and deformation of droplets under the influence of moving wettability steps. A crucial ingredient for simulating the dynamics of droplets is the Cox-Voinov law, which governs the motion of the three-phase contact line. In addition, using the boundary element method allows us to calculate flow fields inside the droplet and not just on its surface and thereby provides a further understanding of the droplet motion.
We first applied our method to sessile droplets and exposed them to a spatially uniform but instantaneous change in wettability. We carefully compared our findings to experimental results from Ref. [41] and found quantitative agreement. The spherical cap model, where the dynamics is governed by the dynamics of the contact line via the Cox-Voinov law, shows small but noticeable differences.
We then applied our numerical scheme to droplets surfing on moving wettability steps. First, using a feedback loop we kept the center of the wettability step at a constant distance or offset from the droplet center. For a range of offsets and step widths, this gives rise to droplets surfing at constant speed on the wettability step into the direction of higher wettability. Inside the droplet a single vortex forms such that the droplet rolls on the substrate. For shallow wettability steps the resulting droplet speed hardly depends on the droplet offset in the range investigated. However, for steeper steps the speed develops a maximum at offsets behind the steepest point of the step so that the droplet experiences regions of lower overall wettability. Additionally, while droplets surfing on steep steps are less symmetric compared to more shallow steps, droplets surfing with the maximum speed are more symmetric compared to slower droplets on the same wettability profile.
In a second study we moved the wettability step at a constant speed with various initial offsets between the step and droplet center. Beyond the maximal droplet speed from the feedback study, steady surfing is not possible. The droplet cannot follow the wettability step and is left behind. Below the maximal surfing speed, droplets always settle to the same offset as observed in the feedback study. However, only the leading positions where the droplet is pushed forward are stable while the lagging positions are unstable in the absence of feedback. We explain this by the curvature of the smooth wettability step: In the leading position the convex variation of the wettability step generates an effective restoring force when the droplet is left behind, while in the lagging position the concave variation cannot stabilize the surfing droplet against small displacements away from the step center.
Our findings demonstrate how a dynamic wettability profile, for example, controlled by structured light can be used to continuously drive droplets forward. We have investigated here the important parameters to optimize the speed of droplets w.r.t. the shape of the wettability profile and its driving speed. Because a single wettability step can be set up to move droplets in any direction, in principle, we provide here a tool to steer droplets along arbitrary paths and with variable speed on a light-switchable substrate. In the future, we plan to investigate droplet dynamics in more complex spatiotemporal wettability patterns motivated by the possibilities of structured light [11].

Conflicts of interest
There are no conflicts to declare.

A Mesh degradation
Over the course of a simulation, mesh fitness M (see Ref. [27] for definition) degrades due to physically necessary changes in the mesh shape. In Fig. 11 we give an Figure 11: The mesh fitness M for a sessile glycerol droplet decreases over time in response to a wettability change at t = 0 s. example of mesh degradation during the initial relaxation of a droplet on a uniform substrate presented in Sec. 3. Notably, M initially decreases roughly linearly while the droplet changes shape, as visible in Fig. 3. However, M decreases further thereafter which implies that, even in the absence of further deformations, numerical error still compounds. Small numerical errors have a compounding effect on M because M is a nonlinear function of vertex positions, and therefore displacements of similar magnitude have a small effect on the initial (effectively optimal) initial mesh for t < 1 s and a disproportionally larger effect on the already deformed mesh for t > 1 s.

B Parameters used for movies in ESI
All movies visualize simulations with material parameters from row 1 of Table 1. Movie M01 corresponds to the last example provided in Sec. 3 and Fig. 5, which means it shows an instantaneous change in uniform wettability from θ eq = 60 • to 120 • . Movie M02 corresponds to the purple triangle at offset −0.2 R 0 in Fig. 8(a). Movie M03 shows a droplet driven by a wettability step from 90 • to 120 • with steepness ∆ y = 1/3, step speed of v s = 2 · 10 −3 R 0 /τ, and an initial lead of s 0 = 2 R 0 . Movie M04 corresponds to the same parameters, except v s = 1 · 10 −3 R 0 /τ. Movie M05 corresponds to the parameters given in Fig. 7. Note that the movies each run at different speeds.