# Nonequilibrium interfacial properties of chemically driven fluids

Yongick Cho  William M. Jacobs<EMAIL_ADDRESS>Department of Chemistry,
Princeton University, Princeton, NJ 08544, USA

###### Abstract

Chemically driven fluids can demix to form condensed droplets that exhibit
phase behaviors not observed at equilibrium. In particular, nonequilibrium
interfacial properties can emerge when the chemical reactions are driven
differentially between the interior and exterior of the phase-separated
droplets. Here, we use a minimal model to study changes in the interfacial
tension between coexisting phases away from equilibrium. Simulations of both
droplet nucleation and interface roughness indicate that the nonequilibrium
interfacial tension can either be increased or decreased relative to its
equilibrium value, depending on whether the driven chemical reactions are
accelerated or decelerated within the droplets. Finally, we show that these
observations can be understood using a predictive theory based on an effective
thermodynamic equilibrium.

##  I Introduction

A “chemically driven” fluid can maintain a nonequilibrium steady-state (NESS)
when provided with a constant external supply of a chemical fuel [1, 2]. By
coupling this fuel source to chemical reactions taking place in the fluid, and
thereby biasing the forward/backward reaction kinetics, it is possible to
drive the fluid away from equilibrium. A NESS is established when the constant
supply of chemical fuel prevents the fluid from relaxing to thermal
equilibrium. For example, in living cells, biomolecules can undergo reactions
involving post-translational modifications, such as changes in conformational
state due to phosphorylation [3, 4, 5]. If the free energy used to drive these
transitions away from equilibrium is derived from a constant supply of ATP,
then the reaction kinetics will break detailed balance, and the intracellular
fluid can be considered to be at a NESS.

Like equilibrium fluids, chemically driven fluids can demix to form coexisting
phases via liquid–liquid phase separation [6, 1]. However, the phase behavior
and phase-separation dynamics of a chemically driven fluid may differ
qualitatively from that of a fluid at equilibrium. For example, coarsening can
be suppressed [7], leading to a monodisperse size distribution of coexisting
droplets at steady state [8, 9], or accelerated [10] relative to equilibrium.
Nonequilibrium droplets can exhibit self-regulatory behaviors such as self-
division [11], self-propulsion [12], and self-organization into complex
internal microstructures that are controlled by the chemical driving forces
[13, 14]. The kinetics of droplet assembly and disassembly can also differ
from equilibrium phase-separation kinetics and can be controlled by tuning the
chemical driving forces [15, 16, 17]. This ability to control the dynamics of
droplet formation suggests that nonequilibrium phase separation may be a
facile means of achieving spatiotemporal regulation within living cells [18,
19].

Understanding how the interfacial tension between coexisting phases changes at
a NESS is essential for describing the phase behavior and dynamics of
chemically driven fluids [15, 20]. In particular, nonequilibrium interfacial
tensions, which imply deviations of the steady-state interfacial properties
from equilibrium, emerge when the driven chemical-reaction kinetics differ
between the interior and exterior of a phase-separated droplet [15]. Such
spatially inhomogeneous reaction kinetics can arise in complex fluids when the
chemical fuel source is uniformly distributed but the catalyst required for
the driven reaction pathway tends to partition into one of the two coexisting
phases. For example, spatially inhomogeneous reaction kinetics can be induced
by the preferential enrichment of enzymes either inside or outside of
biomolecular condensates in living cells [13, 2]. Since enzymatic partitioning
is observed in a wide variety of intracellular condensates as a means of
maintaining chemically specific environments to carry out specialized
biochemical functions [21, 22], it is likely that inhomogeneous reaction
kinetics—and thus nonequilibrium interfacial properties—are a common feature
of nonequilibrium phase-separated droplets in living cells.

In this article, we use a minimal model of a chemically driven fluid to
quantify the effects of driven chemical reactions on the interfacial
properties between coexisting phases at a NESS. We show that simulations of
droplet nucleation and interface roughness lead to compatible inferences of
the nonequilibrium interfacial tension in systems with spatially inhomogeneous
chemical reactions. The rest of this paper is structured as follows: In Sec.
II, we describe the minimal model and consider two chemical-reaction schemes
in which the driven reaction is either accelerated or decelerated inside of
the droplets. In Sec. III, we discuss how a theory based on an “effective
equilibrium” can be used to predict both bulk and interfacial properties of
the coexisting phases. Then in Sec. IV, we show that the nonequilibrium
interfacial tensions inferred from nucleation and roughness simulations agree
with the predictions of our effective-equilibrium theory. Finally, in Sec. V,
we discuss the implications of our results for conducting noninvasive
measurements of interfacial properties in nonequilibrium fluids.

##  II A minimal model for chemically driven fluids

Figure 1: Simulating driven chemical reactions at a phase-separated NESS. (a)
Schematic of an open system showing vapor (left) and liquid (right) phases in
direct coexistence. The equilibrium free-energy difference between the bonding
(B) and inert (I) states in the particle reservoir is $\Delta f_{\text{res}}$.
The effective free-energy differences between the B and I states in the liquid
and vapor phases are $\Delta f_{\text{l}}$ and $\Delta f_{\text{v}}$,
respectively. (b) Kinetic scheme of particle exchange and internal chemical
reactions. In simulations of an open system, each lattice site can
stochastically transition between being empty (E) or being occupied by either
a bonding-state (B) or inert-state (I) particle according to the specified
transition rates. (c) Examples of inhomogeneous chemical-reaction models with
different dependencies on the local potential energy, $u$, experienced by
B-state particles. Orange and blue colors indicate chemical-reaction models I
and II, defined by Eqs. (3) and (4), respectively.

In this section, we present a minimal model of a chemically driven fluid, in
which reactions allow particles to interconvert between two internal
conformational states. In Sec. II.1, we describe the modeling approach first
introduced in Ref. [15]. We then consider two alternative inhomogeneous
chemical-reaction schemes, which result in different reaction kinetics in the
condensed and dilute phases, in Sec. II.2. Implementation details of our
simulations are presented in Sec. II.3.

###  II.1 Modeling driven chemical reactions in an open system

In order to investigate interfacial properties in chemically driven systems,
we employ a lattice model of a fluid consisting of particles in an implicit
solvent [15]. Extending the classical two-dimensional lattice-gas model, we
allow particles to undergo chemical reactions between two internal
conformational states: a bonding (B) state and an inert (I) state. Each
lattice site can be occupied by at most one particle in either internal state.
Short-ranged attractive interactions between B-state particles tend to drive
phase separation, whereas both I-state particles and empty (E) lattice sites,
which represent the implicit solvent, are non-interacting. Specifically, a
particle in the B state engages in nearest-neighbor attractive interactions
with other neighboring B-state particles with pairwise interaction strength
$\epsilon\;(<0)$. These interactions give rise to a local potential energy
$u=\epsilon b\;(\leq 0)$ at each lattice site, where $b$ is the number of
nearest-neighbor lattice sites occupied by B-state particles. I-state
particles, on the other hand, do not interact with nearest-neighbor particles,
and are thus isoenergetic to empty lattice sites, with $u=0$.

Chemical reactions between the B and I conformational states of individual
particles occur via two distinct reaction pathways. An undriven pathway is
governed by equilibrium thermodynamics, while the driven pathway is considered
to be coupled to a constant chemical fuel source. We model these pathways
within the framework of stochastic thermodynamics [23, 24], which treats
reactions as Markovian events taking place at a constant absolute temperature
$T$. Along the undriven pathway, the ratio of the forward and backward
reaction rates is dictated by detailed balance and the internal free-energy
difference between the B and I states in an ideal-gas reservoir, $\Delta
f_{\text{res}}$ (Fig. 1a). Throughout this work, we consider scenarios in
which $\Delta f_{\text{res}}>0$, such that the I state is more stable in the
dilute vapor (v) phase, while attractive interactions are required to
stabilize the B state in the condensed liquid (l) phase. The fugacities
$z_{\text{B}}$ and $z_{\text{I}}$ of B and I-state particles in an equilibrium
ideal gas are thus related to this equilibrium free-energy difference by
$\beta\Delta f_{\text{res}}\equiv-\ln(z_{\text{B}}/z_{\text{I}})$, where
$\beta\equiv(k_{\text{B}}T)^{-1}$. Along the driven pathway, conformational-
state changes are driven out of equilibrium by a chemical-potential difference
$\Delta\mu$ that originates from the constant chemical fuel source. The ratio
of the forward and backward reaction rates on the driven pathway is thus
increased relative to the equilibrium ratio by a factor
$\exp(\beta\Delta\mu)$. The relative flux between the driven and undriven
pathways then determines the steady-state distribution between the
conformational states. In the limit that the flux through the driven pathway
vanishes, we recover the equilibrium distribution. By contrast, if the flux
through the driven pathway dominates, then the system approaches an effective
equilibrium distribution governed by an internal free-energy difference of
$\Delta f_{\text{res}}+\Delta\mu$.

We model the fluid as an open system connected to a particle reservoir, which
allows us to simulate competition between the driven and undriven pathways
under conditions where the reactions are rate-limiting (Fig. 1a). This limit
describes a scenario in which particle transport within the fluid is much
faster than the typical rate of transitions between conformational states
along either pathway. Importantly, this limit is appropriate for studying
interfacial properties that arise purely from the competition between the
driven and undriven pathways in a phase-separated system, which is the focus
of this paper. Open systems also provide advantages for studying phase
coexistence that are analogous to those of the equilibrium grand-canonical
ensemble, such as the elimination of interfaces and reduced finite size
effects [25].

Within the framework of an open system, we implement reactions along the
driven pathway as direct transitions between the B and I states, while
reactions along the undriven pathway occur indirectly via particle exchanges
with the equilibrium particle reservoir (Fig. 1b). Particle insertion and
removal rates depend on the reservoir fugacities, $z_{\text{B}}$ and
$z_{\text{I}}$, the local potential energy, $u$, experienced by a B-state
particle, and the base exchange rate, $D$, between the system and the
reservoir. Direct transitions between the B and I states, corresponding to
driven $\text{I}\rightarrow\text{B}$ and $\text{B}\rightarrow\text{I}$
reactions, proceed with reaction rates $k_{\text{IB}}$ and $k_{\text{BI}}$,
respectively. We use dimensionless reaction rates
$k_{\text{I}\rightarrow\text{B}}\equiv D^{-1}k_{\text{IB}}$ and
$k_{\text{B}\rightarrow\text{I}}\equiv D^{-1}k_{\text{BI}}$ for notational
simplicity throughout the rest of the paper. When the chemical drive
$\Delta\mu$ along the driven pathway is nonzero, the probability of observing
a sequence of events around the single-cycle network in Fig. 1b differs from
that of the time-reversed sequence, breaking detailed balance and leading to a
nonzero net probability current. The chemical drive $\Delta\mu$ is directly
related to the ratio of the probability of observing a forward sequence of
transitions around the cycle in the B-to-I direction relative to that of its
time-reversed sequence [23, 24],

