Droplets on substrates with oscillating wettability

In recent decades novel solid substrates have been designed which change their wettability in response to light or an electrostatic field. Here, we investigate a droplet on substrates with oscillating uniform wettability by varying minimium and maximum contact angles and frequency. To simulate this situation, we use our previous work [Grawitter and Stark, Soft Matter 17, 2454 (2021)], where we implemented the boundary element method in combination with the Cox-Voinov law for the contact-line velocity, to determine the fluid flow inside a droplet. After a transient regime the droplet performs steady oscillations, the amplitude of which decreases with increasing frequency. For slow oscillations our numerical results agree well with the linearized spherical-cap model. They collapse on a master curve when we rescale frequency by a characteristic relaxation time. In contrast, for fast oscillations we observe significant deviations from the master curve. The decay of the susceptibility is weaker and the phase shift between oscillations in wettability and contact angle stays below the predicted $\pi/2$. The reason becomes obvious when studying the combined dynamics of droplet height and contact angle. It reveals non-reciprocal shape changes during one oscillation period even at low frequencies due to the induced fluid flow inside the droplet, which are not captured by the spherical-cap model. Similar periodic non-reciprocal shape changes occur at low frequencies when the droplet is placed on an oscillating nonuniform wettability profile with six-fold symmetry. Such profiles are inspired by the light intensity pattern of Laguerre-Gauss laser modes. Since the non-reciprocal shape changes induce fluid circulation, which is controllable from the outside, our findings envisage to design targeted microfluidic transport of solutes inside the droplet.


