mathfontname=mathsf

# Self-consistent sharp interface theory of active condensate dynamics

Andriy Goychuk<EMAIL_ADDRESS>A.G. and L.D. contributed equally.  Arnold
Sommerfeld Center for Theoretical Physics and Center for NanoScience,
Department of Physics, Ludwig-Maximilians-Universität München, Theresienstraße
37, D-80333 München, Germany  Institute for Medical Engineering and Science,
Massachusetts Institute of Technology, Cambridge, MA 02139, United States
Leonardo Demarchi  A.G. and L.D. contributed equally.  Arnold Sommerfeld
Center for Theoretical Physics and Center for NanoScience, Department of
Physics, Ludwig-Maximilians-Universität München, Theresienstraße 37, D-80333
München, Germany  Sorbonne Université, CNRS, Institut de Biologie Paris-Seine
(IBPS), Laboratoire Jean Perrin (LJP), F-75005 Paris, France  Ivan Maryshev
Arnold Sommerfeld Center for Theoretical Physics and Center for NanoScience,
Department of Physics, Ludwig-Maximilians-Universität München, Theresienstraße
37, D-80333 München, Germany  Erwin Frey<EMAIL_ADDRESS>Arnold Sommerfeld Center
for Theoretical Physics and Center for NanoScience, Department of Physics,
Ludwig-Maximilians-Universität München, Theresienstraße 37, D-80333 München,
Germany  Max Planck School Matter to Life, Hofgartenstraße 8, D-80539 München,
Germany

###### Abstract

Biomolecular condensates help organize the cell cytoplasm and nucleoplasm into
spatial compartments with different chemical compositions. A key feature of
such compositional patterning is the local enrichment of enzymatically active
biomolecules which, after transient binding via molecular interactions,
catalyze reactions among their substrates. Thereby, biomolecular condensates
provide a spatial template for non-uniform concentration profiles of
substrates. In turn, the concentration profiles of substrates, and their
molecular interactions with enzymes, drive enzyme fluxes which can enable
novel non-equilibrium dynamics. To analyze this generic class of systems, with
a current focus on self-propelled droplet motion, we here develop a self-
consistent sharp interface theory. In our theory, we diverge from the usual
bottom-up approach, which involves calculating the dynamics of concentration
profiles based on a given chemical potential gradient. Instead, reminiscent of
control theory, we take the reverse approach by deriving the chemical
potential profile and enzyme fluxes required to maintain a desired condensate
form and dynamics. The chemical potential profile and currents of enzymes come
with a corresponding power dissipation rate, which allows us to derive a
thermodynamic consistency criterion for the passive part of the system (here,
reciprocal enzyme-enzyme interactions). As a first use case of our theory, we
study the role of reciprocal interactions, where the transport of substrates
due to reactions and diffusion is, in part, compensated by redistribution due
to molecular interactions. More generally, our theory applies to mass-
conserved active matter systems with moving phase boundaries.

The formation of biomolecular condensates through liquid-liquid phase
separation has been recognized as a driver of intracellular organization [1,
2, 3, 4, 5, 6, 7]. By tuning the local chemical composition, condensates can
compartmentalize reactions [4], buffer molecules [8], and are involved in
important cellular processes such as gene transcription [9, 10], or midcell
localization during cell division [11]. In contrast to thermal equilibrium
systems, in which the passive constituents of the condensates will completely
segregate and ultimately coarsen into a single droplet through Ostwald
ripening [12, 13, 14, 15], cells contain a myriad of active processes which
supply energy. For example, irreversible chemical reactions that are driven by
a continuous external supply of fuel [16] can enable ‘active droplets’ to
exhibit a wealth of novel phenomena not encountered in thermal equilibrium [6,
3, 17].

Previous studies have considered systems where chemical reactions drive a
continuous turnover of condensate material [17, 18, 19, 20, 21, 22, 23, 24,
25, 26, 27]. The resulting material fluxes lead to multi-droplet coexistence
[21, 22, 25, 26, 27] and droplet division [24]. Here, we build on our previous
work [28] to study a different, mass-conserved, class of systems where enzymes
spontaneously phase separate, or localize to an existing condensate [29].
These enzymes then regulate reactions among other molecules, by transiently
binding substrates and catalyzing their conversion into products via a lowered
activation barrier. For example, in the bacterium _Myxococcus xanthus_ , a
PomXY cluster (moving on the nucleoid) regulates the cycling of PomZ between
two conformations [11, 30, 31, 32]. Conceptually, this setup is analogous to
transcriptional condensates where transcription factors regulate the activity
of RNA Polymerase II, which assembles RNA from individual nucleotides [33, 9].
In turn, attractive electrostatic interactions between the negatively charged
RNA and the positively charged transcription factors favor further condensate
growth. At high RNA concentrations, however, repulsive interactions due to
volume exclusion lead to condensate dissolution [10, 34, 35]. Such feedback
mechanisms between active reactions and passive phase separation are
ubiquitous in the complex intracellular environment, and in membraneless
organelles such as the nucleolus [36].

For this generic class of systems, we have previously shown that the active
chemical turnover and buildup of substrate concentration gradients can, by
driving enzyme fluxes via reciprocal interactions, lead to condensate self-
propulsion, positioning, interrupted coarsening, and division [28]. While we
have previously considered a scenario where enzyme-substrate and enzyme-
product interactions are weak, we here relax this assumption. To that end, we
explicitly account for the effect of enzyme concentration gradients on the
diffusion of substrates and products. Moreover, we significantly extend our
theoretical framework and illustrate how one can determine the velocity of
self-propelling droplets in higher dimensions. We achieve this by developing a
self-consistent sharp interface theory, which proceeds in two steps. At the
outset, we determine the chemical potential and enzyme currents necessary to
keep the droplet in a steady state, maintaining both constant velocity and
shape. In doing so, one can relate the droplet motion to a required rate of
free energy dissipation. By collecting all processes that perform work on the
system into an active force field, we identify the chemical potential with the
free energy of a purely passive sub-component. Following this, we derive a
thermodynamic consistency relation grounded in the fundamental principle that
passive systems reach thermal equilibrium by minimizing their free energy. In
other words, our analysis shows how power dissipation can drive condensate
motion. The physical meaning of the thermodynamic consistency relation is that
the corresponding work is performed by applying an active force field, while
the free energy of the passive condensate material remains constant. Our
theory shows that the condensate motion is driven by the bulk interactions
between the various chemical species and does not require viscous hydrodynamic
coupling. This stands in contrast to previous studies, where liquid droplets
propelled through Marangoni flows in viscous fluids [37, 38, 39], through
active stresses on surfaces [40, 41, 42, 43, 44, 45], or by altering their
wetting properties [46, 47].

The outline of the present article is as follows. In Sec. I, we discuss the
theoretical framework of describing a mass-conserving multi-component mixture
containing enzymes which phase-separate through attractive interactions among
themselves, and also interact with other molecular species such as substrates
and products. Moreover, we give an account of non-equilibrium reactions that
can give rise to inhomogeneous concentration profiles of substrates and
products, thereby inducing an inhomogeneous driving force on the enzymes. As
the central contribution of our work, in Sec. II we show how the core of the
model can be reduced to an implicit description that is independent of the
specific details of the interactions and reactions that we have introduced in
the preceding discussion. In doing so, we derive self-consistency relations
that characterize the enzyme currents, chemical potential profile, and the
resulting droplet velocity for an arbitrary inhomogeneous driving force. We
test our theory and demonstrate that liquid-like droplets are more readily set
in motion than solid-like condensates. In Sec. III, we discuss the impact of
reciprocal interactions. To that end, we first show how reciprocal
interactions can give rise to discontinuities and kinks in these concentration
profiles, how the reactants are redistributed, and how the resulting
concentration profiles can be determined. With these tools in hand, we derive
the conditions to observe droplet self-propulsion in the presence of weakly or
strongly reciprocal interactions, thereby complementing our previous work
[28]. We find that the absence or presence of reciprocal interactions
determines if the onset of self-propulsion is continuous, or discontinuous.

##  I Enzyme-enriched condensates

###  I.1 Thermodynamic currents

Biomolecular condensates can, in principle, consist of many interacting
components [6]. In this study, we investigate regular mixtures consisting of a
solvent with enzyme, substrate, and product molecules. For simplicity, we
focus on the enzyme concentration $c(\mathbf{x},t)$ as the main constituent of
the condensates. We assume that in thermal equilibrium the enzymes can phase-
separate on their own, driven by enzyme-solvent interactions and entropy,
which we encode into a Ginzburg-Landau expansion of the free energy density,

$f_{0}(c)=-\frac{r}{2}\,(c-\tilde{c})^{2}+\frac{u}{4}\,(c-\tilde{c})^{4}+\frac{\kappa}{2}\,|\boldsymbol{\nabla}c|^{2}\,.$
(1)

The control parameter $r$ measures the distance from the critical point and,
for ${r>0}$, indicates attractive enzyme-enzyme interactions that lead to
phase separation. The parameter ${u>0}$ is required to thermodynamically
stabilize the system, so that the free energy density has the form of a
double-well potential. Finally, a positive stiffness parameter, ${\kappa>0}$,
penalizes concentration gradients.

We are interested in a scenario where the concentrations of substrates,
$s(\mathbf{x},t)$, and products, $p(\mathbf{x},t)$, are small and described by
the free energy of an ideal mixture,
${f_{\text{I}}(\varrho)=k_{\text{B}}T\,\varrho\log(\varrho\,\nu)}$, where
${\varrho\in\\{s,p\\}}$, $\nu$ represents the molecular volume, and
$k_{\text{B}}T$ is the thermal energy. By parameterizing enzyme-substrate and
enzyme-product interactions with Flory-Huggins (FH) parameters $\chi_{s}$ and
$\chi_{p}$, respectively, the thermodynamics of the mixture are characterized
by an effective free energy functional,

$\mathcal{F}=\int\\!d^{d}x\,\big{[}f_{0}(c)+f_{\text{I}}(s)+f_{\text{I}}(p)+\chi_{s}\,c\,s+\chi_{p}\,c\,p\big{]}\,,$
(2)

where $d$ refers to the number of spatial dimensions. Note that, in general,
the interaction between different species can be more intricate than the
current bilinear form, e.g., when substrate and product concentrations are
large enough for their virial coefficients to become relevant [10].

When a particle of species ${\varrho\in\\{s,p,c\\}}$ is added or removed at a
specific location $\mathbf{x}$, the system’s free energy incurs thermodynamic
costs given by the corresponding chemical potential
${\mu_{\varrho}=\delta\mathcal{F}/\delta\varrho}$. In the framework of linear
non-equilibrium thermodynamics [48, 49], gradients in chemical potentials act
as thermodynamic forces, corresponding to conservative currents that drive
particle exchange between adjacent points in space. To linear order, one has

$\begin{bmatrix}\mathbf{j}_{s}\\\ \mathbf{j}_{p}\\\
\mathbf{j}_{c}\end{bmatrix}=-\begin{bmatrix}\Lambda(s)&0&0\\\
0&\Lambda(p)&0\\\
0&0&M(c)\end{bmatrix}\cdot\begin{bmatrix}\boldsymbol{\nabla}\mu_{s}\\\
\boldsymbol{\nabla}\mu_{p}\\\ \boldsymbol{\nabla}\mu_{c}\end{bmatrix}\,,$ (3)

where we have disregarded any cross-terms and considered the matrix of Onsager
coefficients (mobilities) to be diagonal. For simplicity, we assume that the
mobilities $\Lambda(s)$ for the substrates and products share the same
functional dependence, which in general is distinct from the mobility $M(c)$
for the enzymes.

In the following, we first study the dynamics of the enzymes, which is
characterized by the continuity equation
${\partial_{t}c=-\boldsymbol{\nabla}\cdot\mathbf{j}_{c}}$. With the above
expression for the enzyme current $\mathbf{j}_{c}$ this takes the form of a
generalized Cahn-Hilliard (CH) model

$\partial_{t}c(\mathbf{x},t)=\boldsymbol{\nabla}\cdot\bigl{[}M(c)\,\boldsymbol{\nabla}\bigl{(}\mu_{0}(c)+\chi_{s}s+\chi_{p}p\bigr{)}\bigr{]}\,.$
(4)

Here, the chemical potential of the enzyme molecules,
${\mu_{0}(c)=-r(c-\tilde{c})+u(c-\tilde{c})^{3}-\kappa\boldsymbol{\nabla}^{2}c}$,
is identical to the CH model [50, 51]. For ${r>0}$, the CH model shows phase
separation into enzyme-rich and enzyme-poor regions, with concentrations
${c_{\pm}=\tilde{c}\pm\sqrt{r/u}}$. In the following, consistently with our
previous work [28], we choose $c_{+}$ as unit of concentration, and
${\epsilon_{0}\coloneqq r\,c_{+}}$ as unit of energy. In addition to the
equilibrium dynamics of the CH model, the above model contains currents that
arise from the interaction with the concentration fields of the substrate and
product molecules. We will next discuss mechanisms that can give rise to non-
uniform substrate and product concentration profiles.

###  I.2 Non-equilibrium reactions

In addition to the thermodynamic currents given by Eq. (3), biological systems
also contain active processes that break detailed balance. As in our prior
work [28], we consider a specific form of activity, where the enzymes $E$
catalyze a chemical reaction that turns substrates $S$ into products $P$ via a
Michaelis-Menten-like kinetics,

$\ce{{\hbox{E + S}}<=>{\hbox{ES}}<=>>{\hbox{EP}}<=>{\hbox{E + P}}}\,.$ (5)

For a thermodynamically consistent formulation of these reaction kinetics, one
needs to take into account the potential barriers that have to be overcome
during each step of the reaction [16, 28]. Here, we focus on a simplified
scenario where the first step $E+S$ ¡=¿ $ES$ is rate-limiting (high potential
barrier), and has no thermodynamically preferred direction (equal forward and
backward rates). In contrast, we assume that the non-rate-limiting second step
$ES$ ¡=¿¿ $EP$ is strongly biased towards products which have lower internal
energy than substrates. The overall reaction then proceeds with a net rate
$k_{1}\,c\,s$ proportional to the concentration ${c\equiv[E]}$ of enzymes and
the concentration ${s\equiv[S]}$ of substrates. In summary, the net reaction
rate is determined by the first reaction step, while free energy is released
in the second reaction step.

To maintain the dynamics and keep the system away from thermodynamic
equilibrium111Alternatively, one could also consider a sufficiently large
domain in which substrates are abundant, so that the system does not reach
thermodynamic equilibrium on the time scales of interest., the free energy
released from converting substrate into product needs to be resupplied by
exchanging product with substrate,

$\ce{{\hbox{P}}+{\hbox{F}}<=>>{\hbox{S}}+{\hbox{W}}}\,,$ (6)

along a separate reaction path which consumes fuel $F$ and releases waste $W$.
We assume that the forward reaction of Eq. (6) is driven by an excessive
abundance of (chemostatted) fuel, while the backward reaction is limited by
negligible concentrations of (chemostatted) waste. The mass action law then
suggests a net reaction rate $k_{2}\,p$ proportional to the concentration
$p\equiv[P]$ of products. In general, the rate constant $k_{2}$ will depend on
the height of the potential barrier for the forward reaction [16, 28] and
could therefore be affected by enzyme-substrate and enzyme-product
interactions. In particular, one needs to specify how the molecules in the
transition state (from product to substrate) interact with the enzymes. We
here assume that the molecules in the transition state still resemble
products, so that the corresponding interaction energy, at the maximum of the
potential barrier, is the same as when product and fuel are still separated222
For instance, this scenario holds when products do not interact with enzymes
and the chemical potential of the transition state remains constant [28]..
Thus, the height of the free energy barrier is set by the difference in the
internal energy of the molecules, an intensive quantity that can be taken
constant, so that $k_{2}$ is concentration-independent.

Taken together, the reaction-diffusion dynamics of substrates and products are
determined by