$\beta\Delta\mu\equiv\ln\left[z_{\text{B}}k_{\text{B}\rightarrow\text{I}}/z_{\text{I}}k_{\text{I}\rightarrow\text{B}}e^{\beta
u}\right].$ (1)

Rearranging Eq. (1) leads to the so-called “local detailed balance” condition
for the direct-transition rates [23, 24],

$k_{\text{B}\rightarrow\text{I}}/k_{\text{I}\rightarrow\text{B}}=\exp(\beta\Delta
f_{\text{res}}+\beta u+\beta\Delta\mu).$ (2)

Thus, Eq. (2) implies that the ratio of the forward and reversed reaction
rates on the driven pathway is increased relative to the ratio on the undriven
pathway by a factor $\exp(\beta\Delta\mu)$, and detailed balance is recovered
only when $\Delta\mu=0$.

###  II.2 Models of inhomogeneous driven chemical reactions

In this paper, we focus on inhomogeneous driven chemical reactions, whose
kinetics differ between the liquid and vapor phases. When the reaction
kinetics are identical, or homogeneous, in both phases, then the steady-state
density distribution of the driven system can be mapped exactly to that of an
equilibrium system [9]. Consequently, no change in the interfacial properties
is observed when the reaction kinetics are homogeneous [15]. By contrast, when
the kinetics differ between the coexisting phases, there is no single
effective-equilibrium model that can describe both phases simultaneously.
Coexisting phases in inhomogeneously driven fluids can only be approximately
described using different effective-equilibrium models, giving rise to
nonequilibrium interfacial properties [15].

Inhomogeneous reaction kinetics are implemented in our model by controlling
the relative fluxes between the driven and undriven pathways in the two
phases. To this end, we tune the ratio of the fluxes by making the base rate
of the driven, direct transitions between the B and I states to be dependent
on a particle’s local environment within the open system. Specifically, we set
the backward rate $k_{\text{I}\rightarrow\text{B}}$ to be dependent on the
local potential energy, $u$. The forward rate,
$k_{\text{B}\rightarrow\text{I}}$, then follows from the local detailed
balance condition, Eq. (2). Here, we consider two chemical reaction schemes in
which $k_{\text{I}\rightarrow\text{B}}$ takes a Metropolis-like form (Fig.
1c):

$\displaystyle k_{\text{I}\rightarrow\text{B}}=$ $\displaystyle
k^{\circ}\min[1,\exp(-\beta\Delta f_{\text{res}}-\beta\Delta\mu-\beta u)]$ (3)
$\displaystyle k_{\text{I}\rightarrow\text{B}}=$ $\displaystyle
k^{\circ}\min[1,\exp(-\beta\Delta f_{\text{res}}-\beta\Delta\mu+\beta u)].$
(4)

The prefactor $k^{\circ}$ sets the relative flux of the driven pathway
compared to the undriven pathway in a dilute vapor phase, where $u\approx 0$.
We refer to the inhomogeneous reaction kinetics defined by Eqs. (3) and (4) as
chemical-reaction models I and II, respectively, throughout the rest of the
paper. In chemical-reaction model I, $k_{\text{I}\rightarrow\text{B}}$
decreases monotonically with respect to $u$, implying that the driven reaction
is faster—all else being equal—in the liquid phase, where the potential energy
tends to be low, than in the vapor phase, where the potential energy tends to
be high. The behavior is the opposite in model II.

###  II.3 Implementation via kinetic Monte Carlo algorithm

To study the consequences of inhomogeneous driven chemical reactions, we
perform kinetic Monte Carlo simulations [26] of our model on a two-dimensional
square lattice with periodic boundary conditions (Fig. 1a). We treat particle
exchanges and direct transitions between B and I states as first-order
Markovian reactions. We then simulate the stochastic evolution of the lattice
with the parameters $z_{\text{B}}$, $z_{\text{I}}$, and $\beta\Delta\mu$ held
constant in time.

All simulation data presented in this paper are obtained at a dimensionless
interaction strength of $\beta\epsilon=-2.95$, which is stronger than the
critical dimensionless interaction strength,
$\beta\epsilon_{\text{c}}={-2\ln(1+\sqrt{2})}$, [27] of the associated
equilibrium lattice-gas model. Two of the three control parameters,
$z_{\text{B}}$, $z_{\text{I}}$, and $\beta\Delta\mu$, must then be specified
in order to define a unique coexistence point. We choose parameters for
simulation by controlling both the extent of the nonequilibrium drive,
$\beta\Delta\mu$, and the average number density of particles in either
internal state in the vapor phase, $\rho_{\text{v}}$. The latter of these
conditions is achieved by tuning the internal free-energy difference $\Delta
f_{\text{res}}$; focusing on $\rho_{\text{v}}$ as opposed to $\Delta
f_{\text{res}}$ makes for an easier connection to experimental settings.
Throughout this work, we consider coexistence points at which the vapor-phase
number density is $\rho_{\text{v}}=0.05$. Finally, the prefactor for the
driven chemical-reaction pathway, $k^{\circ}$, is set to $k^{\circ}=0.1$ for
chemical-reaction model I and $k^{\circ}=1$ for model II. These choices yield
relative timescales for the driven and undriven pathways in which the
interfacial properties clearly deviate from equilibrium, as opposed to the
behavior in the undriven-reaction-dominated ($k^{\circ}\rightarrow 0$) or
driven-reaction-dominated ($k^{\circ}\rightarrow\infty$) limits [15].

##  III Chemical reaction-induced nonequilibrium phase behavior at coexistence

Figure 2: Definitions of nonequilibrium phase coexistence. (a) To define phase
coexistence based on the steady-state density distribution,
$p(\rho_{\text{B}})$, we equate the nonequilibrium potential, $\Phi$, in the
liquid and vapor phases (see text). The liquid and vapor phases are defined by
the location of the top of the barrier, $\rho_{\text{B}}^{*}$, with respect to
$-\ln p(\rho_{\text{B}})$. In this example, the steady-state density
distribution is shown for lattice dimension $L=10$. Insets: Representative
lattice configurations are shown, from left to right, in the vapor phase; at
the top of the barrier, $\rho=\rho_{\text{B}}^{*}$; and in the liquid phase.
Colors indicate particle internal states as in Fig. 1a. (b) To define phase
coexistence based on direct-coexistence simulations, we measure the
probability, $P_{\text{l}}$, that the liquid phase expands to fill the entire
lattice starting from a planar liquid–vapor interface. Coexistence occurs when
$P_{\text{l}}$ equals 50%. In this example, calculations are shown for lattice
dimension $L=16$. Inset: Schematic of the $L\times d$ setup for direct-
coexistence simulations. Simulations are performed within the $L\times L$
shaded window, which tracks the center of the planar liquid–vapor interface.
(c) The two definitions of phase coexistence from panels (a) and (b) match
closely in the macroscopic limit, $L\rightarrow\infty$. Filled circles and
open squares indicate coexistence points determined from the definitions in
panels (a) and (b), respectively. Open diamonds indicate the extrapolated
coexistence points in the macroscopic limit. All data are obtained using
chemical-reaction model I.

In this section, we demonstrate the effect of driven chemical reactions on the
bulk phase behavior of the coexisting phases, which sets the stage for
examining nonequilibrium interfacial properties. In Sec. III.1, we establish
the condition for phase coexistence under nonequilibrium conditions based on
mechanical balance, providing further validation of the approach presented in
Ref. [15]. We then examine how the nonequilibrium phase behavior depends on
the chemical-reaction models defined in Sec. III.2, and we show how these
results can be understood in terms of a theory based on an effective-
equilibrium approximation, which was originally introduced in Ref. [15], in
Sec. III.3.

###  III.1 Determination of phase coexistence

Two phases are in coexistence when there is no net energy or mass flux between
them and they are mechanically balanced [28]. At equilibrium, a variational
principle based on the second law of thermodynamics implies that equal
temperatures, chemical potentials, and pressures between the coexisting phases
lead to energy, mass, and mechanical balance, respectively, and vice versa. In
the grand-canonical ensemble at equilibrium, in which the temperature and
chemical potential are fixed, ensuring that the grand-potential densities are
the same in the two phases satisfies the equal-pressure condition [29].
Although our model is driven out-of-equilibrium, the particle reservoir still
guarantees that the energy and mass conditions are satisfied by providing the
energy and mass input required to keep each phase at steady-state. However,
equating the grand-potential densities in the two phases does not ensure
mechanical balance in fluids that are out of equilibrium due to the general
absence of an equation of state [30, 31].

We therefore examine two closely related definitions of nonequilibrium phase
coexistence based on mechanical balance, which are not guaranteed to coincide
in the case of chemically driven fluids a priori. First, we determine phase
coexistence by equating the steady-state probability of observing the open
system in either the liquid or the vapor phase. We then compare this
definition of coexistence with the result of nonequilibrium direct-coexistence
simulations.

In the steady-state distribution approach, we calculate the steady-state
number-density distribution of the B-state particle, $p(\rho_{\text{B}})$, on
the domain $\rho_{\text{B}}\in[0,1]$. This calculation is performed using an
$L\times L$ lattice. To this end, we employ a form of nonequilibrium umbrella
sampling (NEUS) [32] with respect to a one-dimensional collective variable,
$\rho_{\text{B}}$, that we have developed previously [15]. In this method, we
divide the one-dimensional $\rho_{\text{B}}$ axis into non-overlapping
windows, calculate the steady-state distribution within each window, and
enforce detailed balance between the adjacent windows in order to reconstruct
the entire steady-state distribution. We observe that the system exists with
high probability at either low or high values of $\rho_{\text{B}}$,
corresponding to the vapor and liquid phases, respectively (Fig. 2a). The
phases are separated by a high barrier with respect to $-\ln
p(\rho_{\text{B}})$, as expected for a first-order phase transition [28]. The
height of this barrier scales approximately linearly with the lattice
dimension, $L$, since the barrier reflects the formation of a liquid–vapor
interface on the two-dimensional lattice. Based on the location of the top of
the barrier, $\rho_{\text{B}}^{*}$, we determine the probability of being in
the vapor phase,
$p_{\text{v}}\equiv\int_{0}^{\rho_{\text{B}}^{*}}p(\rho_{\text{B}})d\rho_{\text{B}}$,
versus the liquid phase,
$p_{\text{l}}\equiv\int_{\rho_{\text{B}}^{*}}^{1}p(\rho_{\text{B}})d\rho_{\text{B}}$.
Due to the distinct nature of the phases and the extremely small steady-state
probability near the top of the barrier, the precise definitions of
$p_{\text{v}}$ and $p_{\text{l}}$ are insensitive to the choice of collective
variable and the location of $\rho_{\text{B}}^{*}$. We then define the
nonequilibrium potential difference, $\Delta\Phi$, per lattice site between
the liquid and vapor phases, $\beta\Delta\Phi(L)\equiv
L^{-2}\ln(p_{\text{l}}/p_{\text{v}})$. The potential difference in the
macroscopic limit,
$\Delta\Phi_{\infty}\equiv\lim_{L\rightarrow\infty}\Delta\Phi(L)$, is exactly
equal to the difference between the grand-potential densities of the two
phases in a fluid at equilibrium [33]. We therefore associate
$\Delta\Phi_{\infty}=0$ with nonequilibrium phase coexistence in this
approach. This operational definition of $\Delta\Phi$ is particularly useful,
as it plays the role of a thermodynamic potential for the liquid–vapor phase
transition under nonequilibrium conditions and can be evaluated directly from
NEUS simulations.

