Tuning nucleation kinetics via nonequilibrium chemical reactions

Unlike fluids at thermal equilibrium, biomolecular mixtures in living systems can sustain nonequilibrium steady states, in which active processes modify the conformational states of the constituent molecules. Despite qualitative similarities between liquid-liquid phase separation in these systems, the extent to which the phase-separation kinetics differ remains unclear. Here we show that inhomogeneous chemical reactions can alter the nucleation kinetics of liquid-liquid phase separation in a manner that is consistent with classical nucleation theory, but can only be rationalized by introducing a nonequilibrium interfacial tension. We identify conditions under which nucleation can be accelerated without changing the energetics or supersaturation, thus breaking the correlation between fast nucleation and strong driving forces that is typical of phase separation and self-assembly at thermal equilibrium.

Unlike fluids at thermal equilibrium, biomolecular mixtures in living systems can sustain nonequilibrium steady states, in which active processes modify the conformational states of the constituent molecules. Despite qualitative similarities between liquid-liquid phase separation in these systems, the extent to which the phase-separation kinetics differ remains unclear. Here we show that inhomogeneous chemical reactions can alter the nucleation kinetics of liquid-liquid phase separation in a manner that is consistent with classical nucleation theory, but can only be rationalized by introducing a nonequilibrium interfacial tension. We identify conditions under which nucleation can be accelerated without changing the energetics or supersaturation, thus breaking the correlation between fast nucleation and strong driving forces that is typical of phase separation and self-assembly at thermal equilibrium.
In living systems, phase separation can occur at a nonequilibrium steady state (NESS) as opposed to thermal equilibrium [1]. For example, in active intracellular condensates, biomolecules may be degraded or posttranslationally modified by enzymes that couple conformational changes to the conversion of a chemical fuel, such as ATP, to chemical waste [2]. A variety of mesoscale models have recently been proposed to study the phase behavior of chemically driven fluids [3][4][5][6][7]. A crucial feature of these systems is the potential for inhomogeneous chemical reactions [7,8], meaning that the reactions are driven out of equilibrium to a different extent in each phase of a phase-separating mixture. Such behavior can result from enzymes that preferentially localize to one phase or from chemical fuel gradients that couple to the local density of the phase-separating molecules [2]. Phase separation taking place at an inhomogeneous NESS can exhibit qualitatively different features compared to thermal equilibrium, including suppressed coarsening, monodisperse phase-separated droplet size distributions, and even spontaneous droplet division [9,10].
Inhomogeneous chemical reactions can also affect the kinetics of phase transitions, although the extent to which kinetic pathways at a NESS differ from those at equilibrium is not well understood. Not only do inhomogeneous reactions provide additional control parameters beyond temperature and concentration with which to control a phase transition, but they might also alter the mechanism of phase separation. This possibility contrasts with the typical behavior of equilibrium phaseseparating fluids, in which strong thermodynamic driving forces are typically necessary to initiate homogeneous nucleation at equilibrium unless the system is near a critical point [11,12]. The consequences of this correlation between thermodynamics and nucleation kinetics are well appreciated in the context of molecular self-assembly, especially in cases where strong driving forces are associated with kinetic trapping [13][14][15][16][17]. It has been speculated * wjacobs@princeton.edu that living systems must contend with similar trade-offs in order to harness phase separation for biological functionality [1,18,19].
Here we show that inhomogeneous chemical reactions provide a mechanism to alter the nucleation pathway of a nonequilibrium phase-separating fluid. We first introduce a simulation strategy to compute phase coexistence and nucleation rates in a generic model of liquid-vapor phase separation in systems with inhomogeneous chemical reactions. We then show that the nucleation behavior cannot be described by an equilibrium theory, as the inhomogeneous reactions modify the material properties at the liquid-vapor interface. We rationalize this surprising result using a simple theoretical model and establish how emergent, nonequilibrium interfacial properties can tune the kinetics of phase separation and self-assembly far from equilibrium.
In order to study nonequilibrium phase separation via molecular simulation, we adopt the framework of stochastic thermodynamics [20,21] and consider an open system connected to a particle reservoir at constant chemical potentials (Fig. 1a). For simplicity, we study a twodimensional square lattice model, in which empty lattice sites represent either a vacancy or solvent. To allow for chemical reactions, we assume that the particles have two internal states: a bonding state (B) that promotes phase separation due to attractive nearest neighbor interactions with bond energy < 0, and an inert state (I) that is isoenergetic to an empty lattice site. The internal free-energy difference between the internal states in the reservoir is ∆f res , so that the fugacities of the two states are related by z B /z I = exp(−β∆f res ), where β ≡ (k B T ) −1 . Our model is thus closely related to the equilibrium lattice gas, which exhibits a first-order phase transition between a dilute vapor (v) phase and a condensed liquid (l) phase below a critical temperature [22]. However, unlike the equilibrium lattice gas, the particle transitions between the open system and the reservoir in our model may break detailed balance. As a result, the product of rates for inserting a bonding particle, changing its internal state, and returning it to the reservoir may differ from the product of rates for the reversed se- quence of events by a factor exp(β∆µ), where ∆µ is the chemical potential difference between the chemical fuel and waste used to drive the reactions between internal states in the open system. Nonetheless, we note that performing simulations in an open system has the same advantages as studying equilibrium phase separation in the grand-canonical ensemble, including the elimination of interfaces and a resulting reduction of finite size effects [23,24].
We implement chemical reactions by assuming Markovian transitions, such that internal state changes B I obey local detailed balance, k B→I /k I→B = exp(βu + β∆f res + β∆µ). The local potential energy, u, is equal to the bond energy, , times the number of nearest-neighbor bonding-state particles, and ∆µ is held constant. We obtain homogeneous reactions if we set the backward rate for I → B transitions, k I→B , equal to a constant k • , which represents the ratio between the timescales for internal state changes and particle diffusion. To account for inhomogeneous chemical reactions, we implement I → B transitions such that k I→B depends on the local potential energy. For example, this scheme can describe a unimolecular reaction that may proceed either via an undriven pathway or via a catalyzed pathway that is coupled to the conversion of chemical fuel to waste, with the relative flux through each pathway dependent on the local environment [7]. Specifically, we assume that k I→B takes the Metropolis form We emphasize that the qualitative features of our results do not depend on this specific functional form of k I→B . Instead, the essential feature in the case of inhomogeneous reactions is whether the driven reaction is promoted or suppressed by decreasing the local potential energy. For example, if the reactive flux through the driven pathway is enhanced at low potential energies, then the model can describe either preferentially driven deactivation (B → I) in the liquid phase or preferentially driven activation (I → B) in the vapor phase.
It is useful to characterize the resulting NESS in each phase by introducing an effective internal free-energy difference, β∆f ≡ − log(ρ B /ρ I ) + log exp(−β∆u IB ) I , where ρ B and ρ I are the steady-state number densities of particles in the B and I internal states, respectively, and the second term represents an average of the potential energy change due to converting an I to a B particle at steady state. The explicit dependence of k I→B on the local potential energy in our inhomogeneous reaction model causes ∆f to differ between the liquid and the vapor phases, such that ∆∆f ≡ ∆f l − ∆f v = 0 ( Fig. 1a,d). Although there is no exact mapping between the nonequilibrium lattice gas and an equilibrium model in general, measuring ∆∆f reveals the crucial difference between the inhomogeneous and homogeneous models, the latter of which always has ∆∆f = 0 (see SI). We return to this concept in our theoretical analysis below.
Using kinetic Monte Carlo simulations [25], with particle exchanges as well as internal state switches, we verify the modifications to the lattice gas phase behavior anticipated above. We first identify the conditions for nonequilibrium phase coexistence by equating the total probability of being in the vapor versus the liquid phase at steady state. This is analogous to the equal pressure construction in equilibrium grand-canonical phasecoexistence simulations [23], and further implies that the open system transitions between the liquid and vapor phases with equal forward and backward rates. To this end, we use a form of nonequilibrium umbrella sampling (NEUS) [26] to calculate the steady-state probability as a function of the number density of bonding particles, p(ρ B ). As is characteristic of a first-order phase transition, we observe a barrier with respect to − ln p(ρ B ) that scales with the lattice length L as the system size is increased [27]. Based on the value of the order parameter ρ * B at the top of this barrier, we determine the steady-state probabilities of the vapor and liquid phases, p v ≡ p(ρ B < ρ * B ) and p l ≡ p(ρ B > ρ * B ), respectively (Fig. 1b). We then define the dimensionless thermodynamic driving force between bulk phases to be β∆Φ(L) ≡ L −2 ln (p l /p v ) and associate phase coexistence with ∆Φ = 0 (Fig. 1c,d). We also verify that this criterion is consistent with mechanical equilibrium between the inhomogeneously driven phases by simulating the diffusion of a planar liquid-vapor interface in direct coexistence simulations (see SI). Interestingly, we find that the coexistence point is subject to a fi- nite size effect, in which the deviation of ∆Φ(L) from its value in the thermodynamic limit scales as 1/L. We therefore compute ∆Φ(L) from a series of NEUS simulations and extrapolate to the thermodynamic limit, ∆Φ ≡ ∆Φ(L → ∞). We attribute this effect, which does not appear in the equilibrium lattice gas, to the broken particle-hole symmetry that results from inhomogeneous reactions between the vapor and liquid phases.
We are now able to address the central question of this work: To what extent can equilibrium descriptions of nucleation be applied to phase separation at a NESS? The most widely used theoretical framework for describing nucleation in systems ranging from atomic and molecular fluids to colloidal and biomolecular materials is classical nucleation theory (CNT) [11,12]. In its most general form, equilibrium CNT predicts that nucleation follows a minimum free-energy pathway along a reaction coordinate corresponding to the size of a nucleus of the stable phase. This pathway crosses a free-energy barrier that arises from the competition between the lower thermodynamic potential of the stable phase and the positive interfacial free energy between the nucleus and the bulk metastable phase. CNT predicts that the homogeneous nucleation rate density is the product of a prefactor and a Boltzmann factor corresponding to the height of the barrier, ∆F * ; the prefactor is the product of the monomer number density, ρ 1 ; the speed along the reaction coordinate at the top of the barrier, D * ; and the Zeldovich factor, Γ, that accounts for fluctuations that cross the barrier but return to the metastable state (Fig. 2a). After taking into account the interfacial free energy due to the macroscopic line tension and microscopic nucleus size fluctuations, CNT has been shown to provide a quantitative description of nucleation in the two-dimensional equilibrium lattice gas model [28].
To compute nucleation rates at a NESS, we employ forward-flux sampling (FFS) [29], using the largest cluster size, n, as the reaction coordinate. FFS uses short simulations between milestones along the reaction coordinate to compute the commitment probability to the stable phase, φ(n), and the nucleation rate density, I. The critical nucleus size, n * , is found where φ(n * ) = 1 /2 [30], and the Zeldovich factor can be calculated by fitting φ(n) to an approximately harmonic barrier in the vicinity of n * . We also independently measure the number density of bonding-state monomers in the vapor phase, ρ 1 , and the diffusion coefficient, D * , from nucleus size fluctuations near n * [31]. Finally, we define the supersaturation, S ≡ exp(β∆Φ). We are therefore able to isolate the factor in the CNT rate equation that pertains to the (nonequilibrium) nucleation barrier by computing ln(I/ρ 1 D * Γ) as a function of S, which we control by tuning ∆µ at constant ∆f res .
We first test the prediction of the fundamental nucleation theorem, n * = −∂ ln(I/ρ 1 D * Γ)/∂ ln S + 1, for nucleating a stable liquid phase from a supersaturated vapor phase [12]. This prediction holds as long as the interfacial free energy is independent of the supersaturation, regardless of the functional form of the nucleation barrier. The results of representative simulations shown in the inset of Fig. 2b demonstrate excellent agreement between the critical nucleus sizes obtained from FFS, n * FFS , and the sizes inferred from this theorem, n * nuc. thm. . This provides evidence that the fundamental premise of CNTnamely, that the rate-limiting step coincides with the formation of a critical nucleus of the stable bulk phaseapplies to nucleation at a NESS in the regime β∆Φ 1. The validity of this theorem also provides further support for our definition of the supersaturation.
However, when examining the supersaturation dependence of the apparent nucleation barrier (Fig. 2b), we discover a surprising deviation from the equilibrium lattice gas: Although the interfacial contribution still scales with the perimeter of the two-dimensional nucleus, the inferred line tension, σ, differs from the equilibrium value, σ eq . This deviation only occurs in the case of inhomogeneous reactions, which can be seen by comparing the homogeneous and inhomogeneous results with the equilibrium barrier height in Fig. 2b. These observations indicate that CNT can be extended to describe phase separation at a NESS, but that the apparent nucleation barrier can differ by orders of magnitude from predictions based on equilibrium interfacial properties in the case of nonequilibrium inhomogeneous reactions.
Our simulations confirm that the liquid-vapor interfacial properties are influenced by how far the system is driven out of equilibrium. As a concrete example, we perform simulations corresponding to an isothermal experiment in which the total number of particles is conserved, such that β and ρ v are held constant. We find that the line tension for inhomogeneous reactions deviates farther from the equilibrium value as β∆µ at coexistence is increased (Fig. 3a). We also find that the interfacial properties depend on the relative timescale, k • , between I B transitions and the rate of particle attachment to the nucleus, normalized by its perimeter. To test the sensitivity of the line tension to the ratio of these timescales, we calculate the line tension deviation, ∆σ, while holding ∆f res + ∆µ constant (inset of Fig. 3a). We find that ∆σ is nonzero over a wide range of k • , with the greatest deviation occurring when the timescales are comparable (k • ≈ 1). However, we recover the equilibrium line tension in the limit of either zero reactive flux (k • → 0) or infinitely fast reactions (k • → ∞), as the system reverts either to a true equilibrium or to a NESS in which ∆∆f → 0 (see SI).
We rationalize these results by considering an analytical model that captures the qualitative behavior of the nonequilibrium nucleation kinetics. We make the approximation that particle exchange between the open system and the reservoir relaxes to the steady-state distribution more quickly than the local environment around a particle changes. This approximation provides a good description of both phases at low temperatures arbitrarily far from equilibrium, since in this limit the number densities of particles and vacancies in the vapor and liquid phases, respectively, are extremely low, even though the probability currents associated with futile chemical-reaction cycles can be arbitrarily large. For inhomogeneous models, calculation of the steady-state distribution reveals that ∆∆f = 0 in this ideal limit as expected. To predict the line tension of the liquid-vapor interface at a NESS, we use this "fixed local environment" approximation for the steady-state distribution to calculate an effective free energy for attaching a bonding-state adatom to the interface. We then employ the standard solid-on-solid approximation to calculate the line tension resulting from the effective adatom free energy [32]. The predictions of this theory are shown by solid curves in Fig. 3.
Our key insight from this theory is that nonequilibrium interfacial properties emerge when the bulk phases and the particles at the liquid-vapor interface are described by different effective equilibrium models. The net result is that bonding particles attached to the interface of a critical nucleus may be attracted either more or less strongly, per nearest-neighbor interaction, than in the bulk phases. This simple theory captures both the sign and the approximate functional form of ∆σ with respect to β∆µ coex . It also predicts the nonmonotonic dependence of ∆σ on the relative reaction timescale k • . While the precise form of ∆σ depends on our choice of the potential-energy-dependent backward rate, k I→B , we emphasize that our theory predicts a nonzero ∆σ only if the effective internal free-energy difference, ∆f , is a function of the local potential energy.
Finally, to highlight the control over nucleation rates imparted by inhomogeneous reactions, we compare the nonequilibrium nucleation rate to that of an equilibrium fluid with the same β , S, and ρ v (Fig. 3b). In agreement with our theory, our simulations show that the nucleation rate can be increased by orders of magnitude relative to the corresponding equilibrium system by driving the fluid far from equilibrium (β∆µ 1). The magnitude of this effect is far greater in the inhomogeneous than in the homogeneous model due to the dominant role of the line tension in determining the nucleation barrier, and thus the nucleation rate. Inhomogeneous reactions can therefore break the usual relationship between high supersaturation and fast nucleation, offering a novel way to control nucleation kinetics in nonequilibrium fluids.
In conclusion, we have introduced a strategy for simulating nonequilibrium phase coexistence and phase transformation kinetics within the framework of stochastic thermodynamics. Although the basic premise of CNT applies to the nonequilibrium phase-separating fluids that we have considered, inhomogeneous reactions can result in interfacial tensions that differ substantially from the equilibrium value. Our findings provide further evidence [1,4] that nonequilibrium phase transformations may follow the same scaling laws as equilibrium systems under rather general conditions. Detecting nonequilibrium effects may thus require careful measurements of the interfacial material properties. Our work also reveals a mechanism for decoupling nucleation rates from thermodynamic driving forces at a NESS.
Our results are applicable to a range of experimental systems broadly described as living or active. For example, it has been proposed that stress granules, among other biological condensates, are kept out of equilibrium by inhomogeneous enzyme distributions [2,33,34]. Our prediction of a reduced nonequilibrium surface tension could be tested either in such endogenous systems or using light-activated corelets [35], which have recently been applied to study the nucleation of intracellular condensates [19]. Our model could also be applied to synthetic active polypeptide coacervates [36][37][38] or DNA liquids [39] in which the association/hybridization reactions are engineered to respond to energy input in a manner that is dependent on the local protein/DNA concentration. In both contexts, our results suggest a road map for controlling self-assembly kinetics far from thermal equilibrium.
This work is supported by the National Science Foundation (DMR-2143670). We extend the two-dimensional square lattice-gas model by incorporating two particle internal states: a bonding state (B) and an inert state (I). A particle in the bonding state interacts with nearest-neighbor bonding-state particles with bonding strength < 0. On the other hand, an inert-state particle is isoenergetic to a vacant lattice site and thus does not interact with nearest-neighbor particles.
We consider an open system in contact with a particle reservoir. The fugacities in the reservoir are z B and z I for the B and I internal states, respectively. Utilizing the framework of stochastic thermodynamics, we model the kinetics of particle insertion, removal, and reactions between internal states using Markovian transitions that obey local detailed balance. Particle insertion into an empty (E) lattice site occurs with rates Dz B and Dz I for the B and I internal states, respectively, where D is a rate related to particle exchange between the open system and the reservoir (or the diffusive transport rate in a closed system). Particle removal occurs with rate De βu , where β ≡ (k B T ) −1 and u is the local potential energy that arises from the nearest-neighbor interactions at that lattice site. Lastly, reactions between the B and I states occur with forward and backward rates k BI and k IB . For notational simplicity, we write the dimensionless ratios between the reaction rates and the particle exchange rate as k B→I ≡ D −1 k BI and k I→B ≡ D −1 k IB , respectively, in what follows.
Since the transitions at any particular lattice site involve a single-cycle transition network among E, B, and I states, we can define a nonequilibrium drive ∆µ along the cycle in the B-to-I direction, Nonzero ∆µ results in a nonzero net probability current. Rearranging Eq. (S1) into the ratio of internal reaction rates gives the local detailed balance condition for I B reactions in terms of ∆µ,