$\displaystyle\partial_{t}s$
$\displaystyle=\boldsymbol{\nabla}\cdot\left(D\,\boldsymbol{\nabla}s+\Lambda\,s\,\chi_{s}\boldsymbol{\nabla}c\right)-k_{1}\,c\,s+k_{2}\,p\,,$
(7a) $\displaystyle\partial_{t}p$
$\displaystyle=\boldsymbol{\nabla}\cdot\left(D\,\boldsymbol{\nabla}p+\Lambda\,p\,\chi_{p}\boldsymbol{\nabla}c\right)+k_{1}\,c\,s-k_{2}\,p\,,$
(7b)

where we have assumed that the mobility functions are given by
${\Lambda(s)=\Lambda\,s}$ and ${\Lambda(p)=\Lambda\,p}$ with constant mobility
coefficient $\Lambda$. For thermal Brownian motion of the particles, the
diffusivity is set by the Einstein-Smoluchowsky relation ${D=\Lambda\,k_{B}T}$
[52]. Here, we take the liberty of deviating from this relation, as the system
is considered to be driven away from equilibrium by active processes that
break detailed balance. The decay rate of products sets a characteristic
timescale for the reaction-diffusion dynamics, ${\tau_{0}\coloneqq
k_{2}^{-1}}$, while the corresponding diffusion length sets the length scale,
${l_{0}\coloneqq\sqrt{D/k_{2}}}$.

Figure 1: Concentration profiles of substrates and products for a 3D enzyme
droplet which has radius $R$ (scale bar) and moves with constant speed
${v=4.14\,l_{0}/\tau_{0}}$ through an open domain. The black circle denotes
the boundary of the enzyme-rich region (with concentration $c_{+}$). Colored
arrows illustrate reactive fluxes: Enzymes catalyze the conversion of
substrates into products with rate $k_{1}\,c\,s$, while substrates are
replenished from products with rate $k_{2}\,p$. Upper panel (blue) shows
product concentration profile relative to the far-field, ${\delta
p=p-p_{\infty}}$, while lower panel (red) shows substrate concentration
profile relative to the far-field, ${\delta s=s-s_{\infty}}$, both in units of
the enzyme concentration $c_{+}$ within the droplet. Vector field (thin black
arrows) indicates the excess diffusive currents of substrates and products
(a.u.) relative to the co-moving frame. The remaining parameter values are
${c_{-}=0.5c_{+}}$, ${R=l_{0}}$, ${k_{1}=2k_{2}/c_{+}}$, ${\Lambda=0}$, and
${s+p=c_{+}}$.

In the present work, we focus on the analysis of droplets that exhibit a
traveling steady state which moves with velocity $\mathbf{v}$. These droplets
are characterized by a traveling enzyme concentration profile
${c(\mathbf{x}-\mathbf{v}t)}$, which will be the focus of the next section.
Thus, we substitute the traveling wave ansatz
${s(\mathbf{x},t)=s(\mathbf{x}-\mathbf{v}t)}$ and
${p(\mathbf{x},t)=p(\mathbf{x}-\mathbf{v}t)}$ with the Galilean coordinate
transformation ${\mathbf{z}\coloneqq\mathbf{x}-\mathbf{v}t}$. In the
corresponding co-moving frame, the steady-state concentrations of substrates
and products are determined by

$\displaystyle 0$
$\displaystyle=\boldsymbol{\nabla}\cdot\left(\mathbf{v}\,s+D\,\boldsymbol{\nabla}s+\Lambda\,s\,\chi_{s}\boldsymbol{\nabla}c\right)-k_{1}\,c\,s+k_{2}\,p\,,$
(8a) $\displaystyle 0$
$\displaystyle=\boldsymbol{\nabla}\cdot\left(\mathbf{v}\,p+D\,\boldsymbol{\nabla}p+\Lambda\,p\,\chi_{p}\boldsymbol{\nabla}c\right)+k_{1}\,c\,s-k_{2}\,p\,,$
(8b)

where the gradient $\boldsymbol{\nabla}$ is taken with respect to
$\mathbf{z}$.

Consider now a single round enzyme droplet characterized by a domain
${\mathcal{D}\coloneqq\\{\mathbf{z}:|\mathbf{z}|<R\\}}$ with high enzyme
concentration ${c(\mathbf{z})=c_{+}}$, which is surrounded by an enzyme-
depleted domain (${|\mathbf{z}|\geq R}$) with concentration
${c(\mathbf{z})=c_{-}}$. Given such an enzyme concentration profile, we solved
Eq. (8) by using finite element methods (FEM). To emulate the infinitely large
open space through which the droplet moves, we imposed Dirichlet boundary
conditions at the boundary of the finite simulation domain. For the
corresponding values of the far-field concentrations we used the reactive
equilibria ${s_{\infty}=n/(1+k_{1}c_{\infty}/k_{2})}$ and
${p_{\infty}=n-s_{\infty}}$, where $n$ is the total average concentration of
substrates and products and
${c_{\infty}\coloneqq\lim_{|\mathbf{z}|\to\infty}c(\mathbf{z})}$.

The resulting concentration profiles and fluxes of substrates and products are
shown in Fig. 1. Inside the condensate, where the enzyme concentration is
high, substrates are converted into products. These reactions lead to
concentration gradients which drive diffusive product fluxes out of the
droplet, while substrates are replenished by diffusing into the droplet; see
the vector fields in Fig. 1. The cycle of currents maintaining the non-
equilibrium steady state is closed by reactions which turn products into
substrates outside of the condensate. Finally, as the droplet keeps moving
forward, which corresponds to rearward advection in the co-moving frame, the
extrema of the substrate and product concentration profiles trail behind the
centroid of the condensate. With increasing droplet speed, the maximum of the
product concentration profile and the minimum of the substrate concentration
profile move closer to the trailing edge of the condensate.

##  II Self-consistency relation for driven condensate motion

As a result of enzyme-mediated reactions, the concentration profiles of both
substrates and products are in general non-uniform. These non-uniform
distributions result in enzyme currents driven by enzyme-substrate and enzyme-
product interactions, as described by the generalized CH equation [Eq. (4)].
The associated thermodynamic force is given by

$\mathbf{f}(\mathbf{x},t)\coloneqq-\boldsymbol{\nabla}\big{[}\chi_{s}\,s(\mathbf{x},t)+\chi_{p}\,p(\mathbf{x},t)\big{]}\,.$
(9)

While we are using this specific form here, the theory presented in the
following is independent of the specific physical mechanism that generates the
force field ${\mathbf{f}(\mathbf{x},t)}$, and it applies equally to externally
imposed driving forces or to forces generated from active processes. Given
such an inhomogeneous forcing $\mathbf{f}(\mathbf{x},t)$, this section
outlines how the resulting droplet velocity ${\mathbf{v}}$ can be determined
self-consistently in arbitrary spatial dimensions. This is a cornerstone of
our work and generally applies to any system exhibiting phase separation into
well-defined droplets, irrespective of the underlying free-energy functional
describing the interactions causing phase separation.

_Sharp interface limit._ The central concept underlying our theoretical
analysis is a sharp interface limit, which assumes that the width of the
droplet interface ${w=\sqrt{2\kappa/r}}$ is small compared to all other length
scales relevant to the dynamics. While this width diverges near the critical
point (${r=0}$), it approaches a molecular length scale given by the
characteristic range of protein-protein interactions when phase separation is
sufficiently strong. Therefore, we expect that a sharp interface limit serves
as a good approximation when the condensates (droplets) are much larger in
size than individual molecules. Mathematically, we implement this limit by
taking ${r\rightarrow\infty}$ and ${u\rightarrow\infty}$ while maintaining
finite values for mesoscopic quantities such as the equilibrium enzyme
concentrations, and the chemical potential for the enzymes. This chemical
potential can be written in the form

$\mu_{0}(c)=r\left[-(c-\tilde{c})+\frac{4}{\Delta
c^{2}}(c-\tilde{c})^{3}-\frac{1}{2}w^{2}\boldsymbol{\nabla}^{2}c\right]\,,$
(10)

where ${\Delta c\coloneqq c_{+}-c_{-}}$ is the difference in the equilibrium
enzyme concentrations between the inside and the outside of the droplet. For
the chemical potential to remain finite and thus physical in the asymptotic
limit ${r\to\infty}$, the expression in the square brackets must vanish. This
is precisely the same condition as that of a CH theory, resulting in a
$\tanh$-profile with width $w$, height $\Delta c$, and asymptotic values
${c_{\pm}=\tilde{c}\pm\sqrt{r/u}}$ [50, 51]. In the sharp interface limit we
are considering here, this profile can be considered as piecewise constant.

If there were deviations $\delta c$ around these equilibrium concentrations
$c_{\pm}$, they would quickly relax back to equilibrium because the
corresponding chemical potential,

$\mu_{0}=2r\,\delta c-\frac{1}{2}\,r\,w^{2}\boldsymbol{\nabla}^{2}\delta
c+\mathcal{O}(\delta c^{2})\,,$ (11)

would exhibit large gradients. These chemical potential gradients, with
${|\boldsymbol{\nabla}\mu_{0}|\gg|\chi_{s}\boldsymbol{\nabla}s|+|\chi_{p}\boldsymbol{\nabla}p|}$,
would transiently dominate the dynamics and drive currents that promptly
restore a piecewise constant enzyme concentration profile. As a result, all
variations in the concentration profile are localized at the droplet boundary,
where the competing terms in Eq. (10) are all of the same order and much
larger than the interaction terms with substrates and products, ${r\,\Delta
c\gg|\chi_{s}\,s|+|\chi_{p}\,p|}$.

_Approximation of round droplet shape._ Next, we discuss the geometric shape
of the droplet, which will determine the piecewise constant enzyme
concentration profile in the sharp interface limit. To do so, it is important
to note that the effective surface tension of a droplet in the CH model is
given by [15]

$\gamma=\frac{1}{6}(\Delta c)^{2}\sqrt{2\kappa r}=\frac{1}{6}(\Delta
c)^{2}r\,w\,,$ (12)

with ${\Delta c=c_{+}-c_{-}}$ the density jump at the interface. Since the
interface width $w$ is bounded from below by the molecule diameter, the
effective tension diverges in the limit ${r\rightarrow\infty}$. For this
reason, we restrict ourselves to round droplets with a constant interface
curvature.

###  II.1 Sharp interface theory

As the above considerations show, Eq. (10) is no longer a mathematically well-
defined expression for the chemical potential. Therefore, we take an
alternative approach and consider the chemical potential as a Lagrange
multiplier field that enforces the condition of a piecewise constant
concentration profile: ${c(\mathbf{x})=c_{+}}$ within the domain $\mathcal{D}$
of the droplet and ${c(\mathbf{x})=c_{-}}$ outside. The following analysis
will describe how to obtain an equation that determines the chemical
potential.

_Traveling wave solution._ We first seek a traveling wave solution for the
enzyme dynamics by substituting the ansatz
${c(\mathbf{x},t)=c(\mathbf{x}-\mathbf{v}t)}$ into Eq. (4). One finds that the
steady state of the enzyme concentration profile (in the co-moving frame) is
given by a balance condition between advective terms (left-hand side) and
driving forces (right-hand side)

$-\mathbf{v}\cdot\boldsymbol{\nabla}c(\mathbf{z})=\boldsymbol{\nabla}\cdot\bigl{[}M(c)\,\bigl{(}\boldsymbol{\nabla}\mu_{0}(\mathbf{z})-\mathbf{f}(\mathbf{z})\bigr{)}\bigr{]}\,.$
(13)

As explained in the previous section, there is no explicit form for the
chemical potential $\mu_{0}(\mathbf{z})$ in the sharp interface limit.
Instead, it serves as a Lagrange multiplier field that ensures constant enzyme
densities inside and outside the droplet. This means that the enzyme
concentration profile denoted $c(\mathbf{z})$ and the associated mobility
function, $M(c)$, are assumed to be known. Furthermore, we temporarily assume
knowledge of the droplet velocity $\mathbf{v}$ and will later demonstrate how
to determine it self-consistently. Given the force field
$\mathbf{f}(\mathbf{z})$, and with appropriate boundary conditions which will
be discussed in the next section, one can now use FEM to determine the
chemical potential—the only unknown in the present scheme.

In the sharp interface limit, solving Eq. (13) by FEM requires a careful
analysis since gradients in the enzyme concentration become singular at the
interface. To address this, we rearrange terms in Eq. (13) to express it as:

$\boldsymbol{\nabla}\cdot\mathbf{J}_{\Delta}=0\,,$ (14a) with the current
given by
$\mathbf{J}_{\Delta}=M(c)\bigl{(}-\boldsymbol{\nabla}\mu_{0}+\mathbf{f}\bigr{)}-\mathbf{v}\,\bigl{(}c-c_{\infty}\bigr{)}\,.$
(14b)

Equation (14b) can be viewed as the integral of Eq. (13) and, as explained in
the following, interpreted as the divergence-free current of the enzymes in
the co-moving frame of the droplet. The first term describes the currents that
arise in the laboratory frame due to gradients in the chemical potential of
the enzymes and the given force field $\mathbf{f}(\mathbf{x})$. The second
term is an apparent advective enzyme current, which an observer in the co-
moving frame will see due to the motion relative to the laboratory frame. To
ensure this interpretation we added the integration constant
$\mathbf{v}\,c_{\infty}$ to Eq. (14b). Then, the enzyme currents vanish far
away from the droplet,
${\lim_{|\mathbf{z}|\to\infty}\mathbf{J}_{\Delta}(\mathbf{z})=0}$, where the
enzyme concentration has a low value
${\lim_{|\mathbf{z}|\to\infty}c(\mathbf{z})=c_{\infty}\equiv c_{-}}$. To
approximate these far-field conditions in our numerical implementation, we use
a finite domain with a size much larger than the droplet radius, and impose
no-flux boundary conditions at the domain boundary. This defines the boundary
conditions for the calculation of the chemical potential.

To now apply FEM, we convert Eq. (14a) into an optimization problem (the _weak
form_) by multiplying by a test function $\phi(\mathbf{z})$ and integrating
over the entire domain. In doing so, the singularity in the concentration
profile can be lifted through integration by parts. The numerical solution for
the chemical potential then satisfies

$\int\\!d^{d}z\,\bigl{[}\mathbf{J}_{\Delta}\cdot\boldsymbol{\nabla}\phi\bigr{]}=0\,.$
(15)

Solving this equation with the open-source FEM framework FEniCSx [53, 54, 55,
56], one obtains the Lagrange multiplier field $\mu_{0}(\mathbf{z})$ which
best enforces incompressibility of the enzyme currents [Eq. (14a)] and thus
maintains the enzyme concentration profile. Finally, the corresponding enzyme
currents can be obtained by inserting the numerical solution for the chemical
potential into Eq. (14b). Figure 2(a) shows both the chemical potential
profile and the corresponding enzyme currents for a droplet in a self-
propelling steady state, where the force field
${\mathbf{f}(\mathbf{z})\coloneqq-\boldsymbol{\nabla}\big{[}\chi_{s}\,s(\mathbf{z})+\chi_{p}\,p(\mathbf{z})\big{]}}$,
cf. Eq. (9), is generated by inhomogeneous substrate and product concentration
profiles as discussed in Sec. I.2.

Figure 2: Sharp interface theory for a moving 3D enzyme droplet. (a) Chemical
potential profile of a moving droplet (circle), obtained by solving Eq. (15).
The vector field shows the (mesoscopic) enzyme velocity (a.u.) relative to the
co-moving frame, $\mathbf{J}_{\Delta}/c$, as calculated from Eq. (14b). These
currents are driven by a thermodynamic force, Eq. (9), due to interactions
with substrates and products. The droplet speed, reaction-diffusion dynamics,
as well as substrate and product concentration profiles are taken from Fig. 1.
Remaining parameters: ${\Delta\chi\coloneqq\chi_{p}-\chi_{s}=4r}$ and
${M(c)=Mc}$ with ${M=10D/\epsilon_{0}}$. (b) Chemical potential imbalance
$\Delta\mu_{0}(v)$ as a function of the droplet speed $v$ with stable (filled
circles) and unstable (open circles) steady states. The effective phase space
flow of the droplet speed $v$ driven by the chemical potential imbalance is
shown as arrows. For the indicated stable steady state, the corresponding
enzyme concentration and flow field are shown in panel a), whereas the product
and substrate concentration and flow fields are shown in Fig. 1.

