Pattern formation in a coupled driven diffusive system

TL;DR

This paper introduces a hybrid lattice model to study pattern formation in driven binary mixtures, revealing new phases and behaviors, and providing a continuum description that captures complex phenomena like irregular and perpendicular stripes.

Contribution

The study develops a Field-based Lattice Model combining lattice gas and field theory, uncovering new intermediate phases and deriving coupled PDEs to describe pattern formation in driven systems.

Findings

Identification of an intermediate 'irregular stripes' phase with long-range order

Derivation of coupled PDEs capturing microemulsion and stripe phases

Discovery of novel behaviors like parallel stripes and chaos

Abstract

We investigate pattern formation in a driven mixture of two repulsive particles by introducing a Field-based Lattice Model (FLM), a hybrid model that combines aspects of the driven Widom-Rowlison lattice gas (DWRLG) and its statistical field theory. We find that the FLM effectively captures the bulk behavior of the DWRLG in both low- and high-density phases, suggesting that phase transitions in these models may share a universality class. Under the effect of the drive, the FLM additionally reveals an intermediate regime, not reported in the previous DWRLG studies, characterized by "irregular stripes" with widely fluctuating widths, contrasting with the "regular", well-ordered stripes found at higher densities. In this intermediate phase, the system exhibits long-range order, predominantly perpendicular to the drive direction. To construct a continuum description, we derive two coupled partial differential equations via a gradient expansion of the FLM mean mass-transfer equations, supplemented with additive noise. Designing a numerical solver using the pseudospectral method with dealiasing and stochastic time differencing, we reproduce the low-density microemulsion phase (characterized by a non-zero characteristic wavenumber q*) and perpendicular stripes at high density. We identify the non-zero difference in the characteristic velocities of the fields as a necessary condition for perpendicular stripe formation in the high-density phase. The continuum model also uncovers novel behaviors not previously observed in the FLM, such as stripes aligned parallel to the drive and chaotic patterns. This work highlights how the interplay of external drive, particle interactions, and noise can lead to a rich phenomenology in strongly driven binary mixtures.

Figures

Tables

Table 1. Table 1: Comparison of the charge fluctuation structure factor values (using L = 128 L=128 ), S − ​ ( q ∗ , q y ) S_{-}\left(q^{\ast},q_{y}\right) , at the peak position occurring at q y = 0 q_{y}=0 (third column) and at locations slightly displaced from the peak, q y ​ L / 2 ​ π = 1 , 2 q_{y}\,L/2\pi=1,2 (fourth and fifth columns, respectively) for different values of ρ \rho corresponding to the three regimes: microemulsion (first row), irregular stripes (second row), and regular stripes (third row).

$\rho$	$q^{\ast }$	$S_{-}(q^{\ast },0)$	$S_{-}(q^{\ast },0.05)$	$S_{-}(q^{\ast },0.1)$
$0.39$	$0.1$	$3.55$	$2.47$	$1.27$
$0.60$	$0.1$	$486.55$	$12.90$	$2.26$
$0.80$	$0.54$	$2083.81$	$9.68$	$2.93$

Table 2. Table 2: Comparison of the peak values of the charge fluctuation structure factor (third column) occurring at q y = 0 q_{y}=0 with peak values at q y ​ L / 2 ​ π = 1 , 2 q_{y}\,L/2\pi=1,2 in the irregular ( ρ = 0.4 , δ = 1 \rho=0.4,\delta=1 ) and the regular ( ρ = 0.6 , δ = 1 \rho=0.6,\delta=1 ) stripe regimes. We use g 0 = 0 g_{0}=0 , g − = 0.5 g_{-}=0.5 and λ = 1 \lambda=1 .

$\rho$	$q^{\ast }$	$S_{-}(q^{\ast },0)$	$S_{-}(q^{\ast },0.05)$	$S_{-}(q^{\ast },0.1)$
$0.4$	$0.01$	$181.04$	$3.47$	$0.56$
$0.6$	$0.25$	$826.19$	$0.31$	$0.072$

Equations

$$\sum _{\mathbf{s}}{\rho }_A(\mathbf{s})=\sum _{\mathbf{s}}{\rho }_B(\mathbf{s})=\frac{\rho L^2}{2}.$$