B. Model specification and simulation parameters
In order to conduct numerical simulations, we consider two specific choices for the backward reaction rate, k I→B : where the constant k • sets the relative timescale between the reaction rates and the particle exchange rate with the reservoir. Note that k B→I follows from the local detailed balance condition, Eq. (S2). In the homogeneous model, k I→B is independent of u, while in the inhomogeneous model, k I→B depends on u. The motivation for studying u-independent and u-dependent models is described in detail below. We simulate the stochastic evolution of the system via the kinetic Monte Carlo (kMC) simulation method [25] using the first-order transition rates described above. For convenience, we define intrinsic fugacities λ B and λ I , which are related to the reservoir fugacities by z B ≡ λ B exp(β∆µ) and z I ≡ λ I , so that we can tune the chemical drive while holding the equilibrium free-energy difference between the B and I states constant.
Unless specified otherwise, all simulation data are obtained from 64 × 64 two-dimensional square lattices with periodic boundaries, interaction strength β = −2.95, and relative timescale k • = 10 −1 . The intrinsic fugacities are set such that the total particle density (i.e., the number density of particles in either state) in the vapor phase is ρ v = 0.05 at the coexistence point. The parameters λ B , λ I , and ∆µ are then held constant throughout the course of a simulation.

C. Fixed Local Environment approXimation (FLEX)
To motivate our discussion of inhomogeneous chemical reactions, we introduce a simplified theoretical model that we refer to as the Fixed Local Environment approXimation (FLEX). Within FLEX, we assume that particle exchange between the open system and the reservoir relaxes to the steady state more rapidly than any change in the local configuration, or environment, around a given lattice site. For the two-dimensional square lattice-gas model, the local configuration comprises the four nearest-neighbor lattice sites, which determine the local potential energy u when the lattice site is occupied by a bonding-state particle. Following the convention that W ij indicates the first-order transition rate from state i to state j, we write the transition matrix for a single tagged lattice site with a fixed local environment as where the lattice-site states are ordered (E, B, I). FLEX allows us to predict the distribution among the three lattice-site states at steady state under a constant u determined from the fixed local environment. Solving the master equation at steady state, 0 =ṗ = pW , leads to where p i is the steady-state probability of being in state i, and z B and z I are fugacity-like values that may depend on u, according to the specific functional form of k I→B . However, at equilibrium, z B = z B and z I = z I regardless of k I→B . The apparent internal free-energy difference in the open system, β∆f ≡ − ln(z B /z I ), is related to the free-energy difference in the reservoir, β∆f res ≡ − ln(z B /z I ), by Within FLEX, it is therefore possible to map the non-equilibrium steady state for a single lattice site to an "effective equilibrium" model that has the same steady-state density distribution (but different probability currents). We discuss the implications of this insight in the next section, and we describe the predictions of FLEX for nonequilibrium nucleation kinetics in Sec. IV.