In summary, we have developed an _implicit adiabatic elimination scheme_ :
Whereas standard methods consider how a given chemical potential profile [Eq.
(10)] gradually drives the system toward a state of phase separation,
requiring a dynamical resolution of the concentration profiles on the length
scale of the interface, our approach starts with the premise that the system
undergoes phase separation and treats the corresponding concentration profile
as a fast degree of freedom. This allows us to invert the logical flow, and
instead analyze the chemical potential profile which is required to
adiabatically maintain the concentration profile in a predefined
configuration. As a consequence of the reversal of the logical flow, prompted
by the assumption of a sharp interface, explicit expressions for the chemical
potential and the corresponding free energy of enzyme-enzyme interactions are
abandoned in favor of an implicit description. For the reasoning underlying
our model reduction scheme, the details of the original explicit model are not
important. Instead, the dynamics is universally determined by the force field,
the concentration (order parameter) profile, and the profile velocity. We thus
expect that our approach can also be adapted to Flory-Huggins-like models in
the presence of non-equilibrium processes or to fully non-equilibrium systems
exhibiting moving fronts, such as population expansion or reaction-diffusion
dynamics.

_A thermodynamic consistency relation for the chemical potential yields the
droplet velocity._ So far we have regarded the droplet velocity as a given
quantity. We will now discuss how it can be determined through a condition on
the chemical potential profile required by the thermodynamic consistency of
the steady state.

The sharp interface theory outlined above predicts the chemical potential
field $\mu_{0}(\mathbf{z})$ and the associated enzyme currents in response to
a specified force field $\mathbf{f}(\mathbf{z})$ for any given velocity
$\mathbf{v}$. While the force field $\mathbf{f}(\mathbf{z})$ originates from a
non-equilibrium process or is imposed externally, gradients in the chemical
potential are thermodynamic driving forces. Consequently, these gradients are
linked to physical conditions that could help in determining the droplet
velocity $\mathbf{v}$. To characterize self-propelling and thus polar
droplets, we define the chemical potential imbalance as a function of the
droplet speed,

$\Delta\mu_{0}(v)\coloneqq\int_{\mathcal{D}}\\!d^{d}z\,\bigl{(}\mathbf{\hat{e}}_{v}\cdot\boldsymbol{\nabla}\mu_{0}\bigr{)}\,.$
(16)

This definition is motivated by the fact that it gauges the overall
thermodynamic driving force acting on a droplet: It gives the average chemical
potential gradient within the droplet taking into account rotational symmetry
about the direction of motion ($\mathbf{\hat{e}}_{v}$). The dependence of
$\Delta\mu_{0}(v)$ on the droplet speed is shown in Fig. 2(b), where we use
for specificity a force field $\mathbf{f}(\mathbf{z})$ generated by
inhomogeneous substrate and product concentration profiles, as described in
section I.2.

One may now interpret the function $\Delta\mu_{0}(v)$ as an effective force or
its integral as an effective potential acting on the droplet speed $v$,
closely resembling the dynamics of a one-component nonlinear system. This
indicates that the droplet speed at a steady state is determined by the
condition that the chemical potential imbalance must vanish:

$\Delta\mu_{0}(v)=0\,.$ (17)

In addition, the slope of the chemical potential disequilibrium provides
information about the stability of the steady state: Since positive values of
$\Delta\mu_{0}(v)$ increase the velocity of the droplet, while negative values
tend to slow it down, a negative slope indicates stability. The resulting
phase space flow is represented by the arrows in Fig. 2(b). In the following
section, the above heuristic considerations are substantiated by detailed
thermodynamic arguments.

###  II.2 Thermodynamically consistent chemical potential is a consequence of
energy conservation

We will now illustrate the physical meaning of the thermodynamic consistency
criterion [Eq. (17)], by linking our analysis to non-equilibrium
thermodynamics.

_Power dissipation._ Using the sharp interface theory, we have determined the
chemical potential profile $\mu_{0}(\mathbf{z})$ which is required to maintain
a preset enzyme concentration profile $c(\mathbf{z})$, and a given droplet
velocity $\mathbf{v}$. In turn, in the presence of such a non-uniform chemical
potential333For the Cahn-Hilliard model with a diffuse interface, the
mechanisms giving rise to this chemical potential are encoded in the free
energy density [Eq. (1)]., variations in the concentration profile, e.g., due
to droplet motion, translate into changes in the free energy
${\mathcal{F}_{0}[c]\coloneqq\int\\!d^{d}x\,f_{0}(c)}$ associated with enzyme-
enzyme interactions,

$\delta\mathcal{F}_{0}[c]=\int\\!d^{d}x\,\mu_{0}(\mathbf{x})\,\delta
c(\mathbf{x})\,.$ (18)

Thus, the chemical potential will dissipate free energy with a rate

$\partial_{t}\mathcal{F}_{0}[c]=\int\\!d^{d}x\,\mu_{0}(\mathbf{x})\,\partial_{t}c(\mathbf{x},t)\,.$
(19)

One can rewrite this expression by substituting the traveling wave ansatz,
${c(\mathbf{x},t)=c(\mathbf{x}-\mathbf{v}t)}$, and by partial integration,
which yields

$\partial_{t}\mathcal{F}_{0}[c]=\mathbf{v}\cdot\int\\!d^{d}z\,c(\mathbf{z})\,\boldsymbol{\nabla}\mu_{0}(\mathbf{z})\,.$
(20)

Note that this is simply minus the power dissipated by a flux down a chemical
potential gradient, which can also be derived by using fundamental relations
between power, work and force444 The motion of particles, with net current
${\mathbf{j}(\mathbf{z})=c(\mathbf{z})\,\mathbf{v}(\mathbf{z})}$, in the
presence of a force field $\mathbf{f}(\mathbf{z})$ will, over time, perform
work with rate
$\partial_{t}W=\int\\!d^{d}z\,\mathbf{j}(\mathbf{z})\cdot\,\mathbf{f}(\mathbf{z})\,.$
This general relation between power dissipation and force is independent of
invoking thermodynamic arguments. One recovers Eq. (20) by identifying
$-\boldsymbol{\nabla}\mu_{0}$ as a thermodynamic force which extracts the work
${\delta W=-\delta\mathcal{F}_{0}}$ from the free energy functional. .

Utilizing that the enzyme concentration profile $c(\mathbf{z})$ is piecewise
constant, we partition the aforementioned integral into two domains:
$\mathcal{D}$ within the droplet, and $\Omega\backslash\mathcal{D}$ outside of
it, where $\Omega$ denotes the whole domain. With
${c(\mathbf{z})=c_{+}\,\forall\,\mathbf{z}\in\mathcal{D}}$ and
${c(\mathbf{z})=c_{-}\,\forall\,\mathbf{z}\notin\mathcal{D}}$, one then has

$\partial_{t}\mathcal{F}_{0}=\mathbf{v}\cdot\biggl{[}c_{+}\int_{\mathcal{D}}\\!d^{d}z\,\boldsymbol{\nabla}\mu_{0}+c_{-}\int_{\Omega\backslash\mathcal{D}}\\!d^{d}z\,\boldsymbol{\nabla}\mu_{0}\biggr{]}.$
(21)

To reformulate the second term in the square brackets, we use the Gauss
theorem which implies that the integral of $\boldsymbol{\nabla}\mu_{0}$ over
the whole domain $\Omega$ vanishes

$\int_{\Omega\backslash\mathcal{D}}d^{d}z\,\boldsymbol{\nabla}\mu_{0}+\int_{\mathcal{D}}d^{d}z\,\boldsymbol{\nabla}\mu_{0}=\oint_{\partial\Omega}\\!d^{d-1}\boldsymbol{S}\,\mu_{0}=0\,.$
(22)

Note that the last equality holds because the chemical potential must be
constant in the far field and because the integral of the normal vector over a
closed surface is always zero. Substituting Eq. (22) into Eq. (21), one has

$\partial_{t}\mathcal{F}_{0}=v\,\Delta
c\,\left[\mathbf{\hat{e}}_{v}\cdot\int_{\mathcal{D}}\\!d^{d}z\,\boldsymbol{\nabla}\mu_{0}\right],$
(23)

where we have used $\mathbf{v}=v\,\mathbf{\hat{e}}_{v}$ and $\Delta
c=c_{+}-c_{-}$. By comparing this expression with the definition of the
chemical potential imbalance [Eq. (16)], one finds

$\partial_{t}\mathcal{F}_{0}=v\,\Delta c\,\Delta\mu_{0}(v)\,.$ (24)

Therefore, the heuristically derived thermodynamic consistency relation
${\Delta\mu_{0}(v)=0}$ [Eq. (17)] equates to the minimization of the free
energy $\mathcal{F}_{0}[c]$ associated with enzyme-enzyme interactions. In
other words, in steady state, the currents along the chemical potential
profile $\mu_{0}(\mathbf{z})$ should cease to dissipate power.

Note that, when the internal energy $\mathcal{U}$ of a system component does
not change over time because it is in a (non-equilibrium) steady state, one
can use the thermodynamic relation ${\mathcal{F}=\mathcal{U}-T\,\mathcal{S}}$
to relate Eq. (24) to an entropy production rate,
${\partial_{t}\mathcal{S}=-T^{-1}\partial_{t}\mathcal{F}}$. One can use these
concepts to also calculate the entropy production rate due to flows along the
force field $\mathbf{f}(\mathbf{z})$.

_Thermodynamic housekeeping._ To elucidate the significance of the above
result, we start with a scenario in which the force field is absent,
${\mathbf{f}(\mathbf{z})=0}$. Because there is no net driving force,
thermodynamics requires that the droplet can only be in a stable steady state
when it is at rest, ${v=v^{\star}=0}$. This state is characterized by a
vanishing chemical potential imbalance, $\Delta\mu_{0}(v^{\star})=0$, and by a
vanishing power dissipation. What are then the characteristics of a
hypothetical state with finite velocity? To answer this question, we
determined the chemical potential for a given droplet velocity by solving Eq.
(13); for an illustration see Supplemental Video 1. This shows that for any
state with ${v\neq 0}$, the chemical potential at the trailing edge of the
droplet must be higher than at its leading edge, so that
${v\,\Delta\mu_{0}(v)<0}$, and the power dissipation
$-\partial_{t}\mathcal{F}_{0}$ is positive [Eq. (24)]. Hence, droplet motion
would continuously siphon energy from the enzyme-enzyme interactions. To then,
in turn, adiabatically maintain enzyme-enzyme interactions and the enzyme
concentration profile, one would need to inject energy into the system. Taken
together, energy conservation rules out a self-propelling droplet state in the
absence of an external power supply.

To pursue these thermodynamic arguments further, we next analyze how the power
dissipation changes when the steady state is perturbed. To that end, we first
evaluate

$\left.\frac{\partial}{\partial
v}\,\partial_{t}\mathcal{F}_{0}\right\rvert_{v=v^{\star}}=v^{\star}\Delta
c\,\partial_{v}\Delta\mu_{0}(v^{\star})\,,$ (25)

which vanishes for ${v^{\star}=0}$. Thus, the steady state is marginally
stable at the linear level and one needs the second derivative to determine
its stability:

$\left.\frac{\partial^{2}}{\partial
v^{2}}\,\partial_{t}\mathcal{F}_{0}\right\rvert_{v=v^{\star}=0}=2\Delta
c\,\partial_{v}\Delta\mu_{0}(v^{\star})\biggr{\rvert}_{v^{\star}=0}\,.$ (26)

Based on the arguments laid out in the previous paragraph, the steady state,
which has vanishing power dissipation, is _stable_ when the power dissipation
${-\partial_{t}\mathcal{F}_{0}}$ is _minimal_ and hence Eq. (26) is negative.
This leads us to investigate the implications of a local _maximum_ in the
power dissipation as a function of the velocity.

To that end, we next consider a scenario in which a force field is present.
Specifically, we consider a system where this force field is generated by the
non-uniform concentration profiles of substrates and products [Eq. (9)], which
arise from reactions, diffusion, and advection in the co-moving frame of the
condensate [Eq. (8)]. As in the above analysis, we again solve Eq. (13) for
given droplet velocities and monitor the chemical potential profile; see
Supplemental Video 2. Importantly, since the force field depends on the
droplet velocity, the chemical potential imbalance becomes a non-monotonic
function of the droplet speed $v$; shown in Fig. 2(b) for the setup discussed
here. In particular, there are now several values $v^{\star}$ for the velocity
at which the power dissipation and the chemical potential imbalance vanish,
${\Delta\mu_{0}(v^{\star})=0}$. As shown above in Eq. (24), the power
dissipation $-\partial_{t}\mathcal{F}_{0}$ at ${v^{\star}=0}$ still vanishes
but now corresponds to a _local maximum_ as a function of the droplet velocity
[cf. Fig. 2(b) and Eq. (26)]. To maintain a state with a very small velocity,
one would thus need to take energy out of the system, lest the condensate will
spontaneously accelerate, which suggests that ${v^{\star}=0}$ is an unstable
steady state. Next, we characterize the steady states ${v^{\star}\neq 0}$ and,
without loss of generality, consider the case ${v>0}$. The slope of the
chemical potential imbalance shown in Fig. 2(b) is negative,
${\partial_{v}\Delta\mu_{0}(v^{\star})<0}$. Together with Eq. (25), this
implies that a further increase in the velocity would require constant energy
injection. Conversely, a slight decrease in the velocity would require
constant energy elimination. This is the signature of a stable steady state
with a finite velocity. These insights allow constructing a phase portrait
[c.f. arrows in Fig. 2(b)] reminiscent of dynamical systems theory [57].

###  II.3 Test of the sharp interface theory

To test our sharp interface theory, we compare its predictions, in particular
the thermodynamic consistency condition [Eq. (17)], with numerical simulations
of the full condensate dynamics that take into account the diffuse interface
of the enzyme droplet as described by the generalized CH equation [Eq. (4)].
In both approaches, the enzyme currents are driven by inhomogeneous substrate
and product concentration profiles. In the simulations of the full dynamics,
we determine these concentration profiles by solving the time-dependent
reaction-diffusion equations for substrates and products, which are given in
Sec. I.2, in the laboratory frame [Eq. (7)]. In our sharp interface theory, in
the same way as above, we use the steady-state profiles [Fig. 1] obtained by
solving the corresponding reaction-diffusion equations for substrates and
products in the co-moving frame [Eq. (8)].

To determine the velocities where the thermodynamic consistency condition [Eq.
(17)] is fulfilled, different numerical algorithms can be used. For instance,
one could perform a parametric sweep of the chemical potential imbalance for
different velocity values and directly identify the roots from the resulting
graph, as shown in Fig. 2(b). Here, as a numerically more efficient method, we
use Newton iterations starting from an initial guess for the velocity;
specifically, we set ${v=10.0\,l_{0}/\tau_{0}}$. We find that the
thermodynamic consistency condition correctly predicts the onset of the self-
propulsion of spherical droplets [Fig. 3]. Moreover, as further discussed in
Appendix A, the sharp interface theory reproduces the droplet speed with
reasonable quantitative accuracy. This holds true even for droplets with weak
phase separation [${\Delta c\ll c_{+}}$, Fig. 3] and for those that dissolve
over time due to enzyme-substrate and enzyme-product interactions
[Supplemental Movies 3 and 4 and Appendix B].

Figure 3:  Test of the sharp interface theory for moving enzyme condensates.
The top row shows the predictions of the self-consistent sharp interface
theory, while the bottom row shows the results of the full simulations solving
Eqs. (4) and (7). For this comparison, we consider two different droplet
concentrations, ${c_{-}=0.5\,c_{+}}$ and ${c_{-}=0.1\,c_{+}}$, as indicated in
the graph. The predictions of the sharp interface theory for the phase
boundary of the self-propulsion instability (solid lines) correctly mark the
onset of droplet motion in the full simulations, both for 2D and 3D systems.
Moreover, the theoretically predicted droplet velocities also agree
quantitatively with the simulations. As parameters we used
${\chi_{s}=-0.05r}$, ${\chi_{p}=-0.01r}$, ${\Lambda=0}$, ${w=0.1l_{0}}$,
${R=l_{0}}$, and ${s+p=c_{+}}$. In 2D, we simulated the full dynamics in a
circular domain having radius ${L=5l_{0}}$ for ${c_{-}=0.5c_{+}}$ and radius
${L=7l_{0}}$ for ${c_{-}=0.1c_{+}}$. Full 3D simulations were performed in a
cylindrical domain having radius ${L_{r}=3l_{0}}$ and length ${L_{z}=10l_{0}}$
for ${c_{-}=0.5c_{+}}$ and radius ${L_{r}=4l_{0}}$ and length
${L_{z}=14l_{0}}$ for ${c_{-}=0.1c_{+}}$. The simulation data for
${c_{-}=0.1c_{+}}$ are taken from our previous work, Ref. [28].