Introduction
The shape a liquid droplet forms on a flat surface is determined by the wettability landscape of the surface. On more wettable parts of the surface, the droplet spreads out and on less wettable parts it contracts [1]. In the past two decades researchers have developed substrates the wettability of which can be controlled such that patterns that change in space and/or time can be created [2,3,4]. For example, the wetted substrates become switchable due to a single layer of light-responsive molecules [2,5,6], electro-responsive molecules [7], an array of light-responsive pillars [8,9,10], or more complex nanostructures [11,12]. Regardless of the specific mechanism, controlling the dynamics of a liquid drop by time-varying wettability patterns is not yet fully explored, despite its importance for lab-on-a-chip devices [13] and for targeted deposition of solutes, for exampe, in printing devices [14].
So far, much research has aimed at understanding wetting on substrates with static non-uniform wettability patterns. Early theoretical studies used perturbation methods to analytically estimate the influence of small local wettability gradients on droplets with simple circular or cylindrical shapes [15,16]. Later, experimental and numerical work investigated more complex shapes which occur, for example, when a droplet crosses a static step in wettability [17,18], flows over two neighboring stripes of increased wettability [19], over a checkerboard pattern [20,21], or random spatial fluctuations in wettability [21]. From the perspective of the droplet these patterns also become time-varying if the droplets starts to move, for example, on an inclined substrate [22]. Furthermore, a droplet may trigger a chemical reaction with the substrate and thereby create a wettability gradient, which moves it forward [23]. Taking into account the switchable substrates introduced above, we have shown recently how a chemically inert droplet responds to moving steps in wettability [24]. Thus, the droplet's motion is under external control.
In this article we present a theoretical investigation how a small liquid droplet behaves on a substrate with oscillations in uniform wettability. We also give a brief outlook toward its behavior on a substrate with oscillations of a non-uniform wettability pattern. Generally, small sessile droplets form spherical caps on flat substrates with uniform wettability because that shape minimizes its total surface energy on a flat substrate [25]. However, when the substrate's wettability oscillates, the droplet continually tries to follow but cannot relax to its equilibrium shape. Its continuous motion in turn gives rise to internal flow in the droplet. We are interested in both-the droplet shape during oscillation and the accompanying internal flow. Since at small length scales viscous forces dominate inertia within the droplet, the internal fluid flow is described by the Stokes equations [26]. We solve these using the boundary element method, the implementation of which we have described in a previous work [24]. At the edge of the droplet-substrate interface (the contact line), we use the Cox-Voinov law to calculate the velocity of the contact line [27]. We have previously validated and applied our method to droplets that are steered by moving steps in wettability [24].
The case of droplets exposed to oscillating wettability is distinct from droplets on vertically vibrated substrates, which have been studied in Refs. [28,29,30]. For the droplet such vibrations only play a role in the presence of inertia or external forces. They affect the droplet as a whole, while oscillations in wettability only act via forces at the contact line. Vertical vibrations have been shown to give rise to ripples that travel up the side of the droplet and to generate higher-harmonic deformations of the contact line for water [28,29] as well as for mercury droplets [30]. At large amplitudes of the vibrations the droplet breaks up by ejecting small amounts of liquid at its top [29]. As we study wettability oscillations in the absence of inertia, we will not observe such extreme phenomena in our case.
Our theoretical approach stands alongside two other continuum approaches to dynamic wetting. In the first approach, one uses a thin-film equation to evaluate the droplet dynamics via its height profile [31,32,33,34], which means the contact angle should be small and cannot exceed 90 degrees. Another approach, which we will discuss in detail below, is the spherical cap model. It constraints the shape of the droplet to a spherical cap [27,35] and does not capture fluid flow within the droplet. The spherical cap shape is motivated by the equilibrium shape of a droplet on a substrate with uniform wettability. With our approach, we are able to evaluate the applicability of the spherical cap model and investigate the internal flow field of the droplet.
Our findings add to the microfluidics toolbox [36,37,38] another way to interact with and manipulate droplets by placing them on substrates with oscillating wettability. Specifically, we find that the contact angle oscillations of the droplet decrease with increasing frequency. For slow oscillations this can be well described by the spherical-cap model, which even provides a characteristic time scale to map the oscillations onto a common master curve. However, the master curve is no longer applicable for fast oscillations. A more detailed study of the droplet dynamics in terms of two shape variables, such as contact angle and droplet height, reveals that they oscillate out-of-phase with each other. Thus, the droplet performs a non-reciprocal motion during one oscillation period, which cannot be described by the spherical-cap model. It is due to fluid flow within the droplet, which gives rise to fluid circulation within the droplet.
Our article is structured as follows: In Sect. 2 we reiterate the theoretical basis of our boundary element method applied to dynamic wetting. In Sect. 3 we describe and discuss our findings in detail. First, in Sect. 3.1 we present the basic phenomenology of the droplet oscillations. We analyze them using linearresponse theory and the spherical-cap model, where the contact angle serves as a single shape characteristic. Second, in Sect. 3.2 we look at the coupled dynamics of contact angle and height, reveal the non-reciprocal motion of the droplet, and discuss its implications for the internal fluid flow. Third, in Sect. 3.3 we study a closely-related example of a droplet on a substrate with an oscillating non-uniform wettability pattern. Finally, we conclude in Sect. 4.

Simulation method
The motion of a droplet consisting of an incompressible simple liquid is completely described by the dynamics of its interfaces: the gas-liquid and the solid-liquid interface. The motion of any point s on the interfaces is governed by the fluid velocity field at this point: In the following we summarize how we determine the velocity field both at the droplet surface and its interior.

Stokes flow
We consider droplets in which viscous drag dominates inertia, which is the limit of Stokes or creeping flow [39]. The equations govering the velocity field v of Stokes flow are where µ is viscosity and the second equation is the incompressibility condition which constrains the pressure p. These differential equations can be restated as boundary integral equations [40] using the Oseen tensor O and the associated stress field T: where ∂ D has a corner with inward solid angle α.
where ∂ D is the (time-dependent) surface of the droplet. Because the equation relates velocity v and stress σn, for any surface point r either variable must be prescribed by boundary conditions.