D. Effect of inhomogeneous chemical reactions
We now consider the effect of the internal reaction kinetics on the thermodynamics of a fluid that phase separates into dilute (vapor) and condensed (liquid) phases. Within FLEX, Eq. (S7) shows that under nonequilibrium conditions (∆µ = 0), the effective internal free-energy difference between B and I states in the open system, ∆f , may be different from the free-energy difference in the reservoir, ∆f res . However, as long as k I→B is the same in both the vapor and liquid phases, a single effective equilibrium model provides a common description of the steady-state distribution in both phases. We refer to such a system as having homogeneous internal-state chemical reactions. However, if k I→B is dependent on the local potential energy, which is on average higher in the vapor phase than in the liquid phase, then the effective equilibrium descriptions must be different in the two phases. As a result, a system with a u-dependent k I→B results in inhomogeneous internal-state chemical reactions. The steady-state density distribution in a phase-separated system with inhomogeneous chemical reactions cannot be described by a common effective equilibrium.
Our numerical results indicate that these insights provided by FLEX are also useful for analyzing the full lattice model. Motivated by the FLEX analysis, we quantify the extent of inhomogeneous chemical reactions by estimating the effective free-energy difference between the two particle internal states, ∆f , from simulations of each bulk phase. In order to calculate ∆f for a system at a nonequilibrium steady state (NESS), we start from a grand-canonical model with particle fugacities z B and z I . A lattice configuration is defined by the identities of all the lattice sites, {c(r)}, where c ∈ {E, B, I}. In a translationally symmetric system, we can assume that a tagged particle is located at the origin, r = 0. Including the empty lattice site as an "internal state," E, with z E = 1, we can write the equilibrium probability of the tagged particle being in internal state i as where primed summation is over nearest-neighboring pairs (counting each unique bond once), double primed summation is over the nearest-neighbor sites r of r = 0, and the angle brackets indicate an ensemble average conditioned on the particle at r = 0 being in the indicated internal state. Eq. (S10) can be viewed as a Bennet acceptance ratio in the semi-grand ensemble or a generalization of the Widom insertion method. Finally, we use Eq. (S10) to calculate the effective β∆f ≡ − ln(z B /z I ) in a NESS by substituting p eq with the NESS distribution, p, and averaging over the ensemble of lattice configurations at steady state: In the last step, we have exploited the isoenergeticity of the inert state. Note that Eq. (S11) reduces to Eq. (S7) if we assume a fixed local environment around the tagged lattice site. We can numerically test the effect of applying nonequilibrium drive to the models with homogeneous and inhomogeneous chemical reactions, Eq. (S3), by measuring the effective internal free-energy differences in coexisting liquid and vapor phases, ∆f l and ∆f v , respectively. Fig. S4 shows the free-energy differences calculated according to Eq. (S11) using the NESS distribution obtained from simulations in each phase. The measured values of the effective internal free-energy differences, ∆f l and ∆f v , for nonequilibrium homogeneous systems are identical in the liquid and vapor phases as expected. By contrast, the effective internal free-energy difference in the dilute phase, ∆f v , exhibits a monotonic decrease with respect to the nonequilibrium drive in inhomogeneous systems. The FLEX predictions for ∆f v and ∆f l obtained by evaluating Eq. (S7) at fixed u = 0 and u = 4 , respectively, show the same trends for both nonequilibrium models. It is important to note that the existence of an effective equilibrium that is the same in both phases does not mean that a driven, homogeneous system corresponds to a true equilibrium, because the entropy production rate is always positive due to particle exchange with the reservoir.