In the direct-coexistence approach [34], we consider an open system that is
periodic in the vertical direction but always in contact with the liquid and
vapor phases in the horizontal direction (inset of Fig. 2b). Fluctuations at
the liquid–vapor interface cause the portion of the lattice that is occupied
by the liquid phase to expand or shrink. Starting from an initial
configuration in which the lattice is split equally between the liquid and
vapor phases, we simulate the system until the liquid phase either fills the
lattice entirely or vanishes completely. For simulation efficiency, we only
allow particle exchanges and chemical reactions to take place within a window
of size $L\times L$, which tracks the position of the planar liquid–vapor
interface (inset of Fig. 2b); the width of this window is much larger than the
interface height fluctuations, which ensures that the window boundaries do not
affect the growth or shrinkage of the liquid phase. We then measure the
probability that the lattice is occupied by the liquid phase at the end of the
simulation, $P_{\text{l}}$. Consistent with the notion of mechanical balance
[31], we associate phase coexistence with the condition $P_{\text{l}}=50\%$.

We find that the coexistence points identified by these two approaches match
closely when extrapolated to the macroscopic limit. In both approaches, we
find that the chemical drive at phase coexistence, $\Delta\mu_{\text{coex}}$,
shows a consistent dependence on the lattice dimension, $L$, which allows us
to extrapolate the coexistence condition to the macroscopic limit,
$L\rightarrow\infty$ (Fig. 2c). Interestingly, the $1/L$ scaling that we
observe is reminiscent of effective Coulombic interactions predicted in other
treatments of chemically reactive systems [35, 17]. Even in the case of a
strongly driven system $(|\beta\Delta\mu|\gg 1)$, the two definitions of phase
coexistence yield very similar coexistence points, $\Delta\mu_{\text{coex}}$
(differing by less than $10^{-2}k_{\text{B}}T$), in the macroscopic limit
(Fig. 2c). These differences are typically orders of magnitude smaller than
the chemical drive at coexistence, as well as the range of $\Delta\mu$ values
over which simulate droplet nucleation (see Sec. IV.1). We therefore conclude
that these two coexistence definitions agree well within the range of system
parameters that we consider here, and we utilize the definition
$\beta\Delta\Phi_{\infty}=0$ for operational convenience.

###  III.2 Effective-equilibrium quantification of nonequilibrium phase
behavior

Figure 3: Nonequilibrium phase behavior of coexisting bulk phases. (a) The
phase diagram of the open system depends on the inhomogeneous chemical-
reaction scheme, as illustrated by the differences between chemical-reaction
models I (orange) and II (blue). In both cases, coexistence is calculated for
a vapor phase with a number density of $\rho_{\text{v}}=0.05$. (b) Driven
chemical reactions induce different effective internal free-energy
differences, $\Delta\Delta f\equiv\Delta f_{\text{l}}-\Delta f_{\text{v}}$,
between the coexisting phases (see text). (c) In the inhomogeneous chemical-
reaction schemes considered in this work, the effective internal free-energy
differences arise primarily due to changes in the vapor phase (solid symbols)
as opposed to the liquid phase (open symbols). In all panels, symbols indicate
results obtained from simulations of bulk phases at the coexistence point,
while lines indicate theoretical predictions obtained using FLEX. Shaded
(unshaded) regions in panels (a) and (b) indicate where the liquid (vapor)
phase is predicted to be stable by FLEX. Solid and dotted lines in panel (c)
indicate the vapor and liquid phases, respectively.

Because the relative fluxes depend on the kinetics of the chemical reactions,
reaction models I and II lead to different coexistence lines, as shown in Fig.
3a. In this phase diagram, $\Delta f_{\text{res}}$ decreases with
$\Delta\mu_{\text{coex}}$ along the coexistence line in the case of model I,
while the value of $\Delta f_{\text{res}}$ at coexistence is essentially
independent of the chemical drive in the case of model II. In both cases, the
total number density in the vapor phase is fixed such that
$\rho_{\text{v}}=0.05$ at coexistence.

The differences between the phase diagrams of the two chemical-reaction models
suggest that the bulk properties of the coexisting phases depend on the
details of the reaction models at a given nonequilibrium chemical drive. To
understand how the driven chemical reactions affect the bulk properties of the
coexisting liquid and vapor phases, we map each phase to an effective-
equilibrium model. Here, effective equilibrium implies an equilibrium fluid
with the same B-state interaction strength, $\epsilon$, as in the
nonequilibrium system, but with a potentially different effective internal
free-energy difference between the B and I states, $\Delta f$. This effective
free-energy difference is calculated using a generalized Widom insertion
formula [36, 15],

$\beta\Delta
f\equiv-\ln(\rho_{\text{B}}/\rho_{\text{I}})+\ln\langle\exp(-\beta\Delta
u_{\text{I}\rightarrow\text{B}})\rangle_{\text{I}},$ (5)

where $\rho_{\text{B}}$ and $\rho_{\text{I}}$ are the steady-state number
densities of B and I-state particles, respectively, in a bulk phase; $\Delta
u_{\text{I}\rightarrow\text{B}}$ is the potential energy change due to
changing an I-state particle to the B state; and the angle brackets
$\langle\cdot\rangle_{\text{I}}$ indicate a steady-state average conditioned
on a tagged lattice site being occupied by an I-state particle. In
equilibrium, Eq. 5 reduces to $\beta\Delta f_{\text{eq}}=\beta\Delta
f_{\text{res}}$, where $\Delta f_{\text{eq}}$ is the equilibrium internal
free-energy difference. At a nonzero $\Delta\mu_{\text{coex}}$, we evaluate
Eq. (5) from a steady-state trajectory in which the lattice is occupied
entirely by one phase. In this way, we are able to determine effective free-
energy differences, $\Delta f_{\text{l}}$ and $\Delta f_{\text{v}}$, for the
coexisting liquid and vapor phases, respectively (Fig. 1a). We note that while
this effective equilibrium mapping is not exact in general, it nevertheless
offers useful insights and provides a quantitative foundation for
understanding and predicting nonequilibrium phase behaviors (see Appendix A).

In general, coexisting phases cannot be described by the same effective-
equilibrium model in driven fluids with inhomogeneous chemical reactions. It
is therefore useful to quantify the difference between the effective internal
free-energy differences in the coexisting liquid and vapor phases,
$\Delta\Delta f\equiv\Delta f_{\text{l}}-\Delta f_{\text{v}}$. In Fig. 3b, we
observe that as the system is driven farther away from equilibrium along a
coexistence line, $\Delta\Delta f$ either increases or decreases
monotonically, signaling growing differences between the effective
thermodynamics of the coexisting phases as the magnitude of
$\beta\Delta\mu_{\text{coex}}$ is increased. However, the relationship between
$\beta\Delta\Delta f$ and $\beta\Delta\mu_{\text{coex}}$ depends on the
details of the inhomogeneous chemical-reaction model. We find that
$\beta\Delta\Delta f$ increases with $\beta\Delta\mu_{\text{coex}}$ in the
case of model I, while the relationship is reversed in the case of model II.
The physical origin of this behavior is discussed in the following section.

Unexpectedly, we find that the deviations of $\beta\Delta\Delta f$ from zero
are almost entirely due to the vapor phase for both chemical-reaction models
that we examine here. Fig. 3c shows that $\beta\Delta f_{\text{v}}$ either
decreases or increases monotonically in the case of models I or II,
respectively. By contrast, $\Delta f_{\text{l}}$ exhibits an almost negligible
change at the coexistence conditions. A theoretical explanation for these
observations is discussed below.

###  III.3 FLEX description of nonequilibrium phase behavior

To understand the effect of driven chemical reactions on the nonequilibrium
phase behavior, we utilize the “Fixed Local Environment approXimation" (FLEX)
theoretical framework introduced in Ref. [15] (see Appendix B). The central
assumption of FLEX is that fluctuations of the local potential energy, $u$,
experienced by a B-state particle are small at the level of an individual
lattice site. This assumption is well justified in a bulk phase (see Appendix
A) and turns out to provide qualitatively accurate predictions of interfacial
properties as well. Assuming that we can neglect such fluctuations in the
local “environment” around a lattice site allows us to map the steady-state at
a tagged lattice site to an effective equilibrium at constant $u$. In this
way, FLEX predicts that the effective internal free-energy difference, $\Delta
f(u)$, of a lattice site in an environment where a B-state particle would
experience a local potential energy $u$ is given by (see Appendix B)

$\beta\Delta
f(u)=\beta\Delta{f}_{\text{res}}+\ln\left[\dfrac{1+k_{\text{I}\rightarrow\text{B}}(u)(1+e^{\beta\Delta{f}_{\text{res}}})e^{\beta\Delta\mu}}{1+k_{\text{I}\rightarrow\text{B}}(u)(1+e^{\beta\Delta{f}_{\text{res}}})}\right].$
(6)

This equation implies that the effective internal free-energy difference
between the B and I states at a tagged lattice site is the sum of the
equilibrium internal free-energy difference in the particle reservoir and the
influence of the driven chemical-reaction pathway.

In order to employ FLEX without resorting to simulation, we assume that
$u\approx 4\epsilon$ in the liquid phase, meaning that every lattice is
surrounded by B-state particles, and $u\approx 0$ in the vapor phase, meaning
that particles are sparsely distributed (see Appendix A). Eq. (6) then
predicts that $\Delta\Delta f=\Delta f_{\text{l}}-\Delta f_{\text{v}}=\Delta
f(4\epsilon)-\Delta f(0)$ is an increasing function of the chemical drive when
$k_{\text{I}\rightarrow\text{B}}(4\epsilon)>k_{\text{I}\rightarrow\text{B}}(0)$,
as in model I, while the relationship is the opposite when
$k_{\text{I}\rightarrow\text{B}}(4\epsilon)<k_{\text{I}\rightarrow\text{B}}(0)$,
as in model II. These qualitative predictions agree with the observed
simulation results presented in Fig. 3b. Thus, we attribute the qualitative
differences in $\beta\Delta\Delta f$ at coexistence to the acceleration of the
driven chemical reactions in either the liquid or vapor phase in the case of
chemical-reaction models I or II, respectively (Fig. 1c). Quantitative
predictions of the bulk-phase properties obtained from FLEX also agree well
with the simulation results at phase coexistence (solid lines in Fig. 3; see
Appendix C).