Boundary conditions
On the liquid-gas interface the normal stress balances surface tension forces due to mean curvature κ of the interface, i.e., σn = γκn, where γ is the surface tension of the liquid-gas interface. On the solid liquid interface (at the substrate) two conditions apply: Firstly, the interface cannot deform along its normal e z and therefore v z = 0. Secondly, roughness of the substrate introduces a small amount of slip with slip length λ, which we account for by setting λσn = µv tangential to the interface [41].
As a boundary condition on the contact line, we choose its velocity along the substrate according to the Cox-Voinov law [27,42] v contact = γ where θ eq is the equilibrium contact angle, which defines the wettability of the substrate, and θ dyn is the dynamic or actual contact angle. Only with this separate boundary condition is the problem well-posed because the three-phase contact line is neither clearly part of the gas-liquid interface nor the solid-liquid interface [43,26]. Note also that the Cox-Voinov law excludes the effects of contact angle hysteresis [44,45]. The boundary conditions introduce several material parameters in addition to viscosity µ. For our simulations we choose dimensionless parameters. They correspond, for example, to a droplet with an initial radius R 0 = 100 µm of its circular base area and made of a 90% glycerol and 10% water mixture. This reference system was studied in experiments by de Ruijter et al. [35]. The mixture has µ = 209 mPa · s, kinematic viscosity ν = µ/ρ = 169 mm 2 · s −1 , γ = 65.3 mN · m −1 , λ = 1 nm, and ln(h/λ) = 44. The latter value was observed by de Ruijter et al. by fitting the spherical cap model with Eqs. (10) and (11) mentioned below to their experimental data. Furthermore, we choose the initial contact angle θ dyn to be the time average of θ eq . We calculate and report our data in units of R 0 for length, τ = R 2 0 /ν for time, and F 0 = νµ for force. Thus the remaining dimensionless parameters arẽ γ = γR 0 /F 0 = 0.19 andλ = λ/R 0 = 10 −5 .

Boundary element method
To solve for the velocity field at the droplet's interfaces, we construct a triangular mesh and discretize the integral equation. We then integrate the dynamic of the mesh in time using an adaptive 5th order Runge-Kutta method [46]. The full details of our boundary element method are provided in Ref. [24]. Briefly, to discretize Eq. (3) we divide the droplet surface into polygonal regions, each with a vertex at its center. The polygons are then decomposed into triangles and we integrate separately over each triangle using Gaussian quadrature [47] with 400 sampling points for singular integrands and 9 sampling points for nonsingular integrands. Once we have solved the discretized equation and the surface velocity field is known, we can use the boundary integral equation (3) to evaluate the flow field in the interior of the droplet. Similar numerical approaches have been used to study dewetting of polymer microdroplets [48,49] and for bubbles on a solid surface under the influence of an acoustic field [50].

Droplet on a substrate of oscillating uniform wettability
We consider a droplet on a substrate, where the uniform wettability expressed by the equilibrium contact angle θ eq (t) oscillates with a frequency f between a minimum (θ min eq ) and maximum (θ max eq ) value: For further use below, we note that the wettability oscillation is invariant under time reversal, up to a constant phase shift. We now discuss how the droplet reacts on oscillations in the wettability and compare our numerical results to the outcome from the spherical-cap model. Then, we show that the induced flow in the droplet is nonreciprocal so that it effectively pumps fluid during one oscillation cycle.