II. PHASE COEXISTENCE CALCULATIONS
A. Nonequilibrium Umbrella Sampling (NEUS) We use a form of nonequilibrium umbrella sampling (NEUS) [26] to obtain the steady-state distribution of the bonding-state particle density, ρ B . For this purpose, we divide the entire range of ρ B into non-overlapping boxes and focus on the transition flux between the boxes at steady state. Importantly, transitions are only possible between adjacent boxes in our simulations because each kMC move can insert or remove at most one bonding-state particle. Detailed balance between the boxes always holds under this restriction, regardless of whether the system is driven out of equilibrium, as long as the system is at steady state.
from the left and the right boundaries, respectively, of each box. We then record the average time, t(b), until each of the launched trajectories exits box b; the average time t(ρ B ; b) that a trajectory spends at ρ B within box b; and the probability that a trajectory exits through the i-side boundary, The algorithm is implemented by iteratively obtaining a new ensemble of configurations at each boundary of every box while calculating the steady-state distribution within each box. The initial configurations at each boundary are sampled from brute-force simulations inside of each box, rejecting any kMC events that would allow the system to cross the box boundaries. At each subsequent iteration of the algorithm, we first obtain the transition probabilities, {p ij }, in each box starting from the current ensemble of configurations at the box boundaries. We then compute the steady-state distribution within each box using trajectories launched from the current ensemble of configurations and the calculated {p ij }. We save the configurations from these trajectories that exit the box to form the ensemble for the next iteration of the algorithm. Finally, we average the steady-state distribution and the transition probabilities within each box over the successive iterations and reconstruct the steady-state distribution over the complete range of ρ B as described above.
We apply this algorithm to measure the dimensionless thermodynamic driving force between the bulk phases, β∆Φ ≡ L −2 ln (p l /p v ), as discussed in the main text. Representative results of this algorithm are shown in Fig. S5. The quick decay of ∆Φ during the initial iterations in Fig. S5A indicates that the system rapidly relaxes, after which the ensemble of trajectories remains in the steady state. To verify that NEUS converges to the correct steady-state distribution, we also performed NEUS for an equilibrium lattice gas and confirmed that the results match a cluster expansion [40] up to the third order in z B and z I (Fig. S5B).
When performing systematic calculations at different lattice sizes, we observe a linear dependence of ∆Φ with respect to the inverse of the system size, L −1 , in the case of nonequilibrium inhomogeneous systems (Fig. S5C). By contrast, we do not observe any system-size dependence for the nonequilibrium homogeneous and equilibrium systems (Fig. S5B,D). Thus, in the case of inhomogeneous systems, we calculate the supersaturation in the thermodynamic limit by extrapolating the values obtained from simulations performed in finite systems to the infinite system size (Fig. S5C). We attribute this system-size dependence to the broken particle-hole symmetry between the liquid and the vapor phases induced by the inhomogeneous chemical reactions. We further note that the steady-state distributions for the nonequilibrium homogeneous systems and their equivalent equilibrium systems are identical (Fig. S5D), which is consistent with the prediction of an effective equilibrium for homogeneous systems discussed in Sec. I D.