##  IV Effects of driven chemical reactions on nonequilibrium interfacial
properties

In this section, we demonstrate that driven chemical reactions can induce
changes in the interfacial properties of a nonequilibrium phase-separated
fluid relative to those of an equilibrium fluid. We first infer the
nonequilibrium interfacial tension from simulations of droplet-nucleation
kinetics in Sec. IV.1, extending the results first presented in Ref. [15]. To
further validate the conclusions of Ref. [15], we now perform these
calculations along two independent supersaturation pathways to establish that
the interfacial tension is a material property determined by the coexistence
point. We then show how changes in the interfacial tension that arise in both
chemical-reaction models can be understood within the FLEX framework in Sec.
IV.2. Finally, in Sec. IV.3, we demonstrate that the nonequilibrium
interfacial tension inferred from nucleation simulations is consistent with
the height fluctuations of the liquid–vapor interface at phase coexistence.
The results of this section provide a new, independent test of the
interfacial-tension measurements reported in Ref. [15].

###  IV.1 Inferring a nonequilibrium interfacial tension using classical
nucleation theory

Figure 4: Different supersaturation protocols yield compatible inferences of
the nonequilibrium interfacial tension in simulations of droplet nucleation.
(a) When supersaturating the vapor phase beyond the coexistence point, the
nonequilibrium thermodynamic driving force, $\Delta\Phi_{\infty}$, can be
controlled either by tuning $\Delta\mu$ or $\rho_{\text{v}}$. Shown here are
two such protocols, which vary either $\rho_{\text{v}}$ (squares) or
$\beta\Delta\mu$ (circles). The coexistence point at
$(\beta\Delta\mu_{\text{coex}},\rho_{\text{v}})=(1.87,0.05)$ is indicated by
an empty circle, and dotted lines show FLEX predictions of
equi-$\Delta\Phi_{\infty}$ contours. In this example, the chemical-reaction
kinetics follow model I. (b) The interfacial tension inferred from the
nucleation kinetics is essentially independent of the protocol used to
supersaturate the vapor phase. The effective free-energy barrier, $\beta\Delta
F^{*}$ (see inset and refer to the text for details), is determined at
multiple values of the supersaturation, $\beta\Delta\Phi_{\infty}$. Symbols
correspond to the protocols shown in panel (a). The dotted line is a fit of
the observed nucleation kinetics to classical nucleation theory (CNT), from
which the nonequilibrium interfacial tension, $\sigma_{\text{CNT}}$, is
inferred.

We first infer the nonequilibrium interfacial tension using simulations of
droplet nucleation (Fig. 4a), which we interpret using classical nucleation
theory (CNT) [15]. When applied to equilibrium fluids, CNT predicts that the
nucleation of a liquid droplet from a supersaturated, metastable vapor phase
follows the minimum free-energy pathway along a reaction coordinate
corresponding to the number of particles, $n$, in a cluster. This cluster
represents a “nucleus” of the stable liquid phase (inset of Fig. 4b). The
height of the nucleation barrier on this free-energy landscape, $\Delta
F^{*}$, is determined by a competition between the thermodynamic driving force
favoring the phase transition to the liquid phase and the free-energy cost of
forming the liquid–vapor interface for the liquid droplet. At a NESS, this
thermodynamic driving force must be replaced by an appropriate measure of the
supersaturation with respect to the nonequilibrium coexistence point, as we
discuss below. If $\Delta F^{*}$ is much larger than typical thermal
fluctuations, such that $\beta\Delta F^{*}\gg 1$, then the nucleation rate is
dominated by the rate of crossing this barrier. Under such conditions, CNT
predicts that the nucleation rate density, $J$, follows the phenomenological
law [37]

$J=\rho_{1}D^{*}\Gamma\exp(-\beta\Delta F^{*}),$ (7)

where the height of the free-energy barrier enters in the form of a Boltzmann
factor. Three terms comprise the kinetic prefactor: the number density of
B-state monomers that are not engaged in any attractive interactions,
$\rho_{1}$; the one-dimensional diffusion coefficient along the reaction
coordinate near the top of the nucleation barrier, $D^{*}$; and the Zeldovich
factor, $\Gamma$, which accounts for barrier re-crossing events due to
fluctuations along the reaction coordinate before a nucleation attempt
ultimately succeeds or fails. CNT has been shown to provide a quantitative
description of nucleation kinetics in two-dimensional lattice-gas models at
equilibrium [38].

We compute the nucleation rate density, $J$, starting from a supersaturated
vapor phase using forward-flux sampling (FFS) [39]. While calculating $J$, we
also determine the commitment probability for a successful nucleation event,
$\phi(n)$, along the reaction coordinate. We then compute the Zeldovich factor
by fitting $\phi(n)$ to a harmonic free-energy landscape in the vicinity of
the top of the nucleation barrier where, $n=n^{*}$ [40, 15], at which the
commitment probability equals 50% [41]. The remaining quantities $\rho_{1}$
and $D^{*}$ are independently determined from simulations of the bulk vapor
phase and the evolution of the nucleus size along a trajectory initiated from
the critical nucleus size, $n^{*}$, respectively. From these measurements, we
compute $\beta\Delta F^{*}\equiv-\ln(J/\rho_{1}D^{*}\Gamma)$, which
corresponds to the effective free-energy barrier height for a nonequilibrium
nucleation process.

We then apply CNT to our model of a fluid with driven chemical reactions by
introducing a nonequilibrium interfacial tension, $\sigma_{\text{CNT}}$. To
this end, we assume that the nucleation behavior also admits an effective-
equilibrium description [16]. We therefore fit the observed effective barrier
heights, $\beta\Delta F^{*}$, as a function of the supersaturation,
$\exp(\beta\Delta\Phi_{\infty})$, to an analytical expression for the two-
dimensional lattice-gas free-energy landscape [38] using $\sigma_{\text{CNT}}$
as the sole fitting parameter [15]. This approach assumes that the
thermodynamic driving force for the nucleation process is given by
$\Delta\Phi_{\infty}$, which is consistent with the fundamental CNT assumption
that both the bulk and the interfacial properties of a nucleus are
characteristics of the fluid in the macroscopic limit.

We now consider how the CNT description of the nucleation kinetics depends on
the manner in which a chemically driven fluid is supersaturated. In an
equilibrium lattice-gas model at constant interaction strength, the
thermodynamic driving force can be increased only by increasing the total
particle density in the vapor phase, $\rho_{\text{v}}$. In fluids with driven
chemical reactions, however, the chemical drive, $\Delta\mu$, acts as an
additional degree of freedom and thus provides an alternative means of tuning
the supersaturation. We therefore consider two different supersaturation
protocols (Fig. 4a), in which we either increase $z_{\text{B}}$ and
$z_{\text{I}}$ under fixed $\Delta\mu$ and $\Delta f_{\text{res}}$, or
increase $\Delta\mu$ with $\Delta\mu+\Delta f_{\text{res}}$ fixed [see Eqs.
(3) and (4)]. The former “$\rho_{\text{v}}$ protocol” increases the particle
density in the vapor phase while keeping $\Delta\mu$ and the relative flux
between the reaction pathways constant. By contrast, the latter “$\Delta\mu$
protocol” changes the nonequilibrium potential by tuning the chemical drive
along the driven reaction pathway. In both supersaturation protocols, we keep
the reaction-rate, $k_{\text{I}\rightarrow\text{B}}$, fixed. We also use FLEX
to predict equi-$\Delta\Phi_{\infty}$ contours in the
$(\beta\Delta\mu,\rho_{\text{v}})$ plane (Fig. 4a). Viewing the parameter
space in this way shows that an entire family of supersaturation protocols
could plausibly be defined by monotonically tuning $\beta\Delta\mu$ and
$\rho_{\text{v}}$ away from a coexistence point, where
$\beta\Delta\Phi_{\infty}=0$.

Remarkably, despite the differences between the $\rho_{\text{v}}$ and
$\Delta\mu$ supersaturation protocols, we find that the nucleation kinetics
can be described by CNT to near quantitative accuracy using a common value of
the nonequilibrium interfacial tension, $\sigma_{\text{CNT}}$. This result can
be seen in Fig. 4b, where we plot the effective nucleation barrier as a
function of the supersaturation along the two supersaturation protocols shown
in Fig. 4a. We emphasize that the fitted value of $\sigma_{\text{CNT}}$
differs substantially from the equilibrium interfacial tension in the example
shown in Fig. 4b, consistent with the fact that the fluid is driven away from
equilibrium by inhomogeneous chemical reactions. The good agreement between
the barrier measurements along both protocols, as well as the fit to CNT,
indicates that the effective nucleation barriers along both protocols are
determined by the same value of the nonequilibrium interfacial tension. We
therefore propose that the nonequilibrium interfacial tension is relatively
independent of the supersaturation protocol, as long as the degree of
supersaturation is low ($\beta\Delta\Phi_{\infty}\ll 1$) and the protocols
converge at the same coexistence point. This result agrees with a fundamental
premise of CNT that $\sigma_{\text{CNT}}$ is identical to the interfacial
tension of a macroscopic interface between two coexisting phases [37]. Thus,
for the fitting parameter $\sigma_{\text{CNT}}$ to be considered to be a
nonequilibrium interfacial tension, it is important that its value is
determined by the coexistence point alone. Our simulation results along these
two supersaturation protocols suggest that this is indeed the case.

###  IV.2 Dependence of the nonequilibrium interfacial tension on the chemical
drive at phase coexistence

We find that the nonequilibrium interfacial tension inferred from droplet-
nucleation simulations is a function of the chemical drive at coexistence,
$\Delta\mu_{\text{coex}}$ (Fig. 5a). As the fluid is driven further away from
equilibrium, we find that the interfacial tension determined from nucleation
kinetics, $\sigma_{\text{CNT}}$, tends to deviate further from the equilibrium
value, $\sigma_{\text{eq}}$. However, the manner in which
$\sigma_{\text{CNT}}$ changes as a function of $\Delta\mu_{\text{coex}}$
differs between the two chemical-reaction models. In the case of model I,
$\sigma_{\text{CNT}}$ increases monotonically with respect to
$\Delta\mu_{\text{coex}}$, while it monotonically decreases in the case of
model II. Thus, knowing whether the chemical drive biases the driven reaction
pathway in either the $\text{I}\rightarrow\text{B}$ or the
$\text{B}\rightarrow\text{I}$ direction is insufficient to predict whether the
nonequilibrium interfacial tension will increase or decrease relative to
equilibrium.