$$w_{\mathbf{s}\to {\mathbf{s}}^{\prime }}=\frac{1}{8}\times \{\begin{array}{ll}1+\delta & \mathrm{if}\Delta i=1\\ 1& \mathrm{if}\Delta i=0\\ 1-\delta & \mathrm{if}\Delta i=-1\end{array}$$

$$\Delta {\rho }_A(\mathbf{s})=ϵ{\rho }_A(\mathbf{s})h({\mathbf{s}}^{\prime })/{\rho }_B({\mathbf{s}}^{\prime }),$$

$$/{\rho }_B({\mathbf{s}}^{\prime })\equiv \prod _{{\mathbf{s}}^{\prime \prime }}[1-{\rho }_B({\mathbf{s}}^{\prime \prime })],$$

$${\varphi }_{\pm }(\mathbf{s},t)=\left({\rho }_A(\mathbf{s},t)-\frac{\rho }{2}\right)\pm \left({\rho }_B(\mathbf{s},t)-\frac{\rho }{2}\right),$$

$${\tilde{\varphi }}_{\pm }(\mathbf{q},t)=\sum _{\mathbf{s}}{\varphi }_{\pm }(\mathbf{s},t)e^{i\mathbf{q}\cdot \mathbf{s}}$$

$$S_{\pm }(\mathbf{q},t)=\frac{1}{L^2}⟨{\left|{\tilde{\varphi }}_{\pm }(\mathbf{q},t)\right|}^2⟩.$$

$$\Phi =\frac{S_{-}^{\ast }}{L^2}.$$

$${\Sigma }_{\pm }\equiv L^2\sum _{\mathbf{q}}{S}_{\pm }(\mathbf{q}).$$

$${\rho }^2\le {\Sigma }_{+}\le \rho \text{mbox}and0\le {\Sigma }_{-}\le \rho .$$

$${\bar{\rho }}_{A,B}\left(i\right)\equiv \frac{1}{L}\sum _j{\rho }_{A,B}\left(i,j\right),$$

$${\bar{\varphi }}_{+}\left(i\right)\equiv {\bar{\rho }}_A\left(i\right)+{\bar{\rho }}_B\left(i\right)-\rho ,$$

$${\bar{\varphi }}_{-}\left(i\right)\equiv {\bar{\rho }}_A\left(i\right)-{\bar{\rho }}_B\left(i\right).$$

$${\partial }_t{\varphi }_{\pm }(\mathbf{r},t)+\nabla \cdot {\mathbf{J}}_{\pm }+{\partial }_x{J}_{\pm }^{\delta }=0,$$

$${\partial }_t{\varphi }_{+}$$

$$-\lambda {\nabla }^2{\varphi }_{-}^2-g_{-}{\partial }_x{\varphi }_{-}^2+g_{+}{\partial }_x{\varphi }_{+}^2+{\xi }_{+},$$

$${\partial }_t{\varphi }_{-}$$

$$+\lambda \nabla ({\varphi }_{-}\nabla {\varphi }_{+})+g_0{\partial }_x({\varphi }_{+}{\varphi }_{-})+{\xi }_{-},$$

$$⟨{\xi }_{\pm }(\mathbf{r},t){\xi }_{\pm }({\mathbf{r}}^{\prime },t^{\prime })⟩=-{\sigma }^2{\nabla }^2\delta (\mathbf{r}-{\mathbf{r}}^{\prime })\delta (t-t^{\prime }),$$

$$\{\begin{array}{l}{\varphi }_{-}(\mathbf{r},0)={\sigma }_0\xi (\mathbf{r})\\ {\varphi }_{+}(\mathbf{r},0)=0\end{array}$$

$${\partial }_x\left[\left(g_0{\varphi }_{+}-\Delta v\right){\varphi }_{-}\right].$$

$${\Delta }_t{\rho }_A\left(\mathbf{s},t\right)=-\{{\Delta }_{\mathrm{NN}}{J}_A^{\mathrm{NN}}\left(\mathbf{s},t\right)+\delta {\Delta }_{\mathrm{NN}}{J}_A^{\delta ,\mathrm{NN}}\left(\mathbf{s},t\right)+{\Delta }_{\mathrm{NNN}}{J}_A^{\mathrm{NNN}}\left(\mathbf{s},t\right)+\delta {\Delta }_{\mathrm{NNN}}{J}_A^{\delta ,\mathrm{NNN}}\left(\mathbf{s},t\right)\}$$

$${\Delta }_t{\rho }_A(\mathbf{s},t)\equiv {\rho }_A(\mathbf{s},t+\tau )-{\rho }_A(\mathbf{s},t)$$

$$-{\Delta }_{\mathrm{NN}}{J}_A^{\mathrm{NN}}\left(i,j\right)$$

$$+\frac{1}{8}\sum _{\kappa =\pm 1}\{{\rho }_A(i,j-\kappa )h(i,j)/{\rho }_B^{\updownarrow }(i,j;\kappa )-{\rho }_A(i,j)h(i,j+\kappa )/{\rho }_B^{\updownarrow }(i,j+\kappa ;\kappa )\},$$

$$-{\Delta }_{\mathrm{NN}}{J}_A^{\delta ,\mathrm{NN}}(i,j)$$

$$/{\rho }_B^{\leftrightarrow }(i,j;\kappa )=(1$$

$$(1-{\rho }_B(i,j+\kappa ))(1-{\rho }_B(i,j-\kappa )),$$

$$/{\rho }_B^{\updownarrow }(i,j;\kappa )=(1$$

$$(1-{\rho }_B(i+\kappa ,j))(1-{\rho }_B(i-\kappa ,j))$$

$${\Delta }_t{\varphi }_{\pm }=-\left\{\Delta J_{\pm }+\delta \Delta J_{\pm }^{\delta }\right\}$$

$${\varphi }_{\pm }$$

$$={\varphi }_{\pm }\left(x,y\right)+\ell {\partial }_x{\varphi }_{\pm }\left(x,y\right)+\frac{{\ell }^2}{2!}{\partial }_x^2{\varphi }_{\pm }\left(x,y\right)+...$$

$${\Delta }_t{\varphi }_{\pm }=\tau {\partial }_t{\varphi }_{\pm }+...$$

$${\partial }_t{\varphi }_{+}$$

Full text

I Introduction

While statistical mechanics has excelled in describing systems at equilibrium, a significant portion of the world operates far from equilibrium. Among the most fascinating properties of nonequilibrium systems are the complex spatiotemporal patterns they exhibit [3, 4, 5], including those associated with phase transitions and critical phenomena [6]. There has been growing interest in understanding how such patterns emerge and how they can be used to infer important properties of physical systems, whether active or driven externally.

In the study of systems driven out of equilibrium, significant open questions remain, including how to define phase coexistence [7] or a thermodynamic quantity of the same usefulness as entropy at equilibrium [8]. A driven system is subject to an external force (or “drive”, for simplicity) that maintains it out of equilibrium. Examples include driven diffusive systems in which the drive biases particle diffusion along an axis with periodic boundaries, as in the pioneering studies of driven Ising, Potts, and Blume-Emery-Griffths models (in their lattice-gas representations) [9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19]. Driven diffusive systems exhibit a range of intriguing behaviors [20, 21], including generic long-range correlations at all temperatures [9, 10] and critical properties governed by a nonequilibrium fixed point [22, 23]. At low temperatures, these models typically exhibit phase-separation into stripes which form along the drive direction [9, 10]. Recently, some of the present authors have considered a driven version [1, 2] of the Widom-Rowlinson lattice gas [24] (WRLG). The WRLG is a lattice gas model in which the two species of particles ($A$ and $B$, say) diffuse freely, apart from the constraint of strict exclusion of nearest neighbor (NN) $A$-$B$ pairs. With no interactions except volume exclusion, the only control parameters of the equilibrium WRLG are the overall densities of the two species. Governed by entropy alone, this equilibrium system transitions, as the densities increase, from a homogeneous disordered state to one with phase separation, much like the Ising lattice gas [24]. When both species are driven in the same direction, striking features emerge [1, 2]: In the disordered phase, the $A$-$B$ correlation function develops an intrinsic, preferred wavenumber ($\mathbf{q}^{*}\propto\hat{\mathbf{x}}$, the drive direction)! Correlations grow as the overall density increases, until long-range order appears – in the form of alternating $A$- and $B$-rich stripes, characterized by $\mathbf{q}^{*}$. In stark contrast with the driven Ising case [10], these stripes are always perpendicular to the drive. Further, since $\mathbf{q}^{*}$ is independent of $L$, the system is not scale invariant.

Notably, similar nonequilibrium phenomena, such as pattern formation and phase separation, are found in active matter, where scalar active mixtures have highlighted the importance of nonreciprocal couplings and characteristic length scales far from equilibrium [25, 26, 27, 28, 29]. Although the active matter literature is filled with interesting examples of pattern formation, many such phenomena already appear in simpler externally driven systems. The driven Widom-Rowlinson lattice gas (DWRLG) provides an especially interesting and simple setting for exploring patterns. The model was originally introduced and studied in detail in Ref. [1], where its phase behavior and stripe formation were characterized through simulations. More recently, a coarse-grained field-theoretic description of its disordered phase was formulated using the Doi-Peliti formalism [30, 31] in Ref. [2]. In the latter work, the field equations describe the nonequilibrium dynamics as two coupled driven diffusive systems: the sum and difference of $A$ and $B$ particle densities (referred to as the “density” and “charge” fields in the following), with different characteristic velocities. The field theory is capable of capturing key features observed in the simulations at low particle densities, including the development of a characteristic wavenumber $\mathbf{q}^{*}$, and further highlights the importance of noise in generating this characteristic wavenumber in the disordered phase.

The present study builds upon Refs. [1, 2], bridging the gap between them. More specifically, we seek a minimal continuum description, inspired by the Doi-Peliti formalism, of the ordered (heterogeneous) phase, where open questions remain regarding the structure of the equations of motion and whether they are capable of capturing the stripe orientation observed in lattice simulations. We also revisit the properties of the low-density phase. To this end, we first introduce a new Field-based Lattice Model (FLM) that preserves all the symmetries of the original DWRLG. We find excellent agreement with the phenomenology reported in Refs. [1, 2], strongly suggesting that both models belong to the same universality class. While discrete field theories are commonly employed in the quantum domain [32], their stochastic counterparts are far less explored; our lattice model provides a simple and straightforward example of a stochastic discrete field theory. Next, we perform a gradient expansion (GE) of the FLM master equation, leading to a pair of coupled field equations that, we shall argue, satisfactorily describe both low- and high-density regimes.

The equations we derive are consistent with the Doi-Peliti and Martin-Siggia-Rose-Janssen-De Dominicis (MSRJD) field equations developed for the DWRLG in Ref. [2], but include additional terms necessary for describing a stable ordered phase. In the absence of an external drive, our equations of motion follow the structure of a two-component Cahn-Hilliard-like equation [33]. But a drive breaks rotational symmetry, leading to stripe formation, or traveling waves, similar to phenomena observed in scalar active matter [25, 34]. We demonstrate here that aspects of the phenomenology reported in [25, 34] are already evident in a driven matter context, without requiring non-reciprocal particle interactions or particle activity.

The remainder of the paper is organized as follows: In the next section, we introduce the FLM, providing a brief characterization of the model. In Sec. III, we derive the field equations via a GE, including explicit expressions for the linear coefficients, with additional details outlined in Appendices A and B. In Sec. IV, we integrate these equations numerically using the pseudospectral method combined with stochastic second-order Heun exponential time-differencing (PSSETD2H) [35, 36] and smooth dealiasing techniques [37], with further numerical details provided in Appendices C and D. Both deterministic and stochastic versions of the field equations are analyzed. One of the key numerical findings is that the difference in the characteristic velocity of the charge and density fields plays a crucial role in the nonequilibrium behavior of the model, supporting the conclusions of [2], but in the ordered phase as well. Finally, Sec. V presents our conclusions and outlines potential future directions.

II Field-based Lattice Model

Early studies of Widom-Rowlinson systems were based on two distinct microscopic models, the original one of interacting particles in the continuum [24], and a lattice-gas model [38]. In more recent studies [1, 2], a stochastic field theoretic version was also considered, with the bridge between the two established either phenomenologically or using the Doi-Peliti formalism. All exhibit two distinct regimes: a disordered low-density phase and a high-density ordered phase. Here, we consider a hybrid model 111The rudiments of this model were introduced in [1], where it played the role of an intermediary between simulations of the lattice gas and studies of a phenomenological field theory. In the present study, we investigate a full stochastic version of it., with continuous variables (densities, $\rho_{A,B}\in\left[0,1\right]$) defined on a $L\times L$ lattice 222Generalizations to lattices in higher dimensions are straightforward.. This approach follows the spirit of coarse graining a lattice gas and will be referred to as the Field-based Lattice Model (FLM). Expecting it to be qualitatively similar to the WRLG, we perform simulations to explore its behavior. As we will show, the driven FLM has a richer phenomenology than the DWRLG, in that a third regime appears as the overall density is varied between the homogeneous and striped phases. In this regime, which we call the irregular stripe regime in the following, there is long-range order only in the direction perpendicular to the drive, with well-defined stripes, but of markedly different widths. In the next sections, we present a straightforward gradient expansion of the FLM to arrive at a field-theoretic description with equations of motion for the densities, and compare our results with the field theories studied earlier [1, 2]. We begin with a review of the driven WRLG.

II.1 Model specification

The driven WRLG [1] consists of two species of particles, $A$ and $B$, diffusing on a square lattice of $L^{2}$ sites with periodic boundary conditions. Each site $\mathbf{s}=\left(i,j\right)$ is either empty or occupied by a single particle. In addition, opposite species repel each other: nearest neighbor (NN) $A$-$B$ pairs are forbidden. In the stochastic dynamics, a random particle is chosen and moved to a randomly chosen NN or next NN (NNN) site, provided the destination site is empty and its NN sites are not occupied by a particle of the opposite species. When an external drive $\delta$ is imposed, particle moves in the $+x$ (increasing $i$) direction are favored over moves in the $-x$ (decreasing $i$) direction. Following previous simulation studies, we set the numbers of $A$ and $B$ particles equal, so that we have only two control parameters: the total particle density, $\rho$, and the drive, $\delta$. For $\delta=0$ (the WRLG), the dynamics satisfies detailed balance and the properties of its stationary state can be predicted by considering entropy alone, as every allowed configuration is equally probable. Phase separation occurs when $\rho$ is increased beyond a critical density, $\rho_{c}^{\left(0\right)}\simeq 0.618(1)$, with critical behavior in the Ising universality class [24]. For $\delta>0$ (the DWRLG), detailed balance is violated, the stationary distribution is unknown a priori, while novel phenomena emerge in the nonequilibrium steady-state, for all densities (below, near, and above a $\delta$-dependent $\rho_{c}^{\left(\delta\right)}$, with $\rho_{c}^{\left(\delta\right)}>\rho_{c}^{\left(0\right)}$, in general) [1, 2].

The FLM differs from the DWRLG mainly in having continuous (instead of discrete) variables at each site of the lattice, i.e., $\rho_{A,B}\in\left[0,1\right]$ instead of $0,\pm 1$. Both densities are conserved; here we consider only the case of equal overall densities:

[TABLE]

The evolution is defined through the following set of rules, which are motivated, in part, by the excellent agreement between simulation results for the current-density relation, $J\left(\rho\right)$, and a simple phenomenological form, discussed in more detail in Ref. [1]:

At each iteration, select a site $\mathbf{s}=(i,j)$ at random and one of the eight NN or NNN sites $\mathbf{s}^{\prime}=(i^{\prime},j^{\prime})$, according to the following drive-dependent probabilities

[TABLE]

where $\Delta i=i^{\prime}-i$ and $\delta\in\left[0,1\right]$ is the drive. In other words, hops along (opposite) the direction of the drive are enhanced (suppressed). 2. 2.

Next, select one of the species at random. In case $A$ is chosen, then transfer the following amount of $A$-density from $\mathbf{s}$ to $\mathbf{s}^{\prime}$:

[TABLE]

where $\epsilon$ is a random variable 333In our simulations, $\epsilon$ is chosen from a uniform distribution. From the results, $\epsilon$ serves only to redefine a diffusion length. Note that this multiplicative aspect of the noise may become relevant near the boundary values: $\rho_{A,B}=0,1$. in $\left[0,1\right]$, $\rho_{A}(\mathbf{s})$ is the density of $A$-particles at site $\mathbf{s}$, $h(\mathbf{s}^{\prime})=1-\rho_{A}(\mathbf{s}^{\prime})-\rho_{B}(\mathbf{s}^{\prime})$ is the vacancy density at $\mathbf{s}^{\prime}$, and $\not{\rho}_{B}(\mathbf{s}^{\prime})$ is defined as

[TABLE]

where the product is over all nearest neighbors $\mathbf{s}^{\prime\prime}$ of $\mathbf{s}^{\prime}$, excluding the original site $\mathbf{s}$. Note that this product contains three factors for NN hops, and four in the case of hopping to a NNN. The factor $\not{\rho}_{B}$ models $A$-$B$ repulsion through an estimate of the probability that the nearest neighbors of the target site are free of $B$’s 444The form of Eq. (3) and these suppression factors are motivated in part by the remarkably good agreement between a conjectured formula for the current density and the observed one in DWRLG: Eq. (8) and Fig. 5 in [1].. If instead of $A$, species $B$ is chosen, the same rules apply, simply exchanging $A$ and $B$ in all expressions.

These steps are iterated to generate the evolution. We emphasize that $\Delta\rho_{A,B}(\mathbf{s})\leq h(\mathbf{s}^{\prime})$, so that it is impossible for a site to be “overfilled”. As usual, a Monte Carlo Step (MCS) is defined as $2L^{2}$ iterations, so that, on average, each species at each site is given one attempt to transfer some of its contents.

II.2 Results from simulations and comparisons with the lattice gas

To explore the behavior of this model, we carry out simulations on lattices of size $L$ ranging from $32$ to $256$, varying the overall density, $\rho$, from $0.10$ to $0.95$. Here, we report only results from the extreme cases of $\delta=0$ and $1$, which exhibit all the key phenomenology of the model at intermediate drive values. The main set of studies employs a uniform initial state with $\rho_{A,B}\left(\mathbf{s}\right)=\rho/2$ and, for measurements in the steady-state, up to $100$ realizations per run, for up to $2\cdot 10^{7}$ MCS. With no drive ($\delta=0$), the behavior is essentially the same as the WRLG, displaying a disordered (phase-separated) state at low (high) densities, with the critical density estimated to be $\rho_{c}^{\left(0\right)}\simeq 0.39$. For $\rho$’s below this value, the structure factors are isotropic and display the standard Ornstein-Zernike form. Near criticality, we expect the properties to be in the Ising class [24] and did not pursue this issue. At densities above $\rho_{c}^{\left(0\right)}$, the usual coarsening dynamics is observed, and we verify [43] that our system follows the Lifshitz-Slyozov law, i.e., characteristic length scales associated with the peaks of structure factors grow initially as $t^{1/3}$. For densities closer to unity, the system exhibits arrested coarsening, with the lack of vacancies and $A$-$B$ repulsion leading to considerable frustration. The system begins to lose ergodicity 555Indeed, for $\rho$ as low as just $0.5$, certain configurations are “frozen”. Examples include the checkerboard pattern (all white squares vacant, alternate black ones filled with $A$ or $B$) and similarly alternate rows or columns. Note that such configurations are also present in the WRLG, as $\rho=0,\pm 1$ appears simply as $n=0,\pm 1$. As a result, the particular density $\rho_{g}$ at which this transition occurs may depend on the initial condition. and exhibits glassy behavior. Systems starting with $\rho_{A,B}\left(\mathbf{s}\right)=\rho/2$ quickly freeze into labyrinthine patterns that resemble those of the Turing-Swift-Hohenberg model [45]. A rough estimate of the “glass transition” density is $\rho_{g}\sim 0.9$. Typical configurations from these two (expected) phases and the final glassy regime of the $\delta=0$ system are schematically depicted in Fig. 1(b).

The behavior of the driven FLM (at $\delta=1$ here) is more surprising. Though many of its properties are much the same as those in the DWRLG (i.e., disordered at low densities and ordering into regular stripes at high densities), there are some major differences. In particular, it may be argued that there is an extra regime (as $\rho$ is varied) in the FLM, consisting of a phase in between the analogs of the ordered (high-density) and disordered (low-density) phases of the DWRLG.

At low densities, both the DWRLG and the FLM display a homogeneous, “microemulsion” phase, in which the correlation functions show a periodic structure in the drive direction. For larger $\rho$, both display inhomogeneities. While the DWRLG develops well defined $A$- and $B$-rich stripes when $\rho>\rho_{c}^{\left(0\right)}\sim 0.6$, stripes begin to emerge in the FLM at a lower value, $\rho_{\ell}\sim 0.4$. In addition, the widths of the stripes in the DWRLG are relatively small and (mostly) independent of $\rho$ for maximal drive ($\delta=1$) [1]. By contrast, the stripes in the FLM that first emerge for $\rho>\rho_{\ell}$ appear to be “irregular” in that their widths fluctuate widely. This regime is apparently not accessible to the DWRLG. Such behavior persists up to an “upper” density, $\rho_{u}\sim 0.65$.

As we increase $\rho$ beyond $\rho_{u}$ and all the way up to the glass transition density $\rho_{g}$, we find properties that more closely resemble those of the DWRLG high-density phase: starting from a uniform initial state, the system quickly develops many stripes with thin, regular widths with a well defined characteristic wavenumber $q^{*}$. We will label this the regular stripe regime, $\rho\in\left(\rho_{u},\rho_{g}\right)$. Although investigations with some other initial conditions suggest that there are many long-lived metastable states (i.e., not relaxing to regular stripes within the time limits of our runs), a definitive study of them is challenging and beyond the scope of this article [43]. Instead, our main focus here will be the regular and irregular stripe regimes of the FLM. Although the various regimes encountered in our simulations and numerical integration appear to represent distinct phases, we caution the reader that the term is used somewhat informally. Deferring more systematic studies (including finite-size scaling analyses) required to establish the existence of well defined phases (and associated phase transitions) to future work, our goal is to present qualitative observations in the next subsection.

At the highest densities, excluded-volume considerations mark the behavior of the DWRLG and the driven FLM. In the DWRLG, exclusion of $A$-$B$ NN pairs prohibits the use of random initial configurations. In [1], such densities were studied using single-stripe initial configurations consisting of two blocks (A particles plus vacancies, B particles plus vacancies) separated by empty columns. For the FLM with $\rho>\rho_{g}$, we encounter extremely slow dynamics reminiscent of glass, just as in the $\delta=0$ case. Comparing panels (c) and (d) of Fig. 1, we note that the precipitous drop in $\Phi$, marking the onset of the glassy regime, occurs at nearly the same density in the undriven and maximally-driven cases, suggesting that this kinetic phenomenon is independent of the drive, and occurs well before ordering can develop. To summarize, these different behaviors and regimes for the driven FLM are depicted in Fig. 1(a).

II.3 Brief analysis of various regimes

To highlight the novelty of the irregular stripe regime in the driven FLM, we will present a number of approaches to analyzing the $\delta=1$ system. To best identify the various regimes accessible to this model, we focus here on systems with $L=32,64,128,256$ and three densities, $\rho=0.39,0.60,0.80$, corresponding to the homogeneous, irregular stripe and regular stripe regimes, respectively. As mentioned above, the typical configurations in the first and last regimes are similar to the “microemulsions” and the striped phases of the DWRLG [2]. By contrast, starting with uniform $\rho_{A,B}$ with overall $\rho\in\left(\rho_{\ell},\rho_{u}\right)$, the system quickly (e.g., $\sim 10^{5}$ MCS for $\rho=0.60$ and $L=256$) forms a relatively large number of thin stripes (of regular width, mainly), as if it were settling into regular stripes. Soon thereafter, these stripes suffer a long-wavelength instability, reminiscent of the Eckhaus instability of convection rolls [47], and break up. They subsequently recombine and, by the end of the runs ($\sim 10^{7}$ MCS), settle into relatively few well-defined (irregular) stripes with a wide range of widths. We illustrate this unusual behavior in Fig. 2 with three snapshots of a $\rho=0.60,\ L=256$ system at about $10^{5},10^{6}$, and $10^{7}$ MCS.

There are a number of ways to characterize quantitatively the atypical behavior in the irregular stripe regime. We focus only on the steady-state, deferring an analysis of the interesting evolution process to the future. A standard quantity is the (two-point) correlation function of the fluctuating densities. Instead of deviations of $\rho_{A,B}$ from $\rho/2$, let us consider their sum and difference:

[TABLE]

labeled as the density (DF) and charge (CF) fluctuations: $\phi_{+}$ and $\phi_{-}$, respectively. Translational invariance of the steady-state allows us to exploit their Fourier transforms

[TABLE]

and study the transform of the correlations, i.e., the structure factors (SF) 666Due to the symmetries of the systems we are considering, it can be shown that the cross correlation (associated with $S_{+-}$) vanishes, so that we investigated on these two SFs.

[TABLE]

The average $\left\langle\cdots\right\rangle$ is over the stationary-state ensemble. In practice, we perform both an ensemble average (with up to $100$ independent runs) and time average (of each run, measuring $1000$ times in the last $10^{7}$ MCS).

As in the DWRLG, $S_{-}(\mathbf{q})$ is found to peak at a nonzero wave-vector: $\mathbf{q}^{\ast}\equiv\left(q^{\ast},0\right)$, even in the homogeneous phase, while $S_{+}(\mathbf{q})$ peaks at $2\mathbf{q}^{\ast}$ at high densities (when stripes become more well-established). In fact, the doubling of the frequency for the DF can be explained in simple terms: since the density field marks the interfaces between A- and B-rich regions, a modulation for the CF must be accompanied by a modulation of twice the frequency in the density. Let us denote the value of $S_{-}$ at the peak as $S_{-}^{\ast}\equiv S_{-}(\mathbf{q}^{\ast})$. As in ordinary equilibrium systems (in $d$ dimensions) with transitions to long-range order, we expect this maximum to be $O\left(1\right)$ in the disordered phase, but rise to $O\left(L^{d}\right)$ in the ordered one. Thus, we follow Ref. [1] and define the order parameter as

[TABLE]

Associated with the SFs are sum rules:

[TABLE]

These would be constants in the DWRLG, but, for the FLM, only bounds can be placed on them:

[TABLE]

The two extremes come from a state of complete disorder [$\rho_{A,B}\left(\mathbf{s}\right)=\rho/2$, such as those used as initial states in simulations] and one of complete order [$\rho_{A,B}\left(\mathbf{s}\right)=0$ or $1$ only]. As the crudest measure of order, it is interesting to report that, e.g., $\Sigma_{-}$ is approximately $0.75\%$ of being maximal for $\rho=0.39$, rising to approximately $15\%$ for the irregular stripe regime case of $\rho=0.60$, while reaching about $21\%$ in the high-density $\rho=0.80$ case with regular stripes.

We also demonstrate the different regimes by plotting the $\rho$-dependence of the order parameter, $\Phi\left(\rho\right)$, over a wide range of $\rho$: $\left(0.1,0.9\right)$. In the lower panels of Fig. 1, we illustrate for both $\delta=0$ and $\delta=1$ [Figs. 1(c) and (d), respectively], with $L=64$. In the former, there are clearly two (or three, if the the glassy regime at the highest density is included) regions with very different behavior: $\rho$ below or above $\rho\sim 0.4$. From the associated typical configurations with $L=256$ in (i) and (ii), we identify these regions with the disordered and ordered phases of the $\delta=0$ system, leading to the estimate $\rho_{c}^{\left(0\right)}\simeq 0.39$. By contrast, in Fig. 1(d), there are arguably three (four, including the glassy) regimes for the $\delta=1$ case, from which we estimate $\rho_{\ell}\simeq 0.40$ and $\rho_{u}\simeq 0.65$. We note that the order parameter in the ordered regimes of the driven ($\delta=1$) system is over an order of magnitude larger than in the corresponding regimes of the undriven case. Fig. 3(a) further supports these observations by displaying the peak position of the charge fluctuation structure factor, $q^{*}$, as a function of $\rho$ for $\delta=1$ and $L=128$. A clear transition is observed around $\rho\sim 0.65$, where $q^{*}$ exhibits a sharp jump, indicating a qualitative change in behavior within the ordered phase. Figs. 3(b) and (c) show representative steady-state configurations for four selected densities: $\rho=0.475$ (i), $0.625$ (ii), $0.700$ (iii), and $0.800$ (iv), demonstrating that long-range order has developed along the drive direction.

Another indication of the presence of a new regime between the disordered and regular-stripe regimes is that $\Phi\left(\rho,L\right)$ scales differently with $L$ in three regions. As noted above, $\Phi\sim O\left(L^{-d}\right)$ in the disordered phase, rising to $O\left(1\right)$ in the ordered one for standard statistical systems. For example, for a non-conserved Ising model, the magnetization is nonzero only in the ordered phase (in the thermodynamic limit). In the driven FLM, it can be argued that we encounter anomalous scaling in the intermediate, irregular stripe regime. To illustrate this phenomenon, we display plots of $\Phi\left(\rho,L\right)$ for $\delta=0$ [panels (a)] and $\delta=1$ [panels (b)], scaled by various powers of $L$, in Fig. 4. Note that these are semi-log plots, so that data of a wide range of magnitudes can be shown together. We can discern different scaling properties in the three regimes even for a limited range of $L$’s ($32,64,128$). In the lower panels, $\Phi L^{2}$ (i.e., $S_{-}^{\ast}$ itself) is shown and, noting the data collapse in the low density regime, $\rho<\rho_{\ell}$, for $\delta=1$ (or $\rho<\rho_{c}^{(0)}$ for $\delta=0$), we may conclude that $\Phi$ has the typical scaling behavior of a disordered phase. In the upper panels, we show that the $\Phi$ values collapse best in the high density regime $\rho>\rho_{u}$ for $\delta=1$ (or $\rho>\rho_{c}^{(0)}$ for $\delta=0$), though considerable finite-size corrections may be at play. By contrast, for $\delta=1$, the middle-right panel, corresponding to intermediate densities $\rho_{\ell}<\rho<\rho_{u}$, displays good data collapse for $\Phi L^{p}$, where $p$ is estimated as $p\sim 0.7$ for the best data collapse. Note that no such collapse occurs for $\delta=0$ (or for any other value of $p$ between [math] and $2$), in the middle-left panel, highlighting the absence of a third regime in the equilibrium FLM. The conclusion is that, for the driven FLM, there is a regime where $\Phi$ scales with a power of $L$ that is “intermediate” between the standard values of $-2$ and [math], corresponding to the usual disordered and ordered phases, respectively. This is our irregular stripe regime.

A further argument for the irregular stripe regime comes from more detailed data involving sections of $S_{-}\left(q_{x},q_{y}\right)$. In Fig. 5(a-c), we present an illustration ($\delta=1,L=128$) of how $S_{-}$ varies with $q_{x}$ for the first decade of $q_{y}\,L/2\pi=0,...,10$ in the three regimes. Deep in the homogeneous (microemulsion) regimes, $S_{-}$ roughly follows the Orstein-Zernike form [with two major differences: a discontinuity singularity at $\mathbf{q}\equiv(q_{x},q_{y})=(0,0)$ and a peak at a nonzero, $L$-independent, $\mathbf{q}^{\ast}=\left(q^{\ast},0\right)$]. In other words, the values of $S_{-}\left(q_{x},q_{y}\right)$ are comparable to $S^{\ast}$ (i.e., the peak is quite broad), as seen in Fig. 5(a). In the opposite extreme, deep in the ordered regime, the peak is sharp and narrow (with a value increasing with the system size), as seen in Fig. 5(c). In the irregular stripe regime [Fig. 5(b)], the peak value at $\mathbf{q}^{*}=(q^{*},0)$ falls between the two phases and also falls off more slowly as we increase $q_{y}$.cTo quantitatively appreciate how rapidly $S_{-}$ falls from its maximum values $S_{-}^{\ast}$ as we move away from the ridge, we display the numerical values in Table 1.

The peak positions of the CF SF, which occur at $\left(q^{\ast},0\right)$, are listed in the second column, while the next three columns are associated with the $S_{-}$ values at $\left(q^{\ast},q_{y}\right)$, with $q_{y}\,L/2\pi=0,1$ and $2$. Note that the ratios of the entries in the third and fourth columns are dramatically different in the three cases: $\sim 1.4$, $38$, and $215$ for the microemulsion phase, the irregular stripes, and the regular stripes, respectively. Remarkably, once $q_{y}$ reaches $3\cdot 2\pi/L$, there is hardly any difference between the three rows! We also call attention to the $q^{\ast}$ of the irregular stripe regime, which is the same (for this $L$) as that of the disordered phase, but considerably smaller than that in the regular stripe regime, which may be also observed in Fig. 3 where the stripes at high density are visibly narrower. In Fig. 6, we plot the structure factors for $L=256$. The peak in the low density, disordered phase is well fit by a function of the form $A_{1}q^{2}\left[A_{2}q^{2}+q^{2}(q-q^{*})^{2}\right]^{-1}$, consistent with the behavior observed in the DWRLG [1, 2]. Upon examination of other system sizes, e.g., $L=128$ and $192$, we consistently find $q^{*}=0.101\pm 0.001$ for the disordered regime.

We conjecture that as $\rho$ is increased, the system first has no order, then long-range order just along $\hat{\mathbf{j}}$ (transverse to the drive), and lastly, transitions to a state of full order. A final demonstration of this can be done by binning the column-averaged particle densities:

[TABLE]

which lie in $\left[0,1\right]$. Histograms of these densities are sensitive to ordering only in columns, i.e., a measure of the presence of well-defined stripes (along the direction perpendicular to the drive), but not their widths. Exploiting the equivalents in the DF and CF description, we define

[TABLE]

Since $\bar{\phi}_{+}(i)$ is related to the hole density in a column, it carries crucial information on the location of interfaces between the two species. Our goal is to study histograms, $H\left(\bar{\phi}_{\pm}\right)$, constructed from $1000$ time points sampled every $10^{4}$ MCS in the steady-state, each of which provides $L$ values of $\bar{\phi}_{\pm}$.

When the system is in the homogeneous phase, we expect the histograms to have a simple Gaussian-like peak around [math]. As ordering builds in each column, the peak in $H\left(\bar{\phi}_{-}\right)$ should split symmetrically, much like a histogram of the magnetization in an Ising model when $T$ drops below $T_{c}$. This phenomenon indeed occurs, as illustrated by the blue and orange lines on the left panels of Fig. 7(d) for $\rho=0.39$ and $\rho=0.60$, which correspond to densities in the microemulsion phase (just below stripe formation) and the irregular stripe regime, respectively. In this context, we may roughly estimate $\rho_{\ell}\simeq 0.40$ from the peak split. With further increases of $\rho$, we see that the twin peaks move further apart, but with no obvious signal corresponding to the upper transition density $\rho_{u}$. As shown by the green line on the left panel of Fig. 7(d), the peaks are located close to the extremal values ($\bar{\phi}_{-}=\pm 0.9$) for the case of $\rho=0.80$. On the other hand, the density fluctuation histograms, $H\left(\bar{\phi}_{+}\right)$, seem to be sensitive to the upper transition density $\rho_{u}$ but not to the lower one, $\rho_{\ell}$, as the histograms seen on the right panels of Figs. 7(d) indicate. For $\rho=0.39$ (blue curve), $H\left(\bar{\phi}_{+}\right)$ is still sharply peaked around the origin. Then, we see that $H\left(\bar{\phi}_{+}\right)$ begins to develop considerable asymmetry for $\rho=0.60$ (orange curve), indicating the development of an (asymmetric) bimodal distribution. The histogram is then clearly bimodal at $\rho=0.80$ (green curve), where we see a sharp peak (near $\bar{\phi}_{+}\sim 0.1$) accompanied by a broader distribution of hole-rich columns ($\bar{\phi}_{+}\lesssim-0.2$). The transition in the shape of $H\left(\bar{\phi}_{+}\right)$ can be used to estimate $\rho_{u}$, yielding a consistent result with our previous analyses, $\rho_{u}\simeq 0.65$. Figs. 7(a-c) are the typical configurations of $\bar{\phi}_{\pm}(i)$, with the average DF(CF) in blue(red).

To summarize, we introduced a novel driven diffusive system: a hybrid between the driven Widom-Rowlinson lattice gas [1] and its associated statistical field theory [2]. Following the spirit of coarse graining a lattice gas, this “field-based lattice model” consists of continuous densities defined on a discrete lattice. Though the properties of the undriven version are, as expected, those of the WRLG, unexpected phenomena arise under a drive. While the two are comparable for $\rho\lesssim 0.40$ and share some similarities for $\rho\gtrsim 0.70$, conspicuous differences appear at intermediate densities. We propose the presence of an irregular stripe regime in which, qualitatively, there is long-range order perpendicular to the drive and disorder along the drive direction. To arrive at quantitative conclusions would require extensive simulations beyond the scope of the present study.

III Continuum Description

The out-of-equilibrium transition to a striped regime from a structured disordered regime (microemulsion) of both the FLM and the DWRLG for $\delta>0$ raises the possibility that similar behavior characterizes a broader class of models. Whether a continuum coarse-grained field-theoretic approach can describe such class is unclear, especially for high particle densities. In a previous study [1], the simplest continuum field theory, a linear stochastic partial differential equation (SPDE) based on a FLM-like model, was proposed, yielding a successful phenomenological picture of the low-density phase of the DWLRG. A more systematic derivation of a stochastic field theory, exploiting the Doi-Peliti formalism, was undertaken in Ref. [2], providing good approximations for the noise terms and the drive-dependent quadratic terms in the SPDE. Treating fluctuations perturbatively, the study recovered the chief novel feature of the disordered phase: a charge fluctuation structure factor peak at a nonzero wavelength ($\mathbf{q}^{\ast}\neq 0$), characteristic of a microemulsion. However, approximations necessary within the latter scheme render higher-order terms in the field theory questionable, motivating a different approach for finding a continuum description of the stripe formation at high densities and a full picture of the transition to stripes.

Here, we build on the success of the FLM, developed in Sec. II, examining its continuum limit [integer-valued $\mathbf{s}=\left(i,j\right)$ to continuous $\mathbf{r}=\left(x,y\right)$]. Key details of the calculation are given in Appendices A and B. Specifically, we seek a minimal coarse-grained model of a two-component, coupled diffusive system exhibiting a transition at $\rho_{c}$, such that:

• •

the system is phase-separated for $\rho>\rho^{\mathrm{FT}}_{c}$, exhibiting usual coarsening for $\delta=0$ and stripes perpendicular to the drive for $\delta>0$.

• •

the system is well mixed for $\rho<\rho^{\mathrm{FT}}_{c}$, exhibiting a homogeneous disordered phase for $\delta=0$ and a “microemulsion” with a characteristic length scale along the drive direction for $\delta>0$. We also expect a $\delta$-dependent discontinuity of the SFs at the origin.

Let us begin with a short overview of the procedure: we start with the rules introduced in Sec. II for mass transfer for $\rho_{A,B}(\mathbf{s},t)$ associated with all possible NN and NNN hops to/from site $\mathbf{s}$. Then, we set $\mathbf{r}=\mathbf{s}\ell$, where $\ell$ is a lattice spacing, and write the corresponding difference equations for the updated densities $\rho_{A,B}(\mathbf{s},t+\Delta t)$ due to a single hop occurring over a time interval $\Delta t$. Taking the $\ell,\Delta t\rightarrow 0$ limits and applying a gradient expansion (GE), we cast the difference equations in the form of driven-diffusion continuity equations for the coupled pair of charge and density fluctuations, $\phi_{\pm}(\mathbf{r},t)$:

[TABLE]

where $\mathbf{J}$ and $J^{\delta}$ are the drive-independent and drive-dependent currents, respectively (with $\ell$ and $\Delta t$ absorbed into the coefficients of the various terms in the currents). Seeking minimal expressions for the $J$’s via a GE, we will arrive at a continuum description [see Eq. (III) below] which we study in detail in the next section. Note that the noise terms cannot be derived systematically in this approach. Instead, we will add the simplest form of noise (zero mean, conserved-Gaussian), keeping the variance as a free parameter. To arrive at a minimal description, we limit ourselves to the lowest possible order, in both the fields and gradients, maintaining, in the spirit of Ginzburg-Landau theory, high-order terms only as needed for stability.

Deferring the details of the derivation to Appendix A, we arrive at the following SPDEs,

[TABLE]

where $\xi_{\pm}$ are conserved-Gaussian noises with zero mean, zero cross-correlation, and autocorrelation

[TABLE]

with $\sigma$ an effective noise intensity 777Note that, more generally, the noises acting on $\phi_{\pm}(\mathbf{r},t)$ originate from the noises in the underlying densities $\rho_{A,B}(\mathbf{r},t)$, which are inherently multiplicative: i.e., $\xi_{A}$ vanishes in regions devoid of A particles, and analogously for $\xi_{B}$. If we consider small perturbations from a state with $\phi_{+}\sim\rho$ and $\phi_{-}\sim 0$, the fluctuations are additive to a first approximation (as assumed here).

In the equation above, $D_{\pm}$, $\Gamma_{\pm}$ and $\lambda$ are $\rho$–dependent coefficients, while $g_{\pm,0}$ and $\Delta v$ are proportional to $\delta$. The diffusion constants for the DF and CF are given by $D_{+}$ and $D_{-}$, respectively. In equilibrium, the free energy cost per unit area to form an interface for each corresponding field is $\propto\Gamma_{+,-}$. Moreover, in analogy to Model B dynamics, the charge diffusion coefficient $D_{-}$ changes sign from positive to negative at the transition point and, as such, the fourth-derivative terms are necessary to ensure the stability of the theory as well. We have moved to the co-moving frame of the DF, so that the fields’ characteristic velocities are encoded in a single drive-dependent linear parameter, $\Delta v$. As pointed out in Ref. [2], and as will be shown in our numerical studies, $\Delta v$ plays the key role in generating both the periodic correlations in the low-density phase and stripes in the high-density regime 888A similar feature occurs in the analysis of kinematic waves in drifting crystals [79], where different characteristic velocities also play a central role in creating those periodic structures.. The nonlinear parameters $g_{\pm,0}$ and $\lambda$ account for the drive-independent and drive-dependent contributions to the effective $A$-$B$ interaction, respectively.

The choice of signs in Eq. (III) follows naturally from the GE of the FLM [43], so that the equation is manifestly stable by taking all couplings positive. The linear coefficients are fixed through the GE, but the functional dependence on $\rho$ of the nonlinear parameters is challenging to establish, which is why we shall keep them free. From the functional dependence of $D_{-}$ on $\rho$, we find that this quantity changes sign at $\rho_{c}^{\mathrm{FT}}=0.37\dots$, which is remarkably close to the $\rho_{c}^{\left(0\right)}\simeq 0.39$ result for the critical density estimated from simulations. The exact expression of the linear coefficients obtained via the GE, and key details of the calculation, are outlined in Appendix B. In summary, our formulation in terms of Eq. (III) depends on five free parameters: $\lambda$, $g_{\pm}$, $g_{0}$, and $\sigma$. As noted in Appendix B, the term proportional to $g_{+}$ can be neglected for $\rho>\rho_{c}^{\mathrm{FT}}$. For $\delta=0$, we have $\Delta v=g_{\pm,0}=0$, leaving only the terms proportional to $\lambda$.

Finally, it is interesting to observe that the large-scale properties of the $\delta=0$ system are isotropic, so that there is just one $D,\Gamma,\sigma$ for each field (instead of $D_{x},\,D_{y}$, etc.). Moreover, the fluctuation-dissipation relation constrains the noise correlations to be isotropic, as well. However, it is also known that this relation is violated in driven-diffusive systems and various anisotropies emerge, with the most significant manifestation being a discontinuity singularity of $S\left(\mathbf{q}\right)$ at the origin [20]. Although it is easy to argue qualitatively how the (anisotropic) drive induces such features, it is difficult to obtain quantitative predictions from the microscopic models. Thus, previous studies of the DWRLG [1, 2] tend to rely on phenomenological approaches to arrive at, say, the difference $D_{x}/\sigma_{x}-D_{y}/\sigma_{y}$. In Sec. IV, we show that, despite starting with isotropic choices for parameters like $D$, $\Gamma$, $\lambda$, and $\sigma$, the many surprising phenomena of our driven system 999As a reminder, a linear SPDE can be solved analytically. Not surprisingly, as they correspond to non-interacting fields, having isotropic parameters, the structure factors display only isotropic Ornstein-Zernike forms, despite the presence of $\Delta v$, e.g., the discontinuity singularity, emergence of a peak at $\mathbf{q}\neq\mathbf{0}$, and formation of stripes arise in our numerical solutions due to the interplay of the drive (via $\Delta v$), the interactions (via non-linearities), and the noise.

IV Numerical integration of the Stochastic Field Theory

In this section, we integrate Eq. (III) numerically to study how the possible steady states depend on the various parameters of the model. Our goal in this section is to assess whether the continuum stochastic field theory we propose in Eq. (III) is capable of reproducing the different regimes encountered both in the DWRLG and the FLM, namely the microemulsion structure in the low-density phase, and stripes aligned perpendicular to the drive in the high-density phase. Additionally, we check whether the continuum model also exhibits the irregular stripe regime present in the FLM.

To integrate Eq. (III) numerically in 2D with periodic boundary conditions, we design a numerical SPDE solver using the PSSETD2H and smooth dealiasing. To fix the timestep $\Delta t$, we run simulations at several values of $\Delta t$ and choose the largest value such that there are no appreciable changes in the results; in practice we use $\Delta t=0.5$. Unless otherwise specified, simulations are performed on a domain of size $L=128$, using grid resolutions of $N=200$ points per direction 101010This mesh size is chosen to balance numerical accuracy and computational cost. We have verified that further refining the mesh does not produce any appreciable changes in the results.. We also adopt the following initial conditions for the fields $\phi_{\pm}$:

[TABLE]

where $\sigma_{0}=10^{-2}$ and $\xi(\mathbf{r})$ is a conserved Gaussian noise with zero mean and unit variance.

In the low-density equilibrium regime, $\rho<\rho_{c}^{\mathrm{FT}}=0.37\dots$ and $\delta=0$, for $\lambda\geq 0$, both the CF and DF steady-state structure factors are isotropic and follow the usual Ornstein-Zernike form. For a fixed drive $\delta$ in this regime, exploring the full phenomenology of the model requires scanning a four-dimensional parameter space (or five, if the noise amplitudes are included), which is computationally demanding. Fortunately, we can simplify the task since our main goal is to find the various regimes of the lattice models. Thus, first, at low densities $\rho$, the non-linearities facilitating coarsening dynamics should be irrelevant (as $D_{-}>0$), so we set $\lambda=0$. Then, taking advantage of the perturbative analysis in Ref. [2], we set $g_{0}=4.243$, $g_{-}=0.707$, $g_{+}=2.828$, and $\sigma=0.27$, which should correspond to the microemulsion regime. Fig. 8 shows the steady-state charge fluctuation SF, $S_{-}(\mathbf{q})$, in the resulting microemulsion phase, averaged over $5000$ time points evenly spaced across the final $10^{6}$ time steps, and using $150$ independent realizations. As expected, $S_{-}(\mathbf{q})$ exhibits both a discontinuity at the origin and the emergence of a peak at a nonzero wavevector, $q^{*}=0.190\pm 0.009$, along the drive, closely matching the behavior of the FLM. The peak is well fit by the same fitting form used for the FLM, $A_{1}\,q^{2}\left[A_{2}\,q^{2}+q^{2}(q-q^{*})^{2}\right]^{-1}$. These results are consistent with previous studies [2], where it was also noted that the field theory smoothes the peak region, making it noticeably less sharp than in the lattice-based models.

Another simple case is the high-density equilibrium regime, $\rho>\rho_{c}^{\mathrm{FT}}$ and $\delta=0$. Here, all the drive-dependent couplings vanish and we simply set a positive value for $\lambda$. In this case, a steady-state amplitude is systematically found 111111This is supported by the weakly nonlinear stability analysis of Eq. (III) [43], which can also be used to crudely estimate the amplitude values slightly above the onset, $\rho\gtrsim\rho_{c}^{\mathrm{FT}}$. On the other hand, for $\lambda<0$, Eq. (III) has no bounded solutions and the numerical integration is ill-behaved, thus making necessary the inclusion of higher-order derivative terms. for the charge and density fluctuations. In all the runs performed, we consistently find two types of steady-state configurations: a round domain that resembles a bubble (see Fig. 9), or a single stripe whose orientation varies with realization. In either case, there is a single domain dominated by each species, which is consistent with the results of the FLM in the undriven regime. Furthermore, we verify that our equations of motion follow the Lifshitz-Slyozov law at high densities ($\rho>\rho_{c}^{\mathrm{FT}}$), with the growth of the SF peaks scaling as $t^{1/3}$, so that the characteristics of the coarsening dynamics do not depend on the particular value of $\lambda$. In fact, the cubic saturation terms due to $\lambda$ are more relevant to fix the magnitude of the steady-state amplitude and, as a result, their actual values may depend on the choice of $\lambda$. In the context of the Cahn-Hilliard equation or Model B [54, 55, 56], these couplings are often set to unity, which corresponds to a rescaling of the field values. Therefore, for simplicity, we fix $\lambda=1.0$ for the remainder of this section.

The driven, high-density regime involves a large parameter space. Here, we cannot rely on the coupling values obtained in Ref. [2], as they are specialized for the low-density phase. However, we can reduce the parameter space by discarding the term $g_{+}\phi_{+}^{2}$ in Eq. (III); see the simplifying assumptions in Appendix B. We then proceed by fixing a reasonable value for the noise strengths, $\sigma$. Fig. 9 compares typical runs with and without the conserved noises $\xi_{\pm}$ in (III) for $\rho=0.7$, $\delta=0$; and $\rho=0.7$, $\delta=1.0$, $g_{-}=0.5$, $g_{0}=0$. The linear coefficients are given by the expressions in (B). We find that quantitative properties in the ordered phase are largely insensitive to the noise amplitude within a moderate range 121212This is not true, however, in the low-density regime and at interfaces, as we have verified for the microemulsion phase. In fact, for $\rho<\rho_{c}^{\mathrm{FT}}$, the fields will decay to zero, even if $\delta$ is maximal. Therefore, the noises act as source terms, which are crucial to the development of the microemulsion behavior. In another words, the microemulsion phase cannot be observed without the noise terms in Eqs. (III)., consistent with observations for Models A and B in [54]. Accordingly, we fix the noise intensity to a moderate value, $\sigma=0.01$. At very high noise levels, however, the system displays qualitatively new behaviors, most notably stripe reorientation [see Fig. 9(c,d)], suggesting that the high-noise limit of the model may contain additional phenomena worthy of future investigation.

Therefore, following these simplifications, we are left with two nonlinear couplings to consider: $g_{0}$ and $g_{-}$. To sweep the parameter space, we shall limit ourselves to $g_{0},g_{-}\in[0,1]$, which is a typical range indicated by the different coarse-graining prescriptions in the high-density regime. We also vary $\Delta v$, which depends on the drive and density as shown in Eq. (B). Figure 10 summarizes the four different regimes we encounter via numerical integration:

Stripes perpendicular to the drive emerge already for $g_{-}=g_{0}=0$, solely due to the presence of the term proportional to the difference in the charge and density field characteristic velocities: $\Delta v=v_{-}-v_{+}$ [see Eq. (III)]. However, in this case, the stripes are wavy or vertically modulated due to the absence of the nonlinear, drive-dependent saturation terms. 2. 2.

For $g_{0}\leq g_{0}^{*}$ and $g_{-}>0$, the stripes retain their perpendicular alignment and the vertical modulation disappears. This suggests that the term proportional to $g_{-}$ in Eq. (III) provides stripe stability. This regime has the best agreement between the field equations, the FLM, and the DWRLG phenomenology. 3. 3.

For $g_{0}>g_{0}^{*}$ and $g_{-}\geq 0$, the perpendicular alignment is lost. In this regime, either horizontal stripes emerge or no stripe formation occurs at all, depending on the precise values of $g_{0}$ and $g_{-}$. Large values of $g_{0}$ destabilize vertical stripes, promoting alignment parallel to the drive. 4. 4.

For $0<\Delta v<\Delta v^{*}$, vertical stripes are not formed in any simulation. This supports the conclusion that while $\Delta v>0$ is a necessary condition for stripe formation with the correct alignment, it is not sufficient: $\Delta v$ must exceed a critical threshold $\Delta v^{*}$ to stabilize the vertical stripe pattern.

As an example, for $\rho=0.7$, we can estimate roughly that $\Delta v^{*}=0.08(2)$ (associated, for this density, to $\delta\approx 0.38$) and $g_{0}^{*}=0.02(1)$. Fig. 10 shows typical configurations for the previously mentioned regimes (1)-(4).

Our numerical results suggest that the term proportional to $g_{-}$ in Eq. (III) is necessary to stabilize ordered stripes, analogous to the drive-independent stability term proportional to $\lambda$. Thus, like $\lambda$, we will fix the value of $g_{-}$ in the following. The choice to fix $g_{-}$ is further supported by the observation that, for finite $\Delta v>\Delta v^{*}$, we verified that keeping $0<g_{-}\lesssim 0.5$ does not affect the stripe orientation. Additionally, our simulations show that increasing $g_{0}$ consistently destabilizes perpendicular stripes, moving the system away from the phenomenology of the lattice models. Accordingly, we set $g_{-}=0.5$ and $g_{0}=0$ in what follows 131313Note that it is not strictly necessary to drastically set $g_{0}=0$; any small value $g_{0}\approx 0$ would yield similar results. However, since no particular value is preferred, we choose zero for simplicity and computational efficiency.. In other words, we focus the remainder of our studies on Regime 2 shown in Fig. 10. In this regime, we explore the effect of varying $\rho$ [which will vary the linear terms shown in Eq. (B)] and $\Delta v$ (or equivalently $\delta$), and find that increasing either parameter also increases order along the drive direction. Conversely, when $\rho$ and $\Delta v$ are close to the thresholds $(\rho_{c}^{\mathrm{FT}},\Delta v^{*})$, we find evidence of long-wavelength (Eckhaus-like) instabilities within the perpendicular stripe regime, which act as a mechanism to reduce the number of stripes, similarly to what occurs in the FLM. Based on this, we identify the large and near-threshold parameter regimes as the regular and irregular stripe regimes, respectively. The resemblance goes beyond visual inspection of the configurations: we also detect the clear signatures of the different regimes in the charge SF, with the irregular stripe regime [Fig. 11(a)] showing a broad peak near the origin, whereas the regular stripe regime [Fig. 11(b)] displays a sharp, Dirac delta-like peak accompanied by pronounced higher-order harmonics. However, unlike the FLM irregular stripe regime, we observe weak harmonics around the peak and secondary maxima at $q_{x}>q^{*}$ for $q_{y}\neq 0$, indicating stripes that are more ordered along $x$ yet weakly modulated along 141414In fact, this vertical modulation seems to be a feature of the equations that is already present if $g_{\pm,0}=0$; see, e.g., Fig. 10  (1). $y$. Additionally, in the regular stripe regime, the continuum model further accentuates ordering: the SF is heavily localized at $q_{y}=0$, with modes at any $q_{y}\neq 0$ appearing almost as background noise; the resulting spectrum is strongly peaked and displays deep dips, resembling Bragg peaks in a diffraction grating that are consistent with a highly ordered stripe pattern.

For example, comparing with the results in Table 1, the ratio of peak values $S_{-}(q^{*},q_{y})$ for $q_{y}L/2\pi=0,1$ in the continuum model using $\delta=1$ is $S_{-}\left(q^{\ast},0\right)/S_{-}\left(q^{\ast},0.05\right)\simeq 52$ in the irregular stripe regime ($\rho=0.4$) and $S_{-}\left(q^{\ast},0\right)/S_{-}\left(q^{\ast},0.05\right)\simeq 2700$ in the regular stripe regime ($\rho=0.6$) (nearly an order of magnitude larger than the FLM!). Table 2 lists the peak values for the first three $q_{y}$ modes available in these two cases. In addition, we note that the density at which the regular stripe regime emerges in the continuum model appears to be shifted relative to the FLM: regular stripes are observed already at $\rho=0.6$ for $\delta=1$ (whereas in the FLM, $\rho_{u}\simeq 0.65$). However, note that $\rho$-dependence in the continuum model only refers to the GE results for the linear coefficients [see Eq. (B)]. These assumptions are violated at high densities and we generally expect the location of the transition to be shifted due to the nonlinearities and noise.

These findings are summarized in Fig. 11, which shows representative configurations with the charge fluctuation SF (averaged over $1000$ time points evenly spaced across the last $10^{6}$ time steps and $70$ realizations) for the irregular [Fig. 11(a)] and regular [Fig. 11(b)] stripe regimes. Together, they demonstrate that the continuum model not only reproduces the phenomenology of the FLM but also reveals subtle contrasts, such as the development of tiny harmonics in the intermediate regime 151515We can notice, however, a shoulder at higher $q_{x}$’s in Fig. 5(b) and an amplified jump in SF peak values in the regular stripe regime.

As long-range order builds up along the drive direction, the stripes become not only better defined but also increase in number. We find that the steady-state characteristic modulation of the CF SF, $q^{*}$, increases with $\Delta v$ for $\Delta v>\Delta v^{*}$ at fixed $\rho=0.7,g_{-}=0.5$ and $g_{0}=0$. Notably, this behavior is analogous to and consistent with the results reported in [2], where $q^{*}\sim\Delta v$ was observed in the low-density regime. In Fig. 12(a), we plot the dependence of $q^{*}$ with $\Delta v$, in the regime where perpendicular stripes are well established, i.e., for $\Delta v>\Delta v^{*}$. Panel (b) shows that the CF SF has an increasingly prominent peak along the drive direction as $\Delta v$ grows, signaling the increase of order in that direction, as observed in the FLM. A more detailed analysis of how $q^{*}$ depends on $\Delta v$ is left for future investigation.

Finally, by fixing $\rho$ and varying the drive-dependent parameters $\Delta v$, $g_{0}$, and $g_{-}$, we uncover a wide range of regimes in Eqs. (III), many of which have not been reported in the previous lattice-based studies. These regimes typically occur when $|g_{0}\phi_{+}|\gg\Delta v$. To (partially) understand the origin of these new regimes, we reorganize the drive-dependent terms proportional to $\partial_{x}\phi_{-}$ in Eq. (III) into the combination

[TABLE]

If $|g_{0}\phi_{+}|\gg\Delta v$, we see from Eq. (17) that the flow direction of the CF is strongly dependent on the sign of $\phi_{+}$, facilitating the formation of complex flow patterns.

Fig. 13 shows snapshots taken at long times $t=10^{6}$ that illustrate the variety of new behaviors: regular [Fig. 13(a)] and irregular [Fig. 13(b)] stripes aligned with the drive, spatio-temporal chaos [Fig. 13(c)] and phase separation into bubbles modulated by shock-like envelopes [Fig. 13(d)]. Because these regimes do not appear in either the FLM or the DWRLG, we have restricted our primary analysis to parameter regions that better reflect those models. Nevertheless, Fig. 13 clearly demonstrates that even a minimal set of continuum equations can generate a remarkably rich phenomenology. While a systematic exploration of how each term contributes to the various steady states would be valuable, such an endeavor is computationally demanding and falls beyond the scope of this study. It is worth noting that similar regimes have been observed in other contexts, such as parallel stripes in the KLS model [61] and spatio-temporal chaos in the active non-reciprocal Cahn-Hilliard model [25].

V Discussion and conclusion

We investigate a model of a driven binary mixture (dubbed the “FLM”) which effectively captures the behavior of the driven Widom-Rowlinson lattice gas in both the low- and high-density phases, but includes an additional intermediate regime with irregular stripes. In the high-density, ordered phase, earlier studies [1, 2] report standard coarsening dynamics in the absence of a drive, while the drive leads to the emergence of regular stripes oriented perpendicular to the drive. In Sec. II, we briefly characterized the FLM and observed striking similarities with the DWRLG phenomenology, suggesting that both models may belong to the same universality class within a properly chosen parameter space. Notably, we observe that perpendicular stripes emerge even when the system is initialized with parallel stripes [43], indicating a robust stripe-selection mechanism. This reorientation appears to proceed via an instability that is reminiscent of the Mullins-Sekerka instability of the interfaces [43].

Above the critical density, roughly estimated as $\rho_{c}\simeq 0.39$ for $\delta=0$ and $\rho_{\ell}=0.40$ for $\delta=1$, we identify two additional transition densities in the FLM: (i) For $\delta>0$, the density $\rho_{u}$ ($\simeq 0.65$, for $\delta=1$) marks a transition from short-range ordered (irregular) stripes to long-range ordered (regular) stripes [see, e.g., Figs. 1(a,d)- 5(b,c)]. The irregular stripe region, $\rho_{\ell}<\rho<\rho_{u}$, is characterized by pronounced long-wavelength instabilities, which serve as a mechanism to reduce the number of stripes. (ii) For $\delta\geq 0$, $\rho_{g}\simeq 0.90$; marks the onset of a frozen regime for uniform initial conditions [Figs. 1(a,b)], where coarsening arrests and glassy configurations appear. In Sec. II we offer a number of ways to characterize these two regimes.

In the driven low-density regime, the FLM shows indications of a peak in the CF structure factor [upper panel in Fig. 6(a)] at an $L$-independent characteristic wavevector, $q^{*}=0.101\pm 0.001$. The DF structure factor [lower panel in Fig. 6(a)] also shows a shoulder at $2q^{*}_{-}$, but remains peaked at $2\pi/L$, consistent with the results reported in Refs. [1, 2]. Examining a snapshot of a typical configuration in this regime [Fig. 6(b)], we see that this microemulsion phase emerges from a tendency to form stripes perpendicular to the drive that is frustrated by the stronger effect of the drive in regions having low particle density.

We derive a continuum version of the FLM via a gradient expansion in Sec. III. The analysis yields a system of two (reciprocally) coupled nonlinear stochastic partial differential equations, Eq. (III). These equations are consistent with the field equations developed in Ref. [2] for the low-density phase, but include additional terms that are important for describing the high-density regime, namely the terms proportional to $\lambda$ and $\Gamma_{\pm}$ in Eq. (III). Additionally, the charge fluctuation diffusion constant, $D_{-}$, changes sign at $\rho_{c}^{\mathrm{FT}}=0.37\dots$, marking the onset of the transition from a disordered to an ordered regime. For comparison, the FLM estimates are $\rho_{c}^{(0)}\simeq 0.39$ for the critical density at zero drive ($\delta=0$) and $\rho_{\ell}\simeq 0.40$ as the critical density for (irregular) stripe formation at $\delta=1$.

In Sec. IV, we integrate the field equations numerically. We find that the PSSETD2H method is best suited for this task 161616Interestingly, we found no references in literature that describe these methods in detail, leading us to believe that such a description (and application) could be useful, particularly in the active matter context, since the non-reciprocal Cahn-Hilliard and Active Model B equations are analogous to Eq. (III). Dealiasing techniques are straightforward to implement, effectively minimizing errors from spatial discretization, which is particularly important in the presence of noise 171717In fact, we find that it is necessary to use a much finer mesh in space and time than the one we use in Sec. IV to numerically integrate an equation of type (III) with a noise of moderate intensity, $\sigma\approx 10^{-2}$, without recurring to dealiasing techniques.. The remainder of Sec. IV provides a qualitative overview of the continuum stochastic field theory, with minimal emphasis on quantitative analysis. A central aspect is that noise plays a crucial role in reproducing the features observed in previous studies. We observe not only a discontinuity in the structure factors $S_{\pm}(\mathbf{q})$ at the origin but also a peak in the charge fluctuation structure factor $S_{-}(\mathbf{q})$ at a nonzero wavevector along the drive direction, an indicator of the microemulsion phase.

Turning to the high-density regime, the undriven ($\delta=0$) case displays standard coarsening dynamics consistent with the Lifshitz-Slyozov law in both the lattice model and the field theory. Our main findings for $\delta>0$ are summarized in Figs. 10-11: We identify a regime in the field theory, in which the equations of motion best reproduce the features of the lattice models (see Fig. 11 (2)). In this regime, as both density and drive increase, the continuum model transitions from irregular- to regular-stripe configurations, similarly to the FLM. Additionally, we uncover several novel regimes not previously reported in the lattice studies, including stripes aligned parallel to the drive and turbulent-like behavior (see Fig. 13).

A key finding of Sec. IV is that a nonzero difference in the characteristic velocities of the fields, $\Delta v$, is a necessary condition for the formation of stripes aligned perpendicular to the drive, analogously to what is found for the microemulsion phase behavior in Ref. [2]. When $\Delta v=0$, such stripes do not form, even when the other drive-dependent couplings $g_{0}$ and $g_{-}$ are nonzero. Moreover, the modulation of the CF field (and, consequently, the DF field) along the drive direction, $q^{*}$, increases with the magnitude of $\Delta v$, in agreement with findings from the low-density regime in [2]; see Fig. 12. However, vertically aligned stripes emerge only when $\Delta v$ exceeds a finite threshold.

Finally, it is worth emphasizing that the stochastic continuum equations, Eq. (III), contain just the minimal set of drive-dependent couplings, which, at the mean field level, cannot produce the core phenomenology of the lattice-based models such as the discontinuities and anisotropies of the structure factors. Such phenomena would normally only be accessible through renormalization group analysis of the structure factor corrections due to the noise terms and nonlinearities. Nevertheless, by (numerically) seeking steady-state solutions of the stochastic equations, we effectively access these corrections and find behavior analogous to the nonequilibrium lattice models. Due to both analytical and computational limitations, our study has focused primarily on qualitative, comparative analyses. Nevertheless, we hope that this work provides a useful starting point for developing a more complete field-theoretical framework for this broader class of nonequilibrium systems.

Acknowledgements.

We are grateful to Hugues Chaté and Alexandre Solon for many enlightening discussions, especially on the importance of dealiasing. G. F. O. and R. D. acknowledge support from CNPq, Brazil, under grant No. 161468/2021-7.

Appendix A Gradient expansion I – equations of motion

From the hopping rules stated in Section II, we can write the following equation for the evolution of $\rho_{A}\left(\mathbf{s},t\right)$ (and a similar one for $\rho_{B}$):

[TABLE]

where

[TABLE]

is the change in the field at site $\mathbf{s}$ over a time increment $\tau$, and $\Delta_{\mathrm{NN}}J_{A}^{\mathrm{NN}}$ is the discrete divergence associated with NN hops, etc. Explicitly, the $\Delta_{\mathrm{NN}}J$’s (with $t$ suppressed) are:

[TABLE]

where

[TABLE]

account for the NN $A$-$B$ exclusion for horizontal ($\leftrightarrow$) and vertical ($\updownarrow$) jumps, respectively, and $h(i,j)\equiv 1-\rho_{A}(i,j)-\rho_{B}(i,j)$ for the density of holes at site $\mathbf{s}=(i,j)$. There are analogous expressions for $\Delta_{\mathrm{NNN}}J$’s; see [43] for details. From these, we can construct evolution equations for the deviations (from $\rho$ and [math]) of the density and charge fields (defined in Eq. (5)) $\phi_{\pm}\left(\mathbf{s},t\right)$. Let us write these symbolically as

[TABLE]

Since the $J$’s are polynomials in $\rho_{A,B}$ (up to sixth order), expressions like those in Eq. (A) are polynomial up to $O\left(\phi^{6}\right)$. To arrive at the more familiar form of Eq. (13), we must take the continuum limit and, to simplify the divergence terms, keep only the terms of lowest order in the gradients. The remainder of this Appendix provides one possible approach.

Consider our lattice to have spacing $\ell$ in continuous space ($\mathbf{r}=\left(x,y\right)$), so that the fields above are defined on $\mathbf{r}=\mathbf{s}\ell$. One path to a continuum limit is to let $\ell\rightarrow 0$ and assume that $\phi_{\pm}\left(x,y\right)$ is a smooth function, so that, e.g.,

[TABLE]

Using the same approach for $\Delta_{t}\phi_{\pm}$ and letting the timestep ($\tau$) become infinitesimal, we arrive at

[TABLE]

Combining both sides and keeping only the relevant terms to fourth order in the derivatives and to second order in both fields, we obtain a set of equations similar to the deterministic parts of Eq. (III):

[TABLE]

From the above equation, we see that the undriven currents appearing in Eq. (13), $\mathbf{J}_{\pm}$, read

[TABLE]

while the drive-dependent currents, $J_{\pm}^{\delta}$, are given by

[TABLE]

Note that the form of the terms in Eqs. (A)-(A) is constrained by the $A\leftrightarrow B$ symmetry which implies that $\mathbf{J}_{+}$ ($\mathbf{J}_{-}$) cannot have odd (even) powers of $\phi_{-}$. Let us briefly highlight some characteristics of qualitative importance. First, based on this continuum limit of the FLM, we find that $D_{+}$ is positive for all $\rho$, which implies we can expect the DF to be stable, unless it is driven unstable through coupling to $\phi_{-}$. On the other hand, as expected of systems that display coarsening dynamics, $D_{-}$ changes sign when $\rho$ exceeds a critical value.

Exploiting a Galilean transform, we will work in the co-moving frame of the DF, so that the $v_{+}\phi_{+}$ term is absent in $J_{+}^{\delta}$ and $-v_{-}\phi_{-}$ is replaced by $-(\Delta v)\phi_{-}$, where $\Delta v\equiv v_{-}-v_{+}$. Note that both currents $J_{\pm}^{\delta}$ vanish when $\delta=0$, along with all their individual $v$ and $g$ terms.

Finally, we must add noise terms, $\xi_{\pm}$, which are difficult to derive from first principles, but we know they must conserve the total particle number and have zero mean. Taking the simplest possible approach, let $\xi_{A}$ and $\xi_{B}$ be the noise terms in the equations for $\rho_{A}$ and $\rho_{B}$. To first approximation, these are independent zero-mean conserved Gaussian noises, with variances $\sigma_{A}^{2}$ and $\sigma_{B}^{2}$. By symmetry under exchange of species $A$ and $B$, we must have $\sigma_{A}=\sigma_{B}$. In this case,

[TABLE]

have the same variance $2\sigma^{2}_{A}$, and are uncorrelated, $\langle\xi_{+}\xi_{-}\rangle=0$. Thus, it follows that $\xi_{\pm}$ are zero-mean conserved Gaussian noise terms, uncorrelated in time and with each other, having variance $\sigma^{2}\equiv 2\sigma^{2}_{A}$, i.e.,

[TABLE]

Combining all the ingredients, we arrive at the SPDEs

[TABLE]

For convenience, we collected all the linear/quadratic terms on the first/second line. In Appendix B, we show in some detail how to calculate the exact expressions of the linear coefficients. We also show how a set of assumptions is used to simplify the nonlinear terms in Eq. (25) to arrive at Eq. (III).

Appendix B Gradient expansion II – coefficients

To calculate the exact expression of the linear coefficients, we start by substituting

[TABLE]

in Eq. (18) to write the expressions in Eq. (A) (and the analogous expressions for species B) in terms of the density and charge fluctuations, but this time we only keep the resulting linear piece. Because the change of variables above involves the parameter $\rho$ (which is zeroth order in the fields), and since the $J$’s are sixth-order polynomials in $\rho_{A,B}$, the linear terms will have contributions coming from all possible orders in the fields. As a result, the corresponding linear coefficients are fifth-order polynomials in $\rho$.

The lattice update appearing in the different terms of Eq. (A) can be written in Fourier space in terms of the operator $e^{i\kappa\left(q_{x}+q_{y}\right)}$, such that for any field $f$ on the lattice,

[TABLE]

where $\mathcal{F}^{-1}$ is the inverse discrete Fourier transform.

We take the discrete Fourier transform of the linear part of Eq. (A) and use the above property to express the lattice-updated fields using the exponentials. There will be analogous expressions for NNN hops. We keep working in Fourier space to substitute these expressions in Eq. (18) and write

[TABLE]

where

[TABLE]

All coefficients are polynomials in $\rho$; their explicit forms are omitted for clarity and conciseness. Up to this point, no approximations have been introduced in the linear part: these are exact expressions, with $\mathcal{L}_{\pm}(\mathbf{q})$ denoting the exact linear operators corresponding to the Fourier-space representation of the linear terms in the mass-transfer equation (18). The gradient expansion (GE) of Eq. (27) yields, to lowest order,

[TABLE]

Comparing with expression (27) and collecting similar terms, we find the following coefficients 181818For simplicity, we suppressed all the “physical units” for the coefficients, e.g., $D\propto\ell^{2}\tau^{-1}$, etc. Their detailed forms play no role in the numerical studies.:

[TABLE]

The remaining seven nonlinear coefficients $\lambda_{\pm},\tilde{\lambda}_{\pm}$ and $g_{\pm,0}$ accompanying the quadratic terms in Eq.(25) are still undetermined. Determining their dependence on $\rho$ through the GE is nontrivial, as the expression (18) includes contributions of all orders in the fields $\rho_{A,B}$, not just quadratic. For this reason, we choose to leave the quadratic coefficients as free parameters in the present study. Work is in progress to determine the exact functional dependence of the different nonlinear couplings on $\rho$ [43].

This is a large parameter space (four $\lambda$’s, three $g$’s and $\sigma$), but a significant portion of it simply parametrizes the equilibrium $(\delta=0)$ behavior of the model. Therefore, we can introduce some simplifications to the equilibrium terms in order to better elucidate just the non-equilibrium features of interest. The specific simplifying approximations are

The GE analysis suggests $\tilde{\lambda}_{+}\gg\tilde{\lambda}_{-}\approx 0$ for most realistic values of $\rho$, so we set $\tilde{\lambda}_{-}=0$. 2. 2.

Since stability for $\phi_{+}$ is always provided by $D_{+}>0$, we will drop the term $\lambda_{+}\phi^{2}_{+}$ from our equations 191919Note that such a term originates from a $\phi_{+}^{3}$ potential in the free-energy and, as such, is unstable. Therefore, one must include an additional $\gamma_{+}\,\phi^{3}_{+}$ (a $\phi_{+}^{4}$ potential) in Eq. (25) to ensure the stability of the theory. Since we avoid over-complicating our equations by seeking a minimal description and also considering that the DF is stable to linear-order in $\phi_{+}$, we can safely neglect the term proportional to $\lambda_{+}$.. This choice is additionally supported by a weakly nonlinear stability analysis 202020Based on the approach of Cross et al. [3, 5], the weakly nonlinear stability analysis here consists of assuming that a dominant role in the instability of the homogeneous regime is played by a particular mode of $\phi_{-}$ with the Ansatz $\phi_{-}(\mathbf{r},t)=\alpha(t)\,\mathrm{Re}\,\{e^{i\Psi(t)+i\mathbf{q}\cdot\mathbf{r}}\}$. As $A$-$B$ exclusion tends to favor vacancies between the charged regions, it is reasonable to assume $\beta(t)\,\mathrm{Re}\,\{e^{i2\Psi(t)+i\mathbf{q}\cdot\mathbf{r}}\}$ for $\phi_{+}(\mathbf{r},t)$. In fact, the location of the peaks of the structure factors in the striped regime ($q_{-}=q^{\ast}=q_{+}/2$) displays this type of entraining of $\phi_{+}$ to $\phi_{-}$. Also, a perturbative scheme indicates $\beta\simeq\alpha^{2}$. Details are found in [43]. near the critical $\rho^{\mathrm{FT}}_{c}$. We find that $\phi_{+}$ is of the same order as $\phi_{-}^{2}$, so we drop terms with powers of $\phi_{+}$ larger than 1 when we study the high-density regime. 3. 3.

The GE suggests that the two remaining drive-independent nonlinear couplings, $\lambda_{-}$ and $\tilde{\lambda}_{+}$, are generally of the same order, but of opposite sign: $-\lambda_{-}\sim\tilde{\lambda}_{+}$, with $\tilde{\lambda}_{+}>0$. Thus, we reduce these couplings to a single effective $\lambda$, setting $-\lambda_{-},\,\tilde{\lambda}_{+}=\lambda$.

It is also important to observe that, just as in mean-field theories of other statistical systems, different coarse-graining prescriptions (e.g., the present GE of the FLM or the MSRJD/Doi-Peliti developed in [2]) lead to different values for these coefficients, but they generally agree on their order of magnitude. Using the set of simplifying assumptions in Eq. (25) together with the linear coefficients in Eq. (B), we can write down the final form of the field equations used in this work: Eq. (III).

Appendix C Pseudospectral exponential time-differencing method

In Fourier space, Eq. (III) for $\tilde{\Phi}\equiv\tilde{\Phi}(\mathbf{q},t)=(\widetilde{\phi}_{-}(\mathbf{q},t),\widetilde{\phi}_{+}(\mathbf{q},t))^{T}$ can be written generally as

[TABLE]

where $\mathrm{NL}[\widetilde{\Phi}]$ is the matrix of the nonlinear terms Fourier transform and

[TABLE]

is the matrix of linear operators.

The noise $\eta\equiv\eta(\mathbf{q},t)$ is a vector with complex Gaussian random variable components, with zero mean and correlation

[TABLE]

Note that, since $\Phi(\mathbf{r})$ is a real field, then $\tilde{\Phi}(\mathbf{q})$, the Fourier transform, satisfies the condition $\tilde{\Phi}(-\mathbf{q})=\tilde{\Phi}^{*}$. This should be true for the noise as well.

The pseudospectral method is designed to numerically integrate partial differential equations, where we work directly with Eq. (30) in Fourier space to solve Eq. (III) by choosing an appropriate time-discretizing scheme [67]. The word pseudo refers to the fact that nonlinearities are calculated in direct space by means of the inverse Fourier transform, therefore we do not work fully in reciprocal, or Fourier, space. This method to calculate the nonlinearities introduces aliasing errors [68], which must be dealt with using dealiasing techniques, see Appendix D.

The coefficients of the fields Fourier transform are calculated using the FFTW fast Fourier transform (FFT) routine [69]. Many different conventions are used in literature, such that we advise checking the FFT routine documentation before any implementation. On the $L\times L$ periodic lattice with a $N\times N$ mesh, the physical wave modes $\mathbf{q}$ are related to the index modes by

[TABLE]

where the components $k_{x,y}$ of the index vector $\mathbf{k}$ are

[TABLE]

Using the FFT in a lattice transforms Eq. (30) into a set of $N\times N$ coupled ordinary differential equations, where the Fourier amplitudes of the fields are stored as $N\times N$ complex arrays. We numerically integrate this system of equations in time by using the exponential time-diferencing (ETD) scheme, which is a class of methods that are well adapted to equations of type (30) [35].

The ETD scheme involves a specific change of variables that is similar to the interaction picture in quantum mechanics. Let

[TABLE]

where the exponential of the matrices are understood as the exponential of each entry (as they are diagonal in the charge/density basis). Then, the new variable $\widetilde{U}\equiv\widetilde{U}(\mathbf{q},t)$ satisfies the equation

[TABLE]

We are now in position to use any time-stepping method, such as Euler or Runge-Kutta, to integrate the new equation in time; each will result in a different ETD scheme. Here we use the stochastic second-order Heun (SH) method [36]. It is worth noting that the SH method achieves a balance between precision and simplicity. The deterministic part is approximated with an error of $O(\Delta t^{2})$, while the stochastic part has an error of $O(\Delta t^{1/2})$. Although stochastic versions of the Runge-Kutta (RK) family exist, they are primarily limited to second-order (RK2) methods and are challenging to implement [70].

Reversing the change of variables and denoting $\bullet_{n}$ as the fields at time $t=n\Delta t$, SH involves the calculation of the auxiliary field, $\tilde{A}_{n}\equiv\tilde{A}_{n}(\mathbf{q})$, given by

[TABLE]

where $\zeta_{n}\equiv\zeta_{n}(\mathbf{q})$ are complex Gaussian random variables, with zero mean and correlation

[TABLE]

For discrete $\mathbf{q}$, we further note the following correspondence [36, 54]

[TABLE]

where $\Delta x=L/N$ is the spatial resolution. The fields at the next timestep, $t+\Delta t$, are calculated as

[TABLE]

Observe that, at each iteration, we must repeat the procedure above for every value of $\mathbf{q}$, which is indexed by relation (31).

Appendix D Dealiasing

Although many authors have commented on the importance of dealiasing in pseudospectral methods, e.g., [71, 72, 73], such considerations are rarely mentioned in the physics literature. In particular, dealiasing is extremely important in fluid dynamics simulations, see e.g. [72, 37, 74, 75]. Higher modes (or frequencies) which are not resolved, given our spatial resolution, are aliased to lower modes, generating spurious contributions to their amplitudes. This effect is referred to as “aliasing error” and is especially significant when we must compute nonlinear terms in the equation of interest, or in the case of stochastic equations, since the noise will constantly generate amplitudes for these higher modes. An equation such as (III) will then generally suffer from aliasing errors.

A general procedure to reduce aliasing error is to continuously set to zero amplitudes of modes beyond a certain threshold. Following [37], we define two dealiasing kernels

[TABLE]

where $N$ is the number of points in the lattice and $k=|\mathbf{k}|$ is the modulus of the mode index in Fourier space (see Appendix C). The effect of the first kernel is simply to abruptly set to zero $1/3$ of the amplitudes of the higher modes, while the second kernel sets the amplitudes of the higher modes to zero smoothly. Using the smooth dealiasing kernel, we keep about $12\sim 15\%$ more modes than the sharp dealiasing kernel. For $\rho_{\mathrm{smooth}}$, we use in the exponential the same factor of $36$ that is employed in [37].

In our numerical integrations, we study the effect of both dealiasing methods, but no appreciable differences in the simulations are found. Conversely, if dealiasing is not applied, spurious contributions from higher modes contaminate the small-scale dynamics, appearing either as artifacts in the results or as numerical instabilities. We choose to implement smooth dealiasing in the present work as it has several advantages in comparison with the $1/3$ dealiasing, as discussed in Ref. [37]. The smooth dealiasing is implemented by modifying $\mathrm{NL}[\widetilde{\Phi}_{n}]$ when we compute $\Phi_{n+1}$. Namely, we multiply these non-linear terms by the dealiasing kernel,

[TABLE]

effectively setting to zero the large-$\mathbf{q}$ modes.

If the time-discretizing scheme is a multi-staged method, such as RK2 or RK4, we must implement dealiasing at each stage. As an example, for RK2, two dealiasing operations are required, one for each stage. Therefore, dealiasing is computationally very costly. The method described here is an example of explicit dealiasing. It is worth mentioning that there are other dealiasing techniques, such as implicit [76], random phase shift [77] dealiasing and, more recently, physics-informed machine learning dealiasing [78].