††thanks: ORCID ID: 0000-0001-7989-9271

# Gradient dynamics approach to reactive thin-film hydrodynamics

Florian Voss  Institute of Theoretical Physics, University of Münster,
Wilhelm-Klemm-Str. 9, 48149 Münster, Germany  Uwe Thiele  u.thiele@uni-
muenster.de http://www.uwethiele.de Institute of Theoretical Physics,
University of Münster, Wilhelm-Klemm-Str. 9, 48149 Münster, Germany  Center
for Nonlinear Science (CeNoS), University of Münster, Corrensstr. 2, 48149
Münster, Germany  Center for Multiscale Theory and Computation (CMTC),
University of Münster, Corrensstr. 40, 48149 Münster, Germany

###### Abstract

Wetting and dewetting dynamics of simple and complex liquids is described by
kinetic equations in gradient dynamics form that incorporates the various
coupled dissipative processes in a fully thermodynamically consistent manner.
After briefly reviewing this, we also review how chemical reactions can be
captured by a related gradient dynamics description, assuming detailed
balanced mass action type kinetics. Then, we bring both aspects together and
discuss mesoscopic reactive thin-film hydrodynamics illustrated by two
examples, namely, models for reactive wetting and reactive surfactants. These
models can describe the approach to equilibrium but may also be employed to
study out-of-equilibrium dynamics. In the latter case, one breaks the gradient
dynamics form by chemostatting to obtain active systems. In this way, for
reactive wetting we recover running drops that are driven by chemically
sustained wettability gradients and for drops covered by autocatalytic
reactive surfactants we find complex forms of self-propulsion and self-excited
oscillations.

reactive thin-film hydrodynamics, reactive wetting, reactive surfactants,
gradient dynamics, chemo-mechanical coupling, self-propelled drops

##  1 Introduction

Phenomena involving reactive fluids are found across many scales and range
from the formation of stars to combustion, biofilm growth and many
intracellular processes [1, 2, 3, 4, 5, 6, 7]. Often, occurring chemical
reactions do not only amend existing flows, but give rise to novel phenomena
that result from the interplay of reactions and hydrodynamic transport. In his
seminal paper [8], Turing showed that pattern formation is possible even if
only simple diffusive transport exists in a reactive system. One may then
naturally turn to more intricate hydrodynamic settings, for instance, the
wetting of solid or liquid substrates by (other) liquids [9, 10, 11] and
discuss how chemical reactions might influence the occurring phenomena [12,
13, 14, 15, 16, 17, 18, 19, 20, 21], also see reviews [22, 23, 24].

The statics and dynamics of wetting inherently involves the study of
interacting interfaces and three-phase contact lines which gives rise to a
myriad of fascinating phenomena like drop spreading and dewetting of thin
liquid films [10, 11], contact line motion and instabilities [11, 25], the
lotus effect [26] as well as tears of wine [27] and surfactant-induced
spreading instabilities [11]. One may then expect that the additional presence
of chemical reactions will greatly enrich the spectrum of interfacial effects.
Indeed, examples include chemically driven running droplets [16, 28, 29, 30,
31, 32, 33, 34, 35, 36, 37, 38, 39], and chemically driven Marangoni flows
that cause oscillatory bulk convective motion [40, 41, 42], surface waves [43,
44], fingering instabilities [45], and crawling vesicular aggregates [46].
Such reactive interphase phenomena are highly relevant for technical
applications like soldering [47, 22], but also have many biophysical
applications such as, e.g., intracellular processes related to the wetting
properties of membranes and biomolecular condensates [48, 49, 50].

It is tempting to base theoretical descriptions of such complex dynamic
scenarios, where chemical reactions couple nontrivially to interface-dominated
hydrodynamic transport processes, on a simple addition of mass action chemical
kinetics to existing hydrodynamic descriptions. However, there is a certain
danger that such an ad hoc approach may miss certain cross-influences of the
different aspects that one has added together. A preferential approach could
be to obtain the coupled chemo-hydrodynamic dynamics from a common
thermodynamic framework, namely, a gradient dynamics on an underlying energy
functional as well known from the description of chemically inert systems
described, e.g., by mesoscopic hydrodynamics (thin-film models) [51, 52, 53,
54, 55]. If chemical reactions were incorporated into such an approach,
systems would still be relaxational, i.e., would approach an equilibrium.
Then, specific nonequilibrium conditions could be included in a well-
controlled manner. In this way, couplings that originate in thermodynamic
constraints and that may significantly affect statics and dynamics are still
present, even when the system is permanently driven away from equilibrium, as
is often the case in biological contexts. Familiar examples of reciprocal
pairs of cross-couplings include the Seebeck and Peltier effects
(thermoelectricity) and the Soret and Dufour effects (thermophoresis and
diffusion thermoeffect) [56, 57, 58]. These effects are captured by linear
nonequilibrium thermodynamics and are pairwise related by the Onsager
reciprocity relations [59, 60] as a consequence of microscopic reversibility
(detailed balance). In general, chemical reactions lie outside of the
framework of linear thermodynamics [58]. However, by invoking detailed balance
one may still infer thermodynamically consistent kinetics which then for ideal
systems recover mass action kinetics [61, 62, 63].

Returning specifically to the description of reactive wetting and reactive
surfactants, we may then find chemical reactions that are not only dependent
on the reactant concentrations but also on, e.g., reactant-substrate
interactions, the proximity of the three-phase contact line or the shape of
the liquid-gas interface. In the present work, we review how chemical
reactions can be incorporated into a gradient dynamics model to obtain a
description of reactive thin-film hydrodynamics responsible, e.g., for
reactive wetting. This approach is based on an underlying (free) energy
functional from which thermodynamically fully consistent kinetic equations
describing reactive thin-film hydrodynamics are derived.

The structure of our work is as follows. In Section 2, we review the gradient
dynamics approach to linear nonequilibrium thermodynamics. In Section 3, we
turn to the gradient dynamics description of chemical reactions. We start with
a simple autocatalytic reaction in a spatially homogeneous system, clarify how
mass action kinetics may be re-expressed in terms of a free energy functional
and show that the principle of detailed balance guarantees thermodynamic
consistency. We also discuss how the breaking of detailed balance, e.g. by
chemostatting, can result in persistent nonequilibrium behavior. Section 4
then exemplifies the approach by considering two cases of reactive mesoscopic
thin-film hydrodynamics. Both involve the interplay of chemical reactions and
interfacial effects, e.g., reactant-induced wettability gradients or solutal
Marangoni stresses. Further, we show how breaking the gradient dynamics form
by external chemostats results in active behavior such as the self-propulsion
of drops and the self-excitation of various drop oscillations. We conclude in
Section 5 with a summary and discussion of possible extensions.

##  2 Linear gradient dynamics and applications

Here, we briefly review the gradient dynamics description of systems that are
sufficiently close to thermodynamic equilibrium to allow for a description by
linear nonequilibrium thermodynamics, i.e., the kinetic equations are linear
in the variations of an underlying thermodynamic functional (the thermodynamic
forces). In Section 2.1, we introduce the general kinetic equations for the
dynamics of a $N$ scalar fields and show their thermodynamic consistency. In
Section 2.2 we give two corresponding examples related to reactive thin-film
hydrodynamics.

###  2.1 General Equations

We now consider a spatially extended system, the state of which is fully
determined by $N$ scalar state variables
$\mathbf{u}=\left(u_{1},u_{2},\ldots,u_{N}\right)^{\text{T}}$ which depend on
the spatial coordinates $\vec{x}=(x,y,z)^{\text{T}}$ and time $t$. We assume
that the system relaxes to thermodynamic equilibrium. Depending on the
relevant constraints, different thermodynamic functionals capture this
relaxation process in agreement with the second law of thermodynamics. Here,
we restrict ourselves to closed, isothermal systems. In this case, the
(Helmholtz) free energy functional $F\left[\mathbf{u}\right]$ decreases
monotonically until equilibrium is reached. If the system is sufficiently
close to thermodynamic equilibrium, the general time evolution equations for
the fields $\mathbf{u}$ can be expressed as

$\partial_{t}u_{\alpha}=\nabla\cdot\left[\sum_{\beta=1}^{N}Q^{\mathrm{c}}_{\alpha\beta}\nabla\frac{\delta
F}{\delta
u_{\beta}}\right]-\sum_{\beta=1}^{N}Q^{\mathrm{nc}}_{\alpha\beta}\frac{\delta
F}{\delta u_{\beta}},$ (1)

where $\partial_{t}$ denotes the partial derivative with respect to time $t$
and $\nabla=\left(\partial_{x},\partial_{y},\partial_{z}\right)^{\text{T}}$ is
the spatial gradient operator. The variational derivatives $\delta F/\delta
u_{\beta}$ express thermodynamically conjugate quantities to the fields
$u_{\beta}$ such as the pressure or chemical potentials. In (1), these
variations linearly enter the time evolution equations. Here, the first term
corresponds to transport processes that conserve $\int
u_{\alpha}\text{d}^{3}r$ and that are driven by spatial gradients in pressure
or chemical potentials. The second term is associated with nonconserved
contributions to the dynamics that are related to transitions between the
individual fields and are, for instance, driven by local differences in
chemical potentials. The mobility matrices
$\mathbf{\underline{Q}}^{\mathrm{c}}=\left(Q^{\mathrm{c}}_{\alpha\beta}\right)$
and
$\mathbf{\underline{Q}}^{\mathrm{nc}}=\left(Q^{\mathrm{nc}}_{\alpha\beta}\right)$
relate the energy variations to the fluxes and are positive (semi-)definite
and symmetric, expressing irreversibility of the macroscopic processes and
microscopic reversibility (Onsager relations [59, 60]). We note that the
specific form (1) can be derived from Onsager’s variational principle [64,
54]. Using (1) and the properties of
$\mathbf{\underline{Q}}^{\mathrm{c}},\mathbf{\underline{Q}}^{\mathrm{nc}}$,
one can show that the free energy $F$ decreases monotonically:

$\displaystyle\frac{\text{d}F}{\text{d}t}$
$\displaystyle=\int\left[\sum_{\alpha=1}^{N}\frac{\delta F}{\delta
u_{\alpha}}\frac{\partial u_{\alpha}}{\partial t}\right]\text{d}^{3}r$ (2)
$\displaystyle=\int\left[\sum_{\alpha,\beta=1}^{N}\frac{\delta F}{\delta
n_{\alpha}}\nabla\cdot\left(Q^{\mathrm{c}}_{\alpha\beta}\nabla\frac{\delta
F}{\delta
n_{\beta}}\right)\right]\text{d}^{3}r-\int\left[\sum_{\alpha,\beta=1}^{N}\frac{\delta
F}{\delta n_{\alpha}}Q^{\mathrm{nc}}_{\alpha\beta}\frac{\delta F}{\delta
n_{\beta}}\right]\text{d}^{3}r$ (3)
$\displaystyle=-\int\left[\sum_{\alpha,\beta=1}^{N}Q^{\mathrm{c}}_{\alpha\beta}\left(\nabla\frac{\delta
F}{\delta n_{\alpha}}\right)\cdot\left(\nabla\frac{\delta F}{\delta
n_{\beta}}\right)\right]\text{d}^{3}r-\int\left[\sum_{\alpha,\beta=1}^{N}\frac{\delta
F}{\delta n_{\alpha}}Q^{\mathrm{nc}}_{\alpha\beta}\frac{\delta F}{\delta
n_{\beta}}\right]\text{d}^{3}r$ (4) $\displaystyle\leq 0.$ (5)

Here, integration is performed over the whole system domain. From (2) to (3)
we have used the dynamical equations (1) and from (3) to (4) partial
integration was used on the first term, assuming no-flux or periodic boundary
conditions. The final inequality follows from the positive (semi-)definiteness
of the mobility matrices. Then, if $F$ is bounded from below, it is a Lyapunov
functional to the dynamics (1) which is therefore consistent with (linear)
nonequilibrium thermodynamics.

###  2.2 Examples

To model particular systems, the described gradient dynamics model (1) can be
supplemented by specific choices for the energy functional and the mobility
matrices. In the following, we briefly discuss a few models. For instance, the
diffusion of an ideal gas is described by the diffusion equation

$\partial_{t}n=D\Delta n,$ (6)

where $\Delta=\partial_{x}^{2}+\partial_{y}^{2}+\partial_{z}^{2}$ is the
Laplace operator, $n$ denotes the local particle density (particles per unit
volume) and $D>0$ is the diffusion constant. With the entropic free energy

$F=k_{b}T\int n\left[\ln(n/n_{0})-1\right]\text{d}^{3}r,$ (7)

where $k_{b}$ is the Boltzmann constant, $T$ is the temperature, $n_{0}$ is
some reference density, and with the diffusive mobility

$M(n)=\frac{Dn}{k_{b}T},$ (8)

the diffusion equation (6) can be brought into the gradient dynamics form

$\partial_{t}n=\nabla\cdot\left[M(n)\nabla\frac{\delta F}{\delta n}\right].$
(9)

Note that here $\delta F/\delta n$ corresponds to a chemical potential.

Mesoscopic thin-film equations [54] represent another class of such models. In
the simplest case, the evolution of the local film height $h$ of a thin liquid
film on a homogeneous rigid solid substrate is described by [65, 66, 10]

$\partial_{t}h=-\nabla\cdot\left[Q(h)\nabla\left(\gamma\Delta
h+\Pi(h)\right)\right]$ (10)

with $\nabla=(\partial_{x},\partial_{y})^{\text{T}}$ and
$\Delta=\partial_{x}^{2}+\partial_{y}^{2}$. Here, $\gamma$ is a constant
surface tension and the term $-\gamma\Delta h$ is the Laplace pressure, while
$\Pi(h)$ is a Derjaguin (disjoining) pressure that encodes wettability on the
mesoscopic scale [9, 67, 68]. The mobility factor $Q(h)\geq 0$ typically is a
polynomial in $h$ and depends on the particular model, e.g., the degree of
slip at the substrate [69, 65, 70]. We introduce the free energy functional

$F=\int\left[f(h)+\frac{\gamma}{2}|\nabla h|^{2}\right]\text{d}^{2}r,$ (11)

where the Derjaguin pressure is $\Pi(h)=-f^{\prime}(h)$ with the dash
indicating the univariate derivative. We can then write (10) as [51, 68]

$\partial_{t}h=\nabla\cdot\left[Q(h)\nabla\frac{\delta F}{\delta h}\right].$
(12)

From (12), we can see that $\delta F/\delta h$ corresponds to a pressure (or
its negative). Besides the two shown simple examples, other well-known
examples include the Cahn-Hilliard equation [71, 72], various two- and three-
field thin-film models, e.g., for two-layer films, films of binary mixtures,
surfactant-covered drops and drops on soft and adaptive substrates [73, 74,
75, 76, 52, 77, 78, 55] as well as other soft matter models [79, 64].

##  3 Chemical reactions as nonlinear gradient dynamics

After having treated the gradient dynamics form of systems that are captured
by linear nonequilibrium thermodynamics in Section 2, we now turn to systems
allowing for chemical reactions. As we will see, reactive processes that
follow mass action type kinetics lie outside of the framework of linear
nonequilibrium thermodynamics, i.e., the corresponding reactive fluxes between
chemical species are nonlinear in the variations of the underlying
thermodynamic functional.111We remark that with somewhat ‘unnatural’
mobilities, chemical reactions can be treated by means of linear
nonequilibrium thermodynamics [80, 81, 82] while the resulting fluxes are
still nonlinear in the variations. To motivate this, we first treat ideal
mixtures in Section 3.1 and show that the corresponding mass action type
kinetics can be rewritten employing rates that are nonlinear in the (free)
energy variations. We then show that the general dynamic equations indeed
represent a gradient dynamics when assuming detailed balanced kinetics. In
Section 3.2, we discuss possible deviations from ideal mass action type
kinetics and touch on some of the reactive cross-couplings that follow.

###  3.1 Mass action type kinetics as gradient dynamics

To motivate that chemical processes are insufficiently treated by linear
nonequilibrium thermodynamics, we consider a closed, isothermal box of volume
$V$ at temperature $T$ containing two well-mixed ideal gases, namely species
$N_{1}$ and $N_{2}$. For this setup, the (Helmholtz) free energy appropriately
describes the approach to thermodynamic equilibrium. To be specific, we assume
that the autocatalytic conversion reaction

$2N_{1}+N_{2}\rightleftharpoons 3N_{1}$ (13)

takes place between the two chemical species. For mixtures of ideal gases,
chemical reactions are appropriately described by mass action kinetics [58].
The time evolution of the particle densities $n_{1}$ and $n_{2}$ of the two
gases (particles per unit volume) inside the box is then given by

$\displaystyle\frac{\text{d}n_{1}}{\text{d}t}$
$\displaystyle=r_{\mathrm{f}}n_{1}^{2}n_{2}-r_{\mathrm{b}}n_{1}^{3},$ (14)
$\displaystyle\frac{\text{d}n_{2}}{\text{d}t}$
$\displaystyle=-r_{\mathrm{f}}n_{1}^{2}n_{2}+r_{\mathrm{b}}n_{1}^{3},$ (15)