Figure 5: Inhomogeneous driven chemical-reaction kinetics determine the
nonequilibrium interfacial tension. (a) The deviation of the nonequilibrium
interfacial tension from the equilibrium value, $\Delta\sigma_{\text{CNT}}$,
as a function of the coexistence point at $\rho_{\text{v}}=0.05$. Results are
shown for both inhomogeneous chemical-reaction models I and II. Solid curves
indicate theoretical predictions from FLEX. (b) Increasing the relative
stability of B-state particles in the vapor phase ($\beta\Delta\Delta f>0$)
reduces the nonequilibrium interfacial tension relative to equilibrium, while
decreasing the relative stability $(\beta\Delta\Delta f<0)$ has the opposite
effect. This relationship holds regardless of the inhomogeneous chemical
reaction-kinetics scheme.

Instead, the effective-equilibrium parameter $\Delta\Delta f$ offers a common
explanation of the trends in the nonequilibrium interfacial tension for both
chemical-reaction models (Fig. 5b). This observation implies that the effects
of the driven chemical reactions on the interfacial properties are related to
the differences between the effective-equilibrium behaviors of the coexisting
phases. This relationship can be explained by considering the relative
populations of B and I-state particles at the interface between the coexisting
liquid and vapor phases. We first recall that a low value of $\Delta f$
indicates that the B state is more stable than the I state in terms of its
internal free-energy; thus, decreasing $\Delta f$ increases the probability of
finding a particle in the B state at a given $u$. If $\Delta\Delta f=\Delta
f_{\text{l}}-\Delta f_{\text{v}}$ is positive, then the effective internal
free-energy difference is lower in the vapor phase than in the liquid phase,
meaning that the population of B-state particles is higher in the vapor than
what would be expected on the basis of an effective-equilibrium description of
the liquid phase. A negative $\Delta\Delta f$ implies the opposite behavior.
Because $k_{\text{I}\rightarrow\text{B}}$ monotonically changes with $u$ for
both chemical-reaction models, the sign of $\Delta\Delta f$ also indicates
whether the B-state particle is more or less populated at the liquid–vapor
interface than would be expected on the basis of the liquid phase. Increasing
the population of B-state particles at the interface relative to the liquid
phase, at a positive $\Delta\Delta f$, has the effect of lowering the
effective free-energy cost of forming the interface. As a result, a positive
$\Delta\Delta f$ implies a reduced interfacial tension relative to the
equilibrium value. By the same logic, a negative $\Delta\Delta f$ implies an
increase in the interfacial tension relative to the equilibrium value. These
qualitative arguments agree with the observed changes in the interfacial
tension presented in Fig. 5b.

The effective-equilibrium approach therefore allows us to explain how the
interfacial tension changes along the coexistence line for a given chemical-
reaction model. As noted above, $\Delta\Delta f$ is predominantly determined
by the behavior of the vapor phase in both chemical-reaction models studied
here (Fig. 3c). Thus, in the case of chemical-reaction model I, $\Delta
f_{\text{v}}$ decreases monotonically with respect to
$\Delta\mu_{\text{coex}}$, and the B-state particle population increases, in
the vapor phase. The B-state population at the interface is affected in the
same way, so that the interfacial tension decreases with the applied chemical
drive. In the case of model II, on the other hand, $\Delta f_{\text{v}}$
increases with the applied chemical drive, so that the interfacial tension
shows the opposite dependence. These contrasting trends agree with our
simulation findings in Fig. 5a.

These arguments are also borne out quantitatively within the FLEX framework
(see Appendix D). FLEX assumes that B-state particles at a flat interface
experience a local potential energy of $u=\epsilon$, which lies between
$u\approx 4\epsilon$ and $u\approx 0$ for the liquid and vapor phases,
respectively. Because $k_{\text{I}\rightarrow\text{B}}$ monotonically changes
with $u$ for both chemical-reaction models, Eq. (6) predicts that $\Delta
f(\epsilon)$ lies between $\Delta f_{\text{l}}$ and $\Delta f_{\text{v}}$, so
that a positive (negative) $\Delta\Delta f$ also indicates an increase
(decrease) in the B-state particle population at the interface relative to the
liquid phase. FLEX then predicts that the interfacial tension can be
determined by defining an effective dimensionless interaction energy at the
interface, $\beta\tilde{\epsilon}$, which can differ from the actual
dimensionless interaction energy, $\beta\epsilon$. When $\Delta\Delta f$ is
positive (negative), the relative increase (decrease) in the B-state particles
at the interface means that the effective interaction energy needed to adsorb
B-state particles to the interface is decreased (increased) relative to
$\beta\epsilon$. The nonequilibrium interfacial tension can then be predicted
on the basis of the effective interaction energy, $\beta\tilde{\epsilon}$, at
the interface (see Appendix D). These FLEX predictions are shown as solid
lines in Fig. 5a.

###  IV.3 Inferring a nonequilibrium interfacial tension using direct-
coexistence simulations

We further examine the relationship between driven chemical reactions and the
nonequilibrium interfacial tension by focusing on the roughness of the
liquid–vapor interface. Since our interpretation of the nucleation kinetics
points to nonequilibrium effects at the droplet interface, we hypothesize that
changes in the interfacial properties should be observable in simulations of
liquid and vapor phases in direct coexistence. In order to verify this, we
test whether changes in the interfacial tension match changes in the height
fluctuations of the liquid phase at the interface, which are directly related
to the interfacial tension at equilibrium. If the effective interaction
strength increases, then intuition based on equilibrium phase behavior
suggests that the interface should become smoother due to the greater
effective free-energy cost of introducing curvature at the interface.
Similarly, we anticipate that rougher interfaces will be observed when the
interfacial tension is lower.

In our direct-coexistence simulations, we focus on the location of the planar
liquid–vapor interface in the normal direction. During the simulation, points
on the interface, $\\{h\\}$, defined as the left-most boundary of all rows
that are filled with B-state particles (inset of Fig. 6a), fluctuate as
B-state particles attach and detach from the liquid phase. We then define the
interfacial roughness, $\Delta h$, to be the standard deviation of the heights
$\\{h\\}$,

$\Delta h\equiv\langle(h-\langle h\rangle)^{2}\rangle^{1/2},$ (8)

where the angle brackets denote an average over all $L$ rows of the lattice.
This definition of $\Delta h$ ignores overhanging particles on sides of the
particle rows. We simulate the interface between the coexisting phases on an
$L\times L$ lattice that is equally divided between liquid and vapor phases,
in the same manner as described in Sec. III.1. Due to the system-size
dependence of the coexistence point (Fig. 2c) and the maximum length scale of
fluctuations imposed by the finite lattice dimension $L$ [42], the interfacial
roughness is inherently subject to finite-size effects. We therefore perform
simulations at the $L$-dependent coexistence point of the finite-sized lattice
with $L=64$, as opposed to the coexistence point in the macroscopic limit.

Figure 6: Measurements of interfacial roughness at phase coexistence are
consistent with the nonequilibrium interfacial tension inferred from
nucleation kinetics. (a) Change in the interfacial roughness, $\Delta h$, at
coexistence, compared to the equilibrium value, $\Delta h_{\text{eq}}$.
Symbols show the average of 96 independent simulations of a planar interface
with lattice dimension $L=64$, and error bars indicate standard errors of
mean. Solid lines are theoretical predictions from FLEX. Inset: Schematic of
the interfacial roughness, $\Delta h$, at coexistence. (b) Comparison between
the nonequilibrium interfacial tension inferred from the interfacial
roughness, $\sigma_{\text{rough}}$, and from nucleation kinetics,
$\sigma_{\text{CNT}}$. Orange and blue colors indicate inhomogeneous chemical-
reaction models I and II, respectively.

We interpret the roughness of the liquid–vapor interface using the equilibrium
solid-on-solid (SOS) approximation [42]. In this approximation, B-state
particles are assumed to attach to the interface only at the top of the
particle rows, without overhanging neighboring rows. Assuming that the
interface can be described using the equilibrium SOS approximation, the
observed roughness, $\Delta h$, can be related to the equilibrium value,
$\Delta h_{\text{eq}}$, via

$\dfrac{\Delta h}{\Delta
h_{\text{eq}}}=\dfrac{\sinh(\beta\epsilon/4)}{\sinh(\beta\tilde{\epsilon}/4)},$
(9)

where $\beta\tilde{\epsilon}$ is the effective dimensionless interaction
strength at the interface under the SOS approximation. We can then obtain the
effective interfacial tension, $\sigma_{\text{rough}}(\beta\tilde{\epsilon})$,
by solving for $\beta\tilde{\epsilon}$ using Eq. (9) and applying an
analytical formula for the equilibrium lattice-gas model [see Eq. (17) in
Appendix D].

Analyzing the interfacial roughness within the framework of the SOS model
broadly supports our hypothesis of a chemically driven nonequilibrium
interfacial tension. In Fig. 6a, we show a comparison between the change in
the roughness relative to equilibrium, $(\Delta h-\Delta h_{\text{eq}})/\Delta
h_{\text{eq}}$, and predictions based on FLEX. For these theoretical
predictions, we use the FLEX expression for the effective dimensionless
interaction strength, $\beta\tilde{\epsilon}$ [see Eq. (16) in Appendix D], in
combination with the SOS relation, Eq. (9). The direct-coexistence
measurements and FLEX predictions are consistent in that the interfacial
roughness increases with $\Delta\mu_{\text{coex}}$ in the case of chemical-
reaction model I, while it decreases in the case of model II. These results
are also consistent with the equilibrium intuition that the interfacial
roughness increases when the interfacial tension decreases, and vice versa.

To directly compare the results of inference methods based on nucleation
kinetics and interfacial roughness, we plot the relationship between
$\sigma_{\text{CNT}}$ and $\sigma_{\text{rough}}$ in Fig. 6b. Each point in
this comparison corresponds to a unique nonequilibrium coexistence point,
$\beta\Delta\mu_{\text{coex}}$. Although the dynamic range of
$\sigma_{\text{rough}}$ tends to be smaller than that of
$\sigma_{\text{CNT}}$, which is likely a consequence of the SOS approximation
used in the former measurement, we observe a strong correlation between these
two independent inferences of the nonequilibrium interfacial tension.
Importantly, this correlation holds equally well for both chemical-reaction
models. We therefore conclude that both inference methods yield compatible
measurements of the nonequilibrium interfacial tension, and that this
interfacial property is entirely determined by the nonequilibrium coexistence
conditions.

##  V Discussion

In this paper, we show that driven chemical reactions alter the material
properties of interfaces between coexisting phases in nonequilibrium fluids.
To this end, we consider two models of inhomogeneous chemical-reaction
kinetics, which are either accelerated or decelerated in the condensed liquid
phase relative to the dilute vapor phase, using a minimal lattice model [15].
We find that the interfacial tension, determined either from measurements of
droplet-nucleation kinetics or from measurements of fluctuations at a planar
interface, changes as the system is driven out of equilibrium. Our
interpretation of this effect as a nonequilibrium interfacial tension is
supported by the strong correlation between these two independent
measurements. However, the relationship between the nonequilibrium interfacial
tension and the applied chemical drive depends on the details of the chemical-
reaction model.