Phenomenology
After an initial transient behavior, the droplet settles into a periodic deformation which oscillates with the same frequency as the wettability. For three frequencies separated by a factor of 10, we display the dynamic contact angle in Fig. 1. We observe that the range covered by the dynamic contact angle decreases with increasing frequency. For slow oscillations the contact angle nearly follows the prescribed equilibrium value of the substrate wettability, while for fast oscillations it barely varies leading to an almost steady droplet shape. Note that in our example the droplet does not oscillate about the mean value of θ max eq = 120 • and θ min eq = 60 • . Furthermore, we observe a phase shift between the oscillating equilibrium and dynamic contact angles, which also depends on f . Note, for different combinations of θ max eq Figure 2: Absolute susceptibility |χ| as a function of oscillation frequency f = ω/2π in units of τ −1 for three combinations of θ max eq and θ min eq . The inset shows the corresponding phase shift ∆ϕ.
and θ min eq we provide videos M01-M03 in the Supplementary Material, where a relatively small frequency of To quantify the response of the droplet to the oscillating substrate wettability, we introduce the nonlinear susceptibility χ = |χ| exp(i∆ϕ) with absolute value |χ| and phase shift ∆ϕ. The imaginary unit is indicated by i. We extract |χ| from To calculate the phase shift ∆ϕ( f ), we determine the first Fourier coefficient α dyn of the dynamic contact angle θ dyn , where s → ∞ insures that the oscillations of θ dyn (t) are steady. For the same time intervall we also determine the complex amplitude α eq of the prescribed equilibrium contact angle sand then calculate the phase shift between oscillating wettability and dynamic contact angle from where arg means the phase angle of the complex amplitude.
In Fig. 2 we plot the absolute value |χ| over frequency for different combinations of θ max eq and θ min eq . All three curves show the expected decrease of |χ| with increasing f . Furthermore, for larger difference θ max eq − θ min eq and larger values of the equilibrium contact angles, the curves are shifted to larger frequencies but roughly have the same shape. This suggests by rescaling frequency appropriately, they fall on a master curve. The inset of Fig. 2 plots the corresponding phase shift. For small frequencies of the oscillating wettability the phase shift tends towards zero for the green curve meaning that the dynamic contact angle follows the prescribed wetting angle instantaneously.
To gain more insights and also pursue the idea of the master curve further, we compare our observations to the spherical-cap model [35], where the droplet always keeps the shape of a spherical cap and the dynamics is solely governed by the Cox-Voinov law for the contact line. This will allow us to distinguish phenomena inherent in the contact-line friction from phenomena which are due to the freely deformable liquid-gas interface and the initated fluid flow in the droplet as determined in our boundary element method. For slow temporal variations in the wettability we expect the spherical-cap model to be valid since fluid flow in the droplet is weak, while for larger frequencies deviations should occur. We now go into more detail.
In the spherical-cap model the shape of the droplet is constrained to a single degree of freedom for which it determines a dynamic equation [35]. Here, we follow Ref. [35] and express the model in terms of the dynamic contact angle θ dyn , with To gain some inside and derive a characteristic relaxation time of the spherical-cap model, we linearize it around θ dyn = θ eq in ∆θ = θ dyn − θ eq : The derivative of g is a characteristic relaxation rate τ −1 0 . But, in our case θ eq is a function of time. Nevertheless, to have a constant rate, we calculate the derivative at the mean equilibrium contact angleθ eq = (θ max eq −θ min eq )/2 and obtain (1 − cosθ eq ) 2 (2 + cosθ eq ) 4 (13) Now, we approximate Eq. (12) by using the constant τ 0 instead of the exact derivative of g. Rescaling time by τ 0 , we finally arrive at dθ dyn which is the parameter-free linearized model. Below we will demonstrate that it very nicely fits our computational results for low frequencies and it is the basis for identifying a master curve for |χ(ω = 2π f )|.
Using time rescaled with τ 0 also in the full sphericalcap model, we can rewrite Eqs. (10) and (11) in a nondimensionalized form as Note that here the r.h.s. is independent of the liquidgas surface tension γ lg , the viscosity µ of the fluid, and the Cox-Voinov parameter ln(h/λ), which determines the contact line mobility. All these parameters are subsumed in the relaxation time τ 0 . The only remaining parameters are θ min eq and θ max eq , which also determineθ eq . We now explore the linerarized model. By taking the Fourier transform of Eq. (14), one finds the dynamical susceptibility χ(ω), which quantifies how θ dyn responds to an oscillation in θ eq :θ dyn (ω) = χ(ω)θ eq (ω) with The absolute value |χ| reads and the complex phase ϕ is Note that in the region where the linear model applies, ϕ(ω) is identical to the phase shift ∆ϕ( f = ω/2π) between θ eq (t) and θ dyn (t) introduced above. Because Eqs. (17) and (18) do not depend on any material parameters besides the characteristic relaxation time τ 0 , they are candidates for the master curves for our simulation results in Fig. 2. We simply need to rescale frequency by τ −1 0 calculated for the specific values of θ min eq and θ max eq . In Figs. 3(a) and (b) we display the master curves as dashed lines together with the results from our boundary element method (BEM) using rescaled frequencies. First, in Fig. 3(a) we observe that the BEM results all fall on a common master curve for f < τ −1 0 , which perfectly matches the linear model in this range. For f > τ −1 0 the BEM results do not follow a common master curve and they deviate from the linear model. We observe that the absolute susceptibility enters an algebraic decay with an approximate exponent −0.5 rather than −1 as predicted by the linear model. Second, in Fig. 3(b) we similarly observe that the BEM results approach the linear model for small f , however they start to deviate siginificantly for f τ 0 > 3·10 −1 . The phase shift varies between roughly 0.3π and 0.4π rather than approaching 0.5π as predicted analytically by the linear model in Eq. (18) and indicated as dashed line.
The phase shifts beyond f τ 0 > 3 · 10 −1 apparently depend on θ max eq and θ min eq . To interpret these observations, we distinguish two regimes: a low frequency regime with f τ 0 < 1 and a high frequency regime with f τ 0 > 1. In the low frequency regime, the oscillations are sufficiently slow so that the droplet can adapt its shape and keep it close to a spherical cap. Thus, the linear spherical-cap model is valid. In fact, in the limit of vanishing f it becomes exact as the dynamics becomes quasistatic. According to linear response theory, the imaginary part Imχ, which in our linear model from Eq. (16) reads quantifies dissipation. Because Imχ as well as the phase shift ∆φ from Eq. (18) are linear in ω at small frequencies, this implies that almost no work performed on the droplet through slow wettability oscillations is dissipated by the friction of the contact line.
In the high frequency regime, the droplet deviates far from the linear model since small |χ| means θ dyn barely tracks θ eq (t). But also the full spherical cap model is unable to predict the observed behavior for |χ| and ∆ϕ and, in particular, the deviation from a common master curve. This implies that the droplet shape deviates from a spherical cap when fast wettability oscillations are applied. These deviations occur close to the contact line, which can move quickly without displacing much liquid and thereby bends the free surface, i.e., it locally in-or decreases curvature, which is precluded in the spherical cap model. This explains why θ dyn is more susceptible at large f , meaning it lies above the prediction of the linear model: A small adjustment of the position of the contact line can drastically alter θ dyn and increase in a short time the surface energy of the droplet relative to the equilibrium reference shape for a given θ eq . Thus the work performed on the droplet is not completely dissipated. Therefore, the phase shift angle ∆ϕ is below the value π/2, which is expected for a complete dissipation and which is predicted within the spherical-cap model for large frequencies.
The spherical cap model only includes dissipation at the contact line. It does not account for viscous friction, friction at the substrate interface, or the elasticity of the free surface. All of these are present in the BEM results, however. We now want to focus on these contributions by turning our attention toward the internal flow of the droplet and toward its deformation from the spherical cap shape.