where $r_{\mathrm{f}},r_{\mathrm{b}}>0$ denote the reaction rates in the
forward and backward direction. Thereby, the reaction $2N_{1}+N_{2}\rightarrow
3N_{1}$ is defined as the forward reaction. Note that
$\text{d}(n_{1}+n_{2})/\text{d}t=0$ and the total particle count
$N=(n_{1}+n_{2})V$ in the container is conserved. We consider the associated
chemical potentials of ideal gases [58]

$\mu_{\alpha}=k_{b}T\ln\left(n_{\alpha}/n_{0,\alpha}\right)+\zeta_{\alpha}(T),$
(16)

where $n_{0,\alpha}$ are reference densities and $\zeta_{\alpha}(T)$ are
general dependencies on the temperature $T$. Note that $\zeta_{\alpha}$
represents the chemical potential of species $i$ at the reference density
$n_{0,\alpha}$. For simplicity, we set $\zeta_{\alpha}=0$ in the following.
Using Eq. (16), we rewrite Eqs. (14) and (15) as

$\displaystyle\frac{\text{d}n_{1}}{\text{d}t}$
$\displaystyle=r_{\mathrm{f}}n^{2}_{0,1}n_{0,2}e^{\frac{2\mu_{1}+\mu_{2}}{k_{b}T}}-r_{\mathrm{b}}n^{3}_{0,1}e^{\frac{3\mu_{1}}{k_{b}T}},$
(17) $\displaystyle\frac{\text{d}n_{2}}{\text{d}t}$
$\displaystyle=-r_{\mathrm{f}}n^{2}_{0,1}n_{0,2}e^{\frac{2\mu_{1}+\mu_{2}}{k_{b}T}}+r_{\mathrm{b}}n^{3}_{0,1}e^{\frac{3\mu_{1}}{k_{b}T}}.$
(18)

We know that the gas mixture must approach thermodynamic equilibrium. For
chemical reactions thermodynamic equilibrium is equivalent to the condition
that the affinity $A$ vanishes [58]. The affinity corresponds to the
difference in total chemical energy between the reactants and the products of
a reaction and represents the thermodynamic force that drives the
reaction.222Note that the affinity $A_{j}$ can be related to the reaction
coordinate $\xi_{j}$ which represents the ‘degree of advancement’ of reaction
$j$ (see also [58]). In a well-mixed reactor it is related to the particle
number by $N_{\alpha}=N_{0,\alpha}+\sum_{j}(\nu^{b}_{\alpha j}-\nu^{f}_{\alpha
j})\xi_{j}$ where $N_{0,\alpha}$ are the initial particle counts or
equivalently $\partial\xi_{j}/\partial t=J_{j}$. Differentiating, e.g., the
free energy $F$ with respect to $\xi_{j}$ we obtain $\left(\partial
F/\partial\xi_{j}\right)_{T,V,\xi_{i\neq j}}=\sum_{\alpha}\left(\partial
F/\partial N_{\alpha}\right)_{T,V,N_{\beta\neq\alpha}}\partial
N_{\alpha}/\partial\xi_{j}=\sum_{\alpha}(\nu^{b}_{\alpha j}-\nu^{f}_{\alpha
j})\mu_{\alpha}=A_{j}$, that is, the affinity is the derivative of the
relevant thermodynamic potential with respect to the reaction coordinate [62].
Specifically for reaction (13) the affinity is given by
$A=3\mu_{1}-(2\mu_{1}+\mu_{2})=\mu_{1}-\mu_{2}$. Setting the right hand side
of (17)-(18) to zero for $A=0$, one finds that this condition only holds if
$r=r_{\mathrm{f}}n^{2}_{0,1}n_{0,2}=r_{\mathrm{b}}n^{3}_{0,1}$ is the rate
constant for both the forward and the backward reaction. Only then (17)-(18)
can be reformulated as

$\displaystyle\frac{\text{d}n_{1}}{\text{d}t}$
$\displaystyle=r\left(e^{\frac{2\mu_{1}+\mu_{2}}{k_{b}T}}-e^{\frac{3\mu_{1}}{k_{b}T}}\right),$
(19) $\displaystyle\frac{\text{d}n_{2}}{\text{d}t}$
$\displaystyle=-r\left(e^{\frac{2\mu_{1}+\mu_{2}}{k_{b}T}}-e^{\frac{3\mu_{1}}{k_{b}T}}\right),$
(20)

with a common reaction rate.333Note that we have restricted ourselves to
mixtures of ideal gases in the derivation of (19)-(20). Analogously, one could
consider a reacting ideal solution (or any other ideal system), for which mass
action kinetics are also commonly assumed. The chemical potential of, e.g.,
component $1$ in an ideal solution would read
$\mu_{1}=k_{b}T\ln(\frac{n_{1}}{n})+\zeta_{1}(T)$ with the total particle
density $n_{1}+n_{2}=n$ and where the first term corresponds to the mixing
entropy of $N_{1}$. If one starts with mass action kinetics formulated in
terms of particle densities as in (14)-(15), then the equilibrium constant
$r_{\mathrm{f}}/r_{\mathrm{b}}$ generally becomes a function of $T$ and $n$
[83, 58]. However, this does not change the following general results and the
corresponding mass action kinetics can still be derived from Eqs. (19)-(20) or
(24)-(25). Writing the reactive flux as

$J=J^{f}-J^{b}=r\left(e^{\frac{2\mu_{1}+\mu_{2}}{k_{b}T}}-e^{\frac{3\mu_{1}}{k_{b}T}}\right),$
(21)

at thermodynamic equilibrium we have equal fluxes in both reactive directions,
i.e.,

$J^{f}_{\text{eq}}=J^{b}_{\text{eq}}.$ (22)

Condition (22) is the principle of detailed balance, which follows from
microscopic reversibility [84], and here corresponds to the common reaction
rate $r$ for both reactive directions. While (22) seems rather trivial, as it
is simply the steady state condition for one reaction, in systems with
multiple reactions one may also find steady states for which the reactive
fluxes do not vanish separately, as first discussed by Wegscheider [85]. These
states need to be distinguished from genuine thermodynamic equilibrium given
by detailed balance (22). We note that by introducing the free energy
functional

$F=\int_{V}\left(k_{b}Tn_{1}\left[\ln(n_{1}/n_{0,1})-1\right]+k_{b}Tn_{2}\left[\ln(n_{2}/n_{0,2})-1\right]\right)\text{d}V,$
(23)

we can express the chemical potentials (16) as $\mu_{\alpha}=\frac{\delta
F}{\delta n_{\alpha}}$ and therefore obtain as time evolution equations

$\displaystyle\frac{\text{d}n_{1}}{\text{d}t}$
$\displaystyle=r\left[\exp\left(\frac{2\frac{\delta F}{\delta
n_{1}}+\frac{\delta F}{\delta
n_{2}}}{k_{b}T}\right)-\exp\left(3\frac{\frac{\delta F}{\delta
n_{1}}}{k_{b}T}\right)\right],$ (24)
$\displaystyle\frac{\text{d}n_{2}}{\text{d}t}$
$\displaystyle=-r\left[\exp\left(\frac{2\frac{\delta F}{\delta
n_{1}}+\frac{\delta F}{\delta
n_{2}}}{k_{b}T}\right)-\exp\left(3\frac{\frac{\delta F}{\delta
n_{1}}}{k_{b}T}\right)\right].$ (25)

The more general form of (24) and (25) dates back to Marcelin [61] and De
Donder [62, 86] and suggests applicability to nonideal systems and to more
complicated energy functionals as well as spatially extended systems. While we
have motivated the form (24)-(25) starting from simple mass action kinetics,
we note that the same expression can be derived by considering chemical
reactions as a diffusion process in the (continuous) configuration space of
the reactive complex so that this form indeed remains valid for nonideal
systems [87, 88, 89]. We now show that detailed balanced mass action type
kinetics such as Eqs. (24)-(25) is a gradient dynamics, i.e., that $F$
monotonically decreases. To this end, we treat a more general scenario and
consider a spatially extended system comprising the components
$N_{1},\ldots,N_{Q}$ with the particle densities $n_{1},\ldots,n_{Q}$ in a
closed, thermostatted box of volume $V$. Between the $Q$ components, $R$
chemical reactions may occur. For clarity of notation, we exclusively use
greek letters for field variables and latin letters for reaction indices. Each
reaction $j$ is characterized by its stoichiometric coefficients
$\nu^{f,b}_{\beta j}\geq 0$ for each component $\beta$ in the forward and
backward directions, respectively. Each reaction $j$ can then be summarized as

$\sum_{\beta=1}^{Q}\nu^{f}_{\beta
j}N_{\beta}\leftrightharpoons\sum_{\beta=1}^{Q}\nu^{b}_{\beta j}N_{\beta}.$
(26)

Note that because chemical reactions conserve total mass, the stoichiometric
coefficients must obey the condition

$\sum_{\beta=1}^{Q}m_{\beta}\nu^{f}_{\beta
j}=\sum_{\beta=1}^{Q}m_{\beta}\nu^{b}_{\beta j},$ (27)

where $m_{\beta}$ is the mass per particle for component $\beta$. Eq. (27)
states that the total mass of reactants and products must be identical. For
instance, applying (27) to the autocatalysis reaction (13) implies that the
molecular masses must be equal. In addition to conservation of the total mass
$M=\int_{V}\left[m_{1}n_{1}+m_{2}n_{2}\right]\text{d}V$ the total particle
count $N=\int_{V}\left[n_{1}+n_{2}\right]\text{d}V$ is then also conserved, as
already observed from Eqs. (14)-(15). One should always take care that the
particular choices of molecular masses (which may, e.g., explicitly appear in
some transport coefficients) and of stoichiometric coefficients do not
contradict mass conservation. In analogy to (21), we associate with each
reaction $j$, specified by (26), a detailed balanced mass action type kinetics

$J_{j}=J_{j}^{f}-J_{j}^{b}=r_{j}\left[\exp\left(\frac{\sum_{\beta=1}^{Q}\nu^{f}_{\beta
j}\frac{\delta F}{\delta
n_{\beta}}}{k_{b}T}\right)-\exp\left(\frac{\sum_{\beta=1}^{Q}\nu^{b}_{\beta
j}\frac{\delta F}{\delta n_{\beta}}}{k_{b}T}\right)\right],$ (28)

where $r_{j}>0$ is again the common rate constant for both reaction directions
and $F$ is the free energy (or any other appropriate thermodynamic
functional). Note that (28) can be expanded close to chemical equilibrium
($A=\sum_{\beta=1}^{Q}(\nu^{b}_{\beta j}-\nu^{f}_{\beta j})\frac{\delta
F}{\delta n_{\beta}}\ll k_{b}T$) to obtain a flux that is linear in the
affinity $A$ [58]. To demonstrate applicability to spatially extended systems,
we additionally assume diffusive transport such that the total time evolution
of the particle densities is given by

$\frac{\partial n_{\alpha}}{\partial
t}=\sum_{\beta=1}^{Q}\nabla\cdot\left(L_{\alpha\beta}\nabla\frac{\delta
F}{\delta n_{\beta}}\right)-\sum_{j=1}^{R}(\nu^{f}_{\alpha j}-\nu^{b}_{\alpha
j})J_{j},$ (29)

where the diffusive coefficients $L_{\alpha\beta}$ form a symmetric, positive
(semi)-definite matrix and $\nabla=(\frac{\partial}{\partial
x},\frac{\partial}{\partial y},\frac{\partial}{\partial z})^{\text{T}}$. Note
that (29) conserves the total mass
$M=\int_{V}\left[\sum_{\beta=1}^{Q}m_{\beta}n_{\beta}\right]\text{d}V$ due to
(27). In general, the total number of conserved quantities is determined as
$Q-S$, where $S$ is the dimension of the stoichiometric subspace $\mathcal{S}$
that is spanned by the vectors $\bm{\nu}_{j}=(\nu_{\alpha j})$. Conserved
quantities are then given as linear combinations
$\sum_{\alpha}c_{\alpha}n_{\alpha}$, where $\mathbf{c}=(c_{\alpha})$ are
vectors in the orthogonal complement of $\mathcal{S}$, i.e.,
$\mathbf{c}\cdot\bm{\nu}_{j}=0$ for all $j$ [90]. Since chemical reactions
conserve mass, this implies that particular behavior described for reaction-
diffusion systems with conservation laws [91, 92, 93, 94, 95, 96, 97, 98, 99]
is most likely a better representation of generic behavior than the more
frequently investigated fully open reaction-diffusion systems without any
conservation law.

We next show that (29) constitutes a gradient dynamics. For the free energy
$F$ we have

$\displaystyle\frac{\text{d}F}{\text{d}t}$
$\displaystyle=\int_{V}\left[\sum_{\beta=1}^{Q}\frac{\delta F}{\delta
n_{\beta}}\frac{\partial n_{\beta}}{\partial t}\right]\text{d}^{3}r$ (30)
$\displaystyle=\int_{V}\left[\sum_{\alpha,\beta=1}^{Q}\frac{\delta F}{\delta
n_{\beta}}\nabla\cdot\left(L_{\alpha\beta}\nabla\frac{\delta F}{\delta
n_{\alpha}}\right)\right]\text{d}V-\sum_{j=1}^{R}\int_{V}\left[\sum_{\beta=1}^{Q}\frac{\delta
F}{\delta n_{\beta}}(\nu^{f}_{\beta j}-\nu^{b}_{\beta
j})J_{j}\right]\text{d}^{3}r$ (31)
$\displaystyle=-\int_{V}\left[\sum_{\alpha,\beta=1}^{Q}L_{\alpha\beta}\left(\nabla\frac{\delta
F}{\delta n_{\beta}}\right)\cdot\left(\nabla\frac{\delta F}{\delta
n_{\alpha}}\right)\right]\text{d}V-k_{b}T\sum_{j=1}^{R}\int_{V}\ln\left(J^{f}_{j}/J^{b}_{j}\right)J_{j}\>\text{d}^{3}r$
(32)
$\displaystyle=-\int_{V}\left[\sum_{\alpha,\beta=1}^{Q}L_{\alpha\beta}\left(\nabla\frac{\delta
F}{\delta n_{\beta}}\right)\cdot\left(\nabla\frac{\delta F}{\delta
n_{\alpha}}\right)\right]\text{d}V-k_{b}T\sum_{j=1}^{R}\int_{V}\ln\left(J^{f}_{j}/J^{b}_{j}\right)(J_{j}^{f}-J_{j}^{b})\>\text{d}^{3}r$
(33) $\displaystyle\leq 0.$ (34)

From (30) to (31) we have used the time evolution equations (29). In (32),
partial integration was performed on the first term assuming, e.g., no-flux
boundary conditions. The second term was re-expressed in (32) and (33) using
Eq. (28). We stress that the transformation from (31) to (32) requires a
common rate constant $r_{j}$ for each reaction $j$ and is therefore only
possible for detailed balanced kinetics. The final inequality follows from the
positive definiteness of the transport matrix
$\mathbf{\underline{L}}=(L_{\alpha\beta})$ and from the inequality
$\left[f(b)-f(a)\right]\left[b-a\right]\geq 0$ for any monotonic $f$ if $b\geq
a$. We thus conclude that any bounded $F$ is a Lyapunov functional to the
dynamics (29) with its minimum corresponding to thermodynamic
equilibrium444For an analogous proof for the entropy production, see for
example [100].. There exist several ways to go beyond the described
relaxational dynamics, i.e., to obtain sustained out-of-equilibrium dynamics
(e.g., oscillatory dynamics), as found in many active systems. A common
strategy consists of chemostatting one or several of the reactive components,
i.e., keeping them at a constant chemical potential [101, 102, 63]. This
breaks detailed balance (28), since we can then write the fluxes as

$\displaystyle J_{j}$
$\displaystyle=r_{j}\left[\exp\left(\frac{\sum_{\beta=1}^{Q^{\prime}}\nu^{f}_{\beta
j}\frac{\delta F}{\delta n_{\beta}}+\sum_{\beta=1}^{E}\nu^{f}_{\beta
j}\mu_{\beta,0}}{k_{b}T}\right)-\exp\left(\frac{\sum_{\beta=1}^{Q^{\prime}}\nu^{b}_{\beta
j}\frac{\delta F}{\delta n_{\beta}}+\sum_{\beta=1}^{E}\nu^{b}_{\beta
j}\mu_{\beta,0}}{k_{b}T}\right)\right]$ (35)
$\displaystyle=\tilde{r}^{f}_{j}\exp\left(\frac{\sum_{\beta=1}^{Q^{\prime}}\nu^{f}_{\beta
j}\frac{\delta F}{\delta
n_{\beta}}}{k_{b}T}\right)-\tilde{r}^{b}_{j}\exp\left(\frac{\sum_{\beta=1}^{Q^{\prime}}\nu^{b}_{\beta
j}\frac{\delta F}{\delta n_{\beta}}}{k_{b}T}\right),$ (36)

where $E$ is the number of chemostatted species that are kept at the
respective chemical potentials $\mu_{\beta,0}$ and $Q^{\prime}$ is the number
of nonchemostatted species. The constant chemical potentials can be absorbed
into the effective rates