B. Validation of phase coexistence via direct coexistence simulations
Phase coexistence is established at a NESS when the net flux of particles, the net flux of thermal energy, and the pressure difference between two phases are zero. These conditions are analogous to the equilibrium phase-coexistence criteria of equal chemical potentials, temperatures, and pressures. In our model, coupling to a single particle reservoir regardless of which phase is currently occupying the lattice, along with the local detailed balance condition governing the reaction rates, ensures that the net particle and heat fluxes vanish between the two phases. In order to satisfy mechanical equilibrium, we propose that the nonequilibrium potential difference defined above, ∆Φ, should be set equal to zero. This choice is motivated by analogy to equilibrium statistical mechanics, in which case ∆Φ is equal to the difference between the grand potential densities of the two phases. Because the grand potential is proportional to the pressure at equilibrium, setting ∆Φ equal to zero guarantees mechanical balance at equilibrium. Under nonequilibrium conditions, however, the relation between the nonequilibrium grand potential and the pressure does not hold. Nonetheless, ∆Φ = 0 still implies that the steady-state probabilities of the two phases occupying a given volume are equal at a NESS. In our lattice model, this condition also means that a long trajectory spends an equal amount of time with each phase completely occupying the lattice. We therefore propose that this condition can be used to determine bulk phase coexistence in the thermodynamic limit at a NESS.
We verify that our definition of phase coexistence based on ∆Φ = 0 coincides with mechanical balance between the phases, and thus is a proper extension of the equilibrium concept, by performing stochastic direct coexistence simulations [41]. We simulate two bulk phases in direct contact at a flat interface and allow the system to evolve at steady state until the lattice is fully occupied by either of the bulk phases. Using a slab geometry on a 16 × 64 lattice, we initialize these simulations with half of the lattice in the liquid phase and the other half in the vapor phase. If the bulk phases are in mechanical balance, then the interface should diffuse in either direction without any bias and thus reach either of the absorbing states with equal probability. We therefore measure the probability for a system to reach the liquid phase from the initial condition, P l , to verify that unbiased diffusion coincides with ∆Φ = 0. Fig. S6 shows that in the homogeneous case, P l = 1 /2 at ∆Φ = 0 as expected. Due to the finite-size effects described above, the results of these simulations in the inhomogeneous case depend on both the longitudinal and transverse dimensions of the lattice. From the difference between β∆Φ(L = 16) and β∆Φ(L = ∞) shown in Fig. S5C, we are able to establish that the P l (β∆Φ) curve lies within the bounds shown in Fig. S6. These results are thus consistent with ∆Φ = 0 corresponding to P l = 1 /2 in the inhomogeneous case. We note that the region of uncertainty in β∆Φ shown in Fig. S6 (approximately ±0.01) is much smaller than the magnitude of β∆Φ in all simulations used to test the applicability of classical nucleation theory (β∆Φ ≥ 0.06). Taken together, these results demonstrate that our definition of nonequilibrium phase coexistence appropriately identifies the conditions for mechanical balance in both homogeneous and inhomogeneous nonequilibrium systems.