###  II.4 Analytically solvable limiting case of the sharp interface theory

The sharp interface analysis presented above typically requires a multi-step
numerical evaluation, while closed analytical solutions exist only in special
cases. To illustrate the validity of our numerical approach beyond a
comparison between simulation and theory, we will now discuss one such special
case where an analytic solution can be found. We consider a scenario where the
enzyme flux $\mathbf{J}_{\Delta}$ vanishes outside the condensate. This
applies in three specific cases:  (i) For a 1D system, where the current must
be spatially uniform due to Eq. (14a), the no-flux far-field condition implies
that the current must vanish in the entire domain;  (ii) when the mobility
vanishes in the low-concentration phase, ${M(c_{-})=0}$, it follows from Eq.
(14b) that there can be no flux;  (iii) when the enzymes are completely
depleted in the low-concentration phase, ${c_{-}=0}$, the lack of particles
implies a vanishing net current.  If at least one of these conditions holds,
one can use an approach similar to that described in our earlier work [28]. We
first integrate Eq. (14b) over the domain of the condensate:

$\int_{\mathcal{D}}\\!d^{d}z\,\mathbf{J}_{\Delta}=M(c_{+})\int_{\mathcal{D}}\\!d^{d}z\,\bigl{(}-\boldsymbol{\nabla}\mu_{0}+\mathbf{f}\bigr{)}-V\,\Delta
c\,\mathbf{v}\,,$ (27)

where ${V\coloneqq\int_{\mathcal{D}}\\!d^{d}z}$ is the volume of the droplet.
We will now show that the left-hand side of Eq. (27) is zero. To that end, we
manipulate the expression
${\int\\!dV\,\boldsymbol{\nabla}\cdot(\mathbf{x}\otimes\mathbf{J}_{\Delta})}$
to find

$\int\\!dV\,\mathbf{J}_{\Delta}+\int
dV\,\mathbf{x}\,(\nabla\cdot\mathbf{J}_{\Delta})=\oint
d\mathbf{S}\cdot(\mathbf{x}\otimes\mathbf{J}_{\Delta})\,,$ (28)

where we applied the product rule on the left-hand side and the divergence
theorem on the right-hand side. Taking into account that the enzyme currents
are divergence-free [Eq. (14a)], and assuming that they vanish outside the
condensate and thus along its boundary $\partial\mathcal{D}$, implies that the
average current in the droplet must vanish. This is intuitive, because it
simply means that the center of mass is stationary in the co-moving frame.
Equation (27) can be further simplified by using the thermodynamic consistency
criterion for the chemical potential, Eq. (17). By projecting Eq. (27) on the
axis $\mathbf{\hat{e}}_{v}$ and solving for the droplet speed $v$, one then
finds

$v=\frac{M(c_{+})}{V\,\Delta
c}\int_{\mathcal{D}}\\!d^{d}z\,\mathbf{\hat{e}}_{v}\cdot\mathbf{f}(\mathbf{z})\,.$
(29)

This is also what one would expect if every material point of the condensate
is uniformly transported with the same velocity $\mathbf{v}$. Then, the
average driving force experienced by each material point is balanced by its
viscous friction with the surrounding medium.

However, we would like to stress that this relation, Eq. (29) ceases to be
valid when the enzyme currents are not uniform, such as in systems with more
that one spatial dimension. Specifically, we expect the error in Eq. (29) to
become larger when phase separation is weak, ${\Delta c\ll c_{+}}$. Figure 4
shows that naively using Eq. (29) to calculate the droplet velocity self-
consistently in response to the concentration profiles of substrates and
products, which themselves depend on the droplet velocity, is accurate only
under specific conditions. For instance, it is applicable in scenarios like
the 1D case discussed in Ref. [28], or it can be a good approximation when the
low concentration phase is almost depleted (${c_{-}/c_{+}\approx 0}$).

Interestingly, the simulations of the full dynamics show that self-propulsion
starts much earlier than would be predicted based on Eq. (29), which neglects
the inhomogeneity in the enzyme currents. Thus, liquid-like droplets can move
much faster than one would expect if each of their material points were to be
advected by a uniform velocity (that is, faster than solid-like condensates).
The reason for this is that the net movement of the condensate in terms of its
concentration profile does not require any actual mass transport of the
molecules in the highly concentrated phase over the same distance. This effect
can be attributed to the divergence-free circulation currents that transport
enzymes in the co-moving frame of the condensate: In the high-concentration
phase, enzymes are transported from the leading edge of the droplet to the
trailing edge, which retracts due to the local outflow of enzymes. Conversely,
due to mass conservation, these enzymes return in the low-concentration phase
from the trailing edge to the leading edge, which expands due to the
corresponding influx.

Figure 4:  Naive approach of calculating the droplet velocity. The top row
corresponds to droplets with strong phase separation, ${c_{-}=0.1\,c_{+}}$,
while the bottom row indicates droplets with weak phase separation,
${c_{-}=0.5\,c_{+}}$. Panels show droplet speed (color code) in the full
finite element simulations; data and parameters are identical to Fig. 3. Solid
lines indicate naive theoretical prediction of the phase boundary of the self-
propulsion instability, based on Eq. (29). This approximation gives acceptable
results when phase separation is strong (${c_{-}/c_{+}\sim 0}$, top row), but
fails in geometries whose dimension is larger than one if the enzymes are only
weakly depleted from the low-concentration phase (${c_{-}\sim c_{+}}$, bottom
row).

##  III Reciprocal interactions control type of self-propulsion instability

In our previous work [28], and in Figs. 1-4 we have focused on the limiting
case where interactions are weak and therefore the Flory-Huggins parameters
$\chi_{s,p}$ are small. In this limit, the terms proportional to $\Lambda$ in
the diffusion-reaction equations for the substrates and products, Eq. (7), can
be neglected, corresponding to a non-reciprocal limit. In the following we
will study the more general case ${\Lambda>0}$ and the resulting consequences
of the (partial) restoration of reciprocity. To keep our analysis analytically
tractable555We carry out the analytic calculations with the computer algebra
system Mathematica [58] and provide the corresponding notebooks in Ref. [59].,
we again mostly return to a one-dimensional (1D) system.

###  III.1 Substrate and product concentration profiles

####  III.1.1 Reciprocal interactions induce concentration jumps

A defining feature of taking the sharp interface limit is that the enzyme
concentration profile becomes piecewise constant. Therefore, one can split the
description of the steady state substrate and product concentration profiles
in the co-moving frame [Eq. (8)] into two simpler problems on the domains
inside and outside the droplet,

$\displaystyle 0$
$\displaystyle=\boldsymbol{\nabla}\cdot\left(\mathbf{v}\,s+D\,\boldsymbol{\nabla}s\right)-k_{1}\,c_{\pm}\,s+k_{2}\,p\,,$
(30a) $\displaystyle 0$
$\displaystyle=\boldsymbol{\nabla}\cdot\left(\mathbf{v}\,p+D\,\boldsymbol{\nabla}p\right)+k_{1}\,c_{\pm}\,s-k_{2}\,p\,,$
(30b)

where the reaction rates $k_{1}\,c_{\pm}$ are uniform inside and outside the
droplet, respectively. Importantly, because the gradients in enzyme
concentration vanish in both domains, the terms proportional to the
reciprocity parameter $\Lambda$ are absent in Eq. (30). How then will the
reciprocity parameter factor into the picture? This question will be answered
in the subsequent discussion, elucidating the boundary conditions that
establish a connection between the substrate and product concentration
profiles across the droplet interface.

To ensure particle number conservation of each species, the respective
particle fluxes must be continuous at the droplet interface. Akin to the
domain wall theory used in studies of the totally asymmetric simple exclusion
process [60, 61], this amounts to a balance between changes in diffusion and
mass flux:

$\displaystyle
D\bigl{[}\boldsymbol{\nabla}s|_{\text{in}}-\boldsymbol{\nabla}s|_{\text{out}}\bigr{]}\cdot\mathbf{\hat{e}}_{\perp}$
$\displaystyle=-\mathbf{v}\cdot\mathbf{\hat{e}}_{\perp}\,\bigl{[}s|_{\text{in}}-s|_{\text{out}}\bigr{]}\,,$
(31a) $\displaystyle
D\bigl{[}\boldsymbol{\nabla}p|_{\text{in}}-\boldsymbol{\nabla}p|_{\text{out}}\bigr{]}\cdot\mathbf{\hat{e}}_{\perp}$
$\displaystyle=-\mathbf{v}\cdot\mathbf{\hat{e}}_{\perp}\,\bigl{[}p|_{\text{in}}-p|_{\text{out}}\bigr{]}\,,$
(31b)

where $\mathbf{\hat{e}}_{\perp}$ denotes the unit vector normal to the droplet
interface. This relation can also be formally derived by integrating Eq. (8)
over an infinitesimal line segment which crosses the sharp droplet interface.
Here $\rvert_{\text{in}}$ and $\rvert_{\text{out}}$ indicate the inner and
outer side of the droplet interface, respectively. Thus, the substrate and
product concentration gradients on both sides of the droplet interface will
match when the droplet is at rest (${v=0}$) or when the concentration profiles
are continuous at the interface. However, as we show next, reciprocal
interactions (${\Lambda>0}$) will in general induce concentration jumps at the
interface.

The currents of substrates and products666In analogy to the enzyme currents
[Eq. (14b)], we added integration constants to ensure that the substrate and
product currents vanish far away from the droplet. in Eq. (8),

$\displaystyle\mathbf{J}_{s}$
$\displaystyle=-D\,\boldsymbol{\nabla}s-\Lambda\,s\,\chi_{s}\boldsymbol{\nabla}c-\mathbf{v}\,(s-s_{\infty})\,,$
(32a) $\displaystyle\mathbf{J}_{p}$
$\displaystyle=-D\,\boldsymbol{\nabla}p-\Lambda\,p\,\chi_{p}\boldsymbol{\nabla}c-\mathbf{v}\,(p-p_{\infty})\,,$
(32b)

must not only be continuous to conserve mass but also be finite-valued at the
droplet interface. This is not a given, because the concentration gradient
$\boldsymbol{\nabla}c$ is singular at the interface. To enforce boundedness
for the currents and mesoscopic velocities, we now integrate
$\mathbf{J_{s}}/s$ and $\mathbf{J_{p}}/p$ [cf. Eq. (32)] over an infinitesimal
line segment which crosses the sharp droplet interface. After calculating this
integral and taking the length of the line segment to zero, all except the
first two terms on the right-hand side of Eq. (32) vanish. This leads to the
boundary conditions

$\displaystyle\frac{s|_{\text{in}}}{s|_{\text{out}}}$
$\displaystyle=\exp\biggl{[}-\frac{\Lambda\,\chi_{s}\,\Delta
c}{D}\biggr{]}\,,$ (33a) $\displaystyle\frac{p|_{\text{in}}}{p|_{\text{out}}}$
$\displaystyle=\exp\biggl{[}-\frac{\Lambda\,\chi_{p}\,\Delta
c}{D}\biggr{]}\,.$ (33b)

Thus, for finite $\Lambda$, the substrate and product concentration profiles
in general exhibit a jump at the droplet interface.

When the interactions between enzymes, substrates, and products are
sufficiently weak (${\lvert\Lambda\,\chi_{s,p}\,\Delta c/D\rvert\ll 1}$), the
system reverts to the nonreciprocal limit discussed in our previous work [28].
In this limit, the concentration profiles of substrates and products at the
droplet interface become continuous. In contrast, strong attractive enzyme-
substrate and enzyme-product interactions ($\chi_{s,p}\ll 0$) will lead to an
exponential enrichment of both substrates and products inside of the
condensate. We expect droplet motion to cease once the redistribution of
substrates and products driven by these reciprocal interactions outweighs the
redistribution caused by reactions and diffusion, i.e., when the system is
close enough to thermodynamic equilibrium.

Taken together, conservation of mass dictates that the net flows of substrates
and products at the droplet interface are continuous and finite, and amounts
to Robin boundary conditions specified in Eqs. (31) and (33). This allows us
to revisit our previous analysis [28] for the more general case of
${\Lambda\geq 0}$.

To that end, one could directly solve Eq. (30) in each subdomain (inside and
outside the condensate), connect these solutions at the droplet interface by
using Eqs. (31) and (33), and impose no-flux conditions in the far field.
However, one must also ensure that the solution conserves the total mass
$n_{\infty}$ of substrates and products:

$n_{\infty}=\langle
s+p\rangle=\frac{1}{|\Omega|}\int_{\Omega}\\!d^{d}z\,[s(\mathbf{z})+p(\mathbf{z})]\,,$
(34)

where $|\Omega|$ is the total volume of the integration domain. For an open
domain, where $|\Omega|\to\infty$, this implies that

$\lim_{|\mathbf{z}|\to\pm\infty}[s(\mathbf{z})+p(\mathbf{z})]=n_{\infty}\,,$
(35)

far away from the droplet. At first glance, this additional constraint may
seem trivial because the total concentration of substrates and products,
${n(\mathbf{x})\coloneqq s(\mathbf{x})+p(\mathbf{x})}$, is spatially
uniform777Note that for ${\Lambda=0}$, the total concentration follows a
diffusion equation. for ${\Lambda=0}$. However, this is in general not the
case for ${\Lambda>0}$ and thus requires separate consideration. Specifically,
as we will shown next, a moving droplet can spatially redistribute reactant
mass ${n(\mathbf{x},t)}$ in the laboratory frame. Such mass redistribution is
a key feature and control process of biochemical pattern-forming systems [62,
63, 64].

####  III.1.2 Redistribution of total reactant mass by moving droplets

To analyze the profile of the total concentration of substrates and products,
we add Eqs. (30a) and (30b), which gives

$0=\boldsymbol{\nabla}\cdot\left(n\,\mathbf{v}+D\,\boldsymbol{\nabla}n\right)\,.$
(36)

Because of mass conservation, the reaction terms have cancelled out implying
that, in the co-moving frame, the sum of the advective and diffusive mass
fluxes of the total amount of reactants is divergence-free. Moreover, far away
from the droplet, the total reactant concentration approaches the far-field
value, ${\lim_{|\mathbf{z}|\to\infty}n(\mathbf{z})=n_{\infty}}$, and becomes
homogeneous,
${\lim_{|\mathbf{z}|\to\infty}\boldsymbol{\nabla}n(\mathbf{z})=0}$. Finally,
note that Eq. (36) is valid both inside and outside of the condensate, but not
at the droplet interface where $c(\mathbf{z})$ is singular. Thus, each of
these domains need to be analyzed separately and then connected by the Robin
boundary conditions derived previously. If the substrate and product
concentration fields are smooth and continuous (${\Lambda=0}$), then Eq. (36)
has only ${n(\mathbf{z})=n_{\infty}}$ as solution. In general, however, the
substrate and product concentration fields are neither smooth nor continuous
(${\Lambda>0}$). Since this complicates the analysis considerably, we focus on
1D systems in the following.

In a 1D geometry, the two droplet interfaces define the boundaries between
three spatial subdomains that we label with indices ${i\in\\{-,0,+\\}}$, which
represent the solvent to the left of the droplet ($-$), the droplet ($0$), and
the solvent to the right of the droplet ($+$). Solving Eq. (36), gives for the
total concentration of substrates and products in each subdomain

$n(z)=C_{i}+\Delta C_{i}\,\exp(-vz/D)\,,$ (37)