$\displaystyle\tilde{r}^{f}_{j}$
$\displaystyle=r_{j}\exp\left(\frac{\sum_{\beta=1}^{E}\nu^{f}_{\beta
j}\mu_{\beta,0}}{k_{b}T}\right),$ (37) $\displaystyle\tilde{r}^{b}_{j}$
$\displaystyle=r_{j}\exp\left(\frac{\sum_{\beta=1}^{E}\nu^{b}_{\beta
j}\mu_{\beta,0}}{k_{b}T}\right),$ (38)

which are generally distinct such that the proof (30)-(34) breaks down.
Alternatively, one may directly break detailed balance by introducing
different rates $r^{f}_{j},r^{b}_{j}$ in the forward and backward directions
or by simply assuming irreversible reactions ($r^{b}_{j}=0$ for some $j$).
This may be appropriate, e.g., when treating active protein and enzymatic
reactions [103, 7, 104] or in open systems if the reaction products are
immediately removed. A combination of both strategies can be found in the
construction of the Brusselator model [105, 106], where the presence of
irreversible reactions and assumed constant concentrations of certain
reactants lead to chemical oscillations. Finally, we note that breaking the
principle of detailed balance does not necessarily result in models that allow
for ‘active dynamics’. Even for (well-mixed) open systems, there exist wide
classes of mass action reaction networks with ‘complex balance’ that
guarantees the existence of a Lyapunov-function and a unique equilibrium point
which, however, does not coincide with thermodynamic equilibrium in general
[107, 90, 108].

###  3.2 Reactive cross-couplings in nonideal systems

Having shown that fluxes of the form (28) result in a proper gradient
dynamics, they may now be used to incorporate chemical reactions into models
of systems with more complicated free energies and transport processes
possibly resulting in chemo-mechanical coupling. For nonideal systems, Eq.
(28) then yields mass action type kinetics with reactive fluxes that contain
additional factors from nonideal or even mechanical contributions to the
chemical potentials. A common example comprises the interaction between
chemical reactions and phase separation [101, 109, 63]. There, the free energy
is amended, e.g., by adding Flory-Huggins-type interactions between the
different species [110, 111] as well as interface energies. In consequence,
the reactive fluxes additionally depend on the presence or absence of
interfaces and on the interaction parameters. In the context of reactive
wetting, reactions may be influenced by the proximity of the solid-liquid
interface due to reactant-substrate interactions which may become crucial in
the contact-line region. Examples include effects of reactant-dependent
wettability [16, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39] or of
substrate-mediated condensation [112, 113, 114]. Also, for reactions of
surface-active species (surfactants), the reactive fluxes can depend on the
shape of the liquid-gas interface, leading to an effective chemo-mechanical
coupling. This may be relevant to systems studied, e.g., in [115, 20, 116,
117].

##  4 Applications to reactive thin-film hydrodynamics

Next, we illustrate the gradient dynamics approach from Section 3.1 for two
distinct systems involving wetting phenomena and chemical reactions that
couple to interfacial physics. For simplicity, we restrict ourselves to
shallow drops and thin films of partially wetting liquid on a flat, solid
substrate. In Section 4.1, we treat drops of a suspension (or solution) where
the suspended particles within the drop may adsorb onto the substrate and
particles on the substrate may desorb into the drop. The adsorbed particles in
turn modify the wettability of the substrate, leading to sustained running
drops if the system is driven by chemostatting. In Section 4.2, we treat drops
of a simple liquid covered by two species of insoluble surfactants which are
continuously converted into each other by an autocatalytic reaction.
Additionally, the drop surface is in contact with a chemostatted ‘bath’ that
acts as an external supply of surfactant. Marangoni fluxes then couple the
reactive surface dynamics to the liquid bulk. As the chemostat-induced driving
force is varied, we find, inter alia, Turing patterns of surfactant
distributions, swaying and breathing modes of self-excited drop oscillation as
well as self-propelled drops. In both cases, the gradient dynamics approach
allows for thermodynamically consistent modeling, including cross-couplings
and nonideal energetic contributions.

###  4.1 Reactive wetting – Drops driven by chemically sustained wettability
gradients

As a first example, we consider drops of a partially wetting suspension (or
solution) that are situated on a flat, solid substrate (Fig. 1). The suspended
particles may adsorb from the drop onto the substrate and vice versa, thereby
altering the wettability of the substrate. Here, we treat the case that the
adsorbate renders the substrate less wettable. Similar systems have been
studied experimentally and theoretically in [16, 28, 29, 30, 31, 32, 33, 34,
35, 36, 37, 38, 39] where it has been shown that the induced wettability
gradient leads to self-propelled drop motion. We assume that the occuring
adsorption-desorption process is appropriately described by the conversion
reaction

$A\stackrel{{\scriptstyle r_{0}}}{{\leftrightharpoons}}B,$ (39)

where $A$ and $B$ denote particles that are either suspended in the liquid or
adsorbed onto the substrate. The corresponding particle densities per unit
volume of liquid and per unit substrate area are given by $a(\vec{x},t)$ and
$b(\vec{x},t)$, respectively. Note that for a meaningful representation in the
thin-film geometry, $a(\vec{x},t)$ must correspond to a $z$-averaged quantity.
The reaction itself obeys the principle of detailed balance with the rate
constant $r_{0}$. The local film height of the drop is $h(\vec{x},t)$.

Figure 1: Sketch of the considered geometry. A shallow drop of partially
wetting suspension (or solution) is situated on a flat solid substate. The
particles (or the solute) may adsorb from the drop onto the substrate or
desorb from the substrate into the drop, thus altering the wettability. The
rate constant of this reaction is given by $r_{0}$ in both directions. Due to
the (developing) wettability gradient, the drop moves with the velocity $v$.
The local film height of the drop is $h(\vec{x},t)$, the $z$-averaged local
density of suspended particles per unit liquid volume is $a(\vec{x},t)$ and
the local density of adsorbed particles per unit substrate area is
$b(\vec{x},t)$

With this system, we associate the free energy

$F=\int_{\Omega}\left[f(h,b)+\frac{\gamma}{2}|\nabla
h|^{2}+hg_{a}(a)+g_{b}(b)+h\frac{\sigma_{a}}{2}|\nabla
a|^{2}+\frac{\sigma_{b}}{2}|\nabla b|^{2}\right]\text{d}^{2}r$ (40)

with integration over the system domain $\Omega$. In (40), $f$ denotes the
adsorbate-dependent wetting energy, $\gamma$ denotes the constant surface
tension of the liquid, $g_{a}$ and $g_{b}$ are free energy contributions of
$A$ and $B$ per unit liquid volume and unit substrate area, respectively, and
$\sigma_{a}$, $\sigma_{b}$ denote interfacial stiffnesses that penalize
gradients in $a$ and $b$. First, we compute the liquid pressure and the
chemical potentials as variational derivates of (40). To this end, we observe
that changes in the local particle count of $A$ do not exclusively correspond
to changes in $a(\vec{x},t)$ in the thin-film geometry. This can be seen by
considering the particle count of $A$ in an infinitesimal liquid column of
base area $\text{d}\Omega$, which is given by

$\text{d}A=ah\cdot\text{d}\Omega=\hat{a}\cdot\text{d}\Omega$ (41)

with $\hat{a}=ah$. Because the substrate area is not subject to any dynamics
($\text{d}\Omega=const.$), Eq. (41) expresses a one-to-one correspondence
between the particle count of species $A$ and the field $\hat{a}$, i.e., all
changes in $\text{d}A$ translate to changes in $\hat{a}$ and vice versa.
Therefore, to compute the correct liquid pressure and chemical potentials, the
free energy must be varied with respect to $h$, $b$ and $\hat{a}$ instead of
$a$, as these fields directly relate to respective changes in the local liquid
volume and the local particle counts of $B$ and $A$. The variations are then

$\displaystyle p=\frac{\delta F}{\delta h}$
$\displaystyle=\partial_{h}f-\gamma\Delta
h+g_{a}-ag_{a}^{\prime}+\frac{\sigma_{a}}{2}|\nabla
a|^{2}+\sigma_{a}\frac{a}{h}\nabla\cdot\left(h\nabla a\right),$ (42)
$\displaystyle\mu_{a}=\frac{\delta F}{\delta\hat{a}}$
$\displaystyle=g_{a}^{\prime}-\frac{\sigma_{a}}{h}\nabla\cdot\left(h\nabla
a\right),$ (43) $\displaystyle\mu_{b}=\frac{\delta F}{\delta b}$
$\displaystyle=g_{b}^{\prime}+\partial_{b}f-\sigma_{b}\Delta b,$ (44)

where $g_{a}^{\prime}$ and $g_{b}^{\prime}$ are the univariate derivates of
$g_{a}$ and $g_{b}$, respectively. The dynamical equations then take on the
gradient dynamics form

$\displaystyle\partial_{t}h$
$\displaystyle=\nabla\cdot\left[Q_{hh}\nabla\frac{\delta F}{\delta
h}+Q_{ha}\nabla\frac{\delta F}{\delta\hat{a}}+Q_{hb}\nabla\frac{\delta
F}{\delta b}\right],$ (45) $\displaystyle\partial_{t}\hat{a}$
$\displaystyle=\nabla\cdot\left[Q_{ah}\nabla\frac{\delta F}{\delta
h}+Q_{aa}\nabla\frac{\delta F}{\delta\hat{a}}+Q_{ab}\nabla\frac{\delta
F}{\delta b}\right]+J_{r},$ (46) $\displaystyle\partial_{t}b$
$\displaystyle=\nabla\cdot\left[Q_{bh}\nabla\frac{\delta F}{\delta
h}+Q_{ba}\nabla\frac{\delta F}{\delta\hat{a}}+Q_{bb}\nabla\frac{\delta
F}{\delta b}\right]-J_{r},$ (47)

where the symmetric, positive definite transport matrix
$\mathbf{\underline{Q}}$ is

$\mathbf{\underline{Q}}=\left(\begin{array}[]{ccc}Q_{hh}&Q_{ha}&Q_{hb}\\\
Q_{ah}&Q_{aa}&Q_{ab}\\\
Q_{bh}&Q_{ba}&Q_{bb}\end{array}\right)=\left(\begin{array}[]{ccc}\frac{h^{3}}{3\eta}&\frac{h^{2}\hat{a}}{3\eta}&0\\\
\frac{h^{2}\hat{a}}{3\eta}&\frac{h\hat{a}^{2}}{3\eta}+D_{a}\hat{a}&0\\\
0&0&D_{b}b\end{array}\right)$ (48)

with the dynamic viscosity $\eta$ and the diffusive coefficients $D_{a,b}$.
Note that upper left $2\times 2$-block of (48) directly corresponds to the
transport coefficients of the gradient dynamics formulation for thin films of
mixtures and suspensions in the case without slip at the substrate [75]. In
consequence, only diffusive transport is possible for species $B$ (with
mobility $D_{b}b$). The reactive flux of the adsorption-desorption process is
given by the detailed balanced kinetics

$J_{r}=r(h)\left(\exp\left[\frac{1}{k_{b}T}\frac{\delta F}{\delta
b}\right]-\exp\left[\frac{1}{k_{b}T}\frac{\delta
F}{\delta\hat{a}}\right]\right).$ (49)

We assume that the reaction rate $r(h)$ in (49) is film height-dependent to
impose that adsorption and desorption may only occur where the macroscopic
drop is in contact with the substrate. The rate takes the form

$r(h)=r_{0}\cdot\xi(h)=r_{0}\cdot\frac{1}{2}\left[1+\tanh\left(\frac{h-h_{0}}{\delta
h}\right)\right],$ (50)

where $r_{0}$ is the actual rate constant of the reaction and $\xi(h)$ is a
smooth step-like function. The quantity $h_{0}$ is chosen slightly larger than
the thickness of the equilibrium liquid adsorption layer (precursor film) and
$\Delta h>0$ is a measure for the steepness of the step that is chosen much
smaller than the maximal drop height. Alternatively, the wetting energy may be
adapted in such a way that no solute enters the adsorption layer. We stress
that Eqs. (50) conserves the gradient dynamics form, as (49) still expresses
detailed balanced kinetics.

####  4.1.1 Reduced two-field model with chemostatting

In principle, the free energy (40) can be supplemented by specific choices for
$f,g_{a}$ and $g_{b}$ such that equations (45)-(47) describe the time
evolution of a particular system. However, here we further simplify (45)-(47)
by assuming that species $A$ within the drop is kept at constant and uniform
density $a_{0}$ by means of chemostatting. This may express that it is present
in vast excess of what is needed to cover the substrate. Then, Eqs. (45)-(47)
reduce to the two-field description

$\displaystyle\partial_{t}h$
$\displaystyle=\nabla\cdot\left[\frac{h^{3}}{3\eta}\nabla p\right],$ (51)
$\displaystyle\partial_{t}b$
$\displaystyle=\nabla\cdot\left[D_{b}b\nabla\mu_{b}\right]-J_{r},$ (52)

with the pressure

$\displaystyle p$ $\displaystyle=\partial_{h}f-\gamma\Delta h,$ (53)

and the chemical potential of $B$

$\displaystyle\mu_{b}$
$\displaystyle=g_{b}^{\prime}+\partial_{b}f-\sigma_{b}\Delta b.$ (54)

The reactive flux becomes

$J_{r}=r(h)\left(\exp\left[\frac{\mu_{b}}{k_{b}T}\right]-\exp\left[\frac{\mu_{a,0}}{k_{b}T}\right]\right),$
(55)

with $\mu_{a,0}$ as the chemostatted (uniform and constant) chemical potential
of species $A$ where $r(h)$ is still given by (50). In general, chemostatting
breaks the gradient dynamics form and results in active rather than
thermodynamic behavior. However, in the special case $\mu_{a,0}=0$ there is no
energetic influx associated with chemostatting and dissipation remains
unperturbed. In consequence, the system relaxes to thermodynamic equilibrium
(defined by the separate vanishing of all fluxes) and it can be shown that
Eqs. (51)-(52) then represent a gradient dynamics (see Appendix A.1).

####  4.1.2 Free energy choices

We now supplement the simplified model (51)-(52) with specific choices for $f$
and $g_{b}$. The wetting energy is

$f(h,b)=A_{H}\left(1+\lambda\frac{b}{b_{0}}\right)\cdot\left(-\frac{1}{2h^{2}}+\frac{h_{a}^{3}}{5h^{5}}\right),$
(56)

where $A_{H}$ is the Hamaker constant, $\lambda>0$ is a dimensionless
proportionality factor, $b_{0}$ is some characteristic density of $B$ and
$h_{a}$ denotes the thickness of the equilibrium liquid adsorption layer
(precursor film). The choice (56) models a substrate that is rendered less
wettable by the adsorbate and thus corresponds to the regime of partial
wetting since the spreading coefficient $S=f(h_{a},b)$ [9, 10, 68] is negative
for all $b\geq 0$. Further, we employ the purely entropic local free energy

$g_{b}(b)=k_{b}Tb\left[\ln\left(b/b_{0}\right)-1\right]$ (57)

for species $B$. To study the model numerically, Eqs. (51)-(52) are
nondimensionalized, see Appendix A.3.1. From here on, all symbols denote
dimensionless quantities.

####  4.1.3 Numerical results

Figure 2: (a) Snapshots of the film height profile $h$ and the adsorbate
density profile $b$ at different times in the relaxational case at
$\mu_{a,0}=0$. The drop first moves away from its initial position because
initially at $t=0$ an asymmetric wettability profile is imposed via $b(x,0)$.
Then, the drop slows down due to the equilibration of $b$. As equilibrium is
approached, $b$ becomes uniform inside and outside of the drop with interfaces
in $b$ in the contact line regions. (b) Dependence of the tangents of the
advancing (solid line) and receding (dashed line) contact angles on time.
Starting with an initially symmetric drop, the initial wettability gradient
causes a difference in contact angles, with $\theta_{a}<\theta_{r}$. At the
final equilibrium, $\theta_{a}=\theta_{r}$. (c) The free energy $F$ decreases
monotonically with time since the system is passive. The remaining parameters
are $W=1,\lambda=0.5,r_{0}=0.0025,D_{b}=0.001,h_{0}=2$ and $\delta h=0.3$.
Only part of the computational domain $\Omega=\left[0,2000\right]$ is shown

Figure 3: Space-time plots illustrating the time evolution of (a) the film
height profile $h$ and (b) the adsorbate density profile $b$ in the active
case at $\mu_{a,0}=\ln 6$. The drop first moves due to an initial asymmetric
wettability profile. After a short transient, the drop moves across the
substrate at constant speed while self-sustaining the driving wettability
gradient. (c) Dependence of the tangents of the advancing (solid line) and
receding (dashed line) contact angles on time. After the acceleration phase,
both contact angles remain distinct while the difference very slowly
decreases. (d) Dependence of the drop’s center-of-mass velocity $v$ on time.
After the transient, $v$ slowly decreases due to adsorbate diffusion. (e) The
free energy $F$ first shows a transient. Then, when the drop has reached its
constant speed, $F$ increases linearly. (f) Shown are the corresponding energy
rates associated with dissipation $R_{\text{diss}}$ (dashed line) and
chemostatting $R_{\text{chem}}$ (solid line). Note that $R_{\text{diss}}\geq
0$ per definition. After the transient, the chemostat supplies energy at
constant rate while the system constantly dissipates energy. The remaining
parameters are as in Fig. 2. The computational domain is
$\Omega=\left[0,2000\right]$