A. Forward Flux Sampling (FFS)
We utilize the forward flux sampling (FFS) rare-event simulation method [29] to calculate the nucleation rate starting from the vapor phase. Using the largest cluster of bonding-state particles, n, as the reaction coordinate, we perform FFS using M = 64 milestones from n 0 = 6 to n 63 = 1100, with the spacing between consecutive milestones increasing monotonically from 3 to 50, as we advance the milestones. We first determine the flux across the initial  . S5). We estimate that direct coexistence simulations performed in an infinitely large system would lie within the shaded region in between these curves. milestone, Φ 0 , from a steady-state trajectory in the vapor phase. Likewise, the initial ensemble of configurations at n 0 is obtained by randomly selecting 1000 configurations at n = n 0 from a steady-state trajectory in the vapor phase. We then calculate the probability P (n i+1 |n i ) that a trajectory launched from milestone n i reaches milestone n i+1 before returning to the vapor phase. To this end, we launch trajectories from each milestone n i until we obtain 1000 configurations at n i+1 . We halt the simulation when the probability P reaches unity. Based on these probabilities and the initial flux measurement, the FFS expression for the nucleation rate, I, is (S12)

B. Committing probability and Zeldovich factor
We compute the Zeldovich factor directly from FFS simulations by analyzing the committing probability, φ(n). The quantity φ(n) represents the probability that a system with nucleus size n successfully completes the phase transformation into the stable liquid phase before returning to the metastable vapor phase. We calculate φ(n i ) at each FFS milestone n i based on the milestone probabilities P (n i+1 |n i ), (S13) The critical nucleus size n * is found where φ(n * ) = 1 /2, which is interpolated from the values of φ(n i ). In the diffusive limit, the critical nucleus size coincides with the location of the top of the barrier on the (nonequilibrium) landscape F (n) ≡ −β −1 ln p(n), where p(n) is the steady-state probability of observing a nucleus of size n [30]. In this limit, φ(n) is given by where n = n v and n = n l mark the boundaries of the transition region between the vapor and liquid phases. In the high-barrier limit, the exponential integrands in (S14) are dominated by the barrier height βF (n * ), and we can take saddle point approximation around n = n * , e βF (n) ≈ e βF (n * ) exp βF (n * ) 2 (n − n * ) 2 , (S15) which leads to an approximate form of φ(n): erf[ −βF (n * )/2(n l − n * )] + erf[ −βF (n * )/2(n * − n v )] . (S16) The diffusive-limit condition at φ(n * ) = 1/2 and the liquid-phase boundary condition φ(∞) = 1 further simplify the committing probability into (S17) Note that this approximate form of φ(n) agrees with the vapor phase boundary condition φ(n v ) ≈ 0 as long as n v n * . To compute the Zeldovich factor, we evaluate the second derivative F (n * ) by fitting the committing probabilities at the FFS milestones to Eq. (S17) in the region where 0.25 ≤ φ ≤ 0.75. The Zeldovich factor, is then calculated from the fitted value of F (n * ). Fig. S7A shows that Eq. (S17) works well inside the fitting region.
Furthermore, Fig. S7B shows that the fitted value of the Zeldovich factor and the calculated value determined using the classical nucleation theory (CNT) line tension (see Eq. (S21) below) match one other, so that Eq. (S17) is consistent with CNT.