We can explain the observed deviations in the interfacial properties using a
theoretical framework based on effective-equilibrium models of the coexisting
phases. We interpret our simulation results by mapping each phase in a
nonequilibrium fluid to an effective-equilibrium model whose particle density
distribution is closest to that of the NESS under a constant chemical drive.
By calculating the effective free-energy differences between the internal
states of particles in both phases, we find that we can explain the bulk-phase
properties of driven fluids with different chemical-reaction models. Then, by
applying this effective-equilibrium approach to predict the effective free-
energy cost of forming a liquid–vapor interface, we obtain a consistent
explanation for the observed trends in the nonequilibrium interfacial tensions
with different chemical-reaction models.

The concept of effective equilibrium thus plays an essential role in our
analysis of the effects of driven chemical reactions. Nonetheless, we
emphasize that our simulation approach, which is based on stochastic
thermodynamics, and the underlying lattice model do not make any assumptions
regarding an effective equilibrium. Instead, we utilize the concept of an
effective equilibrium only when interpreting our simulation results or when
generating approximate analytical predictions using FLEX.

Alternatively, the concept of an effective equilibrium can be employed as a
starting point for modeling the thermodynamic and transport properties of
chemically driven fluids. For example, in mesoscale modeling approaches based
on linear irreversible thermodynamics [1, 2], the NESS is obtained starting
from a free-energy functional that implies local equilibration within a fluid
element [9, 16, 17]. In such models, therefore, local effective equilibrium is
assumed to provide a complete description of the system, from which stationary
states and phase-transformation pathways can be computed, rather than an
approximation for analyzing nonequilibrium phenomena as in this paper.
However, this premise may not be justified in nonequilibrium fluids with
driven chemical reactions, since the predictions of the effective equilibrium
may be incorrect at molecular length scales (see Appendix A). Models that
treat fluid elements at a mesoscopic scale may therefore not capture effects
arising as a result of molecular-scale deviations from the effective-
equilibrium assumption. We expect the differences between these approaches
will become more substantial as the system approaches a critical point, since
fluctuations in the local environment may become too large for fluid elements
to be described adequately by a single effective equilibrium.

Although we have studied the effects of driven chemical reactions in a
reaction-limited regime in this paper, we expect that our qualitative
predictions of nonequilibrium interfacial properties will also hold beyond
this regime. Extending our modeling approach to examine the influence of
particle diffusion on nonequilibrium interfacial properties thus represents a
promising direction for future study. Studying nonequilibrium interfacial
properties may also be possible within the context of linear irreversible
thermodynamics [20], although local effective-equilibrium assumptions are
integral to such models as noted above. However, we note that any study of
interfacial properties in mean-field models of chemically driven fluids [9, 2]
will need to disentangle interfacial effects from those arising from changes
in the coexistence conditions. For example, distinguishing between interfacial
and bulk effects is important for comparing nucleation behavior between
stochastic thermodynamics-based [15] and mean-field approaches [16, 17].

Our simulations suggest that measurements of interfacial roughness can serve
as a general method for identifying nonequilibrium effects in driven fluids.
In particular, our analysis indicates that the interfacial tension under
different nonequilibrium coexistence conditions can be inferred quantitatively
regardless of the kinetics of driven chemical reactions. This noninvasive
approach is particularly appealing, as it could be applied directly to living
samples [22] and to “aging” systems whose material properties change over time
[43]. We note that a related strategy has been employed for human cell nuclei
fusion processes, in which the interfacial tension is determined from the
fluctuation–dissipation theorem [44]. This approach is similar to our
proposal, as both strategies are based on the concept of an effective
equilibrium; however, our simulations suggest that our inference approach
could potentially be used under far-from-equilibrium conditions, where
$|\beta\Delta\mu_{\text{coex}}|>1$.

Finally, we note that recently developed experimental settings for studying
chemically driven fluids could provide suitable platforms for testing our
predictions. To realize the inhomogeneous driven chemical reactions that are
essential for observing nonequilibrium interfacial properties, the reaction
kinetics along the driven pathway could be tuned by controlling the
partitioning of enzymes into phase-separated droplets [45] or by engineering
the reactions to be dependent on the local concentration of droplet material
[46, 47]. Investigations of interface roughening could then be used to probe
the nonequilibrium phase behavior of these experimental systems. Since we
anticipate that interface roughening could be observed most easily near the
critical point, it will also be necessary to inspect the critical behavior of
inhomogeneously driven fluids. Such studies present another important
direction for future investigation via theory and simulation.

###### Acknowledgements.

The authors thank Sushant Saryal for helpful discussions and comments on the
manuscript. This work is supported by the National Science Foundation
(DMR-2143670).

## Author Declarations

### Conflict of Interest Statement

The authors have no conflicts to disclose.

##  Appendix A Testing the effective-equilibrium mapping

Figure 7: Chemically driven nonequilibrium fluids cannot be mapped exactly to
effective-equilibrium models. (a–b) In fluids with driven chemical reactions,
the effective internal free-energy difference defined at the level of an
individual lattice site, $\Delta f(u)$, depends on the local potential energy,
$u$, experienced by a B-state particle. The value of $\Delta f(u)$ may differ
from the bulk-phase values, (a) $\Delta f_{\text{l}}$ in the liquid phase and
(b) $\Delta f_{\text{v}}$ in the vapor phase, at coexistence. At equilibrium,
however, lattice configurations obey the Boltzmann distribution, and the two
definitions of the internal free-energy difference coincide. (c) The
distribution of the local environment around I-state particles,
$p_{\text{I}}(u)$, in the coexisting liquid (open symbols) and vapor phases
(solid symbols). The local environment is characterized by the local potential
energy, $u$, that would be experienced if the particle were in the B state.
Gray, orange, and blue colors indicate equilibrium and driven systems with
chemical-reaction models I and II $(\beta\Delta\mu_{\text{coex}}=1.87)$,
respectively. In (b), results for $u=4\epsilon$ are not presented, as such
configurations are extremely rare in the vapor phase.

In this section, we quantitatively assess the mapping between a nonequilibrium
bulk phase and the corresponding effective-equilibrium model. To this end, we
consider an effective internal free-energy difference defined at a microscopic
level with the length scale of an individual lattice site, $\Delta f(u)$, and
compare this with the ensemble-averaged effective internal free-energy
difference, $\Delta f$, defined in Eq. (5). We define $\Delta f(u)$ by mapping
lattice sites at which a B-state particle experiences a local potential energy
$u$ to an effective equilibrium in which the ratio of the number densities for
the B and I-states follows the Boltzmann distribution, i.e.,
$(\rho_{\text{B}}/\rho_{\text{I}})_{u}=\exp[-\beta\Delta f(u)-\beta u]$.
Equivalently, we can write

$\beta\Delta f(u)\equiv-\ln(\rho_{\text{B}}/\rho_{\text{I}})_{u}-\beta u,$
(10)

where the subscript $u$ emphasizes that the number densities are measured at
lattice sites where B-state particles would experience precisely this
potential energy. This definition is closely related to the effective internal
free-energy difference for a bulk phase given by Eq. (5), which may be
considered to be an average of Eq. (10), weighted by the distribution of $u$
across lattice sites in that phase.

Fig. 7a–b show that the internal free-energy differences at the lattice-site
level, $\Delta f(u)$, depend on the local potential energy and are not always
equal to the bulk-phase values, $\Delta f_{\text{l}}$ and $\Delta
f_{\text{v}}$, under nonzero nonequilibrium drive for both chemical-reaction
models. By contrast, the inferred internal free-energy difference is
independent of $u$ at equilibrium. This inconsistency between Eq. (5) and Eq.
(10) demonstrates that the effective-equilibrium mappings for both the liquid
and the vapor phases are not exact, suggesting that exact mappings to
effective-equilibrium models cannot be defined for nonequilibrium fluids in
general.

Despite the inexact nature of the effective-equilibrium approach, the
approximate mapping defined by Eq. (5) nonetheless provides a useful
description of the macroscopic phase behavior because most particles in a
given bulk phase are subject to similar local environments. Fig. 7c shows the
distribution of the local environment around I-state particles,
$p_{\text{I}}(u)$, where $u$ indicates the local potential energy that would
be experienced if the particle were in the B state. This distribution is
sharply peaked at $u=4\epsilon$ in the liquid phase and at $u=0$ in the vapor
phase. Consequently, we are justified in making the approximations $\Delta
f_{\text{l}}\approx\Delta f(4\epsilon)$ and $\Delta f_{\text{v}}\approx\Delta
f(0)$ in our approximate FLEX theory (see Appendix B). Effective equilibrium
is thus a useful conceptual tool as long as the particles in each phase
experience similar local environments.

This effective-equilibrium mapping can be successfully applied to chemically
driven fluids far from a critical point, where the assumption of a constant
local environment in each bulk phase is valid, so that a single effective-
equilibrium model provides a reasonable description of each bulk phase.
However, the concept of an effective equilibrium may no longer be as useful
for fluids near a critical point or when modeling processes such as spinodal
decomposition, where fluctuations cannot be ignored.

##  Appendix B Fixed Local Environment approXimation (FLEX)

In the following sections, we reproduce the description of the “Fixed Local
Environment approXimation" (FLEX) derived in our previous work [15] for the
reader’s convenience. We refer the reader to Ref. [15] for further details. In
FLEX, we assume that fluctuations of of the local configuration, or
“environment”, surrounding a tagged lattice site are negligible. We describe
the local environment in terms of the the potential energy, $u$, experienced
by a B-state particle at the tagged lattice site. In keeping with this small-
fluctuation assumption, we expect that FLEX will be most effective in systems
that are far from a critical point, as discussed in Appendix A.

We use FLEX to map our nonequilibrium model to an effective equilibrium with
identical particle number densities $\tilde{\rho}_{\text{B}}$ and
$\tilde{\rho}_{\text{I}}$. Given a specified $u$, we calculate the steady-
state distribution, $\tilde{\rho}_{i}$ ($i=$ E, B, or I), within a tagged site
site by solving the master equation for the Markovian transition network
presented in Fig. 1b. Given the assumption of a fixed local environment,
$\tilde{\rho}_{\text{B}}$ and $\tilde{\rho}_{\text{I}}$ may regarded as the
number densities of particles in B and I states at the steady-state. The
tagged site is then mapped to an effective-equilibrium model with identical B
and I-state number densities by defining the effective fugacities
$\tilde{z}_{\text{B}}(u)\equiv(\tilde{\rho}_{\text{B}}/\tilde{\rho}_{\text{E}})_{u}\exp(\beta
u)$ and
$\tilde{z}_{\text{I}}(u)\equiv(\tilde{\rho}_{\text{I}}/\tilde{\rho}_{\text{E}})_{u}$,
where the single-site partition function is
$\tilde{\xi}(u)=1+\tilde{z}_{\text{B}}(u)+\tilde{z}_{\text{I}}(u)$.