Deformation and pumping
The shape oscillations of the droplet are accompanied by internal fluid flow. When the wettability increases, the droplet wets more area on the substrate and fluid moves from the top of the droplet through its center to the contact line. When the wettability decreases, the droplet wets less area on the substrate and fluid moves in the opposite direction.
However, this back-and-forth does not cancel out completely and there is a net displacement of fluid after each period, i.e., fluid is pumped within the droplet, as illustrated in Fig. 4. To quantify the net displacement, we place point-like tracer particles in the droplet and track their motion for one full period of oscillation, so that to each starting point r 0 we can assign a displacement where r(t) is the solution of the differential equatioṅ r(t) = v(r, t) with initial condition r(0) = r 0 and v(r, t) is the interior velocity field of the droplet when it is steadily oscillating. In Fig. 4 we display an example of d(r) that is representative for all studied cases. Qualitatively, it shows a circulation of fluid inside the droplet with fluid travelling up through the center of the droplet and down along its free surface in a single toroidal vortex, which covers the interior of the droplet completely. We return to this observation after considering the shape of the droplet. To quantify changes in the droplet shape, one can use, e.g., the dynamic contact angle θ dyn or the droplet height h. In Fig. 5 we plot two oscillation cycles of θ eq (t) in graph (a) together with the corresponding h(t) and its time derivativeḣ(t) in graph (b). We first observe that the height of the droplet adjusts faster to a decrease in wettability than to an increase; meaning, the upward slope of h(t) corresponding to an increase of θ eq (t) is larger than the magnitude of the downward slope. We understand this by studying the mobility of the contact line as a function of the dynamic contact angle. When the droplet is equilibrated, i.e., θ eq = θ dyn and a small change in wettability occurs, i.e. θ eq → θ eq + δθ eq , we calculate the contact line mobility m by linearizing the Cox-Voinov law, Eq. (5). The velocity of the contact line becomes v contact ≈ −mδθ eq with m = γθ 2 eq 3µ ln(h/λ) (21) Apparently, the mobility m can be smaller or larger for droplets with the same θ dyn , depending on θ eq . If the wettability is increasing (smaller θ eq ) m is decreased. Respectively, if the wettability is decreasing (larger θ eq ) m is increased. Thus, when wettability increases over time, contact line mobility decreases and the droplet's shape adjusts more slowly than when wettability decreases over time. This explains the behavior of h(t) in response to θ eq (t).
However, varying rates of deformation are insufficient to explain a net pumping of the liquid because the equations of Stokes flow, Eq. (2), are independent of time. Therefore, a different aspect of the deformation must be responsible for the pumping.
Unlike the spherical cap model our BEM simulations are not constrained to a single degree of freedom. So as a minimal extension for describing the temporal shape variations of the droplet, we investigate the coupled dynamics of two degrees of freedom, contact angle θ dyn and droplet height h. In Fig. 6 we represent the droplet dynamics in the configuration space spanned by these two variables; one oscillation of the droplet corresponds to a closed trajectory. Interestingly, the non-zero area enclosed by the trajectory reveals the dynamics of the droplet as non-reciprocal, meaning under time-reversal the dynamics looks different. In our concrete case the droplet assumes slightly different shapes during increasing and decreasing wettability in the course of one period. The non-reciprocal dynamics is clearly due to the flow field generated inside the droplet. Since the spherical-cap model only has one dynamic variable, its dynamics can only be reciprocal. The dashed line in Fig. 6 shows the model prediction for vanishing frequency. a Note the non-zero area also a Note that due to the discretization of the droplet surface, all θ dyn of the simulated curves are systematically shifted downwards by a Figure 6: Closed trajectories in configuration space projected onto the h-θ dyn plane for various frequencies f (see legend) for θ min eq = 60 • and θ max eq = 120 • . The limiting behavior of the spherical-cap model with f τ 0 1 is indicated by the black dashed line. means that contact angle and height oscillate out-ofphase with each other. From the study of microswimmers at vanishing Reynolds numbers, we know Purcell's scallop theorem [51] which states that non-reciprocal shape changes are needed for microswimmers to move forward. Similarly, the pumping displacement mentioned in the beginning is linked to the non-reciprocal droplet deformation.
To quantify the non-reciprocal shape dynamics of the droplet, we take the area enclosed by the trajectory in configuration space, which, in general, is highdimensional. However, projecting this trajectory into the two-dimensional space spanned by droplet height h and contact angle θ dyn , we can immediately calculate the projected area as where the limit s → ∞ ensures that h(t) and θ dyn (t) perform steady oscillations. In the following, we call the parameter A shape non-reciprocity. The closed trajectories or loops in the configuration space presented in Figure 6 extend further in both dimensions and enclose a larger area A as frequency decreases. An increasing area A means the dynamics of the droplet shape becomes less reciprocal. Figure 7(a) displays the dependence of A on frequency f for all three cases of equilibrium contact angles. We always observe that A decreases for large f , while in the case of θ min eq = 60 • and θ max eq = 120 • the non-reciprocity A has a maximum and decreases toward small f . The latter is expected since in the quasistationary case of vanishing f the motion has to become reciprocal. Notably, even at small frequencies A is relatively large which means even though the dynamic response approaches that of small angle. Otherwise, they should center on the dashed line. We can directly understand the non-zero values for A by considering the basic mechanisms driving the droplet. When wettability changes, it gives rise to uncompensated Young forces [52] at the contact line. The contact line starts to move, which brings θ dyn closer to θ eq and relaxes the Young forces. At the same time, the free surface is bent locally thereby introducing uncompensated surface stresses in the vicinity of the contact line. Those stresses redistribute liquid inside the droplet and eventually affect droplet height h. So, while θ dyn adjusts directly to the changes in wettability, the effect on h is mediated by the initiated flow and, therefore, delayed.
In the quasistatic limit, f → 0, that delay becomes negligible relative to f ; contact angle and height oscillate in synchrony and the droplet's behavior approaches the spherical cap model, Eq. (15). In the limit f → ∞ the surface stresses relax not by redistributing liquid, but because wettability quickly returns to its original value, thereby eliminating the uncompensated Young forces before the contact line moves significantly. Here, the droplet hardly oscillates and thereby approaches a static spherical cap shape. So, in both limits the droplet assumes spherical shapes and only intermediate frequency values f cause a significant deviation from this form.
We now investigate the liquid displacement d(r 0 ) inside the droplet in more depth. To show the link between non-reciprocity and displacement quantitatively, we consider the median pumping speed v c = f 〈d〉 w.r.t. the interior of the droplet, i.e., the 50th percentile value of the spatial distribution of d = |d(r 0 )| multiplied by frequency f . In Fig. 7(b) we observe that v c remains constant for all f with some fluctuations which means it is purely determined by the material properties of the liquid and substrate and not the oscillation frequency. We already mentioned since the spherical-cap model only has a single degree of freedom, it cannot exhibit any circulatory pumping according to the scallop theorem. So, just like the non-reciprocity A, the non-zero median pumping speed v c also shows that the spherical cap model cannot completely describe the droplet's behavior for small f . However, we expect that v c eventually tends to zero since for sufficiently slow motion of the contact line, the droplet will go through a sequence of spherical-cap shapes. Constraints in the simulation time make it impossible to reach the limit f → 0. However, we observe in Fig. 6 and the inset that for f < 10 −2 τ −1 the distance between the two halves of the closed trajectories decreases with decreasing f such that the area A and therefore also v c should ultimately vanish. Note, an analogous behavior was observed, for example, in simulations of the one-armed microswimmer [53]. When its flexible flagellum beats quickly, it bends due to frictional forces and thereby moves nonreciprocally. But when it beats very slowly, it behaves like a rigid rod and thus The equilibrium contact angles varies between θ min eq = 60 • and θ max eq = 120 • .
moves reciprocally and the microswimmer cannot swim forward.
Finally, in Fig. 7(c) we relate the median displacement 〈d〉 directly to the non-reciprocity A and identify a non-linear relation. We checked that the scaling from Sec. 3.1, derived from the spherical cap model, does not produce a master curve for either A or 〈|d|〉 . This corroborates further that the spherical cap model is inapplicable for these quantities.