C. Determination of the nonequilibrium line tension
Ref. [28] extensively tested CNT for the equilibrium two-dimensional square lattice model and validated the CNT expressions for the free-energy landscape along the nucleus-size reaction coordinate, n, and the nucleation rate, I, βF (n) = βσ √ 4πn − nβ∆Φ + 1.25 ln n + d (S19) (S20) In Eq. (S19), d is a constant chosen to equate the gas-phase monomer number density, ρ 1 , and exp[−βF (1)], such that ∆F * ≡ F (n * ) − F (1). Rearranging the expression for the CNT nucleation rate, we obtain an expression for the nucleation barrier height at a given ∆Φ, The critical nucleus size n * follows from the condition F (n * ) = 0, We determine the apparent line tension in nonequilibrium simulations by fitting Eq. (S21) over a range of β∆Φ, using independent measurements of ln(I/ρ 1 D * Γ), and using σ as the sole fitting parameter. We find that we can apply Eq. (S21) to nonequilibrium systems without modifying the functional form of Eq. (S19), although the value of the line tension may differ from the equilibrium value (see Fig. 2B in the main text). In Fig. S8, we demonstrate that the line tension, σ, is the most important variable in determining the nucleation kinetics at a NESS, as expected from equilibrium systems. Within the framework of CNT, there are three independent variables affecting the nucleation kinetics at a fixed supersaturation: ρ 1 , D * , and σ, where the last variable governs both the nucleation barrier and the Zeldovich factor. In Fig. S8, direct comparisons among ρ 1 , D * , Γ, and the apparent barrier height ln(I/ρ 1 D * Γ) at a fixed value of the supersaturation show that the barrier term is indeed the dominant term for both the homogeneous and inhomogeneous nonequilibrium models. Furthermore, in contrast to the barrier term, Γ only shows a weak dependence on ∆µ. This demonstrates that the effect of the nonequilibrium line tension on the nucleation kinetics primarily originates from the nucleation barrier rather than the Zeldovich factor.

IV. THEORETICAL PREDICTIONS FOR THERMODYNAMICS AND KINETICS AT A NESS
In this section, we describe how the approximate theoretical analysis (FLEX) introduced in Sec. I C can be applied to predict the nucleation kinetics at a NESS. We first coarse-grain the two non-interacting isoenergetic states, I and E, into a single inert state in order to map our three-state model onto an effective two-state model. To this end, we sum up the FLEX steady-state distributions for the I and E lattice states in Eq. (S5) and Eq. (S6) assuming a fixed value of the local potential energy, u. We then apply results for the equilibrium lattice-gas model to our effective equilibrium model.

A. FLEX predictions of phase coexistence and entropy production rates
In the equilibrium two-dimensional two-state lattice-gas model, the coexistence condition is given by µ = 2 , where µ is the chemical potential of the particle [22]. This condition reflects the particle-hole symmetry of the two-state lattice-gas model. Within FLEX, we approximate the conditions for phase coexistence by assuming that particle-hole symmetry still applies to the effective equilibrium model. We define the FLEX approximation for the supersaturation, S FLEX , by considering a lattice site surrounded by precisely two bonding state particles, such that the fixed local environment has u = 2 , where z B and z I are defined in Eq. (S5) and Eq. (S6), respectively. Similar to the true equilibrium model, the effective equilibrium approximation given by Eq. (S23) predicts that phase coexistence occurs when the steady-state probability of the bonding and the coarse-grained inert states are equal at a lattice site with u = 2 . FLEX predicts the following expression for the entropy production rate density,Σ FLEX , where j ≡ p B W BI − p I W IB is the net transition flux in the B to I direction. Eq. (S24) indicates thatΣ is always positive unless the system is at equilibrium (β∆µ = 0). Furthermore, the entropy production always depends on u, regardless of the functional form of k I→B , which implies thatΣ should in general differ between the vapor and liquid phases. Within FLEX, we estimate the entropy production rate density in the vapor,Σ v , and in the liquid phase,Σ l , by fixing u = 0 and u = 4 , respectively.  shows that these predictions agree qualitatively with the simulation results. The supersaturation, S, and the steady-state entropy production rate density,Σ, are calculated from NEUS and the simulated trajectory as described in Ref. [21], respectively. The FLEX predictions perfectly match the simulation results for the nonequilibrium homogeneous case, while the inhomogeneous case shows systematic deviations; however, there is a clear linear relation between the predictions and the simulation results even when the system is driven far from equilibrium.