Because the number densities depend on $u$, the liquid and vapor phases, which
are characterized by different average potential energy, are mapped to
distinct effective equilibria. We quantify the difference in the effective
thermodynamic properties of these phases in terms of an effective internal
free-energy difference,
$\beta\Delta\tilde{f}(u)\equiv-\ln(\tilde{z}_{\text{B}}/\tilde{z}_{\text{I}})_{u}$.
Evaluating this quantity with the effective fugacities leads to Eq. (6) in the
main text. The FLEX prediction for the internal free-energy difference
captures the general trend of how each phase deviates from thermal equilibrium
as the system is driven away from equilibrium. This equation also highlights
the importance of controlling the functional dependence of
$k_{\text{I}\rightarrow\text{B}}$ with respect to the local potential energy
in order to realize inhomogeneous chemical reactions.

FLEX illustrates the roles of $\Delta\mu$ and $\Delta f_{\text{res}}$ in
tuning the relative flux of the two competing reaction pathways between the
particle states. The flux through the direct $\text{I}\rightarrow\text{B}$
pathway is
$j_{\text{I}\rightarrow\text{B}}=\tilde{\rho}_{\text{I}}k_{\text{I}\rightarrow\text{B}}$,
while the flux through the indirect
$\text{I}\rightarrow\text{E}\rightarrow\text{B}$ pathway is
$j_{\text{I}\rightarrow\text{E}\rightarrow\text{B}}=\tilde{\rho}_{\text{I}}/(1+e^{\beta\Delta
f_{\text{res}}})$. Thus, the relative flux between the two pathways is

$j_{\text{I}\rightarrow\text{B}}/{j_{\text{I}\rightarrow\text{E}\rightarrow\text{B}}}=k_{\text{I}\rightarrow\text{B}}(1+e^{\beta\Delta
f_{\text{res}}}).$ (11)

FLEX predicts that along the coexistence line, specified by $\Delta
f_{\text{res}}$ and $\Delta\mu$ (Fig. 3a), this relative flux gradually
changes, leading to deviations in the interfacial properties and phase-
transition kinetics relative to an equilibrium fluid.

FLEX makes quantitative predictions of the phase behavior of the open system
at a given nonequilibrium drive, $\Delta\mu$, by assuming that the
particle–hole symmetry of two-dimensional lattice-gas model is a reasonable
approximation at the effective equilibrium. In the equilibrium lattice-gas
model, the liquid and vapor phases coexist at $\mu=2\epsilon$, where $\mu$ is
the chemical potential of a particle, as a result of this particle–hole
symmetry [27]. The equilibrium supersaturation, $S$, is thus approximately
equal to $\exp[\beta(\mu-2\epsilon)]$, which can be interpreted as the ratio
$\rho/(1-\rho)$ at $u=2\epsilon$, where $\rho$ is the number density of
particles. We take this ratio as the definition of supersaturation $\tilde{S}$
in the FLEX-mapped effective equilibrium,

$\tilde{S}\equiv\left[\dfrac{\tilde{\rho}_{\text{B}}}{1-\tilde{\rho}_{\text{B}}}\right]_{u=2\epsilon}=\dfrac{\tilde{z}_{\text{B}}(2\epsilon;\;\Delta\mu)e^{-2\beta\epsilon}}{\tilde{z}_{\text{I}}(2\epsilon;\;\Delta\mu)+1}.$
(12)

We then relate the FLEX supersaturation, $\tilde{S}$, to the “thermodynamic”
driving force of the nonequilibrium phase transition, $\Delta\Phi$, via
$\tilde{S}\equiv\exp(\beta\Delta\Phi)$, and we associate phase coexistence
with $\tilde{S}=1$. The number density of particles in the vapor phase,
$\rho_{\text{v}}$, is taken to be the sum of the B and I-state number
densities in the effective equilibrium at $u=0$,

$\rho_{\text{v}}\equiv\tilde{\rho}_{\text{B}}(0;\;\Delta\mu)+\tilde{\rho}_{\text{I}}(0;\;\Delta\mu).$
(13)

Here we assume that the vapor phase is sparsely populated with particles and
so interactions among them are negligible. Under a specified chemical drive at
phase coexistence, $\Delta\mu_{\text{coex}}$, Eqs. (12) and (13) uniquely
determine the effective fugacities and effective equilibrium of the open
system.

##  Appendix C Application of FLEX to chemical-reaction models

Eqs. (12) and (13) can be simplified to provide a qualitative understanding of
the nonequilibrium phase behavior of the macroscopic liquid and vapor phases
at coexistence. For the kinetic models employed in this paper, the reaction
rate $k_{\text{I}\rightarrow\text{B}}$ plateaus at low potential energies,
either approaching $k^{\circ}$ in model I or $0$ in model II (Fig. 1c). The
plateau allows us to approximately equate
$k_{\text{I}\rightarrow\text{B}}(4\epsilon)$ and
$k_{\text{I}\rightarrow\text{B}}(2\epsilon)$, meaning that the chemical-
reaction kinetics in the liquid phase can be used to determine phase
coexistence. We further assume that the I state is much more stable than the B
state so that $e^{\beta\Delta f_{\text{res}}}\gg 1$. With these
approximations, we can directly evaluate the coexistence conditions for both
chemical-reaction models.

In the case of chemical-reaction model I, $k_{\text{I}\rightarrow\text{B}}$
approaches zero at $u=0$ (Fig. 1b), which leads to identical fugacities in the
reservoir and in the FLEX approximation of the vapor phase. Under this
condition, the FLEX coexistence criterion, $\tilde{S}=1$, leads to a self-
consistent equation for $\Delta f_{\text{res}}$,

$e^{\beta\Delta
f_{\text{res}}}=\rho_{\text{v}}e^{-2\beta\epsilon}\left[\dfrac{1+k_{\text{I}\rightarrow\text{B}}(4\epsilon)(1+e^{\beta\Delta
f_{\text{res}}})}{1+k_{\text{I}\rightarrow\text{B}}(4\epsilon)(1+e^{\beta\Delta
f_{\text{res}}})e^{\beta\Delta\mu}}\right].$ (14)

Substituting Eq. (14) into Eq. (6) yields $\beta\Delta
f_{\text{l}}=\ln(\rho_{\text{v}}e^{-2\beta\epsilon})$, which remains constant
along the coexistence line. When $\rho_{\text{v}}$ is much higher than the
coexistence condition for the equilibrium two-state lattice-gas model,
$\rho_{\text{v}}e^{-2\beta\epsilon}\gg 1$, this prediction for $\Delta
f_{\text{l}}$ agrees well with the FLEX result at equilibrium, $\beta\Delta
f_{\text{res}}=\beta\Delta f_{\text{l}}=\beta\Delta
f_{\text{eq}}=\ln\left[\rho_{\text{v}}(e^{-2\beta\epsilon}+1)-1\right]$.

In the case of chemical-reaction model II, $k_{\text{I}\rightarrow\text{B}}$
becomes small near $u=2\epsilon$, so that the FLEX fugacities become identical
to the reservoir values at this value of $u$. Applying this condition to Eqs.
(12) and (13) leads to a result that is similar to Eq. (14) in the case of
model I,

$\displaystyle\begin{split}e^{\beta\Delta
f_{\text{res}}}&=\rho_{\text{v}}e^{-2\beta\epsilon}\left[\dfrac{1+k_{\text{I}\rightarrow\text{B}}(0)\left\\{1+e^{\beta(\Delta
f_{\text{res}}+\Delta\mu)}\right\\}}{1+k_{\text{I}\rightarrow\text{B}}(0)(1+e^{\beta\Delta
f_{\text{res}}})e^{\beta\Delta\mu}}\right]\\\
&\approx\rho_{\text{v}}e^{-2\beta\epsilon},\end{split}$ (15)

where the approximation on the second line follows from the assumption
$e^{\beta\Delta f_{\text{res}}}\gg 1$. Eq. (15) indicates that $\Delta
f_{\text{res}}$ is nearly independent of $\Delta\mu_{\text{coex}}$, which
agrees with the results shown in Fig. 3a. Evaluating Eq. (6) leads to
$\beta\Delta f_{\text{l}}\approx\beta\Delta f_{\text{res}}$, since
$k_{\text{I}\rightarrow\text{B}}(4\epsilon)\approx 0$ in this model.

##  Appendix D Predicting interfacial properties with FLEX

We now apply FLEX to predict the nonequilibrium interfacial tension. To this
end, we focus on the effective free-energy cost of attaching a B-state
particle to a flat liquid–vapor interface. Since the free-energy cost of
attaching a B-state adatom is $-\beta\epsilon$ in equilibrium, we propose that
an appropriate effective-equilibrium model of a nonequilibrium interface
should reproduce the relation
$\rho_{\text{B}}=\exp(\beta\epsilon)/[1+\exp(\beta\epsilon)]$ at a tagged
lattice site on the interface. We therefore invert this relationship to
extract the effective dimensionless interaction strength,
$\beta\tilde{\epsilon}$, of an adatom at the interface at steady state,

$\beta\tilde{\epsilon}\equiv\ln\left[\dfrac{\tilde{\rho}_{\text{B}}}{1-\tilde{\rho}_{\text{B}}}\right]_{u=\epsilon}=\ln\left[\dfrac{\tilde{z}_{\text{B}}(\epsilon;\;\Delta\mu_{\text{coex}})}{\tilde{z}_{\text{I}}(\epsilon;\;\Delta\mu_{\text{coex}})+1}\right]-\beta\epsilon,$
(16)

where $\tilde{\rho}_{\text{B}}$ is the FLEX approximation of the B-state
number density at coexistence. The notation $[\cdot]_{u=\epsilon}$ indicates
that the tagged lattice site experiences a fixed local environment consisting
of exactly one B-state particle. Evaluating Eq. (9) using this effective
interaction strength predicts the deviation of the interfacial roughness from
the equilibrium value within the FLEX framework.

From the effective dimensionless interaction strength at the interface,
$\beta\tilde{\epsilon}$, we predict the nonequilibrium interfacial tension,
$\sigma(\tilde{\epsilon})$, using the following analytical formula at
equilibrium [48]:

$\sigma(\tilde{\epsilon})=\\!\sqrt{\dfrac{4\tilde{\epsilon}\beta^{-2}}{\pi\chi(\beta)}\\!\\!\int_{\beta_{c}}^{\beta}\\!\\!\\!K^{\prime}\\!\left(\\!\dfrac{8[\cosh(\beta^{\prime}\tilde{\epsilon})\\!-\\!1]}{[\cosh(\beta^{\prime}\tilde{\epsilon})\\!+\\!1]^{2}}\\!\right)\\!\\!\left[\\!\dfrac{\cosh(\beta^{\prime}\tilde{\epsilon})\\!-\\!3}{\sinh(\beta^{\prime}\tilde{\epsilon})}\\!\right]\\!d\beta^{\prime}}\\!,$
(17)