Oscillations of nonuniform wettability profiles
We now extend our investigation to oscillating nonuniform patterns of wettability. Specifically, we choose a pattern where the wettability varies periodically along the contact line of the droplet. As an example, we take the six-fold pattern illustrated in Fig. 8, where the equilibrium contact angle is modulated in space and time: θ eq (r, t) = θ max eq − (θ max eq − θ min eq ) · sin 2 (π f t) · with r = (r cos φ, r sin φ, 0) T and n = 6. The droplet initially sits with its base area centered at r = 0. In an experiment, an equivalent light intensity pattern can be realized with Laguerre-Gauss laser modes [54].
After an initial transient behavior, the droplet settles into steady oscillations analogous to the initial behavior in Sect. 3.1. However, in this case the droplet shape is not axisymmetric but instead it follows the six-fold symmetry of the wettability pattern. A snapshot of the droplet shape is presented in Fig. 8 and the whole dynamics can be seen in a video in Supplementary Material M04.
To study the steady oscillations quantitatively, we introduce the 6th harmonic mode a 6 of the dynamic contact angle along the contact line,  17) and (18), respectively. with k 6 = 6 · 2π/L, the instantaneous length L of the contact line, and the dynamic contact angle θ dyn (s) parameterized by the arc length s of the contact line. In the same way, we calculate the 6th harmonic modeã 6 of the equilibrium contact angle θ eq (s). Both, a 6 and a 6 , are periodic functions in time t when the droplet has reached steady oscillation and so the ratio of their Fourier transforms in time gives again a susceptibility χ( f ). In Fig. 9 we display the absolute value |χ( f )| and the phase shift ∆ϕ given by χ = |χ|e ∆ϕ . The principal behavior is similar to that of an oscillating uniform pattern shown in Fig. 3. However, now the linear model derived from the spherical cap model does not provide a master curve anymore since the droplet shape deviates strongly from the spherical cap. In detail, we observe while |χ| is shifted toward larger f when compared to the linear model (dashed line), ∆ϕ is shifted toward smaller f .
As before in Sect. 3.2, we study the periodic deformation of the droplet by combining two aspects of its shape, a 6 and the droplet height h. For steady oscillations we display closed-loop trajectories for several frequencies f in Fig. 10 and observe that their extent in the h-a 6 plane increases in all directions as f decreases. This is unlike the enclosed area for uniform wettability in Fig. 6, which decreased for f < 10 −3 τ −1 .
As mentioned before in Sect. 3.2, an increase in the enclosed area in configuration space means that the shape dynamics deviates more strongly from a reciprocal motion. While for the oscillating uniform pattern we have the sperical cap model as a reference, which shows reciprocal dynamics and which the droplet should approach for sufficiently small f , such a reference is missing for the nonuniform patterns and the droplet dynamics need not be reciprocal in the limit of small f .