To elucidate some general features of the reduced model (51)-(52), we perform
time simulations of the nondimensional equations (94)-(95) using the finite-
element library oomph-lib [118] with periodic boundary conditions on the one-
dimensional domain $\Omega=\left[0,2000\right]$. Further details on the
numerical methods are given in Appendix B. We consider the passive
(relaxional) case at $\mu_{a,0}=0$ and the active (chemostat-driven) case at
$\mu_{a,0}\neq 0$, while all other parameters are identical. In both cases, we
initialize the simulation with a symmetric drop profile and an asymmetric
concentration profile, i.e., we impose an initial wettability gradient. In the
passive case (see Fig. 2), we find that the drop initially moves away from its
starting position driven by the initial wettability gradient. However, it
slows down subsequently due to the equilibration of $b$ (Fig. 2a). Note that
also the wettability to the right of the drop slowly decreases because $B$
slowly diffuses outside the drop (even across the periodic domain boundaries).
In the final equilibrium state, the concentration of $B$ is uniform inside and
outside the drop with a steep gradient located in the contact line region. The
difference in coverage can be traced back to the reduced chemical potential of
$B$ outside the drop due to the concentration-dependent wetting energy. We
note that the equilibrium value of $b$ inside the drop is determined as ‘bare’
reactive equilibrium, i.e., it must be unity in dimensionless units
(dimensional $b=b_{0}$). The equilibration process is also visible in Fig. 2b
that shows the dependence of the advancing and receding contact angles on
time. They approach the same equilibrium value. As the system is passive, the
free energy $F$ decreases monotonically (Fig. 2c).555While we only show part
of the time evolution of the free energy, we have checked that $F$ always
decreases until the final equilibrium state is reached.

Fig. 3 shows an example for the active case with $\mu_{a,0}>0$. There, we find
continuously moving drops: After a short transient, the drop moves along the
substrate at nearly constant speed maintained by the self-sustained
wettability gradient across its length. Note that the drop’s center-of-mass
velocity $v$ slightly decreases with time (Fig. 3d) due to slow diffusion of
$b$ across the periodic domain boundaries which reduces the wettability in
front of the drop.666A similar effect occurs for no-flux boundary conditions,
where some $B$ diffuses in front of the drop. For boundary conditions which
allow for constant $v$, see [33]. The advancing and receding dynamic contact
angles ($\theta_{a}$ and $\theta_{r}$) differ when the drop moves (Fig. 3c)
with the difference slightly decreasing with time due to the mentioned
diffusion. In contrast to the passive case, after the transient the free
energy increases approximately linearly with time (Fig. 3e). Considering the
rate of energy influx due to the chemostat $R_{\text{chem}}$ and the
dissipation rate $R_{\text{diss}}$ during the sustained drop motion (Fig. 3f,
also see Appendix A.2), both rates are roughly constant since there is a
continuous turnover of suspended particles. Namely, the chemostat replaces
particles in the drop that adsorb onto the substrate. Since the wettability
profile across the drop remains roughly stationary in the comoving frame, the
chemostat must supply energy at a constant rate. Similarly, energy is
dissipated at a constant rate. Therefore, because chemostatting outweighs
dissipation ($R_{\text{chem}}>R_{\text{diss}}$), the free energy increases
linearly, reflecting the continuous coating of the substrate. Strikingly, the
relative difference $(R_{\text{chem}}-R_{\text{diss}})/R_{\text{chem}}$
between the rates of chemostatting and dissipation is less than five percent,
i.e., only a small part of the energy from the chemostat remains in the system
while a much larger part is immediately dissipated.

###  4.2 Reactive surfactants – Sessile drops covered by autocatalytic
surfactants

Our second example are shallow drops of partially wetting liquid that are
covered by two species of insoluble reactive surfactants (Fig. 4). In general,
systems involving reactive surfactants are usually treated starting from a
hydrodynamics perspective where chemical reactions are added in an ad hoc
manner [13, 119, 115, 40, 44, 120]. Here, we employ the described gradient
dynamics approach to describe the interplay between autocatalytic chemical
reactions, Marangoni effects and wetting physics in a thermodynamically
consistent manner, i.e., with an appropriate passive limit. In particular, we
assume that the two species engage in the autocatalytic reaction

$2\Gamma_{1}+\Gamma_{2}\leftrightharpoons 3\Gamma_{1},$ (58)

where $\Gamma_{1}$ and $\Gamma_{2}$ denote the particle densities of the two
species (particles per unit interface area). Additionally, we allow for
surfactant exchange with an external chemostatted ‘bath’.

Figure 4: Sketch of the considered geometry. A shallow drop of a partially
wetting liquid on a flat solid substrate is covered by two species of
insoluble surfactant which engage in an autocatalytic conversion reaction.
Additionally, the two species are in contact with an external bath that acts
as a chemostat and exchanges surfactant with the drop surface
($J_{\text{ex}}$). The film height profile is described by $h$ and the
surfactant density profiles on the film surface by $\Gamma_{1,2}$. For a
gradient dynamics description, the densities have to be projected onto the
substrate plane where the projected densities are
$\tilde{\Gamma}_{1,2}=\xi\Gamma_{1,2}$. Here, $\xi$ is the metric factor of
the film surface

The free energy of the system is given by

$F=\int_{\Omega}\left[f(h)+\xi g(\Gamma_{1},\Gamma_{2})\right]\text{d}^{2}r,$
(59)

where $f(h)$ is the wetting energy and $g(\Gamma_{1},\Gamma_{2})$ is the
energy density of the surfactant-covered free surface. For simplicity, here,
we neglect a possible surfactant dependence of the wetting energy and,
similarly, a dependence of the surface energy on the film height [121]. Note
that we use the exact metric factor $\xi=\sqrt{1+|\nabla h|^{2}}$ and not its
long-wave approximation that was used in Section 4.1. To compute the correct
pressure and chemical potentials from $F$ we need to identify the independent
fields that correspond to the local liquid volume and the local particle
counts of surfactant in a small liquid column. In analogy to the argument in
Section 4.1 for the bulk concentration, the independent fields are the local
film height $h$ and the projected densities
$\tilde{\Gamma}_{\alpha}=\xi\Gamma_{\alpha}$ (particles per unit substrate
area). For further discussion, we refer to [74, 52]. The resulting pressure
and chemical potentials are

$\displaystyle p$ $\displaystyle=\frac{\delta F}{\delta
h}=\partial_{h}f-\nabla\cdot\left[\frac{1}{\xi}\left(g-\Gamma_{1}\partial_{\Gamma_{1}}g-\Gamma_{2}\partial_{\Gamma_{2}}g\right)\nabla
h\right],$ (60) $\displaystyle\mu_{\alpha}$ $\displaystyle=\frac{\delta
F}{\delta\tilde{\Gamma}_{\alpha}}=\partial_{\Gamma_{\alpha}}g,$ (61)

respectively. The coupled time evolution equations are

$\displaystyle\partial_{t}h$
$\displaystyle=\nabla\cdot\left[Q_{hh}\nabla\frac{\delta F}{\delta
h}+Q_{h\Gamma_{1}}\nabla\frac{\delta
F}{\delta\tilde{\Gamma}_{1}}+Q_{h\Gamma_{2}}\nabla\frac{\delta
F}{\delta\tilde{\Gamma}_{2}}\right],$ (62)
$\displaystyle\partial_{t}\tilde{\Gamma}_{1}$
$\displaystyle=\nabla\cdot\left[Q_{\Gamma_{1}h}\nabla\frac{\delta F}{\delta
h}+Q_{\Gamma_{1}\Gamma_{1}}\nabla\frac{\delta
F}{\delta\tilde{\Gamma}_{1}}+Q_{\Gamma_{1}\Gamma_{2}}\nabla\frac{\delta
F}{\delta\tilde{\Gamma}_{2}}\right]+J_{r}+J_{\text{ex},1},$ (63)
$\displaystyle\partial_{t}\tilde{\Gamma}_{2}$
$\displaystyle=\nabla\cdot\left[Q_{\Gamma_{2}h}\nabla\frac{\delta F}{\delta
h}+Q_{\Gamma_{2}\Gamma_{1}}\nabla\frac{\delta
F}{\delta\tilde{\Gamma}_{1}}+Q_{\Gamma_{2}\Gamma_{2}}\nabla\frac{\delta
F}{\delta\tilde{\Gamma}_{2}}\right]-J_{r}+J_{\text{ex},2},$ (64)

where the symmetric positive definite mobility matrix is given by

$\mathbf{\underline{Q}}=\left(\begin{array}[]{ccc}Q_{hh}&Q_{h\Gamma_{1}}&Q_{h\Gamma_{2}}\\\
Q_{\Gamma_{1}h}&Q_{\Gamma_{1}\Gamma_{1}}&Q_{\Gamma_{1}\Gamma_{2}}\\\
Q_{\Gamma_{2}h}&Q_{\Gamma_{2}\Gamma_{1}}&Q_{\Gamma_{2}\Gamma_{2}}\end{array}\right)=\left(\begin{array}[]{ccc}\frac{h^{3}}{3\eta}&\frac{h^{2}\Gamma_{1}}{2\eta}&\frac{h^{2}\Gamma_{2}}{2\eta}\\\
\frac{h^{2}\Gamma_{1}}{2\eta}&\frac{h\Gamma_{1}^{2}}{\eta}+D_{1}\Gamma_{1}&\frac{h\Gamma_{1}\Gamma_{2}}{\eta}\\\
\frac{h^{2}\Gamma_{2}}{2\eta}&\frac{h\Gamma_{1}\Gamma_{2}}{\eta}&\frac{h\Gamma_{2}^{2}}{\eta}+D_{2}\Gamma_{2}\end{array}\right).$
(65)

It represents a straightforward generalization of the mobilities for the
single-surfactant case [74] and contains advective and diffusive
contributions. Here, $D_{\alpha}>0$ are constant diffusive mobilities. For the
reactive flux, we have

$J_{r}=r\left[\exp\left(\frac{2}{k_{b}T}\frac{\delta
F}{\delta\tilde{\Gamma}_{1}}+\frac{1}{k_{b}T}\frac{\delta
F}{\delta\tilde{\Gamma}_{2}}\right)-\exp\left(\frac{3}{k_{b}T}\frac{\delta
F}{\delta\tilde{\Gamma}_{1}}\right)\right],$ (66)

where $r>0$ is the reaction rate. For simplicity, the surfactant exchange with
the external chemostatted bath is assumed to be proportional to the chemical
potential difference between the surfactants on the film and in the bath,
i.e.,

$\displaystyle J_{\text{ex},1}=-\beta_{1}\left(\frac{\delta
F}{\delta\tilde{\Gamma}_{1}}-\mu_{0,1}\right)=-\beta_{1}\frac{\delta
F}{\delta\tilde{\Gamma}_{1}}-\beta\mu,$ (67) $\displaystyle
J_{\text{ex},2}=-\beta_{2}\left(\frac{\delta
F}{\delta\tilde{\Gamma}_{2}}-\mu_{0,2}\right)=-\beta_{2}\frac{\delta
F}{\delta\tilde{\Gamma}_{2}}+\beta\mu.$ (68)

Here, $\beta_{\alpha}>0$ are transition rates for the transfer between the
ambient bath and the film surface and the $\mu_{0,\alpha}$ are the imposed
chemical potentials of the two species in the bath. Here, we restrict our
attention to part of the possible parameter space by assuming the two terms
$\beta_{\alpha}\mu_{0\alpha}$ to be equal in magnitude (given by $\beta\mu$)
and opposite in sign. With $\beta\mu>0$ and the sign choices as in (67)-(68),
$\Gamma_{1}$ tends to transfer from the film surface to the bath, with an
inverse tendency for $\Gamma_{2}$. The parameter $\beta\mu$ then represents
the sustained nonequilibrium driving force of the chemostat. For $\beta\mu=0$
there is no energy flux in or out of the system associated with chemostatting.
In this case, the system is passive and equations (62)-(64) again represents a
gradient dynamics.

Since the dynamical equations are expressed in terms of the projected
densities, the nonconserved terms explicitly depend on the metric factor of
the free surface. This corresponds to a basic form of geometry-induced chemo-
mechanical coupling which, here, leads to an overall slowing-down of the
reaction flux (66) on the convex drop surface. While the metric factor only
weakly deviates from unity in the long-wave limit and can then often be
neglected, on strongly curved surfaces and in bulk fluids chemo-mechanical
couplings can give rise to pattern forming instabilities and sustained
oscillations/motion [122, 123, 124, 125].

####  4.2.1 Free energy choices and further simplifications

Starting from the gradient dynamics form (62)-(64), we now supplement the
model with specific energies and introduce further simplifications. We assume
that the surfactant species only sparsely cover the film surface, implying
purely entropic contributions to the surface free energy, i.e.,

$g_{s}\left(\Gamma_{1},\Gamma_{2}\right)=\gamma^{0}+k_{b}T\,\Gamma_{1}\left[\ln\left(\Gamma_{1}a_{1}^{2}\right)-1\right]+k_{b}T\,\Gamma_{2}\left[\ln\left(\Gamma_{2}a_{2}^{2}\right)-1\right].$
(69)

Here, $\gamma^{0}$ is the surface tension of the bare surface and $a_{1,2}>0$
are characteristic length scales of the two surfactants. Additionally, we use
the simple wetting energy

$f(h)=A_{H}\left(-\frac{1}{2h^{2}}+\frac{h_{a}^{3}}{5h^{5}}\right),$ (70)

where $A_{H}$ is the Hamaker constant and $h_{a}$ is the equilibrium adsorpton
layer height (precursor thickness), corresponding to the minimum of $f(h)$.
Again, we consider the case of partial wetting, i.e., of negative spreading
coefficient $S=f(h_{a})<0$ [9, 10, 68]. The resulting pressure and chemical
potentials are

$\displaystyle p$
$\displaystyle=\partial_{h}f-\nabla\cdot\left[\tfrac{1}{\xi}\left(\gamma^{0}-k_{b}T\,\Gamma_{1}-k_{b}T\,\Gamma_{2}\right)\nabla
h\right],$ (71) $\displaystyle\mu_{\alpha}$
$\displaystyle=k_{b}T\ln(\Gamma_{\alpha}a_{\alpha}^{2}).$ (72)

Finally, we assume that $\xi\approx 1$ in the long-wave limit ($|\nabla h|\ll
1$) and that the changes in surface tension due to the presence of surfactant
are small compared to $\gamma^{0}$. Then the pressure reduces to

$p=\partial_{h}f-\gamma^{0}\Delta h.$ (73)

This approximation does not affect the Maragoni flows but only simplifies the
Laplace pressure. It is commonly used when considering Marangoni effects [65,
11]. The resulting time evolution equations are

$\displaystyle\partial_{t}h=$
$\displaystyle\>\nabla\cdot\left[\frac{h^{3}}{3\eta}\nabla
p\right]+k_{b}T\nabla\cdot\left[\frac{h^{2}}{2\eta}\nabla\Gamma_{1}\right]+k_{b}T\nabla\cdot\left[\frac{h^{2}}{2\eta}\nabla\Gamma_{2}\right],$
(74)
$\displaystyle\begin{split}\partial_{t}\Gamma_{1}=&\>\nabla\cdot\left[\frac{h^{2}\Gamma_{1}}{2\eta}\nabla
p\right]+k_{b}T\nabla\cdot\left[\left(\frac{h\Gamma_{1}}{\eta}+D_{1}\right)\nabla\Gamma_{1}\right]+k_{b}T\nabla\cdot\left[\frac{h\Gamma_{1}}{\eta}\nabla\Gamma_{2}\right]\\\
&+r\left[\left(\Gamma_{1}a_{1}^{2}\right)^{2}\left(\Gamma_{2}a_{2}^{2}\right)-\left(\Gamma_{1}a_{1}^{2}\right)^{3}\right]-\beta_{1}k_{b}T\ln(\Gamma_{1}a_{1}^{2})-\beta\mu,\end{split}$
(75)
$\displaystyle\begin{split}\partial_{t}\Gamma_{2}=&\>\nabla\cdot\left[\frac{h^{2}\Gamma_{2}}{2\eta}\nabla
p\right]+k_{b}T\nabla\cdot\left[\frac{h\Gamma_{2}}{\eta}\nabla\Gamma_{1}\right]+k_{b}T\nabla\cdot\left[\left(\frac{h\Gamma_{2}}{\eta}+D_{2}\right)\nabla\Gamma_{2}\right]\\\
&-r\left[\left(\Gamma_{1}a_{1}^{2}\right)^{2}\left(\Gamma_{2}a_{2}^{2}\right)-\left(\Gamma_{1}a_{1}^{2}\right)^{3}\right]-\beta_{2}k_{b}T\ln(\Gamma_{2}a_{2}^{2})+\beta\mu.\end{split}$
(76)

To numerically study the model, we nondimensionalize (74)-(76). The rescaled
model as well as details on the choice of scales and dimensionless parameters
can be found in Appendix A.3.2. In the following, all symbols denote
dimensionless quantities.