where $C_{i}$ and $\Delta C_{i}$ are integration constants. To determine these
constants, we will now and in the following, without loss of generality,
assume that ${v\geq 0}$. Because the concentration profiles must remain finite
in the far field ${z\to\pm\infty}$, one has ${\Delta C_{-}=0}$ for the solvent
domain left to the droplet. Furthermore, the far-field conditions [Eq. (35)]
imply that ${C_{\pm}=n_{\infty}}$. By summing Eqs. (31a) and (31b), one finds
Robin boundary conditions for the total reactant concentration at the two
droplet interfaces at ${z=\pm R}$, which lead to the conclusion:
${C_{0}=C_{\pm}=n_{\infty}}$. Taken together, the total reactant concentration
is given by

$n(z)=n_{\infty}+\exp(-vz/D)\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color@gray@fill{.5}\times\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\begin{cases}0\,,&z<-R\,,\\\
\Delta C_{0}\,,&|z|\leq R\,,\\\ \Delta C_{+}\,,&z>R\,.\end{cases}$ (38)

The remaining integration constants, $\Delta C_{0}$ and $\Delta C_{+}$, cannot
be determined without taking into account the concentration jumps of the
substrates and products, Eq. (33). Thus, they require resolving the
corresponding concentration profiles, as will be discussed in the following.

####  III.1.3 Inhomogeneous Helmholtz equation determines concentration
profiles

Using mass conservation, ${n(z)=s(z)+p(z)}$, to eliminate the concentration of
products from Eq. (30), we arrive at an inhomogeneous Helmholtz equation with
advection which describes the concentration profile of substrates in each
subdomain,

$0=\partial_{z}\left(v\,s+D\,\partial_{z}s\right)-(k_{1}\,c_{\pm}+k_{2})\,s+k_{2}\,n\,.$
(39)

Recall that $c_{\pm}$ refer to the enzyme concentrations inside and outside
the droplet, respectively. Also note that the total reactant concentration
$n(z)$ is spatially inhomogeneous, as specified by Eq. (38). The boundary
conditions for the substrate profiles at the droplet interfaces are given by
Eq. (31a) and Eq. (33a). These allow us to determine the distribution of
substrates $s(z)$ and the total reactant concentration profile $n(z)$ up to
two constants $\Delta C_{0}$ and $\Delta C_{+}$. To now constrain these two
constants, we turn to the distribution of products, $p(z)=n(z)-s(z)$. In
particular, the product concentration jumps at the two droplet interfaces,
defined by Eq. (33b), specify the remaining two constants $\Delta C_{0}$ and
$\Delta C_{+}$. By solving888Specifically, we used the computer algebra system
Mathematica [58] to solve Eq. (39) within each subdomain; the full expressions
can be found in Ref. [59]. Eq. (39), we determined the concentration profiles
for substrates and products when the droplet was either moving (${v>0}$) or
stationary (${v=0}$).

Our theoretical results for the concentration profiles align closely with our
simulations, as shown in Fig. 5 for a stationary droplet consisting of
enzymes. As we discussed in Sec. I.2, near the center of the droplet, enzyme-
catalyzed reactions increase the concentration of products at the expense of
substrates. These reaction-induced concentration gradients correspond to
particle fluxes which bring substrates towards the condensate and transport
products away. In addition, attractive enzyme-substrate and enzyme-product
interactions increase the concentration of both substrates and products by a
discontinunous concentration jump at the droplet interface [Eq. (33)]. In the
following sections, we will explore the implications of these concentration
profiles on droplet self-propulsion.

Figure 5:  Stationary concentration profiles for enzymes (gray), substrates
(red), and products (blue) of a single droplet in a finite-sized 1D domain
$[-L,L]$. Dots correspond to FEM simulation results, while solid lines are
analytical solutions for the substrate and product concentration profiles in
the sharp interface limit (gray shading indicates droplet). Attractive enzyme-
substrate and enzyme-product interactions, ${\chi_{s}=-0.05r}$,
${\chi_{p}=-0.01r}$, in the presence of reciprocal interactions,
${\Lambda=20D/\epsilon_{0}}$, lead to a concentration jump in substrates and
products. The other parameters are given by ${c_{-}=0.1c_{+}}$,
${w=0.1l_{0}}$, ${R=l_{0}}$, ${L=5l_{0}}$, ${M=100D/\epsilon_{0}}$,
${k_{1}=k_{2}/c_{+}}$, and ${\langle s+p\rangle=c_{+}}$.

###  III.2 Reciprocal interactions can inhibit droplet self-propulsion

Equation (33) demonstrates that substrates and products can be enriched inside
of a droplet simply through attractive interactions with enzymes, without
relying on non-equilibrium reactions. As shown in Fig. 5, this can lead to a
scenario where the concentration of substrates is actually higher in the
condensate than in the surrounding solution, despite being consumed by
enzymatically catalyzed reactions. To illustrate what this means for droplet
self-propulsion, we consider an enzyme condensate which moves (or is pulled)
with a prescribed velocity through an open domain. Moreover, to disentangle
the effects of reciprocal interactions from those of non-equilibrium
reactions, we compare the thermodynamic limit where reactions are absent,
${k_{1}=k_{2}=0}$ and ${\Lambda>0}$, to the nonreciprocal limit, ${k_{1,2}>0}$
and ${\Lambda=0}$, cf. Fig. 6. In the nonreciprocal limit, as shown in our
previous work [28], one finds that products are accumulated and substrates are
depleted at the trailing edge of the droplet. When enzyme-substrate attraction
is stronger than enzyme-product attraction, this can push the droplet towards
its leading edge and thereby drive self-propelled motion. In contrast, in the
thermodynamic limit, products are absent and attractive enzyme-substrate
interactions cause accumulation of substrates at the trailing edge of the
condensate. Thus, the droplet will be pulled back towards its trailing edge,
and condensate motion cannot be sustained. Based on these competing
mechanisms, we expect that reciprocal interactions can potentially inhibit
reaction-induced droplet self-propulsion.

Figure 6:  Concentration profiles for substrates (red), and products (blue)
for a single droplet (gray) which is moving at a prescribed velocity through
an open domain. In the nonreciprocal case, substrates are depleted at the
trailing edge of the condensate, while products are enriched there. In
contrast, in the thermodynamic limit, substrates are enriched at the trailing
edge while products are not produced and are thereby absent from the system.
Parameters unless specified otherwise: ${\chi_{s}=-0.05r}$,
${\chi_{p}=-0.01r}$, ${\Lambda=4D/\epsilon_{0}}$, ${c_{-}=0.1c_{+}}$,
${R=l_{0}}$, ${k_{1}=k_{2}/c_{+}}$, and ${\langle s+p\rangle=c_{+}}$.

Having determined the substrate and product concentration profiles, we now
quantify the effect of reciprocal interactions on the self-propulsion
instability. As we consider a 1D system, one could directly use Eq. (29) to
obtain the droplet velocity in response to any inhomogeneous distribution of
substrates and products, as we demonstrated recently [28],

$v=-\frac{Mc_{+}}{2R\Delta c}\big{[}\chi_{s}\,\Delta s(v)+\chi_{p}\,\Delta
p(v)\big{]}\,.$ (40)

Here $\Delta s(v)$ and $\Delta p(v)$ are the substrate and product
concentration differences, respectively, between the right and the left
boundary of the moving droplet. However, from this approach it does not
immediately become clear which of the self-consistent solutions are the stable
attractors for the dynamics. To close this gap in understanding, we follow a
slightly different route by exploiting the framework presented in Sec. II,
which not only yields the fixed points but also informs about their stability.

To that end, we compute the chemical potential imbalance $\Delta\mu_{0}(v)$,
as defined by Eq. (16) and illustrated in Fig. 2(b) for a 3D droplet. For a 1D
droplet, as discussed in Sec. II.4, our theory simplifies considerably because
the divergence-free enzyme currents must vanish in the far-field and thus in
the entire 1D domain, $J_{\Delta}=0$. By projecting Eq. (27) on the axis
$\mathbf{\hat{e}}_{v}$, substituting the definition of the chemical potential
imbalance [Eq. (16)] and solving for the latter, one finds

$\Delta\mu_{0}(v)=\int_{\mathcal{D}}\\!dz\,\bigl{(}\mathbf{\hat{e}}_{v}\cdot\mathbf{f}\bigr{)}-\frac{V\,\Delta
c}{M(c_{+})}\,v\,,$ (41)

where ${\Delta c=c_{+}-c_{-}}$ and ${V=2R}$ is the droplet volume. Finally,
substituting the thermodynamic force field [Eq. (9)] shows that the chemical
potential imbalance for a 1D droplet is given by

$\Delta\mu_{0}(v)=-\chi_{s}\,\Delta s(v)-\chi_{p}\,\Delta p(v)-\frac{2R\Delta
c}{Mc_{+}}\,v\,,$ (42)

where we have assumed ${M(c_{+})=Mc_{+}}$ with constant $M$. Note that Eqs.
(42) and (40) are equivalent when the chemical potential imbalance vanishes,
which corresponds to the fixed points of the dynamics. As discussed in Sec.
II.2, unstable (stable) fixed points are characterized by a positive
(negative) slope of the chemical potential imbalance as a function of the
velocity [Fig. 2(b)]. Thus, the chemical potential imbalance contains all
information necessary to analyze the phase space flow.

Note that there is a subtlety when evaluating the concentration differences of
the substrates $\Delta s(v)$ and the products $\Delta p(v)$ between the two
droplet interfaces. Since we are in the sharp interface limit, it is not
entirely clear when evaluating the integrals

$\displaystyle\Delta s(v)$
$\displaystyle\coloneqq\int_{\mathcal{D}}\\!dz\,\partial_{z}s(z)\,,$ (43a)
$\displaystyle\Delta p(v)$
$\displaystyle\coloneqq\int_{\mathcal{D}}\\!dz\,\partial_{z}p(z)\,,$ (43b)

which value we should choose for the lower and upper limits. These could be
the values at the inner or outer boundaries or an interpolation between the
two. The difference between these choices arises from the discontinuities in
the concentration profiles which, taking the substrate concentrations at the
right droplet interface as an example, imply ${\lim_{\varepsilon\to
0}s(R-\varepsilon)\neq\lim_{\varepsilon\to 0}s(R+\varepsilon)}$. Currently, we
cannot resolve which is the correct choice solely based on theoretical
arguments. Hence, we proceeded heuristically and compared the various choices
with simulation data for the droplet speed as a function of the enzyme
mobility and the reaction rates. It turned out that these simulation data are
only compatible with choosing a droplet domain that excludes its boundary,
${\mathcal{D}=\\{z:|z|<R\\}}$, which for the enzyme concentration profile
implies ${c(z)=c_{+}\,\forall\,|z|<R}$ else $c_{-}$. In this case, one takes
the concentration values at the inner sides of the droplet interfaces when
evaluating Eq. (43). This ambiguity could be resolved in future studies by
using an enzyme concentration profile with a smooth interface. Notwithstanding
these limitations of the sharp interface approximation, as we show below, it
quantitatively predicts the droplet velocity in agreement with our simulation
results.

###  III.3 Discussion and test of the sharp interface theory

Given the limitations of the sharp interface approximation just discussed, it
is crucial to verify its predictions against full FEM simulations, which
feature dynamic and smooth concentration profiles of substrates and products
[Eq. (7)] as well as enzymes [Eq. (4)]. In this section, after quantitatively
comparing theory and simulations, we also discuss central features of the
bifurcation which leads to the onset of droplet self-propulsion.

To predict the speed of the biomolecular condensates with the sharp interface
theory, we proceed as outlined in Sec. II.3 for the nonreciprocal case. First,
we determine the chemical potential imbalance $\Delta\mu_{0}(v)$ [Eq. (42)] as
a function of the droplet velocity. Then, we select the maximal droplet
velocity $v^{\star}$ which satisfies the thermodynamic consistency criterion
${\Delta\mu_{0}(v^{\star})=0}$ [Eq. (17)]. For larger droplet velocities, the
chemical potential imbalance is a monotonically decreasing and strictly
negative function [Eq. (42)],

$\Delta\mu_{0}(v)\approx-\frac{2R\Delta c}{Mc_{+}}\,(v-v^{\star})\,.$ (44)

This implies ${\partial_{v}\Delta\mu_{0}(v^{\star})<0}$ and, following the
discussion in Sec. (II.2), that higher droplet velocities ${v>v^{\star}}$
would require additional energy input ${-\partial_{t}\mathcal{F}_{0}>0}$ [Eq.
(24)]. Therefore, we conclude that the velocity $v^{\star}$ must correspond to
a stable steady state. A more detailed discussion of the bifurcation diagram
is deferred to a later point of this section.

As shown in Fig. 7, we find good agreement of the velocities and shape of the
phase diagram predicted by the sharp interface theory with our full
simulations. As in the nonreciprocal case, see Appendix A, the predictions
become less accurate for large reaction rates, likely because the diffusion
length ${l_{\pm}=\sqrt{D/(k_{1}c_{\pm}+k_{2})}}$ becomes shorter compared to
the width of the interface. This parameter regime conflicts with the sharp
interface limit, which requires that the interface width is small compared to
all other relevant length scales.

Figure 7: Speed of a self-propelled 1D droplet as a function of the enzyme
mobility $M$ and the reaction rate $k_{1}$. The first panel shows the
predicted droplet speed $v^{\star}$ of the stable steady state (blue color
scale), while the second panel shows the unstable branch $v^{\circ}$ (red
color scale), from the self-consistent sharp interface theory [Eqs. (42) and
(17)]. The third panel shows the measured droplet speed $v_{\text{sim}}$ (blue
color scale) in the full simulations [Eqs. (4) and (7)]. The fourth panel
shows the percentile error of the prediction,
${\Delta\coloneqq(v^{\star}-v_{\text{sim}})/v_{\text{sim}}}$ (blue-orange
color scale). Parameters: ${\chi_{s}=-0.05\,r}$, ${\chi_{p}=-0.01\,r}$,
${\Lambda=4D/\epsilon_{0}}$, ${c_{-}=0.1\,c_{+}}$, ${R=l_{0}}$, and ${\langle
s+p\rangle=c_{+}}$. The additional parameters required for the full FEM
simulations are given by ${w=0.1\,l_{0}}$, and ${L=30\,l_{0}}$.

Strikingly, unlike in the nonreciprocal case (${\Lambda=0}$), in which the
droplet velocity increased continuously from ${v=0}$ upon reaching a critical
enzyme mobility or reaction rate [Fig. 3], a sudden jump in the droplet speed
from ${v=0}$ to a finite value is observed for ${\Lambda>0}$ [Fig. 7]. This
points towards a subcritical bifurcation. The bifurcation diagram of a
subcritical transition is characterized by the stable branch ${v=v^{\star}}$,
the trivial branch ${v=0}$, which can change its stability depending on the
parameters, and an unstable branch ${v=v^{\circ}}$ with ${0\leq v^{\circ}\leq
v^{\star}}$. The binodal phase boundary is defined by the appearance of the
stable branch ${v=v^{\star}}$, which is shown in the first panel of Fig. 7. In
contrast, the spinodal phase boundary is characterized by the disappearance of
the unstable branch ${v=v^{\circ}}$ at the intersection ${v^{\circ}=0}$. This
unstable branch is present within the parameter space examined in Fig. 7, as
shown in the second panel, i.e., we are in the binodal regime. To then observe
droplet self-propulsion in this bistable regime between the binodal and the
spinodal line, one needs to provide a sufficiently large perturbation. To
provide such a perturbation in our simulations, we initialized the
concentration profiles of substrates and products with added noise.

In the previous paragraph, we mapped out and discussed the phase diagram of
the droplet velocity as a function of the enzyme mobility and reaction rates.
Next, we aim to better understand the role of the mobilities of substrates and
products, ${\Lambda(s)=\Lambda s}$ and ${\Lambda(p)=\Lambda p}$. To that end,
we now keep the reaction rates fixed while varying the enzyme mobility $M$
[for reference, a slice along the dashed arrows in the first two panels of
Fig. 7 corresponds to the bifurcation diagram depicted in the highlighted
panel in Fig. 8] and the reciprocity parameter $\Lambda$.