Conclusions
We have studied liquid droplets on substrates with oscillating wettability focussing on their shape and internal fluid flow. When starting the wettability oscillations, the droplets go through a transient period, where the mean contact angle relaxes, and then settles into steady oscillations w.r.t. shape and contact angle. The amplitude of the contact-angle oscillations decreases with increasing frequency and, of course, increases with the amplitude of the wettability oscillations. At small frequencies the amplitude and phase shift of the oscillations follow the linearized and parameter-free spherical-cap model, while they deviate from it for larger frequencies, where the free droplet surface deviates noticeably from a spherical cap. As a result, amplitude and phase shift fall onto a master curve for small frequencies when we scale frequency with the decay time of the linearized spherical-cap model. Upon further analysis of the droplet shape and the interior flow field, we have found that also for small frequencies the spherical cap model cannot account for the droplet's behavior completely. The droplet shape deforms non-reciprocally which becomes evident when tracking its contact angle and height over time. The non-reciprocal dynamics of the droplet shape gives rise to a time-dependent internal flow field which displaces point-like tracer particles over the course of a full period of oscillation rather than returning them to their original position. Notably, the volume median of the displacement per period, the pumping speed, is constant w.r.t. frequency in the investigated frequency range and only depends on material properties of the droplet and substrate. Most importantly, the circulatory pumping of fluid inside the droplet is captured only by solving the full equations of Stokes flow, which we did using the boundary element method.
Repeating the same analysis for a droplet which is controlled by an oscillating nonuniform pattern of wet-tability with six-fold symmetry, we observed a very similar behavior. The droplet settles into steady oscillations, the amplitude of which decreases when the oscillation frequency increases. Furthermore, its shape deforms non-reciprocally, which becomes more pronounced at smaller oscillation frequencies.
Because the circulation is stimulated from outside by oscillations in wettability, it provides a mechanism for controlling the transport and possibly mixing of solutes inside the droplet. Internal transport through external stimuli has been used, for example, to precisely deposit a solute during the evaporation of droplets with lightresponsive surfactants [55]. Note, however, that depending on the contact angle evaporation contributes to and disturbs the flow field close to a moving contact line [56]. By designing specfic spatio-temporal wetttability patterns, such as moving steps in wettability [24], we envisage to similarly control the precise transport of solutes through substrates with switchable wettability.

Conflicts of Interest
There are no conflicts to declare.

A ESI
We provide four example movies. They show a droplet on a substrate, the wettability of which oscillates with frequency f = 10 −3 . The first three movies show a droplet on a substrate with uniform wettability. Movie M01 corresponds to an oscillation with θ min eq = 60 • and θ max eq = 120 • , movie M02 corresponds to θ min eq = 45 • and θ max eq = 90 • , and movie M03 corresponds to θ min eq = 30 • and θ max eq = 60 • . The final movie M04 corresponds to the nonuniform wettability profile given by Eq. (23) with θ min eq = 60 • and θ max eq = 120 • .