####  4.2.2 Numerical results

Figure 5: (a) Partial bifurcation diagram for drop states in dependence of the chemostat driving strength $\beta\mu$. The solution measure corresponds to the norm $||h||_{2}$. Shown are linearly stable branches (solid lines) and sections of the unstable branches (dashed lines) that directly relate to them. The horizontal black line represents the base state - the simple resting drop with homogeneous surfactant concentration. Blue lines correspond to steady drops with localized surfactant Turing-like patterns, orange lines are branches of steadily traveling drops and purple ones represent standing wave-like drop oscillations. Labeled filled circles denote selected states shown in Fig. 6. Panel (b) magnifies the region where a supercritical drift-pitchfork bifurcations occurs resulting in the emergence of the first two branches of traveling drops (left- and right traveling, otherwise identical states). The corresponding drop velocities $v$ are given in panel (c). The remaining parameters are $W=3,\delta=1,D_{1}=1.4,D_{2}=0.1,r=0.8,\beta_{1}=1$ and $\beta_{2}=0.5$ Figure 6: (a) Selected profiles of steady (P1 to P3) and travelling (T1) drop states. Film height $h$ and surfactant $\Gamma_{1}$ are given as solid and dashed lines, respectively, at the points marked by corresponding letters in Fig. 5a. The steadily traveling drop is shown in the comoving frame with speed $v$, the arrow indicates the direction of travel. (b), (c) Space-time plots of standing wave-like drop-surfactant oscillations at the points H1 and H2 in Fig. 5a. Shown are the film height $h$ and the surfactant $\Gamma_{1}$ over a time-span of two periods. The periods are (a) $T\approx 10.58$ and (b) $T\approx 20.71$ Figure 7: Space-time plots of standing wave-like drop-surfactant oscillations corresponding to (a) a stable breathing mode and (b) an unstable swaying mode. The parameters are (a) $W=4,\delta=1,D_{1}=2,D_{2}=0.15,r=1.1,\beta_{1}=6,\beta_{2}$ and $\beta\mu=5.3992$ for a drop volume $V\approx 816$ and domain size $L=100$, and (b) $W=1.504,\delta=1,D_{1}=4,D_{2}=0.3,r=1.2,\beta_{1}=6,\beta_{2}=1$ and $\beta\mu=5.7868$ for $V=720$ and $L=100$. Both modes are shown over a time-span of two periods. The periods are (a) $T\approx 6.04$ and (b) $T\approx 9.28$ Figure 8: (a) Space-time plot for a stable modulated traveling wave-like ‘crawling’ drop-surfactant state in the laboratory frame at $W=7,D_{1}=1.3,D_{2}=0.1,r=0.8,\beta_{1}=1.4,\beta_{2}=0.7$, $\beta\mu=2.2$, $V\approx 816$, $L=100$. (b) The L2-norm $||\psi||_{2}=\sqrt{1/L\int_{0}^{L}\left(h^{2}+p^{2}+\Gamma_{1}^{2}+\Gamma_{2}^{2}\right)\>\text{d}x}$ as a function of time. Because the L2-norm is invariant w.r.t. state translation, the observed modulated traveling wave results in time-periodic behavior Table 1: Summary of the symbol scheme used to differentiate bifurcation classes Symbol | Base branch | Bifurcation type
---|---|---
| steady state | pitchfork
$\Diamond$ | steady state | Hopf
standing wave | torus
$\blacktriangle$ | steady state | drift-pitchfork

We now show a few typical features of chemo-mechanical drop dynamics described
by Eqs. (74)-(76). We focus on two-dimensional drops, i.e., one-dimensional
substrates, and analyze the nondimensional model (107)-(109) given in Appendix
A.3.2. In our analysis we use numerical path continuation, employing the
continuation package pde2path [126], supplemented by time simulations obtained
with the finite-elements library oomph-lib [118]. Further details are given in
Appendix B. For the bifurcation diagram we use as solution measure the
$L_{2}$-norm of the film height

$||h||_{2}=\sqrt{\frac{1}{L}\int_{-L/2}^{L/2}h^{2}\>\text{d}x}$ (77)

with the domain size $L$. For time-periodic states we use the time-averaged
norm $1/T\int_{0}^{T}||h||_{2}\>\text{d}t$, where $T$ is the temporal period.
We use the chemostat driving strength $\beta\mu$ as main control parameter and
investigate how the drop behavior changes as the system is driven further away
from thermodynamic equilibrium. Most results are obtained for a periodic
domain of size $L=100$ and a drop volume of
$V=\int_{-L/2}^{L/2}h\>\text{d}x\approx 816$, corresponding to small drops
with a maximal height of about 20 times the adsorption layer height. In the
bifurcation diagram we indicate different bifurcations using the specific
symbols listed in Table 1.

Fig. 5a presents the resulting bifurcation diagram in dependence of the
driving strength $\beta\mu$. Only a selected part of the rich bifurcation
diagram is shown, namely, linearly stable resting, traveling and oscillating
drop states and relevant sections of directly related unstable branches. For
clarity, we do not show less relevant branches of unstable states. Selected
corresponding drop and surfactant profiles are given in Fig. 6.

The horizontal line in Fig. 5a corresponds to the base state – a simple
resting sessile drop with homogeneous surfactant density. It is the only shown
state that already exists in the passive limit ($\beta\mu=0$). As $\beta\mu$
is increased, the simple drop becomes unstable in a double pitchfork
bifurcation where two branches of unstable steady drop states with localized
Turing-like surfactant patterns emerge subcritically (blue lines). The first
patterned branch (once unstable) is then rendered linearly stable in a saddle-
node bifurcation and is shortly after again destabilized at a drift-pitchfork
bifurcation (filled triangle) [127, 128]. There, a pair of branches of
steadily traveling drops emerges (orange line) that connect the first
patterned branch to the second one (doubly unstable). The corresponding drop
velocities are given in Fig. 5b. Note that the traveling drops are mostly
stable states but become unstable in a saddle-node bifurcation shortly before
they connect to the second patterned branch in a drift-pitchfork bifurcation,
see Fig. 5b. Beyond this bifurcation, the second patterned branch stabilizes
in a saddle-node bifurcation, remains linearly stable for an extended
$\beta\mu$-range, before it becomes unstable at a Hopf bifurcation, and is
here not further discussed. There, a branch of time-periodic states emerges
that correspond to standing wave-like oscillations.777Due to the relaxed
numerical tolerance and mesh resolution necessary for the continuation of this
branch, the end of this branch does not fully match the location of the Hopf
point of the patterned steady branch. However, we have checked that the
profiles and periods match at bifurcation and have therefore shifted the
L2-norm of the branch by $0.00165$ for presentational purposes. The stable
oscillatory state (H1) is illustrated in Fig. 6b. The oscillation in the
surfactant pattern is relatively weak and the one of the drop shape is not
visible on the scale of the figure, indicating a weak coupling between
chemical reactions and bulk liquid. Remarkably, there exists a range of
driving strengths where the simple drop branch, the second patterned branch
and the branch of steadily traveling drops are simultaneously linearly stable.

Time simulations at larger values of $\beta\mu$ also result in time-periodic
states with noticeable surfactant-bulk coupling (cf. state H2 in Fig. 5a,
space-time plots in Fig. 6c) leading to left-right symmetric drop
oscillations, i.e., a breathing mode. Following the corresponding branch by
continuation, we find that it destabilizes in a saddle-node bifurcation at a
smaller value of $\beta\mu$ and in a torus bifurcation at a larger $\beta\mu$.
The emerging linearly stable quasi-periodic drop-surfactant oscillations are
also found in time simulations initialized just beyond this bifurcation (not
shown).

Despite its already rather complex structure, the bifurcation diagram in Fig.
5a contains general features that are also found at stronger surfactant-bulk
coupling, however, within much more involved bifurcation structures (not
shown). Seemingly universal key elements are the emergence of patterned
surfactant branches in double pitchfork bifurcations from the simple drop
branch, the pairwise interconnection of these patterned branches by branches
of steadily traveling drops as ‘rung states’ and the virtually exclusive
existence of quasi-periodic states at large driving strengths. Furthermore, in
general, we find two modes of drop oscillation, breathing and swaying drops,
which are illustrated in Fig. 7. For breathing modes (Fig. 7a), the film
height profile and the surfactant distributions are left-right symmetric at
all times. Swaying modes (Fig. 7b) are generally antisymmetric, i.e., film
height and surfactant profiles that are offset by half a temporal period are
related by reflection. Both modes normally emerge in Hopf bifurcations.
Lastly, in time simulations we also find stable modulated traveling drops that
represent a more complicated form of transport (Fig. 8). Such states may for
instance emerge from drift-pitchfork bifurcations of standing wave-like states
or from Hopf bifurcations of steadily traveling drop states. In time
simulations, the drop begins to drift after an initial spontaneous breaking of
the left-right symmetry and it oscillates asymmetrically while moving. In the
case of Fig. 8a, the drop periodically forms a large protrusion in the
direction of propagation and subsequently contracts again to a cap-like shape,
resulting in net movement of the whole drop. The periodicity in the comoving
frame is best appreciated in the dependence of the L2-norm
$||\psi||_{2}=\sqrt{1/L\int_{0}^{L}\left(h^{2}+p^{2}+\Gamma_{1}^{2}+\Gamma_{2}^{2}\right)\>\text{d}x}$
on time (Fig. 8b). This measure is minimal and maximal when the drop’s
protrusion is at maximal length and when the drop is nearly cap-like (fully
contracted), respectively.

##  5 Conclusion

We have introduced a way to write reactive mesoscopic hydrodynamics (reactive
thin-film models) in gradient dynamics form. On the one hand, this
thermodynamically consistent description allows to consider passive processes
where the system ultimately approaches thermodynamic equilibrium. On the other
hand, it allows for a controlled introduction of persistent nonequilibrium
driving by making the dynamical couplings nonreciprocal (breaking the Onsager
relations of transport coefficients and introducing indefinite transport
matrices [129, 130], breaking detailed balance of chemical reactions [103,
104]) or by chemostatting certain reactant species [101, 102, 63]. One could
also render the interactions nonreciprocal by introducing driving force terms
that cannot be derived from an energy functional [131, 132, 133, 134, 135].

To introduce the general gradient dynamics form we have first briefly reviewed
how thin-film equations without chemical reactions can be written as gradient
dynamics in the framework of linear nonequilibrium thermodynamics. Next, we
have reviewed how chemical reactions governed by detailed balance are brought
in gradient dynamics form, however, this time employing nonlinear
nonequilibrium thermodynamics. Then, we have combined the linear approach for
transport processes (and possibly simple transfer processes between phases)
with the nonlinear approach for chemical reactions. The resulting rather
general model model may be employed to describe various phenomena that involve
the physics of wetting for simple or complex liquids and chemical reactions
that couple to interface dynamics.

In particular, we have illustrated the gradient dynamics form for chemical
reactions first for a single reaction with mass action kinetics in a spatially
homogeneous system. By re-expressing the reaction rates via the chemical
potentials and introducing the principle of detailed balance from simple
thermodynamic considerations, we have motivated a well-known form of mass
action type kinetics. The approach corresponds to a gradient dynamics with
nonlinear relations between the variations of the underlying energy functional
and the thermodynamic fluxes. Active and sustained out-of-equilibrium dynamics
can then be obtained by breaking the detailed balanced structure of the
kinetics. This may, for instance, be achieved by chemostatting certain
reactant species or by directly assuming different reaction rates for the
forward and backward reactions (e.g., for irreversible reactions).

The gradient dynamics form retains its thermodynamic consistency even for
nonideal systems. This naturally leads to modifications of standard mass
action reaction kinetics when nonideal contributions to the energy functional
are present. In this way, a wide variety of cross-couplings between chemical
reactions and fluid-mechanical phenomena may be captured. Here, we have
focused on the combination of chemical reactions and wetting physics. We have
demonstrated the usefulness of the approach at the example of two mesoscopic
hydrodynamic models for shallow drops of partially wetting liquid on
homogeneous rigid solid substrates – describing reactive wetting and reactive
surfactants.

First, we have considered active wetting: drops of solutions or suspensions
where the suspended particles may adsorb from the drop onto the substrate or
desorb from the substrate into the drop, thereby altering the substrate’s
wettability. In this case, one finds that the adsorption-desorption process is
influenced by the proximity of the contact line due to interactions between
the adsorbed particles and the substrate. At thermodynamic equilibrium, this
leads to an additional interface between regions with different but
respectively homogeneous adsorbate coverage that is located in the contact
line region. When driving the system, e.g., by chemostatting the solute within
the drop, the drop may perpetually move across the substrate while rendering
the underlying subtrate less wettable. Such drops driven by chemically
sustained wettability gradients had already been experimentally and
theoretically studied [16, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39].
Here, we recover these results with a model based on a gradient dynamics form
where controlled chemostatting allows to quantify the ‘distance from a passive
situation’.

Second, we have considered reactive surfactants: sessile drops of simple
liquids are covered by two species of insoluble surfactant which engage in an
autocatalytic conversion reaction. Bringing the surfactants on the free
surface additionally in contact with a chemostat, i.e., a vast ‘surfactant
bath’, a through-flow of energy and mass is introduced. Systems involving
reactive surfactants had already been modeled in an ad hoc manner, e.g., for
liquid film flows in [13, 44]. Here, employing the developed gradient dynamics
description, we have found the emergence of a geometry-induced chemo-
mechanical coupling resulting from the interdependence of the autocatalytic
reaction and the free surface geometry that is known to result in chemo-
mechanical instabilities on strongly curved surfaces [122, 123]. Note,
however, that this effect vanishes when employing a long-wave approximation as
a final step of the derivation, as we have done here. However, the reactive
surfactant dynamics still couples to the liquid bulk within the drop via
solutal Marangoni forces. Focusing on the chemostatted out-of-equilibrium
case, we have used numerical path continuation and time simulations to
determine the rather complex bifurcation structure that shows how increasingly
complex modes of motion emerge with increasing chemostatic driving strength.
Described examples include steady localized Turing-like patterns of
surfactants, self-excited breathing and swaying modes of standing wave-like
drop oscillations as well as stationary and time-periodic self-propelled drop
motion. The latter represents a modulated traveling wave-like behavior and the
observed periodic emergence and disappearance of advancing protrusions might
be seen as a primitive model of cell propulsion where time-periodic protrusion
dynamics can also be observed [136]. Further, it is remarkable that the
bifurcation structure shows features that are quite similar to the snaking of
localized states described for the Swift-Hohenberg equation [137] and its
massconserving equivalent, the phase-field-crystal (PFC) model [138]. However,
in contrast to these cases, here, the localization is caused by the finite
drop extension that limits the region where Turing patterns can develop. This
should be further scrutinized in the future.

The presented general framework for the gradient dynamics-based description of
reactive thin-film hydrodynamics can be easily adapted to many related
situations and systems by incorporating other contributions to the energy such
as interaction terms between surfactants (e.g., to account for surfactant
phase transitions [139]) or by adapting the mobilities to account for other
transport channels. In the case of our second example, one could introduce
film height-dependent rate constants [mobilities] in the nonconserved
[conserved] part of the dynamics (62)-(64) to suppress surfactant reactions
[diffusion] in the ultrathin adsorption layer outside the drop. Finally,
although we have restricted ourselves to mesoscopic hydrodynamic models for
the description of sessile drops, we stress that the presented gradient
dynamics approach may also be employed in many other systems where chemical
reactions and fluid flows interact.

## Statements and Declarations

### Acknowledgements

FV and UT thank Yutaka Sumino for valuable discussions on reactive wetting, FV
further thanks Jan Diekmann and Christopher Henkel for helpful advice
regarding the numerics shown in Section 4.1.3.

### Author contributions

Conceptualization: FV, UT; Methodology: FV, UT; Formal analysis: FV; Software:
FV; Investigation: FV; Visualization: FV; Interpretation: FV, UT; Writing -
original draft preparation: FV, UT; Writing - review and editing: FV, UT;
Supervision: UT; Funding acquisition: UT; Project administration: UT;

### Competing interests

The authors declare no competing interests.

### Data Availability Statement

The numerical data underlying this study and the python scripts used to
present them will be made publicly available at the data repository Zenodo.

##  Appendix A Appendix

###  A.1 Gradient dynamics form of the reduced two-field model for
$\mu_{a,0}=0$

Here we show that the reduced model (51)-(52) is a gradient dynamics for
$\mu_{a,0}=0$ as stated in Section 4.1.1. First, we introduce the (reduced)
free energy

$F=\int_{\Omega}\left[f(h,b)+\frac{\gamma}{2}|\nabla
h|^{2}+g_{b}(b)+\frac{\sigma_{b}}{2}|\nabla b|^{2}\right]\text{d}^{2}r$ (78)

and use it to obtain the liquid pressure and the chemical potential of $B$ as

$\displaystyle p=\frac{\delta F}{\delta h}$
$\displaystyle=\partial_{h}f-\gamma\Delta h,$ (79)
$\displaystyle\mu_{b}=\frac{\delta F}{\delta b}$
$\displaystyle=g_{b}^{\prime}+\partial_{b}f-\sigma_{b}\Delta b.$ (80)