Figure 8: Droplet speed $v^{\star}$ (blue color scale) predicted by the sharp
interface theory, as a function of the enzyme mobility $M$ and the mobility
$\Lambda$ of substrates and products in response to a gradient in enzymes.
Bifurcation diagrams on the right correspond to slices along the dashed arrows
in the left panel, for several selected values of $\Lambda$. Highlighted panel
corresponds to slice along the dashed arrow in Fig. 7. Solid lines indicate
stable branches, whereas dashed lines indicate branches of unstable fixed
points. Parameters: ${\chi_{s}=-0.05\,r}$, ${\chi_{p}=-0.01\,r}$,
${c_{-}=0.1c_{+}}$, ${R=l_{0}}$, ${k_{1}=k_{2}/c_{+}}$, and ${\langle
s+p\rangle=c_{+}}$.

First, we observe that the transition to self-propulsion is supercritical
(continuous) for small substrate and product mobilities ${\Lambda=0}$, as
shown in the first panel on the right side of Fig. 8. For general
${\Lambda>0}$, condensate self-propulsion sets in through a subcritical,
imperfect pitchfork bifurcation. Moreover, consistent with our arguments in
the previous section, we observe that an increase in $\Lambda$ reduces the
size of the parameter region where droplet self-propulsion can be observed.
Specifically, our findings indicate the existence of an upper bound for
$\Lambda$, below which droplets display self-propulsion. By using the
Einstein-Smoluchowsky relation ${D=\Lambda\,k_{B}T}$ to relate the mobility of
substrates and products to their diffusion coefficient in the case of thermal
Brownian motion [52], the vertical axis in Fig. 8 can be read as
${\Lambda\epsilon_{0}/D=\epsilon_{0}/(k_{B}T)}$. This suggests that for the
parameters used in Fig. 8, the energy scale for the enzyme-substrate and
enzyme-product interactions must be no more than ${5k_{B}T}$ to observe
droplet self-propulsion. The results in Fig. 8 also show that the mobility of
substrates and products must be small compared to the mobility of enzymes.

Given these observations, what are necessary conditions to observe droplet
self-propulsion for large substrate and product mobilities ${\Lambda}$ and for
strong enzyme-substrate and enzyme-product interactions? To answer this
question, in the following we use the sharp interface theory to further
elucidate how ${\Lambda}$ affects droplet self-propulsion and, moreover, when
the transition is subcritical (discontinuous) or supercritical (continuous).

###  III.4 Analysis of the role of reciprocity

Before we continue with our analysis, we first summarize a few key results
that will be used in the following. In our previous work [28] we established
the mobility of the enzymes $M$ as a key control parameter for observing
droplet self-propulsion. We have found that in the absence of reciprocal
interactions (${\Lambda=0}$), self-propulsion always sets in if the mobility
is sufficiently large, cf. Fig. 3 and Fig. 8. This can be understood from Eq.
(42), which we rewrite as

$\Delta\mu_{0}(v)=\Delta\mu^{\star}_{0}(v)-\frac{2R\Delta c}{Mc_{+}}\,v\,.$
(45)

This chemical potential imbalance has an upper bound given by the limit
$M\to\infty$, defining the maximal chemical potential imbalance

$\Delta\mu^{\star}_{0}(v)\coloneqq\lim_{M\rightarrow\infty}\Delta\mu_{0}(v)=-\chi_{s}\Delta
s(v)-\chi_{p}\Delta p(v)\,.$ (46)

As discussed in Sec. II.2, steady states of the condensate dynamics are
defined by the roots of Eq. (45). Moreover, the thermodynamic housekeeping
analysis in Sec. II.2 showed that stable steady states are characterized by a
negative slope ${\partial_{v}\Delta\mu_{0}(v^{\star})<0}$ whereas unstable
states have a positive slope ${\partial_{v}\Delta\mu_{0}(v^{\star})>0}$.
Geometrically, this means that—at stable steady states—the curve
$\Delta\mu^{\star}_{0}(v)$ intersects the line ${(2R\Delta c)/(Mc_{+})\,v}$
from above. Finally, as discussed in the previous section, the largest root of
Eq. (45) must correspond to a stable steady state. Therefore, one can map out
the bifurcation diagram simply by determining the roots.

The above criteria to observe condensate self-propulsion can only be fulfilled
if the maximal chemical potential imbalance $\Delta\mu^{\star}_{0}(v)$ [Eq.
(46)] has a band of positive (negative) values for positive (negative)
velocities. In that case, as mentioned at the beginning in this section, Eq.
(45) is guaranteed to have non-trivial roots for ${v\neq 0}$ if the enzyme
mobility is sufficiently high,

$M>\frac{2R\Delta c}{c_{+}}\,\min_{v}\frac{v}{\Delta\mu^{\star}_{0}(v)}>0\,.$
(47)

In the following, we exploit these simple rules to improve our understanding
of the onset of condensate self-propulsion. Because the system does not have a
preferred spatial direction, without loss of generality, we focus on positive
droplet velocities ${v\geq 0}$.

####  III.4.1 Weakly reciprocal droplets

First, we consider a scenario where the reciprocity parameter $\Lambda$ is
small. To that end, we expand the maximal chemical potential imbalance to
first order in $\Lambda\chi_{s,p}\Delta c/D$,

$\Delta\mu^{\star}_{0}(v)\approx\Delta\mu^{\star}_{0}(v)\bigr{\rvert}_{\Lambda=0}+\partial_{\Lambda}\Delta\mu^{\star}_{0}(v)\bigr{\rvert}_{\Lambda=0}\,\Lambda\,,$
(48)

which implicitly assumes weak enzyme-substrate and enzyme-product
interactions. Based on the discussion in the previous section, droplet self-
propulsion is only possible if the maximal chemical potential imbalance
$\Delta\mu^{\star}_{0}(v)$ has a band of positive values for ${v>0}$. To gain
a clearer understanding of when this condition holds, we will next discuss the
characteristic features of the two leading-order contributions to the series
expansion in Eq. (48).

Figure 9: Maximal chemical potential imbalance $\Delta\mu_{0}^{\star}(v)$ [Eq.
(46)] as a function of the droplet speed $v$, for different values of the
parameter $\Lambda$. In the limiting case ${M\to\infty}$, steady states are
defined by the roots of the maximal chemical potential imbalance. Arrows
indicate how the chemical potential imbalance affects the droplet speed $v$.
Stable (filled circles) and unstable (open circles) steady states are
indicated for ${\Lambda=4D/\epsilon_{0}}$. Color code of curves is the same as
in Fig. 8 and in Fig. 7. The dashed red line shows how the maximal chemical
potential imbalance is affected by the parameter $\Lambda$. Parameters:
${\chi_{s}=-0.05\,r}$, ${\chi_{p}=-0.01\,r}$, ${c_{-}=0.1\,c_{+}}$,
${R=l_{0}}$, ${k_{1}=k_{2}/c_{+}}$, and ${\langle s+p\rangle=c_{+}}$.

_Leading order: nonreciprocal limit._ The first term in Eq. (48),
${\Delta\mu^{\star}_{0}(v)\rvert_{\Lambda=0}}$, recapitulates the limit of
nonreciprocal droplets [black curve in Fig. 9], which we have discussed in our
previous work [28]. In this limit, an example of the substrate and product
concentration profiles, ${s(z)\rvert_{\Lambda=0}}$ and
${p(z)\rvert_{\Lambda=0}}$, is depicted in the left panel of Fig. 6. As these
concentration profiles show, reactions and diffusion lead to depletion of
substrates and enrichment of products at the trailing edge of the moving
condensate, in comparison to its leading edge. Thus, considering Eq. (46),
${\Delta\mu^{\star}_{0}(v)\rvert_{\Lambda=0}}$ must be strictly positive for
${v>0}$ when enzyme-substrate attraction dominates over enzyme-product
interactions. As we discussed in the previous section, the existence of a band
of positive values ${\Delta\mu^{\star}_{0}(v)\rvert_{\Lambda=0}>0}$ for
${v>0}$ fulfills the necessary condition for droplet motion. Moreover, the
positive slope of the curve ${\Delta\mu^{\star}_{0}(v)\rvert_{\Lambda=0}>0}$
at the fixed point ${v=0}$ is indicative of an unstable state. Hence, we
conclude that the leading-order term in the series expansion of Eq. (48)
promotes droplet self-propulsion.

It is important to note that this leading-order term is bound from above for
positive droplet velocities. This can be deduced from the fact that
${\Delta\mu^{\star}_{0}(v)\rvert_{\Lambda=0}}$ must vanish in the following
two limiting cases. For ${v=0}$, the chemical potential imbalance vanishes
because the system is symmetric under parity transformations. In the limit
${v\to\infty}$, the substrate and product concentration profiles remain flat,
as the time required to traverse one droplet diameter, $2R/v$, is much shorter
than the timescales of reactions, $\tau_{0}$, and diffusion, $R^{2}/D$, so
that the densities equilibrate and become uniform. Since the leading-order
term ${\Delta\mu^{\star}_{0}(v)\rvert_{\Lambda=0}}$ is strictly positive for
positive velocities ${v>0}$, it follows that it must have a maximum for some
finite value ${v>0}$ [black curve in Fig. 9]. As a consequence, the correction
term in the series expansion of Eq. (48) can exceed this first-order
contribution.

_Reciprocal correction._ The second term in Eq. (48),
${\partial_{\Lambda}\Delta\mu^{\star}_{0}(v)\rvert_{\Lambda=0}}$, is a
correction that accounts for weak reciprocity [dashed red curve in Fig. 9].
The corresponding changes in the substrate and product999Note that in the
equilibrium scenario depicted in Fig. 6, only substrates are present in the
system. concentration profiles can be understood from the right panel of Fig.
6: Their attractive interactions with enzymes, ${\chi_{s,p}<0}$, will favor
influx at the leading edge but hinder outflux of substrates and products at
the trailing edge of the droplet. Hence, substrates and products will
accumulate at the trailing edge of the condensate which acts as a moving
potential barrier101010Repulsive attractions would instead cause accumulation
at the leading edge.. These concentration profiles, considering Eq. (46), lead
to ${\partial_{\Lambda}\Delta\mu^{\star}_{0}(v)\rvert_{\Lambda=0}}$ being
strictly negative and decreasing monotonically with increasing droplet speed,
as more substrates and products accumulate at the trailing edge of the
condensate. Therefore, with increasing $\Lambda$, the correction term can
cause the maximal chemical potential imbalance [Eq. (48)] to become strictly
negative for ${v>0}$, leaving only a single stable fixed point ${v=0}$ [solid
red curve in Fig. 9]. In that case, droplet self-propulsion is completely
suppressed.

_Bistability and criticality._ In section III.3, we discussed that for
${\Lambda>0}$, self-propulsion typically sets in through a subcritical
bifurcation when varying the enzyme mobility $M$, and that there is a bistable
region. In the following, we use our framework to further elucidate these
phenomena. Recall that the condition for a stable steady state is that the
curve $\Delta\mu^{\star}_{0}(v)$ cuts the line ${(2R\Delta c)/(Mc_{+})\,v}$
from above. To observe a supercritical bifurcation when varying the enzyme
mobility $M$, the fixed point ${v=0}$ must lose its stability concomitantly
with the emergence of a new branch of stable fixed points for ${v\geq 0}$. An
important feature of a supercritical transition is that the velocity does not
suddenly jump to a finite value when varying the enzyme mobility $M$ for
example, but instead continuously increases from ${v=0}$. Geometrically, this
is only possible if the maximal chemical potential imbalance
$\Delta\mu^{\star}_{0}(v)$ grows but curves downward,
${\partial_{v}^{2}\Delta\mu^{\star}_{0}(v)\rvert_{v=0}\leq 0}$, when the
velocity is increased starting from the fixed point ${v=0}$. Using
Mathematica, we found that the leading order term in the series expansion in
Eq. (48) always has vanishing curvature,
${\partial_{v}^{2}\Delta\mu^{\star}_{0}(v)\rvert_{\Lambda=0,v=0}=0}$. While we
were not able to completely map out the correction term, we found that it
curves upward as a function of the droplet velocity, for small droplet radii
$R$ and for the parameters studied in Fig. 10. This further supports the
conclusion that for ${\Lambda>0}$, in a wide range of parameters, the onset of
condensate self-propulsion is subcritical.

For sufficiently high enzyme mobility $M$, a subcritical bifurcation occurs
when the fixed point ${v=0}$ becomes unstable. This is only possible if the
maximal chemical potential imbalance $\Delta\mu^{\star}_{0}(v)$ has a positive
slope at the fixed point ${v=0}$. Using the series expansion in Eq. (48), this
leads to

$\Lambda<-\frac{\partial_{v}\Delta\mu^{\star}_{0}(v)\bigr{\rvert}_{\Lambda=0,v=0}}{\partial_{v}\partial_{\Lambda}\Delta\mu^{\star}_{0}(v)\bigr{\rvert}_{\Lambda=0,v=0}}\,,$
(49)

where we have used the fact that
${\partial_{\Lambda}\Delta\mu^{\star}_{0}(v)}\leq 0$. In Fig. 10, the dark
blue region shows when this criterion is fulfilled for different values of the
reciprocity parameter $\Lambda$ and the condensate radius $R$. Even if the
slope criterion Eq. (49) is not fulfilled, the system can still be bistable as
we explain next. More specifically, to permit condensate self-propulsion for
sufficiently high enzyme mobility, it suffices if the maximal chemical
potential imbalance $\Delta\mu^{\star}_{0}(v)$ has a positive maximum [cf.
light blue curve in Fig. 9],

$\max_{v>0}\Delta\mu^{\star}_{0}(v)>0\,,$ (50)

for positive velocities. In Fig. 10, the light blue region shows when
criterion Eq. (50) is fulfilled for different values of the reciprocity
parameter $\Lambda$ and the condensate radius $R$. The inequality (50) is
automatically satisfied when the maximum of the leading-order term in Eq. (48)
is larger than the saturation value of the correction, which leads to the
weaker constraint

$\Lambda\lesssim-\frac{\max_{v}\Delta\mu^{\star}_{0}(v)\bigr{\rvert}_{\Lambda=0}}{\lim_{v\rightarrow\infty}\partial_{\Lambda}\Delta\mu^{\star}_{0}(v)\bigr{\rvert}_{\Lambda=0}}\,.$
(51)

Here, we have used the fact that
${\partial_{\Lambda}\Delta\mu^{\star}_{0}(v)\leq 0}$. This approximation gives
a lower bound for the critical value of $\Lambda$ below which self-propulsion
is possible [black dotted line in Fig. 10].

In summary, we have explored how the magnitude of the reciprocal interactions,
quantified by $\Lambda$, controls the onset of condensate self-propulsion. In
the bistable region, it requires a sufficiently large perturbation of the
system to excite the condensates into a self-propelling state. Importantly,
all of our results so far suggest that the reciprocity parameter $\Lambda$,
and hence enzyme-substrate and enzyme-product interactions, must be weak to
observe condensate motion. How can we relax this constraint?

Figure 10: Phase diagram illustrating how droplet self-propulsion depends on
the condensate radius $R$ and on the mobility of substrates and products in
response to a gradient in enzymes, ${\Lambda(s)=\Lambda s}$ and
${\Lambda(p)=\Lambda p}$, in the sharp interface approximation. The light blue
region shows the parameter regime where droplet self-propulsion is bistable
[Eq. (50)]. The corresponding phase boundary can be crudely approximated by
Eq. (51), as shown by the black dotted line. The dark blue region shows the
spinodal parameter regime where only self-propelling condensate states are
stable [Eq. (49)]. For simplicity, we consider the limit $M\to\infty$.
Parameters: ${\chi_{s}=-0.05\,r}$, ${\chi_{p}=-0.01\,r}$,
${c_{-}=0.1\,c_{+}}$, ${k_{1}=k_{2}/c_{+}}$, and ${\langle s+p\rangle=c_{+}}$.

####  III.4.2 Droplets in the strongly reciprocal regime

