# Droplets on substrates with oscillating wettability

Josua Grawitter and Holger Stark

[always]halhttp://hal.archives-ouvertes.fr/eprint HALhal
[always]semanticscholarhttps://api.semanticscholar.org/CorpusID:eprint
SEMANTICSCHOLARsemanticscholar

Droplets on substrates with oscillating wettability†

Josua Grawitter111Technische Universität Berlin, Institut für Theoretische
Physik, Straße des 17. Juni 135, 10623 Berlin, Germany.,222E-mail:
<EMAIL_ADDRESS>and Holger Stark111Technische Universität
Berlin, Institut für Theoretische Physik, Straße des 17. Juni 135, 10623
Berlin, Germany.,333E-mail<EMAIL_ADDRESS>July 30, 2021

††footnotetext: †Electronic Supplementary Information (ESI) available: Their
descriptions are collected in Appendix A.

Abstract
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.

##  1 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 checker-board 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
linear-response 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.

##  2 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 $\bm{s}$ on the interfaces is
governed by the fluid velocity field at this point:

$\bm{\dot{s}}=\bm{v}(\bm{s})\leavevmode\nobreak\ .$ (1)

In the following we summarize how we determine the velocity field both at the
droplet surface and its interior.

###  2.1 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 $\bm{v}$ of Stokes flow are

$\displaystyle\mu\nabla^{2}\bm{v}$ $\displaystyle=\nabla p,$
$\displaystyle\nabla\cdot\bm{v}$ $\displaystyle=0\leavevmode\nobreak\ ,$ (2)

where $\mu$ 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
$\mathsf{O}$ and the associated stress field $\mathsf{T}$:

$c(\bm{r})\bm{v}(\bm{r})=\oint\limits_{\partial
D}\mathsf{O}(\bm{r-r^{\prime}})\bm{\sigma
n}(r^{\prime})\,\mathrm{d}^{2}\bm{r^{\prime}}\\\ -\oint\limits_{\partial
D}\bm{v}(\bm{r^{\prime}})\cdot\mathsf{T}(\bm{r-r^{\prime}})\bm{n}(\bm{r^{\prime}})\,\mathrm{d}^{2}\bm{r^{\prime}}$
(3)

with

$c(\bm{r})=\begin{cases}1&\text{for }\bm{r}\in D\setminus\partial D\\\
\frac{1}{2}&\text{for }\bm{r}\in\partial D\text{, where $\partial D$ is
smooth}\\\ \frac{\alpha}{4\pi}&\text{for }\bm{r}\in\partial D\text{, where
$\partial D$ has a corner}\\\ &\text{with inward {solid} angle
$\alpha$.}\end{cases}$ (4)

where $\partial D$ is the (time-dependent) surface of the droplet. Because the
equation relates velocity $\bm{v}$ and stress $\bm{\sigma}\bm{n}$, for any
surface point $\bm{r}$ either variable must be prescribed by boundary
conditions.

###  2.2 Boundary conditions