where $K^{\prime}$ is the elliptic integral of the first kind,
$\chi(\beta)=[1-\sinh^{-4}({\beta\tilde{\epsilon}/2})]^{1/8}$, and
$\beta_{\text{c}}$ is the inverse critical temperature given by
$\beta_{\text{c}}|\tilde{\epsilon}|=2\ln(1+\sqrt{2})$. This prediction
captures the general trend of the interfacial tension under nonequilibrium
conditions with different chemical-reaction kinetics (Fig. 5a).

## References

  * 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).
  * Zwicker [2022] D. Zwicker, “The intertwined physics of active chemical reactions and phase separation,” Curr. Opin. Coll. Inter. Sci. 61, 101606 (2022).
  * Monahan  _et al._ [2017] Z. Monahan, V. H. Ryan, A. M. Janke, K. A. Burke, S. N. Rhoads, G. H. Zerze, R. O’Meally, G. L. Dignon, A. E. Conicella, W. Zheng, R. B. Best, R. N. Cole, J. Mittal, F. Shewmaker, and N. L. Fawzi, “Phosphorylation of the FUS low-complexity domain disrupts phase separation, aggregation, and toxicity,” EMBO J. 36, 2951 (2017).
  * Kim  _et al._ [2019] T. H. Kim, B. Tsang, R. M. Vernon, N. Sonenberg, L. E. Kay, and J. D. Forman-Kay, “Phospho-dependent phase separation of FMRP and CAPRIN1 recapitulates regulation of translation and deadenylation,” Science 365, 825 (2019).
  * Nosella and Forman-Kay [2021] M. L. Nosella and J. D. Forman-Kay, “Phosphorylation-dependent regulation of messenger RNA transcription, processing and translation within biomolecular condensates,” Curr. Opin. Cell Biol. 69, 30 (2021).
  * Berry, Brangwynne, and Haataja [2018] J. Berry, C. P. Brangwynne, and M. Haataja, “Physical principles of intracellular organization via active and passive phase transitions,” Rep. Prog. Phys. 81, 046601 (2018).
  * Zwicker, Hyman, and Jülicher [2015] D. Zwicker, A. A. Hyman, and F. Jülicher, “Suppression of Ostwald ripening in active emulsions,” Phys. Rev. E 92, 012317 (2015).
  * 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).
  * Kirschbaum and Zwicker [2021] J. Kirschbaum and D. Zwicker, “Controlling biomolecular condensates via chemical reactions,” J. R. Soc. Interface 18, 20210255 (2021).
  * Tena-Solsona  _et al._ [2021] M. Tena-Solsona, J. Janssen, C. Wanzke, F. Schnitter, H. Park, B. Rieß, J. M. Gibbs, C. A. Weber, and J. Boekhoven, “Accelerated ripening in chemically fueled emulsions,” Chem. Sys. Chem. 3, e2000034 (2021).
  * 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,” Nature Phys. 13, 408 (2017).
  * 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).
  * Bauermann  _et al._ [2022] J. Bauermann, S. Laha, P. M. McCall, F. Jülicher, and C. A. Weber, “Chemical kinetics and mass action in coexisting phases,” J. Am. Chem. Soc. 144, 19294 (2022).
  * Donau  _et al._ [2022] C. Donau, F. Späth, M. Stasi, A. M. Bergmann, and J. Boekhoven, “Phase transitions in chemically fueled, multiphase complex coacervate droplets,” Angew. Chem. Int. Ed. 61, e202211905 (2022).
  * Cho and Jacobs [2023] Y. Cho and W. M. Jacobs, “Tuning nucleation kinetics via nonequilibrium chemical reactions,” Phys. Rev. Lett. 130, 128203 (2023).
  * Cates and Nardini [2023] M. E. Cates and C. Nardini, “Classical nucleation theory for active fluid phase separation,” Phys. Rev. Lett. 130, 098203 (2023).
  * Ziethen, Kirschbaum, and Zwicker [2023] N. Ziethen, J. Kirschbaum, and D. Zwicker, “Nucleation of chemically active droplets,” Phys. Rev. Lett. 130, 248201 (2023).
  * Nott  _et al._ [2015] T. J. Nott, E. Petsalaki, P. Farber, D. Jervis, E. Fussner, A. Plochowietz, T. D. Craggs, D. P. Bazett-Jones, T. Pawson, J. D. Forman-Kay, and A. J. Baldwin, “Phase transition of a disordered nuage protein generates environmentally responsive membraneless organelles,” Mol. Cell 57, 936 (2015).
  * Söding  _et al._ [2020] J. Söding, D. Zwicker, S. Sohrabi-Jahromi, M. Boehning, and J. Kirschbaum, “Mechanisms for active regulation of biomolecular condensates,” Trends Cell Biol. 30, 4 (2020).
  * Besse  _et al._ [2023] M. Besse, G. Fausti, M. E. Cates, B. Delamotte, and C. Nardini, “Interface roughening in nonequilibrium phase-separated systems,” Phys. Rev. Lett. 130, 187102 (2023).
  * O’Flynn and Mittag [2021] B. G. O’Flynn and T. Mittag, “The role of liquid–liquid phase separation in regulating enzyme activity,” Curr. Opin. Cell Biol. 69, 70 (2021).
  * Schuster  _et al._ [2021] B. S. Schuster, R. M. Regy, E. M. Dolan, A. Kanchi Ranganath, N. Jovic, S. D. Khare, Z. Shi, and J. Mittal, “Biomolecular condensates: Sequence determinants of phase separation, microstructural organization, enzymatic activity, and material properties,” J. Phys. Chem. B 125, 3441 (2021).
  * Seifert [2012] U. Seifert, “Stochastic thermodynamics, fluctuation theorems and molecular machines,” Rep. Prog. Phys. 75, 126001 (2012).
  * Van den Broeck and Esposito [2015] C. Van den Broeck and M. Esposito, “Ensemble and trajectory thermodynamics: A brief introduction,” Physica A 418, 6 (2015).
  * Wilding [1995] N. B. Wilding, “Critical-point and coexistence-curve properties of the Lennard-Jones fluid: A finite-size scaling study,” Phys. Rev. E 52, 602 (1995).
  * Gillespie [2007] D. T. Gillespie, “Stochastic simulation of chemical kinetics,” Annu. Rev. Phys. Chem 58, 35 (2007).
  * Pathria [1996] R. K. Pathria,  _Statistical mechanics_ , 2nd ed. (Butterworth-Heinemann, Oxford, 1996).
  * Chandler [1987] D. Chandler,  _Introduction to modern statistical mechanics_ (Oxford University Press, New York, 1987).
  * Hansen and McDonald [2013] J.-P. Hansen and I. R. McDonald,  _Theory of simple liquids: with applications to soft matter_ , 4th ed. (Academic press, Oxford, 2013).
  * Solon  _et al._ [2015] A. P. Solon, Y. Fily, A. Baskaran, M. E. Cates, Y. Kafri, M. Kardar, and J. Tailleur, “Pressure is not a state function for generic active fluids,” Nature Phys. 11, 673 (2015).
  * Omar  _et al._ [2023] A. K. Omar, H. Row, S. A. Mallory, and J. F. Brady, “Mechanical theory of nonequilibrium coexistence and motility-induced phase separation,” Proc. Natl. Acad. Sci. USA 120, e2219900120 (2023).
  * Warmflash, Bhimalapuram, and Dinner [2007] A. Warmflash, P. Bhimalapuram, and A. R. Dinner, “Umbrella sampling for nonequilibrium processes,” J. Chem. Phys. 127, 154112 (2007).
  * Wilding and Bruce [1992] N. B. Wilding and A. D. Bruce, “Density fluctuations and field mixing in the critical fluid,” J. Phys.: Condens. Matter 4, 3087 (1992).
  * Noya, Vega, and de Miguel [2008] E. G. Noya, C. Vega, and E. de Miguel, “Determination of the melting point of hard spheres from direct coexistence simulation methods,” J. Chem. Phys. 128, 154507 (2008).
  * Christensen, Elder, and Fogedby [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).
  * Allen and Tildesley [2017] M. P. Allen and D. J. Tildesley,  _Computer simulation of liquids_ , 2nd ed. (Oxford University Press, New York, 2017).
  * Oxtoby [1992] D. W. Oxtoby, “Homogeneous nucleation: Theory and experiment,” J. Phys.: Condens. Matter 4, 7627 (1992).
  * Ryu and Cai [2010] S. Ryu and W. Cai, “Validity of classical nucleation theory for Ising models,” Phys. Rev. E 81, 030601 (2010).
  * Allen, Valeriani, and ten Wolde [2009] R. J. Allen, C. Valeriani, and P. R. ten Wolde, “Forward flux sampling for rare event simulations,” J. Phys.: Condens. Matter 21, 463102 (2009).
  * Chaudhury and Makarov [2010] S. Chaudhury and D. Makarov, “A harmonic transition state approximation for the duration of reactive events in complex molecular arrangements,” J. Chem. Phys. 133, 034118 (2010).
  * Hummer [2004] G. Hummer, “From transition paths to transition states and rate coefficients,” J. Chem. Phys. 120, 516 (2004).
  * Saito [1996] Y. Saito,  _Statistical physics of crystal growth_ (World Scientific, Singapore, 1996).
  * Jawerth  _et al._ [2020] L. Jawerth, E. Fischer-Friedrich, S. Saha, J. Wang, T. Franzmann, X. Zhang, J. Sachweh, M. Ruer, M. Ijavi, S. Saha, J. Mahamid, A. A. Hyman, and F. Jülicher, “Protein condensates as aging Maxwell fluids,” Science 370, 1317 (2020).
  * Caragine, Haley, and Zidovska [2018] C. M. Caragine, S. C. Haley, and A. Zidovska, “Surface fuctuations and coalescence of nucleolar droplets in the human cell nucleus,” Phys. Rev. Lett. 121, 148101 (2018).
  * Saleh, Jeon, and Liedl [2020] O. A. Saleh, B.-j. Jeon, and T. Liedl, “Enzymatic degradation of liquid droplets of DNA is modulated near the phase boundary,” Proc. Natl. Acad. Sci. U.S.A. 117, 16160 (2020).
  * Späth  _et al._ [2021] F. Späth, C. Donau, A. M. Bergmann, M. Kränzlein, C. V. Synatschke, B. Rieger, and J. Boekhoven, “Molecular design of chemically fueled peptide–polyelectrolyte coacervate-based assemblies,” J. Am. Chem. Soc. 143, 4782 (2021).
  * Nakashima  _et al._ [2021] K. K. Nakashima, M. H. van Haren, A. A. M. André, I. Robu, and E. Spruijt, “Active coacervate droplets are protocells that grow and resist Ostwald ripening,” Nat. Comm. 12, 3819 (2021).
  * Shneidman, Jackson, and Beatty [1999] V. A. Shneidman, K. A. Jackson, and K. M. Beatty, “On the applicability of the classical nucleation theory in an Ising system,” J. Chem. Phys. 111, 6932 (1999).