So far, we have shown that droplet self-propulsion can be inhibited if the
reciprocity parameter $\Lambda$ is sufficiently large. Then, is it possible at
all to observe droplet self-propulsion when
${\Lambda\,\epsilon_{0}/D\rightarrow\infty}$? This would correspond to a limit
at which enzyme-substrate and enzyme-product interactions are strong, so that
the effects due to reciprocity become very important. To understand the
particular significance of this case, we again use the Einstein-Smoluchowsky
relation ${D=\Lambda\,k_{B}T}$ to relate the mobility of substrates and
products to their diffusion coefficient in the case of thermal Brownian motion
[52]. The horizontal axis in Fig. 8 can then be interpreted as
${M\epsilon_{0}/D=(M/\Lambda)\,\epsilon_{0}/(k_{B}T)}$ and the vertical axis
as ${\Lambda\epsilon_{0}/D=\epsilon_{0}/(k_{B}T)}$. Thus, for the parameters
discussed so far, weak enzyme-substrate and enzyme-product interactions, with
energy scale ${\epsilon_{0}\sim k_{B}T}$, would necessitate a large ratio of
mobilities ${(M/\Lambda)}$ to induce droplet self-propulsion. This can be the
case, for example, if the enzymes correspond to small proteins (such as
transcription factors and components of the transcription machinery) and the
products correspond to large RNA molecules111111To treat this scenario, one
would need to relax the assumption of substrates and products having equal
diffusion coefficients. More specifically, the produced RNA polymer would have
much smaller diffusion coefficient than the substrate nucleotides.,121212If
the enzymes are proteins that cleave RNA polymers, then the products would
have higher diffusion coefficient than the substrates..

To explore if droplet self-propulsion is possible for a small ratio of
mobilities ${(M/\Lambda)}$, we will now study the extreme limit
${\Lambda\,\epsilon_{0}/D\rightarrow\infty}$. So far, we have focused on a
scenario where both substrates and products show attractive interactions with
enzymes, with ${\chi_{s}<\chi_{p}<0}$. As shown in Appendix C, and in
agreement with our results in Sec. III.3 and Sec. III.4.1, in this case the
state ${v=0}$ will always remain stable in the limit
${\Lambda\,\epsilon_{0}/D\rightarrow\infty}$. Hence, we will now modify the
above assumptions, keeping the enzyme-substrate interactions attractive,
${\chi_{s}<0}$, but now assuming that enzyme-product interactions are
repulsive, ${\chi_{p}>0}$; both are assumed to be weaker than enzyme-enzyme
interactions so that the sharp interface approximation remains valid.
Moreover, we consider condensates where the concentration of enzymes is
enriched inside the droplet, ${\Delta c>0}$. With these assumptions, one can
determine the concentration jumps for substrates and products at the droplet
interface, Eq. (33), in the limit
${\Lambda\,\epsilon_{0}/D\rightarrow\infty}$:

$\displaystyle\frac{s|_{\text{in}}}{s|_{\text{out}}}$
$\displaystyle=\exp\biggl{[}-\frac{\Lambda\,\chi_{s}\,\Delta
c}{D}\biggr{]}\to\infty\,,$ (52a)
$\displaystyle\frac{p|_{\text{in}}}{p|_{\text{out}}}$
$\displaystyle=\exp\biggl{[}-\frac{\Lambda\,\chi_{p}\,\Delta c}{D}\biggr{]}\to
0\,.$ (52b)

This limit implies ${p|_{\text{in}}=0}$ and ${s|_{\text{out}}=0}$, which
considerably simplifies the further analysis. Next, similar as in Sec.
III.4.1, we study the stability of the fixed point at ${v=0}$. A non-moving
droplet state, ${v=0}$, is unstable when the slope of the chemical potential
imbalance is positive,

$\partial_{v}\Delta\mu_{0}(v)\rvert_{v=0}=\partial_{v}\Delta\mu_{0}^{\star}(v)\rvert_{v=0}-\frac{2R\Delta
c}{Mc_{+}}>0\,,$ (53)

where the first equality follows from Eq. (45). After solving for the
substrate and product concentration profiles, as discussed in Sec. III.1 but
now with the simplified boundary conditions Eq. (52), and substituting these
profiles in Eq. (46), one has131313We carried out the analytic calculations
with the computer algebra system Mathematica [58] and provide the
corresponding notebooks in Ref. [59].

$\partial_{v}\Delta\mu_{0}^{\star}(v)\rvert_{v=0}=\frac{2R\,n_{\infty}\chi_{s}}{D(l_{0}^{2}-l_{+}^{2})}\color[rgb]{.5,.5,.5}\definecolor[named]{pgfstrokecolor}{rgb}{.5,.5,.5}\pgfsys@color@gray@stroke{.5}\pgfsys@color<EMAIL_ADDRESS>\times\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\left[l_{-}l_{+}\coth(R/l_{+})+l_{+}^{2}-l_{0}^{2}\right]\,,$
(54)

where ${l_{\pm}=\sqrt{D/(k_{1}c_{\pm}+k_{2})}}$ are the diffusion lengths
inside and outside the condensate, respectively. For inequality (53) to be
fulfilled, the maximal chemical potential imbalance must have a positive slope
[Eq. (54)]. Note that, since ${l_{\pm}<l_{0}=\sqrt{D/k_{2}}}$ and
$\chi_{s}<0$, the factor in the first line of Eq. (54) is always negative.
Therefore, a first necessary criterion for the emergence of droplet self-
propulsion follows by requiring that the second line of Eq. (54) is also
negative, so that the maximal chemical potential imbalance has a positive
slope. This suggests that droplets can exhibit self-propulsion if the droplet
size $R$ is sufficiently large. Because ${\coth(R/l_{+})>1}$ and
${\chi_{s}<0}$, Eq. (54) has the following bound:

$\partial_{v}\Delta\mu_{0}^{\star}(v)\rvert_{v=0}\leq\frac{2R\,n_{\infty}\chi_{s}}{D(l_{0}^{2}-l_{+}^{2})}\left[l_{-}l_{+}+l_{+}^{2}-l_{0}^{2}\right],$
(55)

leading to the additional criterion ${l_{0}^{2}>l_{-}l_{+}+l_{+}^{2}}$. This
condition can be rewritten by substituting the expressions for the diffusion
lengths ${l_{\pm}=\sqrt{D/(k_{1}c_{\pm}+k_{2})}}$, leading to

$\left(\frac{k_{1}c_{+}}{k_{2}}\right)^{2}>\frac{1+k_{1}c_{+}/k_{2}}{1+k_{1}c_{-}/k_{2}}\geq
1+\frac{k_{1}c_{+}}{k_{2}}\,,$ (56)

where the last inequality approximately becomes an equality in the limit of
strong phase separation, ${c_{-}/c_{+}\approx 0}$. In this limit, one can
solve the quadratic inequality (56) exactly, leading to the criterion
${k_{1}c_{+}/k_{2}>(\sqrt{5}+1)/2}$ for the reaction rates. Thus, the turnover
rate of substrates must be sufficiently fast for droplets to show self-
propulsion.

Finally, substituting Eq. (54) into the inequality (53) and solving for the
enzyme mobility leads to a lower enzyme mobility bound,

$\frac{M}{\Lambda}\geq-\frac{k_{B}T}{\chi_{s}n_{\infty}}\frac{\Delta
c}{c_{+}}\left[1-\frac{l_{-}l_{+}\coth(R/l_{+})}{l_{0}^{2}-l_{+}^{2}}\right]^{-1}\,,$
(57)

above which condensates will begin to self-propel; note that we have again
used the Einstein-Smoluchowsky relation ${D=\Lambda\,k_{B}T}$. Thus, the
critical ratio between the mobility of enzymes and the mobility of substrates
and products can be lowered by decreasing the enzyme’s propensity to phase
separate, $\Delta c$, or by increasing the interaction strength, $\chi_{s}$.
Importantly, the critical mobility ratio ${M/\Lambda\propto-
k_{B}T/(\chi_{s}n_{\infty})}$ can become very small in the limit of strong
enzyme-substrate interactions ${\chi_{s}\ll-k_{B}T/n_{\infty}}$. To summarize
this section, we have found that, in general, droplet self-propulsion requires
sufficiently fast reaction rates, sufficiently large condensates, and
sufficiently large enzyme mobility.

##  IV Discussion

We have studied the non-equilibrium dynamics of biomolecular condensates
(droplets) containing a high concentration of enzymes. Through their catalytic
role in biochemical reactions, these enzymes establish a spatial framework
that governs the arrangement of substrates and products. In turn, the
resulting substrate and product concentration profiles drive the enzyme flows
through molecular interactions. This interplay between non-equilibrium
reactions and reciprocal biomolecular interactions can give rise to novel
phenomena, such as droplet self-propulsion, which we have here studied in
depth.

We have identified the criteria for droplet self-propulsion in the two
opposing limits of weak or strong enzyme-substrate and enzyme-product
interactions. These limits correspond to small or large values of the
reciprocity parameter $\Lambda$, which represents the magnitude of the Onsager
coefficient (mobility) defining the extent of the currents induced by
gradients in the chemical potentials of substrates and products. In the limit
of weak interactions, our first analysis indicated that droplet self-
propulsion requires a very small ratio of mobilities ${\Lambda/M\ll 1}$ [Fig.
8], where $M$ is the mobility of enzymes. This could, for example, apply to
experiments where self-propelling droplets leave a long low-pH trail by
producing acid molecules [65]. However, there is also experimental evidence of
self-propelling urease-containing droplets [66, 67] which seem to disobey the
above restriction. Motivated by these experiments, we used our theory to
investigate a scenario where the reciprocity parameter $\Lambda$ is large or,
in other words, enzyme-substrate and enzyme-product interactions are strong.
Our analysis revealed that, in this case droplet self-propulsion is possible
for repulsive enzyme-product interactions, $\chi_{p}>0$, if the reaction rates
are sufficiently fast and the condensate is sufficiently large. These results
qualitatively agree with Ref. [67], where motion was only observed for
droplets whose radius exceeded $20\text{\,}\mathrm{\SIUnitSymbolMicro m}$.

Interestingly, we have found that liquid-like droplets generally move faster
than solid-like condensates. The reason for this finding is that droplets can
translocate simply by increasing solubility at one interface and decreasing
solubility at the opposing interface [68]. This will cause the condensate side
with higher solubility to move inward (dissolve) and the side with lower
solubility to move outward (grow) and can thereby move the condensate without
a net transport of mass. In contrast, the motion of a solid-like condensate
implies a transport of its entire mass, and will hence be limited by the
mobility of its molecular components.

_Adiabatic elimination scheme._ To gain these insights, we have developed a
self-consistent sharp-interface theory that describes the dynamics of active
phase-separated systems. Using an implicit adiabatic elimination scheme—based
on the assumption that the concentration profile of the enzymes is always
maintained adiabatically in a steady state—we inferred the chemical potential
and currents of the enzymes necessary to maintain a state with a given droplet
shape and velocity. Finally, we determined which value of the droplet speed is
thermodynamically consistent, by calculating the rate of power dissipation.
The power dissipation reveals if a given droplet state can be realized without
supplying additional energy or removing excess energy from the system. This
thermodynamic housekeeping analysis suggests that increasing the droplet speed
$|v|$ compared to a stable steady state would require energy influx, whereas
decreasing the droplet speed $|v|$ would require energy outflux.

_Droplet shape._ We restricted our analysis to a scenario where the
condensates remain in a round shape. This is generally a good approximation
near the onset of rotational symmetry breaking and self-propulsion, or if the
surface tension of the condensates is large; the latter case is realized for
strong enzyme-enzyme interactions. However, the approximation of a round
droplet shape will need to be relaxed if the enzyme-enzyme interactions are
sufficiently weak compared to the enzyme-substrate and enzyme-product
interactions. To then ensure that the droplet surface tension remains a
finite-valued material parameter, the interface stiffness of the droplet will
need to be inversely proportional to the strength of enzyme-enzyme
interactions. Then, droplets could deform into non-spherical shapes with non-
uniform interface curvature, thereby leading to Laplace pressure variations
that affect the enzyme currents. Such a generalization would require the
derivation of conditions for the droplet shape, which will be addressed in
future studies. This would aid future studies in investigating changes and
instabilities in droplet shape [24].

_Analogy to fluid mechanics._ At the heart of our theoretical analysis, we
studied enzyme droplets in a moving steady state by taking a sharp interface
limit. To ensure stationarity of the concentration profiles in the co-moving
frame, we constructed a chemical potential such that the enzyme currents were
forced to be divergence-free. This closely resembles the strategy used for
incompressible fluids, where an implicitly defined finite pressure field
enforces the fluid incompressibility condition. To make this analogy more
explicit, we will now compare the overdamped dynamics of the enzyme currents
[shown as vector field in Fig. 2(a)] to the inviscid hydrodynamics described
by Euler’s equations [69]. For an inviscid fluid with zero shear viscosity,
the fluid velocity $\mathbf{u}$ is determined by a balance between inertial
forces and local driving by applying an external force $\mathbf{g}$, with a
pressure field $p$ enforcing fluid incompressibility:

$\displaystyle\partial_{t}(\rho\,\mathbf{u})$
$\displaystyle=-\boldsymbol{\nabla}p+\mathbf{g}\,,$ (58a)
$\displaystyle\boldsymbol{\nabla}\cdot\mathbf{u}$ $\displaystyle=0\,,$ (58b)

where $\rho$ is the fluid density and fluid convection is neglected. In
analogy, on the right-hand side of Eq. (14b), the chemical potential $\mu_{0}$
plays the role of an effective pressure field that enforces incompressibility
as described by Eq. (14a), while the remaining terms are analogous to an
applied force. In contrast to Euler’s equations, however, the current on the
left-hand side of Eq. (14b), when divided by the enzyme mobility, corresponds
to drag friction forces instead of inertial forces, and is therefore
irreversible.

Finally, we note that the description of our model in terms of currents can
naturally incorporate hydrodynamic interactions. We expect that our self-
consistent sharp interface theory can be readily applied to such a scenario.
By choosing the mobility function $M(c)$ to model a solid where only the dense
phase is mobile, and the force field $\mathbf{f}$ to model interactions with a
neutral or charged solute, for example, one could compare the diffusiophoresis
[70] of solid colloids to that of liquid droplets. Even in the absence of
hydrodynamic interactions, our results already suggest that liquid-like
condensates can be transported much faster than solid particles, because such
a transport only requires moving the concentration profile instead of every
molecule contained in the condensate.

_Reciprocal interactions._ Using the theoretical framework developed here, we
have re-examined the phenomenon of self-propulsion of droplets, which was
studied in our previous work [28] in the non-reciprocal limit (${\Lambda=0}$),
for cases where reciprocity is restored. We found that for ${\Lambda>0}$
interactions with enzymes generally lead to discontinuities in the substrate
and product concentration profiles. As reciprocity moves the system closer to
thermodynamic equilibrium, additional conditions must be met to observe the
self-propulsion of the droplets:  (i) For weak enzyme-substrate and enzyme-
product interactions, the reciprocity parameter must lie below a certain
threshold, Eq. (51), to observe self-propulsion. In general, the onset of
self-propulsion will then be discontinuous (supercritical) as one varies, for
example, the enzyme mobility $M$. A second, lower threshold determines if the
system is bistable (binodal regime), or if the onset of self-propulsion is
spontaneous (spinodal regime). The onset of droplet self-propulsion only
becomes continuous (supercritical) in the nonreciprocal limit (${\Lambda=0}$).
(ii) In principle, droplet self-propulsion is also possible for large
reciprocity parameters, ${\Lambda\to\infty}$, if enzyme-product interactions
are repulsive, ${\chi_{p}>0}$. However, the size of the droplet must then be
considerable, and the non-equilibrium reaction rates must be sufficiently
fast.  Similar to the non-reciprocal limit in Ref. [28], as briefly discussed
in Appendix D, in the reciprocal case one can also observe droplet
positioning, elongation, and division.

_Fisher waves and other model systems._ The theory of moving droplet
interfaces developed here was in part inspired by the analysis of Fisher waves
[71], which model the expansion of growing populations. However, in contrast
to the model studied in the present work, which was derived based on
thermodynamic arguments and incorporated phase separation _a priori_ ,
population dynamics are arguably far from thermal equilibrium and described by
phenomenological models. It would therefore be interesting to test if the
ideas underlying the framework presented in this work can be adapted to front
propagation in such far-from-equilibrium systems. If that is the case, then
the derivation of balance conditions for the arising pseudo-chemical
potentials could pave the way to borrow powerful tools from non-equilibrium
thermodynamics for gaining new insights into ecological systems. However, we
also note that a central aspect of our analysis was the conservation law for
the total mass of enzymes, while population dynamics models lack such
conservation laws. Therefore, we expect that the present framework can be more
readily adopted for the analysis of non-reciprocal mass-conserved models [72,
73] or mass-conserved reaction-diffusion systems with advection [74].