We can then re-express (51)-(52) as

$\displaystyle\partial_{t}h$
$\displaystyle=\nabla\cdot\left[Q_{hh}\nabla\frac{\delta F}{\delta h}\right],$
(81) $\displaystyle\partial_{t}b$
$\displaystyle=\nabla\cdot\left[Q_{bb}\nabla\frac{\delta F}{\delta
b}\right]-J_{r},$ (82)

where

$\mathbf{\underline{Q}}=\left(\begin{array}[]{cc}Q_{hh}&0\\\
0&Q_{bb}\end{array}\right)=\left(\begin{array}[]{cc}\frac{h^{3}}{3\eta}&0\\\
0&D_{b}b\end{array}\right)$ (83)

is again a symmetric, positive (semi-)definite matrix. Therefore, the
transport processes only lead to a reduction of the free energy. For the
chemical reaction with $\mu_{a,0}=0$, we have the reactive flux

$J_{r}=r(h)\left[\exp\left(\frac{\delta F/\delta
b}{k_{b}T}\right)-1\right]=r(h)\left[X-1\right],$ (84)

where $X=\exp\left(\frac{\delta F/\delta b}{k_{b}T}\right)$. We now show that
this flux is purely dissipative. To this end, we consider the change in $F$
due to reactive processes

$\displaystyle\left(\frac{\text{d}F}{\text{d}t}\right)_{\text{react}}$
$\displaystyle=\int_{\Omega}\frac{\delta F}{\delta
b}\left(-J_{r}\right)\text{d}^{2}r$ (85)
$\displaystyle=-k_{b}T\int_{\Omega}r(h)\ln X\left(X-1\right)\text{d}^{2}r$
(86) $\displaystyle\leq 0.$ (87)

Here, we have only used the definitions of $J_{r}$ and $X$. The final
inequality follows from the inequality $\ln X\cdot(X-1)\geq 0$. We note that
this inequality does not hold for the more general expression $\ln
X\cdot(X-X_{0})$, i.e., for $\mu_{a,0}\neq 0$. Therefore, as thermodynamics
suggests, Eqs. (51)-(52) are (only) a gradient dynamics for $\mu_{a,0}=0$.

###  A.2 Computation of energetic influx and dissipation for $\mu_{a,0}\neq 0$

In Section 4.1.3 we have discussed the rates of dissipation and of energetic
influx from the chemostat for $\mu_{a,0}\neq 0$ in the reduced model
(51)-(52). Here, we briefly show how to compute these quantities by comparison
with the passive case $\mu_{a,0}=0$. To this end, we compute the total change
in free energy due to reactions for $\mu_{a,0}\neq 0$

$\displaystyle\left(\frac{\text{d}F}{\text{d}t}\right)_{\text{react}}$
$\displaystyle=\int_{\Omega}\frac{\delta F}{\delta
b}\left(-J_{r}\right)\text{d}^{2}r$ (88)
$\displaystyle=-k_{b}T\int_{\Omega}r(h)\ln
X\left[X-\exp\left(\frac{\mu_{a,0}}{k_{b}T}\right)\right]\text{d}^{2}r$ (89)
$\displaystyle=\underbrace{-k_{b}T\int_{\Omega}r(h)\ln
X\left[X-1\right]\text{d}^{2}r}_{\leq
0}+\underbrace{k_{b}T\int_{\Omega}r(h)\ln
X\left[\exp\left(\frac{\mu_{a,0}}{k_{b}T}\right)-1\right]\text{d}^{2}r}_{=R_{\text{chem}}}.$
(90)

Here, we have again used $X=\exp\left(\frac{\delta F/_{\text{red}}\delta
b}{k_{b}T}\right)$ and the definition (84) for $J_{r}$. In line (90), we have
separated the purely dissipative contributions which we know from the passive
case in Appendix A.1. We then read off

$R_{\text{chem}}=k_{b}T\int_{\Omega}r(h)\ln
X\left[\exp\left(\frac{\mu_{a,0}}{k_{b}T}-1\right)\right]\text{d}^{2}r.$ (91)

Note that $R_{\text{chem}}$ may be positive or negative, i.e., the chemostat
may be an energetic source or sink. If one chooses $\mu_{a,0}>0$, however, we
usually find $R_{\text{chem}}>0$ after brief initial transients since then
there is an energetic cost associated with replacing particles of species $a$
that adsorb onto the substrate. We now obtain the dissipation rate
$R_{\text{diss}}$ from $\text{d}F/\text{d}t=-R_{\text{diss}}+R_{\text{chem}}$,
where $\text{d}F/\text{d}t$ is the total derivative of $F$. Note that we have
$R_{\text{diss}}\geq 0$.

###  A.3 Details on nondimensionalization

####  A.3.1 Drops driven by chemically sustained wettability gradients

We here briefly describe the derivation of the nondimensional version of Eqs.
(51)-(52). For the nondimensionalization, we introduce the scales

$(x,y)=L(\tilde{x},\tilde{y})\hskip 14.22636pth=l\tilde{h}\hskip
14.22636ptt=\tau\tilde{t}\hskip 14.22636ptb=b_{0}\tilde{b},$ (92)

with the scaling factors

$L=l\sqrt{\frac{\gamma}{k_{b}Tb_{0}}}\hskip 14.22636ptl=h_{a}\hskip
14.22636pt\tau=\frac{L^{2}\eta}{h_{a}k_{b}Tb_{0}}.$ (93)

Dimensionless quantities are denoted by tildes. The dimensionless time
evolution equations then take the form

$\displaystyle\partial_{\tilde{t}}\tilde{h}$
$\displaystyle=\tilde{\nabla}\cdot\left[\frac{\tilde{h}^{3}}{3}\tilde{\nabla}\tilde{p}\right],$
(94) $\displaystyle\partial_{\tilde{t}}\tilde{b}$
$\displaystyle=\tilde{\nabla}\cdot\left[\tilde{D}_{b}\tilde{b}\tilde{\nabla}\tilde{\mu}_{b}\right]-\tilde{J}_{r},$
(95)

with the nondimensional pressure and chemical potentials given as

$\displaystyle\tilde{p}$
$\displaystyle=W\cdot\left(1+\lambda\tilde{b}\right)\cdot\left(\frac{1}{\tilde{h}^{3}}-\frac{1}{\tilde{h}^{6}}\right)-\tilde{\Delta}\tilde{h},$
(96) $\displaystyle\tilde{\mu}_{b}$ $\displaystyle=\ln\tilde{b}+\lambda
W\left(-\frac{1}{2\tilde{h}^{2}}+\frac{1}{5\tilde{h}^{5}}\right)-\tilde{\sigma}_{b}\tilde{\Delta}\tilde{b}.$
(97)

and with

$\tilde{J}_{r}=\tilde{r}(\tilde{h})\left[\exp\left(\tilde{\mu}_{b}\right)-\exp\left(\tilde{\mu}_{a,0}\right)\right],$
(98)

where the reaction rate is given by

$\tilde{r}(\tilde{h})=\tilde{r}_{0}\xi(\tilde{h}).$ (99)

The remaining dimensionless quantities are

$\displaystyle W$ $\displaystyle=\frac{A_{H}}{h_{a}^{2}b_{0}k_{b}T},$ (100)
$\displaystyle\tilde{\sigma}_{b}$
$\displaystyle=\frac{b_{0}}{L^{2}k_{b}T}\sigma_{b},$ (101)
$\displaystyle\tilde{D}_{b}$ $\displaystyle=\frac{\tau k_{b}T}{L^{2}}D_{b},$
(102) $\displaystyle\tilde{r}_{0}$ $\displaystyle=\frac{\tau}{b_{0}}r_{0},$
(103) $\displaystyle\tilde{\mu}_{a,0}$
$\displaystyle=\frac{1}{k_{b}T}\mu_{a,0},$ (104)

where the wettability parameter $W$ describes the ratio of wetting energy to
entropic substrate adsorbate contributions. In the main text, tildes denoting
nondimensional quantities are omitted.

####  A.3.2 Sessile drops covered by insoluble autocatalytic surfactants

We here briefly summarize the nondimensionalization of Eqs. (74)-(76). To this
end, we introduce the scales

$t=\tau\tilde{t},\hskip 28.45274pt(x,y)=L(\tilde{x},\tilde{y}),\hskip
28.45274pth=l\tilde{h},\hskip
28.45274pt\Gamma_{\alpha}=a_{1}a_{2}\tilde{\Gamma}_{\alpha}\hskip
28.45274pt(f,g_{s})=\kappa(\tilde{f},\tilde{g}_{s}),$ (105)

where dimensionless quantities are denoted by tildes. The respective scales
are chosen as

$\tau=\frac{L^{2}\eta}{\kappa l},\hskip
28.45274ptL=\sqrt{\frac{\gamma^{0}}{\kappa}}l,\hskip 28.45274ptl=h_{a},\hskip
28.45274pt\kappa=\frac{k_{b}T}{a_{1}a_{2}}.$ (106)

Note that for this particular choice of scales the long-wave approximation
$L\gg l$ simplifies to $\gamma^{0}\gg\tfrac{k_{b}T}{a_{1}a_{2}}$. Therefore,
neglecting surfactant contributions to the capillary pressure, as done in
(73), is consistent with the present scaling. The nondimensional equations are
now given by

$\displaystyle\partial_{\tilde{t}}\tilde{h}=$
$\displaystyle\>\tilde{\nabla}\cdot\left[\frac{\tilde{h}^{3}}{3}\tilde{\nabla}\left[W\left(\frac{1}{\tilde{h}^{3}}-\frac{1}{\tilde{h}^{6}}\right)-\tilde{\Delta}\tilde{h}\right]\right]+\tilde{\nabla}\cdot\left[\frac{\tilde{h}^{2}}{2}\tilde{\nabla}\tilde{\Gamma}_{1}\right]+\tilde{\nabla}\cdot\left[\frac{\tilde{h}^{2}}{2}\tilde{\nabla}\tilde{\Gamma}_{2}\right],$
(107)
$\displaystyle\begin{split}\partial_{\tilde{t}}\tilde{\Gamma}_{1}=&\>\tilde{\nabla}\cdot\left[\frac{\tilde{h}^{2}\tilde{\Gamma}_{1}}{2}\tilde{\nabla}\left[W\left(\frac{1}{\tilde{h}^{3}}-\frac{1}{\tilde{h}^{6}}\right)-\tilde{\Delta}\tilde{h}\right]\right]+\tilde{\nabla}\cdot\left[\left(\tilde{h}\tilde{\Gamma}_{1}+\tilde{D}_{1}\right)\tilde{\nabla}\tilde{\Gamma}_{1}\right]+\tilde{\nabla}\cdot\left[\tilde{h}\tilde{\Gamma}_{1}\tilde{\nabla}\tilde{\Gamma}_{2}\right]\\\
&+\tilde{r}\left[\delta\tilde{\Gamma}_{1}^{2}\tilde{\Gamma}_{2}-\delta^{3}\tilde{\Gamma}_{1}^{3}\right]-\tilde{\beta}_{1}\ln(\tilde{\Gamma}_{1}\delta)-\tilde{\beta\mu},\end{split}$
(108)
$\displaystyle\begin{split}\partial_{\tilde{t}}\tilde{\Gamma}_{2}=&\>\tilde{\nabla}\cdot\left[\frac{\tilde{h}^{2}\tilde{\Gamma}_{2}}{2}\tilde{\nabla}\left[W\left(\frac{1}{\tilde{h}^{3}}-\frac{1}{\tilde{h}^{6}}\right)-\tilde{\Delta}\tilde{h}\right]\right]+\tilde{\nabla}\cdot\left[\tilde{h}\tilde{\Gamma}_{2}\tilde{\nabla}\tilde{\Gamma}_{1}\right]+\tilde{\nabla}\cdot\left[\left(\tilde{h}\tilde{\Gamma}_{2}+\tilde{D}_{2}\right)\tilde{\Gamma}_{2}\right]\\\
&-\tilde{r}\left[\delta\tilde{\Gamma}_{1}^{2}\tilde{\Gamma}_{2}-\delta^{3}\tilde{\Gamma}_{1}^{3}\right]-\tilde{\beta}_{2}\ln(\tilde{\Gamma}_{2}\delta^{-1})+\tilde{\beta\mu}.\end{split}$
(109)

Here, the dimensionless wettability parameter

$W=\frac{A_{H}a_{1}a_{2}}{h_{a}^{2}k_{b}T}$ (110)

expresses the ratio of wetting energy to surfactant energy. The remaining
dimensionless quantities are

$\displaystyle\delta$ $\displaystyle=\frac{a_{1}}{a_{2}},$ (111)
$\displaystyle\tilde{D}_{\alpha}$ $\displaystyle=\frac{\eta
a_{1}a_{2}}{h_{a}}D_{\alpha},$ (112) $\displaystyle\tilde{r}$
$\displaystyle=\tau
a_{1}a_{2}r=\frac{\gamma^{0}h_{a}\eta(a_{1}a_{2})^{3}}{(k_{b}T)^{2}}r,$ (113)
$\displaystyle\tilde{\beta}_{\alpha}$ $\displaystyle=\tau
a_{1}a_{2}k_{b}T\beta_{\alpha}=\frac{\gamma^{0}h_{a}\eta(a_{1}a_{2})^{3}}{k_{b}T}\beta_{\alpha},$
(114) $\displaystyle\tilde{\beta\mu}$ $\displaystyle=\tau
a_{1}a_{2}\beta\mu=\frac{\gamma^{0}h_{a}\eta(a_{1}a_{2})^{3}}{(k_{b}T)^{2}}\beta\mu.$
(115)

In the main text, we omit tildes denoting dimensionless quantities.

##  Appendix B Numerical methods

Here, we briefly describe the numerical methods that were employed to obtain
the results presented in Sections 4.1.3 and 4.2.2.

###  B.1 Time simulations

For time simulations, we use the finite-element library oomph-lib [118] on a
1D periodic domain with spatial adaptivity. Typical numbers of elements are
$n_{x}=5000$ (Section 4.1.3) and $n_{x}=1000$ (Section 4.2.2). Time-stepping
is employed using a second-order backward differentiation scheme (BDF[2]) and
is usually adaptive. An exception is the data shown in Figs. 3a,b which was
obtained with a constant timestep of $\text{d}t=5$ to produce equidistant
spacing in the figures.

###  B.2 Path continuation

For numerical continuation, we discretize the partial differential equations
(PDEs) using the finite-element method on a 1D periodic domain and employ the
continuation package pde2path [126] which uses pseudo-arclength continuation
with a predictor-corrector method to follow solution branches. In particular,
pde2path allows for the treatment of PDEs with continuous symmetries such as
volume conservation or translational invariance, which are both present in
Eqs. (74)-(76) with periodic boundaries. We use static meshes with a typical
number of points $n_{x}=300$ and $n_{t}=40$ in space and time, respectively.
For the continuation of time-periodic solutions, we do not employ branch-
switching at Hopf bifurcations for technical reasons. Instead, we interpolate
time-simulation data onto the static mesh and use it as an initial guess to
start continuation from a point on the branch of time-periodic states.