On the liquid-gas interface the normal stress balances surface tension forces
due to mean curvature $\kappa$ of the interface, i.e., $\bm{\sigma
n}=\gamma\kappa\bm{n}$, where $\gamma$ 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 $\bm{e}_{z}$ and
therefore $v_{z}=0$. Secondly, roughness of the substrate introduces a small
amount of slip with slip length $\lambda$, which we account for by setting
$\lambda\bm{\sigma n}=\mu\bm{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_{\mathrm{contact}}=\frac{\gamma}{9\mu\ln(h/\lambda)}(\theta_{\mathrm{dyn}}^{3}-\theta_{\mathrm{eq}}^{3})\,,$
(5)

where $\theta_{\mathrm{eq}}$ is the equilibrium contact angle, which defines
the wettability of the substrate, and $\theta_{\mathrm{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 $\mu$. For our simulations we choose dimensionless parameters. They
correspond, for example, to a droplet with an initial radius
$R_{0}=100\,\mu\mathrm{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
$\mu=209\,\mathrm{mPa}\\!\cdot\\!\mathrm{s}$, kinematic viscosity
$\nu=\mu/\rho=169\,\mathrm{mm}^{2}\\!\cdot\\!\mathrm{s}^{-1}$,
$\gamma=65.3\,\mathrm{mN}\\!\cdot\\!\mathrm{m}^{-1}$,
$\lambda=1\,\mathrm{nm}$, and $\ln(h/\lambda)=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 $\theta_{\mathrm{dyn}}$ to be the time
average of $\theta_{\mathrm{eq}}$. We calculate and report our data in units
of $R_{0}$ for length, $\tau=R_{0}^{2}/\nu$ for time, and $F_{0}=\nu\mu$ for
force. Thus the remaining dimensionless parameters are $\tilde{\gamma}=\gamma
R_{0}/F_{0}=0.19$ and $\tilde{\lambda}=\lambda/R_{0}=10^{-5}$.

###  2.3 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].

##  3 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 $\theta_{\mathrm{eq}}(t)$ oscillates with a
frequency $f$ between a minimum ($\theta_{\mathrm{eq}}^{\mathrm{min}}$) and
maximum ($\theta_{\mathrm{eq}}^{\mathrm{max}}$) value:

$\theta_{\mathrm{eq}}(t)=\theta_{\mathrm{eq}}^{\mathrm{min}}+(\theta_{\mathrm{eq}}^{\mathrm{max}}-\theta_{\mathrm{eq}}^{\mathrm{min}})\cdot\sin^{2}\left(\pi
ft\right)$ (6)

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 non-reciprocal so that
it effectively pumps fluid during one oscillation cycle.

###  3.1 Phenomenology

Figure 1:  Contact angle oscillations in response to oscillations in
wettability for three frequencies $f$ with
$\theta_{\mathrm{eq}}^{\mathrm{max}}=120^{\circ}$ and
$\theta_{\mathrm{eq}}^{\mathrm{min}}=60^{\circ}$.

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 $\theta_{\mathrm{eq}}^{\mathrm{max}}=120^{\circ}$ and
$\theta_{\mathrm{eq}}^{\mathrm{min}}=60^{\circ}$. Furthermore, we observe a
phase shift between the oscillating equilibrium and dynamic contact angles,
which also depends on $f$. Note, for different combinations of
$\theta_{\mathrm{eq}}^{\mathrm{max}}$ and
$\theta_{\mathrm{eq}}^{\mathrm{min}}$ we provide videos M01–M03 in the
Supplementary Material, where a relatively small frequency of
$f=10^{-3}\,\tau^{-1}$ is used (for details see Appendix A).

To quantify the response of the droplet to the oscillating substrate
wettability, we introduce the nonlinear susceptibility
$\chi=|\chi|\exp(\mathrm{i}\Delta\varphi)$ with absolute value $|\chi|$ and
phase shift $\Delta\varphi$. The imaginary unit is indicated by $\mathrm{i}$.
We extract $|\chi|$ from

$|\chi(f)|=\frac{\max_{t}\\{\theta_{\mathrm{dyn}}(t)\\}-\min_{t}\\{\theta_{\mathrm{dyn}}(t)\\}}{\theta_{\mathrm{eq}}^{\mathrm{max}}-\theta_{\mathrm{eq}}^{\mathrm{min}}}\leavevmode\nobreak\
.$ (7)

To calculate the phase shift $\Delta\varphi(f)$, we determine the first
Fourier coefficient $\alpha_{\mathrm{dyn}}$ of the dynamic contact angle
$\theta_{\mathrm{dyn}}$,

$\alpha_{\mathrm{dyn}}(f=T^{-1})=\lim\limits_{s\to\infty}\frac{1}{T}\int_{s}^{s+T}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\theta_{\mathrm{dyn}}(t)\,\mathrm{e}^{\mathrm{i}\,2\pi
t/T}\,\mathrm{d}t\,,$ (8)

where $s\to\infty$ insures that the oscillations of $\theta_{\mathrm{dyn}}(t)$
are steady. For the same time intervall we also determine the complex
amplitude $\alpha_{\mathrm{eq}}$ of the prescribed equilibrium contact angle
sand then calculate the phase shift between oscillating wettability and
dynamic contact angle from

$\Delta\varphi(f)=\arg[\alpha_{\mathrm{dyn}}(f)]-\arg[\alpha_{\mathrm{eq}}(f)]\,,$
(9)

where $\arg$ means the phase angle of the complex amplitude.

Figure 2:  Absolute susceptibility $|\chi|$ as a function of oscillation
frequency $f=\omega/2\pi$ in units of $\tau^{-1}$ for three combinations of
$\theta_{\mathrm{eq}}^{\mathrm{max}}$ and
$\theta_{\mathrm{eq}}^{\mathrm{min}}$. The inset shows the corresponding phase
shift $\Delta\varphi$.

In Fig. 2 we plot the absolute value $|\chi|$ over frequency for different
combinations of $\theta_{\mathrm{eq}}^{\mathrm{max}}$ and
$\theta_{\mathrm{eq}}^{\mathrm{min}}$. All three curves show the expected
decrease of $|\chi|$ with increasing $f$. Furthermore, for larger difference
$\theta_{\mathrm{eq}}^{\mathrm{max}}-\theta_{\mathrm{eq}}^{\mathrm{min}}$ 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
$\theta_{\mathrm{dyn}}$,

$\frac{\mathrm{d}\theta_{\mathrm{dyn}}}{\mathrm{d}t}=-g(\theta_{\mathrm{dyn}},\theta_{\mathrm{eq}})$
(10)

with

$g(\theta_{\mathrm{dyn}},\theta_{\mathrm{eq}})=\frac{\gamma_{\text{lg}}}{9\mu\ln(h/\lambda)}\sqrt[3]{\frac{\pi}{3V}}(\theta_{\mathrm{dyn}}^{3}-\theta_{\mathrm{eq}}^{3})\\\
\times\sqrt[3]{(1-\cos\theta_{\mathrm{dyn}})^{2}\;(2+\cos\theta_{\mathrm{dyn}})^{4}}$
(11)

To gain some inside and derive a characteristic relaxation time of the
spherical-cap model, we linearize it around
$\theta_{\mathrm{dyn}}=\theta_{\mathrm{eq}}$ in
$\Delta\theta=\theta_{\mathrm{dyn}}-\theta_{\mathrm{eq}}$:

$\frac{\mathrm{d}\theta_{\mathrm{dyn}}}{\mathrm{d}t}=-\left.\frac{\partial
g}{\partial\theta_{\mathrm{dyn}}}\right|_{\theta_{\mathrm{dyn}}=\theta_{\mathrm{eq}}}(\theta_{\mathrm{dyn}}-\theta_{\mathrm{eq}})\,.$
(12)

The derivative of $g$ is a characteristic relaxation rate $\tau_{0}^{-1}$.
But, in our case $\theta_{\mathrm{eq}}$ is a function of time. Nevertheless,
to have a constant rate, we calculate the derivative at the mean equilibrium
contact angle
$\bar{\theta}_{\mathrm{eq}}=(\theta_{\mathrm{eq}}^{\mathrm{max}}-\theta_{\mathrm{eq}}^{\mathrm{min}})/2$
and obtain

$\tau_{0}^{-1}(\bar{\theta}_{\mathrm{eq}})=\frac{\gamma_{\text{lg}}\bar{\theta}_{\mathrm{eq}}^{2}}{3\mu\ln(h/\lambda)}\sqrt[3]{\frac{\pi}{3V}(1-\cos\bar{\theta}_{\mathrm{eq}})^{2}\;(2+\cos\bar{\theta}_{\mathrm{eq}})^{4}}$
(13)

Now, we approximate Eq. (12) by using the constant $\tau_{0}$ instead of the
exact derivative of $g$. Rescaling time by $\tau_{0}$, we finally arrive at

$\frac{\mathrm{d}\theta_{\mathrm{dyn}}}{\mathrm{d}(t/\tau_{0})}\approx-[\theta_{\mathrm{dyn}}-\theta_{\mathrm{eq}}(t)]\,,$
(14)

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 $|\chi(\omega=2\pi f)|$.

Using time rescaled with $\tau_{0}$ also in the full spherical-cap model, we
can rewrite Eqs. (10) and (11) in a non-dimensionalized form as

$\frac{\mathrm{d}\theta_{\mathrm{dyn}}}{\mathrm{d}(t/\tau_{0})}=-\frac{1}{3\bar{\theta}_{\mathrm{eq}}^{2}}(\theta_{\mathrm{dyn}}^{3}-\theta_{\mathrm{eq}}(t)^{3})\\\
\times\sqrt[3]{\frac{(1-\cos\theta_{\mathrm{dyn}})^{2}\;(2+\cos\theta_{\mathrm{dyn}})^{4}}{(1-\cos\bar{\theta}_{\mathrm{eq}})^{2}\;(2+\cos\bar{\theta}_{\mathrm{eq}})^{4}}}$
(15)

Note that here the r.h.s. is independent of the liquid-gas surface tension
$\gamma_{\mathrm{lg}}$, the viscosity $\mu$ of the fluid, and the Cox-Voinov
parameter $\ln(h/\lambda)$, which determines the contact line mobility. All
these parameters are subsumed in the relaxation time $\tau_{0}$. The only
remaining parameters are $\theta_{\mathrm{eq}}^{\mathrm{min}}$ and
$\theta_{\mathrm{eq}}^{\mathrm{max}}$, which also determine
$\bar{\theta}_{\mathrm{eq}}$.

We now explore the linerarized model. By taking the Fourier transform of Eq.
(14), one finds the dynamical susceptibility $\chi(\omega)$, which quantifies
how $\theta_{\mathrm{dyn}}$ responds to an oscillation in
$\theta_{\mathrm{eq}}$:
$\hat{\theta}_{\mathrm{dyn}}(\omega)=\chi(\omega)\hat{\theta}_{\mathrm{eq}}(\omega)$
with

$\chi(\omega)=\frac{1}{1-\mathrm{i}\,\tau_{0}\omega}\,.$ (16)

The absolute value $|\chi|$ reads

$|\chi(\omega)|=\sqrt{\frac{1}{1+(\tau_{0}\omega)^{2}}}$ (17)

and the complex phase $\varphi$ is

$\varphi(\omega)=\arctan(\tau_{0}\omega)\,.$ (18)

Note that in the region where the linear model applies, $\varphi(\omega)$ is
identical to the phase shift $\Delta\varphi(f=\omega/2\pi)$ between
$\theta_{\mathrm{eq}}(t)$ and $\theta_{\mathrm{dyn}}(t)$ introduced above.

Figure 3:  Absolute susceptibility $|\chi|$ (a) and phase shift
$\Delta\varphi$ (b) as a function of oscillation frequency $f=\omega/2\pi$ in
units of $f_{0}=\tau_{0}^{-1}$ determined in the simulations for three
combinations of $\theta_{\mathrm{eq}}^{\mathrm{max}}$ and
$\theta_{\mathrm{eq}}^{\mathrm{min}}$ (coloured crosses with solid lines to
guide the eye, see legend). The dashed grey lines indicate the prediction of
the linear model in Eqs. (17) and (18), respectively. The dotted black line in
(a) shows a $f^{-0.5}$ scaling.

Because Eqs. (17) and (18) do not depend on any material parameters besides
the characteristic relaxation time $\tau_{0}$, they are candidates for the
master curves for our simulation results in Fig. 2. We simply need to rescale
frequency by $\tau_{0}^{-1}$ calculated for the specific values of
$\theta_{\mathrm{eq}}^{\mathrm{min}}$ and
$\theta_{\mathrm{eq}}^{\mathrm{max}}$. 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<\tau_{0}^{-1}$, which perfectly matches the linear model in this range. For
$f>\tau_{0}^{-1}$ 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\tau_{0}>3\cdot 10^{-1}$. The phase
shift varies between roughly $0.3\pi$ and $0.4\pi$ rather than approaching
$0.5\pi$ as predicted analytically by the linear model in Eq. (18) and
indicated as dashed line. The phase shifts beyond $f\tau_{0}>3\cdot 10^{-1}$
apparently depend on $\theta_{\mathrm{eq}}^{\mathrm{max}}$ and
$\theta_{\mathrm{eq}}^{\mathrm{min}}$.

To interpret these observations, we distinguish two regimes: a low frequency
regime with $f\tau_{0}<1$ and a high frequency regime with $f\tau_{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 $\mathrm{Im}\chi$, which in our linear
model from Eq. (16) reads

$\mathrm{Im}\chi(\omega)=\frac{\tau_{0}\omega}{1+(\tau_{0}\omega)^{2}}\,,$
(19)

quantifies dissipation. Because $\mathrm{Im}\chi$ as well as the phase shift
$\Delta\phi$ from Eq. (18) are linear in $\omega$ 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 $|\chi|$ means $\theta_{\mathrm{dyn}}$ barely tracks
$\theta_{\mathrm{eq}}(t)$. But also the full spherical cap model is unable to
predict the observed behavior for $|\chi|$ and $\Delta\varphi$ 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 $\theta_{\mathrm{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 $\theta_{\mathrm{dyn}}$ and increase in a short time the surface energy
of the droplet relative to the equilibrium reference shape for a given
$\theta_{\mathrm{eq}}$. Thus the work performed on the droplet is not
completely dissipated. Therefore, the phase shift angle $\Delta\varphi$ is
below the value $\pi/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.

###  3.2 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.

Figure 4: Displacements $\bm{d}(\bm{r})$ of tracers after one full period in a
cross section of a droplet on a substrate with oscillating wettability with
$\theta_{\mathrm{eq}}^{\mathrm{min}}=60^{\circ}$,
$\theta_{\mathrm{eq}}^{\mathrm{max}}=120^{\circ}$, and $f=10^{-3}\,\tau^{-1}$.

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 $\bm{r}_{0}$ we can
assign a displacement

$\bm{d}(\bm{r}_{0})=\bm{r}(T)-\bm{r}_{0}\,,$ (20)

where $\bm{r}(t)$ is the solution of the differential equation
$\dot{\bm{r}}(t)=\bm{v}(\bm{r},t)$ with initial condition
$\bm{r}(0)=\bm{r}_{0}$ and $\bm{v}(\bm{r},t)$ is the interior velocity field
of the droplet when it is steadily oscillating. In Fig. 4 we display an
example of $\bm{d}(\bm{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.

Figure 5: Droplet shape changes for $f=10^{-3}\,\tau^{-1}$,
$\theta_{\mathrm{eq}}^{\mathrm{max}}=60^{\circ}$ and
$\theta_{\mathrm{eq}}^{\mathrm{min}}=30^{\circ}$. (a) Wettability of the
substrate as a function of time characterized by the equilibrium contact angle
$\theta_{\mathrm{eq}}$. (b) Deformation of the droplet quantified by the
height $h(t)$ and rate of deformation quantified by $\dot{h}(t)$ plotted
versus time. Note that $\dot{h}(t)$ indicates unequal upward and downward
slopes for the height dynamics.

To quantify changes in the droplet shape, one can use, _e.g._ , the dynamic
contact angle $\theta_{\mathrm{dyn}}$ or the droplet height $h$. In Fig. 5 we
plot two oscillation cycles of $\theta_{\text{eq}}(t)$ in graph (a) together
with the corresponding $h(t)$ and its time derivative $\dot{h}(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 $\theta_{\text{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._ ,
$\theta_{\mathrm{eq}}=\theta_{\mathrm{dyn}}$ and a small change in wettability
occurs, _i.e._
$\theta_{\mathrm{eq}}\to\theta_{\mathrm{eq}}+\delta\theta_{\mathrm{eq}}$, we
calculate the contact line mobility $m$ by linearizing the Cox-Voinov law, Eq.
(5). The velocity of the contact line becomes

$v_{\mathrm{contact}}\approx-m\delta\theta_{\mathrm{eq}}\quad\text{with}\quad
m=\frac{\gamma\theta_{\mathrm{eq}}^{2}}{3\mu\ln(h/\lambda)}\,$ (21)

Apparently, the mobility $m$ can be smaller or larger for droplets with the
same $\theta_{\mathrm{dyn}}$, depending on $\theta_{\mathrm{eq}}$. If the
wettability is increasing (smaller $\theta_{\mathrm{eq}}$) $m$ is decreased.
Respectively, if the wettability is decreasing (larger $\theta_{\mathrm{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 $\theta_{\text{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 $\theta_{\mathrm{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.aaaNote that due to the discretization of the droplet surface, all
$\theta_{\mathrm{dyn}}$ of the simulated curves are systematically shifted
downwards by a small angle. Otherwise, they should center on the dashed line.
Note the non-zero area also means that contact angle and height oscillate out-
of-phase 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.

Figure 6:  Closed trajectories in configuration space projected onto the
$h$-$\theta_{\mathrm{dyn}}$ plane for various frequencies $f$ (see legend) for
$\theta_{\mathrm{eq}}^{\mathrm{min}}=60^{\circ}$ and
$\theta_{\mathrm{eq}}^{\mathrm{max}}=120^{\circ}$. The limiting behavior of
the spherical-cap model with $f\tau_{0}\ll 1$ is indicated by the black dashed
line.

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 high-
dimensional. However, projecting this trajectory into the two-dimensional
space spanned by droplet height $h$ and contact angle $\theta_{\mathrm{dyn}}$,
we can immediately calculate the projected area as

$A=\lim\limits_{s\to\infty}\oint_{s}^{s+T}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!h(t)\,\dot{\theta}_{\mathrm{dyn}}(t)\,\mathrm{d}t$
(22)

where the limit $s\to\infty$ ensures that $h(t)$ and
$\theta_{\mathrm{dyn}}(t)$ perform steady oscillations. In the following, we
call the parameter $A$ shape non-reciprocity.

Figure 7:  Shape non-reciprocity $A$ (a) and pumping velocity
$v_{\mathrm{c}}=f\langle d\rangle_{\mathcal{D}}$ (b) as a function of
oscillation frequency $f$ determined for the three combinations of
$\theta_{\mathrm{eq}}^{\mathrm{min}}$ and
$\theta_{\mathrm{eq}}^{\mathrm{max}}$ (colored points with dashed lines to
guide the eye, see legend). (c) Displacement $\langle d\rangle_{\mathcal{D}}$
plotted directly against $A$ for the same combinations, with an arrow
indicating the direction of increasing $f$.

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
$\theta_{\mathrm{eq}}^{\mathrm{min}}=60^{\circ}$ and
$\theta_{\mathrm{eq}}^{\mathrm{max}}=120^{\circ}$ 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 the spherical cap model, there
is still a significant difference between our BEM dynamics and the spherical
cap model.

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 $\theta_{\mathrm{dyn}}$ closer to
$\theta_{\mathrm{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
$\theta_{\mathrm{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\to 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\to\infty$ 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 $\bm{d}(\bm{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\langle d\rangle_{\mathcal{D}}$ w.r.t. the interior of the droplet,
_i.e._ , the 50th percentile value of the spatial distribution of
$d=|\bm{d}(\bm{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_{\mathrm{c}}$ also shows
that the spherical cap model cannot completely describe the droplet’s behavior
for small $f$. However, we expect that $v_{\mathrm{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\rightarrow 0$. However, we
observe in Fig. 6 and the inset that for $f<10^{-2}\,\tau^{-1}$ the distance
between the two halves of the closed trajectories decreases with decreasing
$f$ such that the area $A$ and therefore also $v_{\mathrm{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 moves reciprocally and the microswimmer cannot swim forward.

Finally, in Fig. 7(c) we relate the median displacement $\langle
d\rangle_{\mathcal{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
$\langle|\bm{d}|\rangle_{\mathcal{D}}$. This corroborates further that the
spherical cap model is inapplicable for these quantities.

###  3.3 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:

$\theta_{\mathrm{eq}}(\bm{r},t)=\theta_{\mathrm{eq}}^{\mathrm{max}}-(\theta_{\mathrm{eq}}^{\mathrm{max}}-\theta_{\mathrm{eq}}^{\mathrm{min}})\cdot\sin^{2}\left(\pi
ft\right)\cdot\\\
\sin^{2}\left(\frac{n}{2}\phi\right)\left[1-\exp{\left(-10\frac{|\bm{r}|^{2}}{R_{0}^{2}}\right)}\right]$
(23)

with $\bm{r}=(r\cos\phi,r\sin\phi,0)^{\mathrm{T}}$ and $n=6$. The droplet
initially sits with its base area centered at $\bm{r}=0$. In an experiment, an
equivalent light intensity pattern can be realized with Laguerre-Gauss laser
modes [54].

Figure 8: Snapshot of a droplet on a substrate with a nonuniform wettability
pattern with six-fold symmetry. The equilibrium contact angles varies between
$\theta_{\mathrm{eq}}^{\mathrm{min}}=60^{\circ}$ and
$\theta_{\mathrm{eq}}^{\mathrm{max}}=120^{\circ}$.

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,

$a_{6}=\frac{1}{L}\int_{0}^{L}\theta_{\mathrm{dyn}}(s)\cos(\mathrm{i}\,k_{6}s)\,\mathrm{d}s$
(24)

with $k_{6}=6\cdot 2\pi/L$, the instantaneous length $L$ of the contact line,
and the dynamic contact angle $\theta_{\mathrm{dyn}}(s)$ parameterized by the
arc length $s$ of the contact line. In the same way, we calculate the 6th
harmonic mode $\tilde{a}_{6}$ of the equilibrium contact angle
$\theta_{\mathrm{eq}}(s)$. Both, $a_{6}$ and $\tilde{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
$\chi(f)$. In Fig. 9 we display the absolute value $|\chi(f)|$ and the phase
shift $\Delta\varphi$ given by $\chi=|\chi|\mathrm{e}^{\Delta\varphi}$. 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 $|\chi|$ is
shifted toward larger $f$ when compared to the linear model (dashed line),
$\Delta\varphi$ is shifted toward smaller $f$.

Figure 9:  Absolute susceptibility $|\chi|$ and phase shift $\Delta\varphi$
(inset) as a function of oscillation frequency $f$ in units of
$f_{0}=\tau_{0}^{-1}$ determined in the simulations for two combinations of
$\theta_{\mathrm{eq}}^{\mathrm{max}}$ and
$\theta_{\mathrm{eq}}^{\mathrm{min}}$ (coloured crosses and dots with solid
lines to guide the eye, see legend). The dashed grey lines indicate the
prediction of the linear model in Eqs. (17) and (18), respectively.

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}\,\tau^{-1}$.

Figure 10: Closed trajectories in configuration space projected onto the
$h$-$a_{6}$ plane for various frequencies $f$ (see legend) for
$\theta_{\mathrm{eq}}^{\mathrm{min}}=60^{\circ}$ and
$\theta_{\mathrm{eq}}^{\mathrm{max}}=120^{\circ}$.

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$.

##  4 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 wettability 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 light-responsive 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.

##  Appendix 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
$\theta_{\mathrm{eq}}^{\mathrm{min}}=60^{\circ}$ and
$\theta_{\mathrm{eq}}^{\mathrm{max}}=120^{\circ}$, movie M02 corresponds to
$\theta_{\mathrm{eq}}^{\mathrm{min}}=45^{\circ}$ and
$\theta_{\mathrm{eq}}^{\mathrm{max}}=90^{\circ}$, and movie M03 corresponds to
$\theta_{\mathrm{eq}}^{\mathrm{min}}=30^{\circ}$ and
$\theta_{\mathrm{eq}}^{\mathrm{max}}=60^{\circ}$. The final movie M04
corresponds to the nonuniform wettability profile given by Eq. (23) with
$\theta_{\mathrm{eq}}^{\mathrm{min}}=60^{\circ}$ and
$\theta_{\mathrm{eq}}^{\mathrm{max}}=120^{\circ}$.

## Acknowledgements

We acknowledge financial support from DFG (German Research Foundation) via
Collaborative Research Center 910.

fielddoidefault0openaccess fielddoidefault0openaccess
fielddoidefault0openaccess fieldurldefault0openaccess
fieldurldefault0openaccess fieldurldefault0openaccess
fielddoidefault0openaccess fielddoidefault0openaccess
fieldurldefault0openaccess fieldurldefault0openaccess
fielddoidefault0openaccess fielddoidefault0openaccess
fieldurldefault0openaccess

## References

  * [1] Daniel Bonn et al.  “Wetting and spreading”  In _Rev. Mod. Phys._ 81, 2009, pp. 739  DOI: 10.1103/RevModPhys.81.739
  * [2] Kunihiro Ichimura, Sangkeun Oh and Masaru Nakagawa  “Light-Driven Motion of Liquids on a Photoresponsive Surface”  In _Science_ 288, 2000, pp. 1624  DOI: 10.1126/science.288.5471.1624
  * [3] Takahiro Seki  “A Wide Array of Photoinduced Motions in Molecular and Macromolecular Assemblies at Interfaces”  In _Bull. Chem. Soc. Jpn._ 91.7, 2018, pp. 1026–1057  DOI: 10.1246/bcsj.20180076
  * [4] Robert Malinowski, Ivan P. Parkin and Giorgio Volpe  “Advances towards programmable droplet transport on solid surfaces and its applications”  In _Chem. Soc. Rev._ 49, 2020, pp. 7879  DOI: 10.1039/D0CS00268B
  * [5] N. Delorme, J.-F. Bardeau, A. Bulou and F. Poncin-Epaillard  “Azobenzene-Containing Monolayer with Photoswitchable Wettability”  In _Langmuir_ 21, 2005, pp. 12278  DOI: 10.1021/la051517x
  * [6] Bruce C. Bunker  “Reversible switching of interfacial interactions”  In _Mater. Sci. Eng. R_ 62, 2008, pp. 157  DOI: 10.1016/j.mser.2008.06.001
  * [7] Joerg Lahann et al.  “A Reversibly Switching Surface”  In _Science_ 299, 2003, pp. 371  DOI: 10.1126/science.1078933
  * [8] Wuhui Jiang et al.  “Photo-switched wettability on an electrostatic self-assembly azobenzene monolayer”  In _Chem. Commun._ , 2005, pp. 3550  DOI: 10.1039/b504479k
  * [9] Federica Pirani et al.  “Light-Driven Reversible Shaping of Individual Azopolymeric Micro-Pillars”  In _Sci. Rep._ 6, 2016, pp. 31702  DOI: 10.1038/srep31702
  * [10] Stefano Luigi Oscurato, Fabio Borbone, Pasqualino Maddalena and Antonio Ambrosio  “Light-Driven Wettability Tailoring of Azopolymer Surfaces with Reconfigured Three-Dimensional Posts”  In _ACS Appl. Mater. Interfaces_ 9.35, 2017, pp. 30133–30142  DOI: 10.1021/acsami.7b08025
  * [11] Ho Sun Lim et al.  “Photoreversibly Switchable Superhydrophobic Surface with Erasable and Rewritable Pattern”  In _J. Am. Chem. Soc._ 128, 2006, pp. 14458  DOI: 10.1021/ja0655901
  * [12] Ho Sun Lim et al.  “UV-Driven Reversible Switching of a Roselike Vanadium Oxide Film between Superhydrophobicity and Superhydrophilicity”  In _J. Am. Chem. Soc._ 129, 2007, pp. 4128  DOI: 10.1021/ja0692579
  * [13] Ehsan Samiei, Maryam Tabrizian and Mina Hoorfar  “A review of digital microfluidics as portable platforms for lab-on a-chip applications”  In _Lab Chip_ 16 The Royal Society of Chemistry, 2016, pp. 2376–2396  DOI: 10.1039/C6LC00387G
  * [14] Dileep Mampallil and Huseyin Burak Eral  “A review on suppression and utilization of the coffee-ring effect”  In _Advances in Colloid and Interface Science_ 252, 2018, pp. 38–54  DOI: 10.1016/j.cis.2017.12.008
  * [15] H.. Greenspan  “On the motion of a small viscous droplet that wets a surface”  In _J. Fluid. Mech._ 84, 1978, pp. 125  DOI: 10.1017/S0022112078000075
  * [16] F. Brochard  “Motions of Droplets on Solid Surfaces Induced by Chemical or Thermal Gradients”  In _Langmuir_ 5, 1989, pp. 432  DOI: 10.1021/la00086a025
  * [17] K.. Glasner  “A boundary integral formulation of quasi-steady fluid wetting”  In _J. Comp. Phys._ 207, 2005, pp. 529  DOI: 10.1016/j.jcp.2005.01.022
  * [18] A. Moosavi, M. Rauscher and S. Dietrich  “Size dependent motion of nanodroplets on chemical steps”  In _J. Chem. Phys._ 129.4, 2008, pp. 044706  DOI: 10.1063/1.2955860
  * [19] Siegfried Dietrich, M.. Popescu and M. Rauscher  “Wetting on structured substrates”  In _J. Phys.: Condens. Matter_ 17, 2005, pp. S577  DOI: 10.1088/0953-8984/17/9/017
  * [20] Ondřej Kašpar et al.  “Confinement of water droplets on rectangular micro/nano-arrayed surfaces”  In _Lab Chip_ 16, 2016, pp. 2487  DOI: 10.1039/c6lc00622a
  * [21] Nikos Savva, Danny Groves and Serafim Kalliadasis  “Droplet dynamics on chemically heterogeneous substrates”  In _J. Fluid Mech._ 859, 2019, pp. 321  DOI: 10.1017/jfm.2018.758
  * [22] S. Engelnkemper et al.  “Continuation for Thin Film Hydrodynamics and Related Scalar Problems”  In _Computational Modelling of Bifurcations and Instabilities in Fluid Dynamics_ Cham: Springer International Publishing, 2019, pp. 459–501  DOI: 10.1007/978-3-319-91494-7_13
  * [23] Uwe Thiele, Karin John and Markus Bär  “Dynamical Model for Chemically Driven Running Droplets”  In _Phys. Rev. Lett._ 93, 2004, pp. 027802  DOI: 10.1103/PhysRevLett.93.027802
  * [24] Josua Grawitter and Holger Stark  “Steering droplets on substrates using moving steps in wettability”  In _Soft Matter_ 17, 2021, pp. 2454  DOI: 10.1039/D0SM02082F
  * [25] Pierre-Gilles Gennes  “Wetting: statics and dynamics”  In _Rev. Mod. Phys._ 57, 1985, pp. 827  DOI: 10.1103/RevModPhys.57.827
  * [26] Chun Huh and L.. Scriven  “Hydrodynamic Model of Steady Movement of a Solid / Liquid / Fluid Contact Line”  In _J. Colloid Interface Sci._ 35, 1971, pp. 85  DOI: 10.1016/0021-9797(71)90188-3
  * [27] O.. Voinov  “Hydrodynamics of Wetting”  In _Fluid Dyn._ 11, 1976, pp. 714  DOI: 10.1007/BF01012963
  * [28] X. Noblin, A. Buguin and F. Brochard-Wyart  “Vibrations of sessile drops”  In _Eur. Phys. J. Spec. Top._ 166, 2009, pp. 7–10  DOI: 10.1140/epjst/e2009-00869-y
  * [29] Seiichi Sudo et al.  “The Dynamic Behavior of Liquid Droplets on Vibrating Plate”  In _J. Jpn. Soc. Exp. Mech._ 10, 2010, pp. s38–s45  DOI: 10.11395/jjsem.10.s38
  * [30] Tanu Singla et al.  “Dynamics of a vertically vibrating mercury drop”  In _AIP Adv._ 9.4, 2019, pp. 045204  DOI: 10.1063/1.5088043
  * [31] Gary F. Teletzke, H. Davis and L.E. Scriven  “Wetting hydrodynamics”  In _Rev. Phys. appl._ 23, 1988, pp. 989–1007  DOI: 10.1051/rphysap:01988002306098900
  * [32] Alexander Oron, Stephen H. Davis and S. Bankoff  “Long-scale evolution of thin liquid films”  In _Rev. Mod. Phys._ 69, 1997, pp. 931  DOI: 10.1103/RevModPhys.69.931
  * [33] M.. Popescu, G. Oshanin, S. Dietrich and A.-M. Cazabat  “Precursor films in wetting phenomena”  In _J. Phys.: Condens. Matter_ 24, 2012, pp. 243102  DOI: 10.1088/0953-8984/24/24/243102
  * [34] Uwe Thiele  “Recent advances in and future challenges for mesoscopic hydrodynamicmodelling of complex wetting”  In _Colloids Surf. A_ 553, 2018, pp. 487  DOI: 10.1016/j.colsurfa.2018.05.049
  * [35] M.. Ruijter et al.  “Contact Angle Relaxation during the Spreading of Partially Wetting Drops”  In _Langmuir_ 13, 1997, pp. 7293  DOI: 10.1021/la970825v
  * [36] Todd M. Squires and Stephen R. Quake  “Microfluidics: Fluid physics at the nanoliter scale”  In _Rev. Mod. Phys._ 77, 2005, pp. 977  DOI: 10.1103/RevModPhys.77.977
  * [37] Damien Baigl  “Photo-actuation of liquids for light-driven microfluidics: state of the art and perspectives”  In _Lab Chip_ 12, 2012, pp. 3637  DOI: 10.1039/c2lc40596b
  * [38] Qingwen Dai, Wei Huang, Xiaolei Wang and M.M. Khonsari  “Directional interfacial motion of liquids: Fundamentals, evaluations, and manipulation strategies”  In _Tribol. Int._ 154, 2021, pp. 106749  DOI: 10.1016/j.triboint.2020.106749
  * [39] Sangtae Kim and Sepp J. Karrila  “Microhydrodynamics”  Mineola/NY: Dover Publications, 2005
  * [40] C. Pozrikidis  “Boundary integral and singularity methods for linearized viscous flow”  Cambridge: Cambridge University Press, 1992
  * [41] Silvia Jiménez Bolaños and Bogdan Vernescu  “Derivation of the Navier slip and slip length for viscous flows over a rough boundary”  In _Phys. Fluids_ 29, 2017, pp. 057103  DOI: 10.1063/1.4982899
  * [42] R.. Cox  “The dynamics of the spreading of liquids on a solid surface”  In _J. Fluid Mech._ 168, 1986, pp. 169  DOI: 10.1017/S0022112086000332
  * [43] H.. Moffatt  “Viscous and resistive eddies near a sharp corner”  In _J. Fluid Mech._ 18, 1964, pp. 1–18  DOI: 10.1017/S0022112064000015
  * [44] H.. Eral, D…. ’t Mannetje and J.. Oh  “Contact angle hysteresis: a review of fundamentals and applications”  In _Colloid Polym. Sci._ 291, 2013, pp. 247  DOI: 10.1007/s00396-012-2796-6
  * [45] Jacco H. Snoeijer and Bruno Andreotti  “Moving Contact Lines: Scales, Regimes, and Dynamical ransitions”  In _Annu. Rev. Fluid Mech._ 45, 2013, pp. 269  DOI: 10.1146/annurev-fluid-011212-140734
  * [46] Ch. Tsitouras  “Runge–Kutta pairs of order 5(4) satisfying only the first column simplifying assumption”  In _Comput. Math. Appl._ 62, 2011, pp. 770  DOI: 10.1016/j.camwa.2011.06.002
  * [47] John T. Katsikadelis  “The Boundary Element Method for Engineers and Scientists”  Elsevier, 2016
  * [48] Joshua D. McGraw et al.  “Slip-mediated dewetting of polymer microdroplets”  In _Proc. Natl. Acad. Sci. U.S.A._ 113, 2016, pp. 1168  DOI: 10.1073/pnas.1513565113
  * [49] Tak Shing Chan et al.  “Morphological evolution of microscopic dewetting droplets with slip”  In _J. Fluid Mech._ 828, 2017, pp. 271  DOI: 10.1017/jfm.2017.515
  * [50] Yu.. Pityuk, O.. Abramova, N.. Gumerov and I.. Akhatov  “Boundary Element Modeling of Dynamics of a Bubble in Contact with a Solid Surface at Low Reynolds Numbers”  In _Math. Models Comput. Simul._ 10, 2018, pp. 209  DOI: 10.1134/S2070048218020102
  * [51] E.. Purcell  “Life at low Reynolds number”  In _Am. J. Phys._ 45, 1977, pp. 3  DOI: 10.1119/1.10903
  * [52] Thomas Young  “An essay on the cohesion of fluids”  In _Philos. Trans. R. Soc._ 95, 1805, pp. 65–87  DOI: 10.1098/rstl.1805.0005
  * [53] E. Gauger and Holger Stark  “Numerical study of an artificial microswimmer”  In _Phys. Rev. E_ 74, 2006, pp. 021907  DOI: 10.1103/PhysRevE.74.021907
  * [54] L. Allen, M.. Beijersbergen, R… Spreeuw and J.. Woerdman  “Orbital angular momentum of light and the transformation of Laguerre-Gaussian laser modes”  In _Phys. Rev. A_ 45, 1992, pp. 8185  DOI: 10.1103/PhysRevA.45.8185
  * [55] Subramanyan Namboodiri Varanakkottu et al.  “Light-Directed Particle Patterning by Evaporative Optical Marangoni Assembly”  In _Nano Lett._ 16, 2016, pp. 644  DOI: 10.1021/acs.nanolett.5b04377
  * [56] Hanneke Gelderblom, Oscar Bloemen and Jacco H. Snoeijer  “Stokes flow near the contact line of an evaporating drop”  In _J. Fluid Mech._ 709, 2012, pp. 69  DOI: 10.1017/jfm.2012.321