_Conclusion and outlook._ The interplay of reciprocal biomolecular
interactions that cause phase separation, coupled with out-of-equilibrium
chemical reactions, is a widespread organizational motif in living cells. By
developing a versatile sharp interface theory we have gained insight into the
specific conditions under which this motif can lead to the self-propulsion of
condensates. The mechanism underlying condensate motion can be understood as a
gradient in solubility which leads to translocation of droplets [67]. More
specifically, attractive interactions with substrates in the surrounding
solution and repulsion with products in the condensate, locally increase the
solubility of the enzymes. This leads to the dissolution of the trailing edge
of the droplet, which is enriched in products and depleted in substrates, and
growth of the leading edge of the droplet.

The ability of biomolecular condensates to migrate along concentration
gradients, which was recently recognized theoretically [68, 28, 75] and
demonstrated experimentally for pH [65, 66, 67] or salt [76] gradients, is
reminiscent of the diffusiophoresis of colloids [70]. For colloidal
diffusiophoresis, however, hydrodynamic shear stresses play an essential role
because interactions between the colloidal surface and the inhomogeneously
distributed solute lead to an effective surface tension gradient [70]. While
we have here neglected fluid viscosity, its consideration in future studies
would make it possible to further elucidate the parallels and differences
between droplet and colloid motion. Moreover, accounting for fluid mechanics
would provide a more accurate description of experimental systems in which
condensates move due to chemical gradients such as salt or pH variations.

The response to an applied gradient in concentration, or mechanical (that is,
viscoelastic) properties of the surrounding medium is a unifying theme among
soft and living matter. In addition to chemical gradients which drive
diffusiophoresis and cellular chemotaxis, for example, it was shown that
droplets [77, 78] and cells [79] can also migrate along stiffness gradients.
Given the success of phase field models in describing cell migration [80, 81,
82, 83, 84, 85], we hypothesize that some of the ideas developed in the
present manuscript could also apply to cellular dynamics.

###### Acknowledgements.

We thank Dominik Schumacher and Lotte Søgaard-Andersen for helpful
discussions. We acknowledge financial support by the German Research
Foundation (DFG) through TRR 174 (Project ID No. 269423233) and SFB1032
(Project ID No. 201269156) and the Excellence Cluster ORIGINS under Germany’s
Excellence Strategy (EXC-2094-390783311). During his time at the Ludwig-
Maximilians-Universität München, AG was supported by a DFG fellowship through
the Graduate School of Quantitative Biosciences Munich (QBM). During his time
at the Massachusetts Institute of Technology, AG was supported by the National
Science Foundation (NSF) through grant number 2044895 and by an EMBO
Postdoctoral Fellowship (ALTF 259-2022). During his time at Sorbonne
University, LD has received funding from the European Union’s Horizon 2020
research and innovation programme under the Marie Skłodowska-Curie grant
agreement No. 860949. IM has received funding from the European Union’s
Framework Programme for Research and Innovation Horizon 2020 under the Marie
Skłodowska-Curie Grant Agreement No. 754388 (LMU Research Fellows) and from
LMU excellent, funded by the Federal Ministry of Education and Research (BMBF)
and the Free State of Bavaria under the Excellence Strategy of the German
Federal Government and the Länder. EF was supported by the European Union
(ERC, CellGeom, project number 101097810) and the Chan-Zuckerberg Initiative
(CZI).

## References

  * Brangwynne  _et al._ [2009] C. P. Brangwynne, C. R. Eckmann, D. S. Courson, A. Rybarska, C. Hoege, J. Gharakhani, F. Jülicher, and A. A. Hyman, Germline P granules are liquid droplets that localize by controlled dissolution/condensation, Science 324, 1729 (2009).
  * Hyman  _et al._ [2014] A. A. Hyman, C. A. Weber, and F. Jülicher, Liquid-liquid phase separation in biology, Annu. Rev. Cell Dev. Biol. 30, 39 (2014).
  * 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 (2017).
  * Lyon  _et al._ [2021] A. S. Lyon, W. B. Peeples, and M. K. Rosen, A framework for understanding the functions of biomolecular condensates across scales, Nat. Rev. Mol. Cell Biol. 22, 215 (2021).
  * Alberti  _et al._ [2019] S. Alberti, A. Gladfelter, and T. Mittag, Considerations and challenges in studying liquid-liquid phase separation and biomolecular condensates, Cell 176, 419 (2019).
  * Shin and Brangwynne [2017] Y. Shin and C. P. Brangwynne, Liquid phase condensation in cell physiology and disease, Science 357, eaaf4382 (2017).
  * Choi  _et al._ [2020] J.-M. Choi, A. S. Holehouse, and R. V. Pappu, Physical principles underlying the complex biology of intracellular phase transitions, Annu. Rev. Biophys. 49, 107 (2020).
  * Klosin  _et al._ [2020] A. Klosin, F. Oltsch, T. Harmon, A. Honigmann, F. Jülicher, A. A. Hyman, and C. Zechner, Phase separation provides a mechanism to reduce noise in cells, Science 367, 464 (2020).
  * Shrinivas  _et al._ [2019] K. Shrinivas, B. R. Sabari, E. L. Coffey, I. A. Klein, A. Boija, A. V. Zamudio, J. Schuijers, N. M. Hannett, P. A. Sharp, R. A. Young, and A. K. Chakraborty, Enhancer features that drive formation of transcriptional condensates, Mol. Cell 75, 549 (2019).
  * Henninger  _et al._ [2021] J. E. Henninger, O. Oksuz, K. Shrinivas, I. Sagi, G. LeRoy, M. M. Zheng, J. O. Andrews, A. V. Zamudio, C. Lazaris, N. M. Hannett, T. I. Lee, P. A. Sharp, I. I. Cissé, A. K. Chakraborty, and R. A. Young, RNA-mediated feedback control of transcriptional condensates, Cell 184, 207 (2021).
  * Schumacher  _et al._ [2017] D. Schumacher, S. Bergeler, A. Harms, J. Vonck, S. Huneke-Vogt, E. Frey, and L. Søgaard-Andersen, The PomXYZ proteins self-organize on the bacterial nucleoid to stimulate cell division, Dev. Cell 41, 299 (2017).
  * Ostwald [1900] W. Ostwald, About the supposed isomerism of the red and yellow mercury oxide and the surface tension of solid bodies, Z. Phys. Chem. 34U, 495 (1900).
  * Lifshitz and Slyozov [1961] I. M. Lifshitz and V. V. Slyozov, The kinetics of precipitation from supersaturated solid solutions, J. Phys. Chem. Solids 19, 35 (1961).
  * Wagner [1961] C. Wagner, Theorie der Alterung von Niederschlägen durch Umlösen (Ostwald-Reifung), Z. Elektrochem. 65, 581 (1961).
  * Bray [1994] A. J. Bray, Theory of phase-ordering kinetics, Adv. Phys. 43, 357 (1994).
  * Zwicker [2022] D. Zwicker, The intertwined physics of active chemical reactions and phase separation, Curr. Opin. Colloid Interface Sci. 61, 101606 (2022).
  * Weber  _et al._ [2019] C. A. Weber, D. Zwicker, F. Jülicher, and C. F. Lee, Physics of active emulsions, Rep. Prog. Phys. 82, 064601 (2019).
  * Glotzer  _et al._ [1994] S. C. Glotzer, D. Stauffer, and N. Jan, Monte Carlo simulations of phase separation in chemically reactive binary mixtures, Phys. Rev. Lett. 72, 4109 (1994).
  * Glotzer  _et al._ [1995] S. C. Glotzer, E. A. Di Marzio, and M. Muthukumar, Reaction-controlled morphology of phase-separating mixtures, Phys. Rev. Lett. 74, 2034 (1995).
  * Christensen  _et al._ [1996] J. J. Christensen, K. Elder, and H. C. Fogedby, Phase segregation dynamics of a chemically reactive binary mixture, Phys. Rev. E 54, R2212 (1996).
  * Carati and Lefever [1997] D. Carati and R. Lefever, Chemical freezing of phase separation in immiscible binary mixtures, Phys. Rev. E 56, 3127 (1997).
  * Zwicker  _et al._ [2015] D. Zwicker, A. A. Hyman, and F. Jülicher, Suppression of Ostwald ripening in active emulsions, Phys. Rev. E 92, 012317 (2015).
  * Lamorgese and Mauri [2016] A. Lamorgese and R. Mauri, Spinodal decomposition of chemically reactive binary mixtures, Phys. Rev. E 94, 022605 (2016).
  * Zwicker  _et al._ [2017] D. Zwicker, R. Seyboldt, C. A. Weber, A. A. Hyman, and F. Jülicher, Growth and division of active droplets provides a model for protocells, Nat. Phys. 13, 408 (2017).
  * Wurtz and Lee [2018] J. D. Wurtz and C. F. Lee, Chemical-reaction-controlled phase separated drops: Formation, size selection, and coarsening, Phys. Rev. Lett. 120, 078102 (2018).
  * Li and Cates [2020] Y. I. Li and M. E. Cates, Non-equilibrium phase separation with reactions: a canonical model and its behaviour, J. Stat. Mech. 2020, 053206 (2020).
  * Kirschbaum and Zwicker [2021] J. Kirschbaum and D. Zwicker, Controlling biomolecular condensates via chemical reactions, J. R. Soc. Interface 18, 20210255 (2021).
  * 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).
  * [29] See Supplemental Material at [URL will be inserted by publisher] for further details and an additional analysis of the model, which includes Refs. .
  * Bergeler and Frey [2018] S. Bergeler and E. Frey, Regulation of Pom cluster dynamics in Myxococcus xanthus, PLOS Comp. Biol. 14, e1006358 (2018).
  * Kober  _et al._ [2019] M. Kober, S. Bergeler, and E. Frey, Can a flux-based mechanism explain protein cluster positioning in a three-dimensional cell geometry?, Biophys. J. 117, 420 (2019).
  * Hanauer  _et al._ [2021] C. Hanauer, S. Bergeler, E. Frey, and C. P. Broedersz, Theory of active intracellular transport by DNA relaying, Phys. Rev. Lett. 127, 138101 (2021).
  * Hnisz  _et al._ [2017] D. Hnisz, K. Shrinivas, R. A. Young, A. K. Chakraborty, and P. A. Sharp, A phase separation model for transcriptional control, Cell 169, 13–23 (2017).
  * Schede  _et al._ [2023] H. H. Schede, P. Natarajan, A. K. Chakraborty, and K. Shrinivas, A model for organization and regulation of nuclear condensates by gene activity, Nat. Commun. 14, 4152 (2023).
  * Natarajan  _et al._ [2023] P. Natarajan, K. Shrinivas, and A. K. Chakraborty, A model for cis-regulation of transcriptional condensates and gene expression by proximal lncrnas, Biophys. J. 122, 2757–2772 (2023).
  * Lafontaine  _et al._ [2021] D. L. J. Lafontaine, J. A. Riback, R. Bascetin, and C. P. Brangwynne, The nucleolus as a multiphase liquid condensate, Nat. Rev. Mol. Cell Biol. 22, 165–182 (2021).
  * 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).
  * Dwivedi  _et al._ [2022] P. Dwivedi, D. Pillai, and R. Mangal, Self-propelled swimming droplets, Curr. Opin. Colloid Interface Sci. 61, 101614 (2022).
  * Michelin [2023] S. Michelin, Self-propulsion of chemically active droplets, Annu. Rev. Fluid Mech. 55, 77–101 (2023).
  * Joanny and Ramaswamy [2012] J.-F. Joanny and S. Ramaswamy, A drop of active matter, J. Fluid Mech. 705, 46 (2012).
  * Khoromskaia and Alexander [2015] D. Khoromskaia and G. P. Alexander, Motility of active fluid drops on surfaces, Phys. Rev. E 92, 062311 (2015).
  * Whitfield and Hawkins [2016] C. A. Whitfield and R. J. Hawkins, Instabilities, motion and deformation of active fluid droplets, New J. Phys. 18, 123016 (2016).
  * Kree  _et al._ [2017] R. Kree, P. S. Burada, and A. Zippelius, From active stresses and forces to self-propulsion of droplets, J. Fluid Mech. 821, 595 (2017).
  * Yoshinaga [2019] N. Yoshinaga, Self-propulsion of an active polar drop, J. Chem. Phys. 150, 184904 (2019).
  * Trinschek  _et al._ [2020] S. Trinschek, F. Stegemerten, K. John, and U. Thiele, Thin-film modeling of resting and moving active droplets, Phys. Rev. E 101, 062802 (2020).
  * 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).
  * 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 (2005).
  * De Groot and Mazur [2013] S. R. De Groot and P. Mazur,  _Non-Equilibrium Thermodynamics_, Dover Books on Physics (Dover Publications, 2013).
  * Balian [2007] R. Balian,  _From Microphysics to Macrophysics_ (Springer, Berlin, Heidelberg, 2007).
  * Cahn and Hilliard [1958] J. W. Cahn and J. E. Hilliard, Free energy of a nonuniform system. I. Interfacial free energy, The Journal of Chemical Physics 28, 258–267 (1958).
  * Cahn [1961] J. W. Cahn, On spinodal decomposition, Acta Metall. 9, 795 (1961).
  * Frey and Kroy [2005] E. Frey and K. Kroy, Brownian motion: a paradigm of soft matter and biological physics, Ann. Phys. 517, 20 (2005).
  * Logg  _et al._ [2012] A. Logg, K.-A. Mardal, G. N. Wells, _et al._ ,  _Automated Solution of Differential Equations by the Finite Element Method_, edited by A. Logg, K.-A. Mardal, and G. N. Wells (Springer, Berlin, Heidelberg, 2012).
  * Alnæs  _et al._ [2014] M. S. Alnæs, A. Logg, K. B. Ølgaard, M. E. Rognes, and G. N. Wells, Unified form language: A domain-specific language for weak formulations of partial differential equations, ACM Trans. Math. Softw. 40, 10.1145/2566630 (2014).
  * Scroggs  _et al._ [2022a] M. W. Scroggs, J. S. Dokken, C. N. Richardson, and G. N. Wells, Construction of arbitrary order finite element degree-of-freedom maps on polygonal and polyhedral cell meshes, ACM Trans. Math. Softw. 48, 10.1145/3524456 (2022a).
  * Scroggs  _et al._ [2022b] M. W. Scroggs, I. A. Baratta, C. N. Richardson, and G. N. Wells, Basix: a runtime finite element basis evaluation library, J. Open Source Softw. 7, 3982 (2022b).
  * Strogatz [2018] S. H. Strogatz,  _Nonlinear Dynamics and Chaos_ , second edition ed. (CRC Press, Boca Raton, 2018).
  * [58] Wolfram Research, Inc., Mathematica, Version 12.2, Champaign, IL, 2020.
  * Goychuk  _et al._ [2023] A. Goychuk, L. Demarchi, I. Maryshev, and E. Frey, Source code - The role of non-reciprocity in the fluid-like dynamics of enzyme-enriched droplets (2023).
  * Kolomeisky  _et al._ [1998] A. B. Kolomeisky, G. M. Schütz, E. B. Kolomeisky, and J. P. Straley, Phase diagram of one-dimensional driven lattice gases with open boundaries, Journal of Physics A: Mathematical and General 31, 6911 (1998).
  * Parmeggiani  _et al._ [2003] A. Parmeggiani, T. Franosch, and E. Frey, Phase coexistence in driven one-dimensional transport, Phys. Rev. Lett. 90, 086601 (2003).
  * Halatek  _et al._ [2018] J. Halatek, F. Brauns, and E. Frey, Self-organization principles of intracellular pattern formation, Philos. Trans. R. Soc. B: Biological Sciences 373, 20170107 (2018).
  * Halatek and Frey [2018] J. Halatek and E. Frey, Rethinking pattern formation in reaction–diffusion systems, Nat. Phys. 14, 507–514 (2018).
  * 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).
  * Hanczyc [2011] M. M. Hanczyc, Metabolism and motility in prebiotic structures, Phil. Trans. R. Soc. B 366, 2885–2893 (2011).