## References

  * Girichidis et al. [2020] P. Girichidis, S. S. R. Offner, A. G. Kritsuk, R. S. Klessen, P. Hennebelle, J. M. D. Kruijssen, M. G. H. Krause, S. C. O. Glover, and M. Padovani.  Physical processes in star formation.  _Space Sci. Rev._ , 216, 2020.  doi: 10.1007/s11214-020-00693-8.
  * De Wit [2020] A. De Wit.  Chemo-hydrodynamic patterns and instabilities.  _Annu. Rev. Fluid Mech._ , 52:531–555, 2020.  doi: 10.1146/annurev-fluid-010719-060349.
  * Liberman [2010] M. A. Liberman.  _Introduction to Physics and Chemistry of Combustion_.  Springer Berlin Heidelberg, Berlin, Heidelberg, 2010.
  * Donlan [2002] R. M. Donlan.  Biofilms: Microbial life on surfaces.  _Emerg. Infect. Dis._ , 8:881–890, 2002.  doi: 10.3201/eid0809.020063.
  * Trinschek et al. [2017] S. Trinschek, K. John, S. Lecuyer, and U. Thiele.  Continuous vs. arrested spreading of biofilms at solid-gas interfaces \- the role of surface forces.  _Phys. Rev. Lett._ , 119:078003, 2017.  doi: 10.1103/PhysRevLett.119.078003.
  * Nakashima et al. [2019] K. K. Nakashima, M. A. Vibhute, and E. Spruijt.  Biomolecular chemistry in liquid phase separated compartments.  _Front. Mol. Biosci._ , 6, 2019.  doi: 10.3389/fmolb.2019.00021.
  * Burkart et al. [2022] T. Burkart, M. C. Wigbers, L. Wurthner, and E. Frey.  Control of protein-based pattern formation via guiding cues.  _Nat. Rev. Phys._ , 2022.  doi: 10.1038/s42254-022-00461-3.
  * Turing [1952] A. M. Turing.  The chemical basis of morphogenesis.  _Philos. Trans. R. Soc. Lond. Ser. B-Biol. Sci._ , 237:37–72, 1952.  doi: 10.1098/rstb.1952.0012.
  * de Gennes [1985] P. G. de Gennes.  Wetting: Statics and dynamics.  _Rev. Mod. Phys._ , 57:827–863, 1985.  doi: 10.1103/RevModPhys.57.827.
  * Bonn et al. [2009] D. Bonn, J. Eggers, J. Indekeu, J. Meunier, and E. Rolley.  Wetting and spreading.  _Rev. Mod. Phys._ , 81:739–805, 2009.  doi: 10.1103/RevModPhys.81.739.
  * Craster and Matar [2009] R. V. Craster and O. K. Matar.  Dynamics and stability of thin liquid films.  _Rev. Mod. Phys._ , 81:1131–1198, 2009.  doi: 10.1103/RevModPhys.81.1131.
  * Pismen [1984] L. M. Pismen.  Composition and flow patterns due to chemo-Marangoni instability in liquid films.  _J. Colloid Interface Sci._ , 102:237–247, 1984.  doi: 10.1016/0021-9797(84)90216-9.
  * Dagan and Pismen [1984] Z. Dagan and L. M. Pismen.  Marangoni waves induced by a multistable chemical reaction on thin liquid films.  _J. Colloid Interface Sci._ , 99:215–225, 1984.  doi: 10.1016/0021-9797(84)90103-6.
  * Bain and Whitesides [1989] C. D. Bain and G. M. Whitesides.  A study by contact-angle of the acid-base behavior of monolayers containing $\omega$-mercaptocarboxylic acids adsorbed on gold: An example of reactive spreading.  _Langmuir_ , 5:1370–1378, 1989.  doi: 10.1021/la00090a019.
  * Meinkohn and Mikhailov [1993] D. Meinkohn and A. Mikhailov.  Motion of thin wetting liquid films in the presence of surface chemical reactions.  _Physica A_ , 198:25–45, 1993.  doi: 10.1016/0378-4371(93)90182-4.
  * Bain et al. [1994] C. D. Bain, G. D. Burnetthall, and R. R. Montgomerie.  Rapid motion of liquid-drops.  _Nature_ , 372:414–415, 1994.  doi: 10.1038/372414a0.
  * de Gennes [1997] P. G. de Gennes.  Forced wetting by a reactive fluid.  _Europhys. Lett._ , 39:407–412, 1997.
  * Zheng et al. [1998] D. W. Zheng, W. Wen, and K. N. Tu.  Reactive wetting- and dewetting-induced diffusion-limited aggregation.  _Phys. Rev. E_ , 57:R3719–R3722, 1998.
  * Saiz et al. [2000] E. Saiz, R. M. Cannon, and A. P. Tomsia.  Reactive spreading: Adsorption, ridging and compound formation.  _Acta Materialia_ , 48:4449–4462, 2000.
  * van Nierop et al. [2006] E. A. van Nierop, A. Ajdari, and H. A. Stone.  Reactive spreading and recoil of oil on water.  _Phys. Fluids_ , 18:038105, 2006.  doi: 10.1063/1.2187068.
  * Hanczyc [2011] M. M. Hanczyc.  Metabolism and motility in prebiotic structures.  _Philos. Trans. R. Soc. B-Biol. Sci._ , 366:2885–2893, 2011\.  doi: 10.1098/rstb.2011.0141.
  * Kumar and Prabhu [2007] G. Kumar and K. N. Prabhu.  Review of non-reactive and reactive wetting of liquids on surfaces.  _Adv. Colloid Interface Sci._ , 133:61–89, 2007.  ISSN 0001-8686.  doi: 10.1016/j.cis.2007.04.009.
  * Shioi et al. [2010] A. Shioi, T. Ban, and Y. Morimune.  Autonomously moving colloidal objects that resemble living matter.  _Entropy_ , 12:2308–2332, 2010.  doi: 10.3390/e12112308.
  * Eustathopoulos and Voytovych [2016] N. Eustathopoulos and R. Voytovych.  The role of reactivity in wetting by liquid metals: A review.  _J. Mater. Sci._ , 51:425–437, 2016.  doi: 10.1007/s10853-015-9331-3.
  * Snoeijer and Andreotti [2013] J. H. Snoeijer and B. Andreotti.  Moving contact lines: Scales, regimes, and dynamical transitions.  _Annu. Rev. Fluid Mech._ , 45:269–292, 2013.  doi: 10.1146/annurev-fluid-011212-140734.
  * Koch et al. [2008] K. Koch, B. Bhushan, and W. Barthlott.  Diversity of structure, morphology and wetting of plant surfaces.  _Soft Matter_ , 4:1943–1963, 2008.  doi: 10.1039/b804854a.
  * Neogi [1985] P. Neogi.  Tears-of-wine and related phenomena.  _J. Colloid Interface Sci._ , 105:94–101, 1985.  doi: 10.1016/0021-9797(85)90350-9.
  * Domingues Dos Santos and Ondarçuhu [1995] F. Domingues Dos Santos and T. Ondarçuhu.  Free-running droplets.  _Phys. Rev. Lett._ , 75:2972–2975, 1995.
  * Brochard-Wyart and de Gennes [1995] F. Brochard-Wyart and P.-G. de Gennes.  Spontaneous motion of a reactive droplet.  _C. R. Acad. Sci. Ser. II_ , 321:285–288, 1995.
  * de Gennes [1998] P. G. de Gennes.  The dynamics of reactive wetting on solid surfaces.  _Physica A_ , 249:196–205, 1998.
  * Lee and Laibinis [2000] S. W. Lee and P. E. Laibinis.  Directed movement of liquids on patterned surfaces using noncovalent molecular adsorption.  _J. Am. Chem. Soc._ , 122:5395–5396, 2000.
  * Lee et al. [2002] S. W. Lee, D. Y. Kwok, and P. E. Laibinis.  Chemical influences on adsorption-mediated self-propelled drop movement.  _Phys. Rev. E_ , 65:051602, 2002.  doi: 10.1103/PhysRevE.65.051602.
  * Thiele et al. [2004] U. Thiele, K. John, and M. Bär.  Dynamical model for chemically driven running droplets.  _Phys. Rev. Lett._ , 93:027802, 2004.
  * Sumino et al. [2005] Y. Sumino, H. Kitahata, K. Yoshikawa, M. Nagayama, S. M. Nomura, N. Magome, and Y. Mori.  Chemosensitive running droplet.  _Phys. Rev. E_ , 72:041603, 2005.  doi: 10.1103/PhysRevE.72.041603.
  * John et al. [2005] K. John, M. Bär, and U. Thiele.  Self-propelled running droplets on solid substrates driven by chemical reactions.  _Eur. Phys. J. E_ , 18:183–199, 2005.  doi: 10.1140/epje/i2005-10039-1.
  * Sumino et al. [2006] Y. Sumino, N. Magome, and K. Yoshikawa.  An oil droplet that spontaneously climbs up stairs.  _Prog. Theor. Phys. Suppl._ , 161:348–351, 2006.  doi: 10.1143/ptps.161.348.
  * Kitahata [2006] H. Kitahata.  Spontaneous motion of a droplet coupled with chemical reaction.  _Prog. Theor. Phys. Suppl._ , 161:220–223, 2006.
  * Yao et al. [2012] X. Yao, H. Bai, J. Ju, D. Zhou, J. Li, H. Zhang, B. Yang, and L. Jiang.  Running droplet of interfacial chemical reaction flow.  _Soft Matter_ , 8:5988–5991, 2012.  doi: 10.1039/c2sm25153a.
  * Arscott [2016] S. Arscott.  Dynamic chemically driven dewetting, spreading, and self-running of sessile droplets on crystalline silicon.  _Langmuir_ , 32:12611–12622, 2016.  doi: 10.1021/acs.langmuir.6b03287.
  * Diewald et al. [1996] M. Diewald, K. Matthiessen, S. C. Müller, and H. R. Brand.  Oscillatory hydrodynamic flow due to concentration dependence of surface tension.  _Phys. Rev. Lett._ , 77:4466–4469, 1996.
  * Maass et al. [2016] C. C. Maass, C. Krüger, S. Herminghaus, and C. Bahr.  Swimming droplets.  _Annu. Rev. Condens. Matter Phys._ , 7:171–193, 2016.  doi: 10.1146/annurev-conmatphys-031115-011517.
  * Michelin [2023] S. Michelin.  Self-propulsion of chemically active droplets.  _Annu. Rev. Fluid Mech._ , 55:77–101, 2023.  doi: 10.1146/annurev-fluid-120720-012204.
  * Pismen [1997] L. M. Pismen.  Interaction of reaction-diffusion fronts and Marangoni flow on the interface of a deep fluid.  _Phys. Rev. Lett._ , 78:382–385, 1997.
  * Pereira et al. [2007] A. Pereira, P. M. J. Trevelyan, U. Thiele, and S. Kalliadasis.  Dynamics of a horizontal thin liquid film in the presence of reactive surfactants.  _Phys. Fluids_ , 19:112102, 2007.  doi: 10.1063/1.2775938.
  * Riolfo et al. [2012] L. A. Riolfo, J. Carballido-Landeira, C. O. Bounds, J. A. Pojman, S. Kalliadasis, and A. De Wit.  Experimental reaction-driven liquid film fingering instability.  _Chem. Phys. Lett._ , 534:13–18, 2012.  doi: 10.1016/j.cplett.2012.03.024.
  * Nakada et al. [2016] M. Nakada, Y. Fujikami, M. Kawaguchi, D. Yamamoto, and A. Shioi.  A molecular assembly that crawls on a solid substrate with a metabolic-like process.  _Mol. Syst. Des. Eng._ , 1:208–215, 2016.  doi: 10.1039/c5–m˝e00012–b˝.
  * Braun et al. [1995] R. J. Braun, B. T. Murray, W. J. Boettinger, and G. B. McFadden.  Lubrication theory for reactive spreading of a thin drop.  _Phys. Fluids_ , 7:1797–1810, 1995.
  * Banani et al. [2017] S. F. Banani, H. O. Lee, A. A. Hyman, and M. K. Rosen.  Biomolecular condensates: organizers of cellular biochemistry.  _Nat. Rev. Mol. Cell Biol._ , 18:285–298, 2017.  doi: 10.1038/nrm.2017.7.
  * Gouveia et al. [2022] B. Gouveia, Y. Kim, J. W. Shaevitz, S. Petry, H. A. Stone, and C. P. Brangwynne.  Capillary forces generated by biomolecular condensates.  _Nature_ , 609:255–264, 2022.  doi: 10.1038/s41586-022-05138-6.
  * Mangiarotti et al. [2023] A. Mangiarotti, N. N. Chen, Z. L. Zhao, R. Lipowsky, and R. Dimova.  Wetting and complex remodeling of membranes by biomolecular condensates.  _Nat. Commun._ , 14:2809, 2023.  doi: 10.1038/s41467-023-37955-2.
  * Mitlin [1993] V. S. Mitlin.  Dewetting of solid surface: Analogy with spinodal decomposition.  _J. Colloid Interface Sci._ , 156:491–497, 1993.  doi: 10.1006/jcis.1993.1142.
  * Thiele et al. [2016] U. Thiele, A. J. Archer, and L. M. Pismen.  Gradient dynamics models for liquid films with soluble surfactant.  _Phys. Rev. Fluids_ , 1:083903, 2016.  doi: 10.1103/PhysRevFluids.1.083903.
  * Hartmann et al. [2023] S. Hartmann, C. Diddens, M. Jalaal, and U. Thiele.  Sessile drop evaporation in a gap - crossover between diffusion-limited and phase transition-limited regime.  _J. Fluid Mech._ , 960:A32, 2023.  doi: 10.1017/jfm.2023.176.
  * Thiele [2018] U. Thiele.  Recent advances in and future challenges for mesoscopic hydrodynamic modelling of complex wetting.  _Colloid Surf. A_ , 553:487–495, 2018.  doi: 10.1016/j.colsurfa.2018.05.049.
  * Hartmann et al. [2024] S. Hartmann, J. Diekmann, D. Greve, and U. Thiele.  Drops on polymer brushes – advances in thin-film modelling of adaptive substrates.  _Langmuir_ , x:x, 2024.  doi: 10.1021/acs.langmuir.3c03313.
  * Miller [1960] D. G. Miller.  Thermodynamics of irreversible processes. the experimental verification of the onsager reciprocal relations.  _Chem. Rev._ , 60:15–37, 1960.  doi: 10.1021/cr60203a003.
  * Rahman and Saghir [2014] M.A. Rahman and M.Z. Saghir.  Thermodiffusion or soret effect: Historical review.  _Int. J. Heat Mass Transf._ , 73:693–705, 2014.  doi: 10.1016/j.ijheatmasstransfer.2014.02.057.
  * de Groot and Mazur [1984] S. R. de Groot and P. Mazur.  _Non-equilibrium Thermodynamics_.  Dover publications, New York, 1984.
  * Onsager [1931a] L. Onsager.  Reciprocal relations in irreversible processes. I.  _Phys. Rev._ , 37:405–426, 1931a.  doi: 10.1103/PhysRev.37.405.
  * Onsager [1931b] L. Onsager.  Reciprocal relations in irreversible processes. II.  _Phys. Rev._ , 38:2265–2279, 1931b.  doi: 10.1103/PhysRev.38.2265.
  * Marcelin [1914] M. R. Marcelin.  _Contribution a l’etude de la cinétique physico-chimique_.  PhD thesis, Faculté des Sciences de Paris, 1914.
  * De Donder [1936] T. De Donder.  _Thermodynamic Theory of Affinity: A Book of Principles_.  Oxford University Press, Oxford, 1936.
  * Zwicker [2022] D. Zwicker.  The intertwined physics of active chemical reactions and phase separation.  _Curr. Opin. Colloid Interface Sci._ , 61:101606, 2022.  doi: 10.1016/j.cocis.2022.101606.
  * Doi [2011] Masao Doi.  Onsager’s variational principle in soft matter.  _J. Phys.: Condens. Matter_ , 23:284118, 2011.  doi: 10.1088/0953-8984/23/28/284118.
  * Oron et al. [1997] A. Oron, S. H. Davis, and S. G. Bankoff.  Long-scale evolution of thin liquid films.  _Rev. Mod. Phys._ , 69:931–980, 1997.  doi: 10.1103/RevModPhys.69.931.
  * Thiele [2007] U. Thiele.  Structure formation in thin liquid films.  In S. Kalliadasis and U. Thiele, editors, _Thin Films of Soft Matter_ , pages 25–93. Springer Vienna, Vienna, 2007.  doi: 10.1007/978-3-211-69808-2“˙2.
  * Starov and Velarde [2009] V. M. Starov and M. G. Velarde.  Surface forces and wetting phenomena.  _J. Phys.-Condens. Matter_ , 21:464121, 2009.  doi: 10.1088/0953-8984/21/46/464121.
  * Thiele [2010] U. Thiele.  Thin film evolution equations from (evaporating) dewetting liquid layers to epitaxial growth.  _J. Phys.: Condens. Matter_ , 22:084019, 2010.  doi: 10.1088/0953-8984/22/8/084019.
  * Miksis and Davis [1994] M. J. Miksis and S. H. Davis.  Slip over rough and coated surfaces.  _J. Fluid Mech._ , 273:125–139, 1994.  doi: 10.1017/s0022112094001874.
  * Münch et al. [2005] A. Münch, B. Wagner, and T. P. Witelski.  Lubrication models with small to large slip lengths.  _J. Eng. Math._ , 53:359–383, 2005.  doi: 10.1007/s10665-005-9020-3.
  * Cahn and Hilliard [1958] J. W. Cahn and J. E. Hilliard.  Free energy of a nonuniform system. 1. Interfacial free energy.  _J. Chem. Phys._ , 28:258–267, 1958.  doi: 10.1063/1.1744102.
  * Bray [1994] A. J. Bray.  Theory of phase-ordering kinetics.  _Adv. Phys._ , 43:357–459, 1994.  doi: 10.1080/00018739400101505.
  * Pototsky et al. [2005] A. Pototsky, M. Bestehorn, D. Merkt, and U. Thiele.  Morphology changes in the evolution of liquid two-layer films.  _J. Chem. Phys._ , 122:224711, 2005.  doi: 10.1063/1.1927512.
  * Thiele et al. [2012] U. Thiele, A. J. Archer, and M. Plapp.  Thermodynamically consistent description of the hydrodynamics of free surfaces covered by insoluble surfactants of high concentration.  _Phys. Fluids_ , 24:102107, 2012.  doi: 10.1063/1.4758476.
  * Thiele et al. [2013a] U. Thiele, D. V. Todorova, and H. Lopez.  Gradient dynamics description for films of mixtures and suspensions: dewetting triggered by coupled film height and concentration fluctuations.  _Phys. Rev. Lett._ , 111:117801, 2013a.  doi: 10.1103/PhysRevLett.111.117801.
  * Huth et al. [2015] R. Huth, S. Jachalski, G. Kitavtsev, and D. Peschka.  Gradient flow perspective on thin-film bilayer flows.  _J. Eng. Math._ , 94:43–61, 2015.  doi: 10.1007/s10665-014-9698-1.
  * Thiele and Hartmann [2020] U. Thiele and S. Hartmann.  Gradient dynamics model for drops spreading on polymer brushes.  _Eur. Phys. J.-Spec. Top._ , 229:1819–1832, 2020.  doi: 10.1140/epjst/e2020-900231-2.
  * Henkel et al. [2021] C. Henkel, J. H. Snoeijer, and U. Thiele.  Gradient-dynamics model for liquid drops on elastic substrates.  _Soft Matter_ , 17:10359–10375, 2021.  doi: 10.1039/D1SM01032H.
  * Archer and Rauscher [2004] A. J. Archer and M. Rauscher.  Dynamical density functional theory for interacting brownian particles: Stochastic or deterministic?  _J. Phys. A-Math. Gen._ , 37:9325–9333, 2004.  doi: 10.1088/0305-4470/37/40/001.
  * Grmela and Öttinger [1997] M. Grmela and H. C. Öttinger.  Dynamics and thermodynamics of complex fluids. i. development of a general formalism.  _Phys. Rev. E_ , 56:6620–6632, 1997.  doi: 10.1103/PhysRevE.56.6620.
  * Öttinger and Grmela [1997] H. C. Öttinger and M. Grmela.  Dynamics and thermodynamics of complex fluids. ii. illustrations of a general formalism.  _Phys. Rev. E_ , 56:6633–6655, 1997.  doi: 10.1103/PhysRevE.56.6633.
  * Mielke [2011] A. Mielke.  A gradient structure for reaction-diffusion systems and for energy-drift-diffusion systems.  _Nonlinearity_ , 24:1329–1346, 2011.  doi: 10.1088/0951-7715/24/4/016.
  * Koudriavtsev et al. [2001] A. Koudriavtsev, R. F. Jameson, and W. Linert.  _The Law of Mass Action_.  Springer Berlin Heidelberg, 2001.  doi: 10.1007/978-3-642-56770-4.
  * Gorban [2014] A.N. Gorban.  Detailed balance in micro- and macrokinetics and micro-distinguishability of macro-processes.  _Results Phys._ , 4:142–147, 2014.  doi: https://doi.org/10.1016/j.rinp.2014.09.002.
  * Wegscheider [1902] R. Wegscheider.  Über simultane Gleichgewichte und die Beziehungen zwischen Thermodynamik und Reaktionskinetik homogener Systeme.  _Z. Phys. Chem._ , 39U:257–303, 1902.  doi: doi:10.1515/zpch-1902-3919.
  * Van Rysselberghe [1958] P. Van Rysselberghe.  Reaction rates and affinities.  _J. Chem. Phys._ , 29:640–642, 1958.  doi: 10.1063/1.1744552.
  * Pagonabarraga et al. [1997] I. Pagonabarraga, A. PerezMadrid, and J. M. Rubi.  Fluctuating hydrodynamics approach to chemical reactions.  _Physica A_ , 237:205–219, 1997.  doi: 10.1016/S0378-4371(96)00377-9.
  * Hanggi et al. [1990] P. Hanggi, P. Talkner, and M. Borkovec.  Reaction-rate theory - 50 years after kramers.  _Rev. Mod. Phys._ , 62:251–341, 1990.  doi: 10.1103/RevModPhys.62.251.
  * Peters [2017] B. Peters.  _Reaction rate theory and rare events_.  Elsevier B.V., Amsterdam, 2017.
  * Feinberg [2019] M. Feinberg.  _Foundations of Chemical Reaction Network Theory_.  Springer International Publishing, 2019.  doi: 10.1007/978-3-030-03858-8.
  * John and Bär [2005] K. John and M. Bär.  Travelling lipid domains in a dynamic model for protein-induced pattern formation in biomembranes.  _Phys. Biol._ , 2:123–132, 2005.  doi: 10.1088/1478-3975/2/2/005.
  * Ishihara et al. [2007] S. Ishihara, M. Otsuji, and A. Mochizuki.  Transient and steady state of mass-conserved reaction-diffusion systems.  _Phys. Rev. E_ , 75:015203, 2007.  doi: 10.1103/PhysRevE.75.015203.
  * Morita and Ogawa [2010] Y. Morita and T. Ogawa.  Stability and bifurcation of nonconstant solutions to a reaction-diffusion system with conservation of mass.  _Nonlinearity_ , 23:1387–1411, 2010.  doi: 10.1088/0951-7715/23/6/007.
  * Alonso and Bär [2010] S. Alonso and M. Bär.  Phase separation and bistability in a three-dimensional model for protein domain formation at biomembranes.  _Phys. Biol._ , 7:046012, 2010.  doi: 10.1088/1478-3975/7/4/046012.
  * Kuwamura and Morita [2015] M. Kuwamura and Y. Morita.  Perturbations and dynamics of reaction-diffusion systems with mass conservation.  _Phys. Rev. E_ , 92:012908, 2015.  doi: 10.1103/PhysRevE.92.012908.
  * Bernitt et al. [2017] E. Bernitt, H. G. Dobereiner, N. S. Gov, and A. Yochelis.  Fronts and waves of actin polymerization in a bistability-based mechanism of circular dorsal ruffles.  _Nat. Commun._ , 8:15863, 2017.  doi: 10.1038/ncomms15863.
  * Halatek and Frey [2018] J. Halatek and E. Frey.  Rethinking pattern formation in reaction-diffusion systems.  _Nature Phys._ , 14:507–514, 2018.  doi: 10.1038/s41567-017-0040-5.
  * Brauns et al. [2020] F. Brauns, J. Halatek, and E. Frey.  Phase-space geometry of mass-conserving reaction-diffusion dynamics.  _Phys. Rev. X_ , 10:041036, 2020.  doi: 10.1103/PhysRevX.10.041036.
  * Brauns et al. [2021] F. Brauns, H. Weyer, J. Halatek, J. Yoon, and E. Frey.  Wavelength selection by interrupted coarsening in reaction-diffusion systems.  _Phys. Rev. Lett._ , 126:104101, 2021.  doi: 10.1103/physrevlett.126.104101.
  * Gaspard [2022] P. Gaspard.  _The Statistical Mechanics of Irreversible Phenomena_.  Cambridge University Press, 2022.  doi: 10.1017/9781108563055.
  * Kirschbaum and Zwicker [2021] J. Kirschbaum and D. Zwicker.  Controlling biomolecular condensates via chemical reactions.  _J. R. Soc. Interface_ , 18, 2021.  doi: 10.1098/rsif.2021.0255.
  * Cotton et al. [2022] M. W. Cotton, R. Golestanian, and J. Agudo-Canalejo.  Catalysis-induced phase separation and autoregulation of enzymatic activity.  _Phys. Rev. Lett._ , 129:158101, 2022.  doi: 10.1103/PhysRevLett.129.158101.
  * Demarchi et al. [2023] L. Demarchi, A. Goychuk, I. Maryshev, and E. Frey.  Enzyme-enriched condensates show self-propulsion, positioning, and coexistence.  _Phys. Rev. Lett._ , 130:128401, 2023.  doi: 10.1103/PhysRevLett.130.128401.
  * Golestanian et al. [2005] R. Golestanian, T. B. Liverpool, and A. Ajdari.  Propulsion of a molecular machine by asymmetric distribution of reaction products.  _Phys. Rev. Lett._ , 94:220801, 2005.  doi: 10.1103/physrevlett.94.220801.
  * Prigogine and Lefever [1968] I. Prigogine and R. Lefever.  Symmetry breaking instabilities in dissipative systems. II.  _J. Chem. Phys._ , 48:1695–1700, 1968.  doi: 10.1063/1.1668896.
  * Nicolis [1999] G. Nicolis.  _Introduction to Nonlinear Science_.  Cambridge University Press, Cambridge, 1999.  doi: 10.1017/CBO9781139170802.
  * Horn and Jackson [1972] F. Horn and R. Jackson.  General mass action kinetics.  _Arch. Ration. Mech. Anal._ , 47:81–116, 1972.  doi: 10.1007/BF00251225.
  * Rao and Esposito [2016] R. Rao and M. Esposito.  Nonequilibrium thermodynamics of chemical reaction networks: Wisdom from stochastic thermodynamics.  _Phys. Rev. X_ , 6:041064, 2016.  doi: 10.1103/physrevx.6.041064.
  * Ziethen et al. [2023] N. Ziethen, J. Kirschbaum, and D. Zwicker.  Nucleation of chemically active droplets.  _Phys. Rev. Lett._ , 130:248201, 2023.  doi: 10.1103/PhysRevLett.130.248201.
  * Flory [1942] P. J. Flory.  Thermodynamics of High Polymer Solutions.  _J. Chem. Phys._ , 10:51–61, 1942.  doi: 10.1063/1.1723621.
  * Huggins [1941] M. L. Huggins.  Solutions of Long Chain Compounds.  _J. Chem. Phys._ , 9:440–440, 1941.  doi: 10.1063/1.1750930.
  * Spratte et al. [1994] K. Spratte, L. F. Chi, and H. Riegler.  Physisorption instabilities during dynamic Langmuir wetting.  _Europhys. Lett._ , 25:211–217, 1994.  doi: 10.1209/0295-5075/25/3/010.
  * Köpf et al. [2010] M. H. Köpf, S. V. Gurevich, R. Friedrich, and L. Chi.  Pattern formation in monolayer transfer systems with substrate-mediated condensation.  _Langmuir_ , 26:10444–10447, 2010.  doi: 10.1021/la101900z.
  * Wilczek et al. [2015] M. Wilczek, W. B. H. Tewes, S. V. Gurevich, M. H. Köpf, L. Chi, and U. Thiele.  Modelling pattern formation in dip-coating experiments.  _Math. Model. Nat. Phenom._ , 10:44–60, 2015.  doi: 10.1051/mmnp/201510402.
  * Rednikov et al. [1994] A. Y. Rednikov, Y. S. Ryazantsev, and M. G. Velarde.  Drop motion with surfactant transfer in a homogeneous surrounding.  _Phys. Fluids_ , 6:451–468, 1994.
  * Sumino and Yoshikawa [2014] Y. Sumino and K. Yoshikawa.  Amoeba-like motion of an oil droplet.  _Eur. Phys. J.-Spec. Top._ , 223:1345–1352, 2014.  doi: 10.1140/epjst/e2014-02194-x.
  * Nanzai et al. [2016] B. Nanzai, M. Kato, and M. Igawa.  Spontaneous motion of various oil droplets in aqueous solution of trimethyl alkyl ammonium with different carbon chain lengths.  _Colloid Surf. A-Physicochem. Eng. Asp._ , 504:154–160, 2016.  doi: 10.1016/j.colsurfa.2016.04.063.
  * Heil and Hazel [2006] M. Heil and A. L. Hazel.  Oomph-lib - an object-oriented multi-physics finite-element library.  In H.-J. Bungartz and M. Schäfer, editors, _Fluid-Structure Interaction: Modelling, Simulation, Optimisation_ , pages 19–49. Springer, Berlin, Heidelberg, 2006.  doi: 10.1007/3-540-34596-5˙2.
  * Buyevich et al. [1993] Yu. A. Buyevich, L. M. Rabinovich, and A. V. Vyazmin.  Chemo-Marangoni convection I. Linear analysis and criteria of instability.  _J. Colloid Interface Sci._ , 157:202–210, 1993.  doi: 10.1006/jcis.1993.1177.
  * Chaithanya et al. [2022] K. V. S. Chaithanya, S. A. Shenoy, and P. Dayal.  Hydrodynamics of a confined active Belousov-Zhabotinsky droplet.  _Phys. Rev. E_ , 106:065103, 2022.  doi: 10.1103/PhysRevE.106.065103.
  * Thiele et al. [2018] U. Thiele, J. H. Snoeijer, S. Trinschek, and K. John.  Equilibrium contact angle and adsorption layer properties with surfactants.  _Langmuir_ , 34:7210–7221, 2018.  doi: 10.1021/acs.langmuir.8b00513.  Also see Erratum: Langmuir, 35, 4788-4789 (2019), doi:10.1021/acs.langmuir.9b00616.
  * Mietke et al. [2019] A. Mietke, V. Jemseena, K. V. Kumar, I. F. Sbalzarini, and F. Jülicher.  Minimal model of cellular symmetry breaking.  _Phys. Rev. Lett._ , 123:188101, 2019.  doi: 10.1103/physrevlett.123.188101.
  * Würthner et al. [2023] L. Würthner, A. Goychuk, and E. Frey.  Geometry-induced patterns through mechanochemical coupling.  _Phys. Rev. E_ , 108:014404, 2023.  doi: 10.1103/physreve.108.014404.
  * Bois et al. [2011] J. S. Bois, F. Jülicher, and S. W. Grill.  Pattern formation in active fluids.  _Phys. Rev. Lett._ , 106:028103, 2011.  doi: 10.1103/PhysRevLett.106.028103.
  * Gross et al. [2017] P. Gross, K. V. Kumar, and S. W. Grill.  How active mechanics and regulatory biochemistry combine to form patterns in development.  _Annu. Rev. Biophys._ , 46:337–356, 2017.  doi: 10.1146/annurev-biophys-070816-033602.
  * Uecker et al. [2014] H. Uecker, D. Wetzel, and J. D. M. Rademacher.  pde2path - a Matlab package for continuation and bifurcation in 2D elliptic systems.  _Numer. Math.-Theory Methods Appl._ , 7:58–106, 2014.  doi: 10.4208/nmtma.2014.1231nm.
  * Goldstein et al. [1991] R. E. Goldstein, G. H. Gunaratne, L. Gil, and P. Coullet.  Hydrodynamic and interfacial patterns with broken space-time symmetry.  _Phys. Rev. A_ , 43:6700–6721, 1991.  doi: 10.1103/PhysRevA.43.6700.
  * Fauve et al. [1991] S. Fauve, S. Douady, and O. Thual.  Drift instabilities of cellular-patterns.  _J. Phys. II_ , 1:311–322, 1991.  doi: 10.1051/jp2:1991170.
  * Frohoff-Hülsmann et al. [2023] T. Frohoff-Hülsmann, M. P. Holl, E. Knobloch, S. V. Gurevich, and U. Thiele.  Stationary broken parity states in active matter models.  _Phys. Rev. E_ , 107:064210, 2023.  doi: 10.1103/PhysRevE.107.064210.
  * Greve et al. [2024] Daniel Greve, Giorgio Lovato, Tobias Frohoff-Hülsmann, and Uwe Thiele.  Maxwell construction for a nonreciprocal cahn-hilliard model.  _arXiv:2402.08634_ , 2024.
  * Trinschek et al. [2020] S. Trinschek, F. Stegemerten, K. John, and U. Thiele.  Thin-film modelling of resting and moving active droplets.  _Phys. Rev. E_ , 101:062802, 2020.  doi: 10.1103/PhysRevE.101.062802.
  * Saha et al. [2020] S. Saha, J. Agudo-Canalejo, and R. Golestanian.  Scalar active mixtures: The non-reciprocal Cahn-Hilliard model.  _Phys. Rev. X_ , 10:041009, 2020.  doi: 10.1103/PhysRevX.10.041009.
  * Frohoff-Hülsmann et al. [2021] T. Frohoff-Hülsmann, J. Wrembel, and U. Thiele.  Suppression of coarsening and emergence of oscillatory behavior in a Cahn-Hilliard model with nonvariational coupling.  _Phys. Rev. E_ , 103:042602, 2021.  doi: 10.1103/PhysRevE.103.042602.
  * Stegemerten et al. [2022] F. Stegemerten, K. John, and U. Thiele.  Symmetry-breaking, motion and bistability of active drops through polarization-surface coupling.  _Soft Matter_ , 18:5823–5832, 2022.  doi: 10.1039/D2SM00648K.
  * Frohoff-Hülsmann and Thiele [2023] T. Frohoff-Hülsmann and U. Thiele.  Nonreciprocal cahn-hilliard model emerges as a universal amplitude equation.  _Phys. Rev. Lett._ , 131:107201, 2023.  doi: 10.1103/PhysRevLett.131.107201.
  * Amiri et al. [2023] B. Amiri, J. C. J. Heyn, C. Schreiber, J. O. Rädler, and M. Falcke.  On multistability and constitutive relations of cell motion on fibronectin lanes.  _Biophys. J._ , 122:753–766, 2023.  doi: 10.1016/j.bpj.2023.02.001.
  * Burke and Knobloch [2006] J. Burke and E. Knobloch.  Localized states in the generalized Swift-Hohenberg equation.  _Phys. Rev. E_ , 73:056211, 2006.  doi: 10.1103/PhysRevE.73.056211.
  * Thiele et al. [2013b] U. Thiele, A. J. Archer, M. J. Robbins, H. Gomez, and E. Knobloch.  Localized states in the conserved Swift-Hohenberg equation with cubic nonlinearity.  _Phys. Rev. E_ , 87:042915, 2013b.  doi: 10.1103/PhysRevE.87.042915.
  * Köpf et al. [2012] M. H. Köpf, S. V. Gurevich, R. Friedrich, and U. Thiele.  Substrate-mediated pattern formation in monolayer transfer: a reduced model.  _New J. Phys._ , 14:023016, 2012.  doi: 10.1088/1367-2630/14/2/023016.