B. FLEX prediction of the nonequilibrium line tension
We use FLEX to predict the line tension contribution to the CNT rate equation, Eq. (S21), by considering the attachment of a single bonding-state particle to a flat liquid-vapor interface at steady state. We perform this calculation at the FLEX-predicted coexistence point, as we expect the line tension to be independent of the supersaturation. We further assume that the interface can be described using the restricted solid-on-solid (RSOS) model [32]. At coexistence, our approximation for the supersaturation, S FLEX , implies that the formation of a kink at the interface (Fig. S10A) incurs no (effective) free-energy penalty in the effective equilibrium model. We can therefore predict the line tension from the (effective) free-energy cost of attaching a single bonding-state adatom to the interface in the RSOS model. We take the negative of this free-energy cost as the new interaction strength β˜ based on the particle-hole symmetry in the effective equilibrium model and predict the line tension using the equilibrium expression [42] where K is the elliptic integral of the first kind, χ(β) = 1 − sinh −4 (β˜ /2) 1/8 , and β c is the (inverse) critical temperature given by β c˜ = 2 ln(1 + √ 2). We find that this RSOS approximation qualitatively explains the decreasing trend of the line tension with respect to the nonequilibrium drive for the inhomogeneous model (see Fig. 3a in the main text). By contrast, we do not observe any change in the predicted line tension in the homogeneous case because the interfacial properties are described by an effective equilibrium model that is common to both phases.
Based on this analysis, we postulate that if the inferred nonequilibrium line tension differs from the equilibrium value, then the coexisting phases at a NESS must be thermodynamically inhomogeneous and thus described by different effective equilibrium models. Both FLEX and our simulation results support this postulated relationship between the thermodynamic inhomogeneity and the interfacial properties. In Fig. S10B, the degree of inhomogeneity, ∆∆f = ∆f l − ∆f v , is calculated from simulation trajectories obtained in each phase at steady state using Eq. (S11), while the FLEX predictions are calculating using Eq. (S7) and assuming that u = 4 for ∆f l and u = 0 for ∆f v . We coarse-grain the inert particle (blue) and unoccupied (empty) lattice sites and thus map the upper configuration to the lower configuration with coarse-grained inert lattice sites (gray). From the steady-state distribution calculated using FLEX, we obtain an effective bonding energy, β˜ , for an isolated adatom as shown. (B) The relationship between the thermodynamic inhomogeneity, ∆∆f , at coexistence and the deviation of the line tension from equilibrium, ∆σ.
(C) The dependence of ∆∆f on the relative timescale, k • , at the same conditions shown in the inset of Fig. 3a in the main text. Solid lines show the FLEX predictions, and marks report simulation results. Data are shown for nonequilibrium homogeneous (blue) and inhomogeneous (orange) models.
Both simulation and theory are consistent with our prediction that a deviation in the line tension (∆σ = 0) implies a nonzero ∆∆f . At the same time, simulation and theory both show that the converse does not necessarily hold, as nonzero values of ∆∆f may not result in nonzero values of ∆σ. FLEX suggests that this latter relationship is dependent on the precise functional form of k I→B . For the inhomogeneous model, we find that the maximum of β∆∆f occurs when k • ≈ 1, which is also when the line tension deviates furthest from the equilibrium value ( Fig. S10C; see also the inset of Fig. 3a). This observation further supports our hypothesis that changes in the interfacial properties are only possible when the two coexisting nonequilibrium phases do not share a common effective equilibrium description. We note that, however, that the two limits k • → 0 and k • → ∞ do not correspond to the same steady-state distribution. In the limit k • → 0, k I→B also approaches 0 regardless of its functional form, and the system reverts back to a true equilibrium with the change of variable λ B → λ B exp(β∆µ), so that z B = z B and z I = z I (Eq. (S5) and Eq. (S6)). In the limit k I→B → ∞, however, the fugacities in the system and the reservoir are not identical unless the system is at equilibrium (β∆µ = 0).

C. FLEX prediction of the nonequilibrium nucleation kinetics
We can derive approximate expressions for the various factors governing the nucleation kinetics, ρ 1 , D * , Γ, and β∆F * , within the FLEX framework. For the monomer density in the vapor phase, ρ 1 , we assume that the bondingstate particles are sparsely distributed and thus approximate the local potential energy u as zero for all lattice sites. Then ρ 1 and the total particle density, ρ v , are approximated as where the steady-state distributions p B , p I , and p E are given by Eq. (S5) and Eq. (S6). We approximate the diffusion coefficient at the top of the nucleation barrier, D * , as the rate of attaching a bondingstate adatom to a circular nucleus of size n * , which is given by Eq. (S22). We assume that the adatom interacts only with the critical nucleus and that the local environment can therefore be described by u = . Under this condition, the mean time to insert a bonding-state adatom into an unoccupied lattice site,T , is given bỹ We approximate the attachment as a first-order transition and define the bonding-state adatom attachment rate, w + , to be the inverse of the mean time, Finally, we approximate the diffusion coefficient as the product of the perimeter of a circular critical nucleus, √ 4πn * , and the bonding-state adatom attachment rate per lattice site, w + : Given the nonequilibrium line tension σ estimated from the RSOS approximation described above, the barrier height ∆F * = F (n * ) − F (1) is calculated from Eq. (S19), and the Zeldovich factor Γ is given by where the critical nucleus size n * is estimated from Eq. (S22). Comparisons between these FLEX predictions and simulation results for both equilibrium and nonequilibrium systems are shown in Fig. S11. Overall, we find qualitative agreement for all four factors in the CNT rate equation. Importantly, FLEX qualitatively predicts the enhanced kinetics for nonequilibrium systems shown in Fig. 3B in the main text: Homogeneous systems show substantial enhancement only for the diffusion coefficient, while the most substantial contribution to the enhanced kinetics in the inhomogeneous case result from the apparent nucleation barrier, as indicated in Fig. S8A.