Driven Widom-Rowlinson lattice gas

# Driven Widom-Rowlinson lattice gas

## Abstract

In the Widom-Rowlinson lattice gas, two particle species (A, B) diffuse freely via particle-hole exchange, subject to both on-site exclusion and prohibition of A-B nearest-neighbor pairs. As an athermal system, the overall densities are the only control parameters. As the densities increase, an entropically driven phase transition occurs, leading to ordered states with A- and B-rich domains separated by hole-rich interfaces. Using Monte Carlo simulations, we analyze the effect of imposing a drive on this system, biasing particle moves along one direction. Our study parallels that for a driven Ising lattice gas – the Katz-Lebowitz-Spohn (KLS) model, which displays atypical collective behavior, e.g., structure factors with discontinuity singularities and ordered states with domains only parallel to the drive. Here, other novel features emerge, including structure factors with kink singularities (best fitted to ), maxima at non-vanishing wavevector values, oscillating correlation functions, and ordering into multiple striped domains perpendicular to the drive, with a preferred wavelength depending on density and drive intensity. Moreover, the (hole-rich) interfaces between the domains are statistically rough (whether driven or not), in sharp contrast with those in the KLS model, in which the drive suppresses interfacial roughness. Defining a novel order parameter (to account for the emergence of multistripe states), we map out the phase diagram in the density-drive plane and present preliminary evidence for a critical phase in this driven lattice gas.

## I Introduction

Driven lattice gases have played a central role in far-from-equilibrium thermodynamics and statistical mechanics, in particular in the study of phase transitions schmittmann-zia (); marro () and steady-state thermodynamics (SST) ST (). In this context, most studies have focussed on the lattice gas with attractive nearest-neighbor (NN) interactions pioneered by Katz, Lebowitz, and Spohn KLS () (KLS). In equilibrium, this model is the lattice gas version of the Ising model LeeYang (). Under a drive favoring particle motion along one axis (and suppressing motion in the opposite sense), the KLS model exhibits a variety of remarkable collective behaviors schmittmann-zia (); marro (), such as generic long-range correlations LRC (); Onuki79 (); Onuki81 (); Kirkpatrick16 (), ‘negative response’ ZiaPM (), non-Ising critical properties JSLC (), anisotropic phase separation KLS (), and anomalous interfacial fluctuations leung88 (); InterfaceCorr (). More recently, it has been extensively studied in connection with the difficulties in formulating a zeroth law of thermodynamics for systems driven into steady states far from equilibrium pradhan (); sst1 (); sst3 (). Another model of interest is the driven lattice gas with NN exclusion (NNE lattice gas) drnne (), which exhibits a complex pattern of phase ordering and jamming, and has also proved useful in efforts to define a nonequilibrium chemical potential sst2 (). The NNE lattice gas is an athermal model exhibiting a phase transition to sublattice ordering. While athermal models are convenient in that they involve only a single effective parameter (the dimensionless chemical potential ), the NNE lattice gas does not exhibit macroscopically distinguishable phases. An athermal model that does exhibit phase separation is the lattice Widom-Rowlinson (WR) model or Widom-Rowlinson lattice gas (WRLG) wrlg (). It is likely the simplest model with purely repulsive interactions to undergo phase separation. This study is devoted to a driven WRLG.

The equilibrium WRLG was introduced in wrlg (), as a discrete-space version of the original (continuous-space) Widom-Rowlinson model WR70 (); RW82 (). Each site may be vacant, or occupied by a particle belonging to species A or species B. Nearest-neighbor A-B pairs are prohibited, as is multiple occupancy of a site. (The latter restriction departs from the original, continuous-space WR model, in which particles belonging to the same species do not interact at all.) In wrlg () it is shown that the WRLG (with equal densities of A and B particles, ) exhibits a continuous phase transition between a disordered, low-density phase and coexisting A- and B-rich phases at a critical density . Best estimates for the critical density are 0.618(1) and 0.3542(1) on the square and simple cubic lattices, respectively. (Figures in parentheses denote statistical uncertainties.) Results on scaling behavior wrlg () support an Ising-like critical point, as expected on grounds of symmetry. Subsequent studies of the WRLG have focussed on multi-species versions Nielaba-Lebowitz (). Here we study the driven WRLG: as in the KLS model, and the driven NNE lattice gas, particle motion is favored along one direction, and suppressed in the opposite sense, as if an electric or other external field were acting on the particles. To model a nonequilibrium steady state, it is essential that the system have periodic boundaries along the drive direction. (For impenetrable boundaries the system corresponds to an equilibrium fluid subject to a simple linear gravitational potential.)

In the NNE lattice gas above the critical density, the sublattice-ordered phases are equivalent under a unit translation in the or direction, whereas the symmetry connecting the A- and B-rich phases in the WRLG is global exchange of A and B particles. In the driven NNE lattice gas with first- and second-neighbor hopping (the model studied by Szolnoki and Szabo, szolnoki-szabo02 ()), particles are free to hop between sublattices, so there is global sublattice ordering rather than spatial separation into two phases; in this case no stripe pattern is observed.

In the single-component lattice gas with finite repulsive NN interactions, (equivalent in equilibrium to the antiferromagnetic Ising model), the drive destroys sublattice ordering szabo94 (); szabo96 (). For infinite NN repulsion, i.e., the nearest-neighbor exclusion (NNE) model, phase separation under a drive is possible, as well as formation of dynamically arrested regions, if the dynamics allows NN hopping only drnne (). Including hopping to second neighbors changes the phase behavior: instead of separation into high- and low-density regions, there is a continuous, Isinglike transition to a phase with sublattice ordering, as in equilibrium szolnoki-szabo02 (). In the present work, the dynamics includes both NN and second-neighbor hopping.

There are two principal motivations for studying the driven WRLG. First, aside from the examples already cited, not many driven systems exhibiting phase separation have been examined in detail; the lack of attractive interactions leads to a very different sort of phase ordering under a drive, as compared with the KLS model. In this context, we should mention that the WRLG is a spin-1 or 3-state system, a class that includes well known systems such as the Potts Potts () and Blume-Emery-Griffiths BEG () models. Some driven 3-state models have been investigated previously Threestate (), but none display the novel properties of the system we present here. Second, it is valuable to have further examples of model systems for testing SST.

There are two main foci in the present work: anomalous correlations (even in the disordered phase) and remarkably rich varieties of phase segregated order. Our principal source of information is Monte Carlo simulation of the model on the square lattice, which we complement with simple mean-field theory approaches. The first issue is best displayed in terms of the structure factor , with the subscripts denoting wavevector components along, and perpendicular to, the drive. In addition to displaying a discontinuity singularity at the origin [] as observed in similar driven diffusive systems Huang93 (); Zia95 (); Korniss97 (); Praestgaard00 (); Vidigal07 (); Greulich08 (), exhibits a kink singularity, best fitted to ! This novel feature is accompanied by a maximum of located at some non-vanishing value, associated with periodic correlations of composition fluctuations along the drive direction. Let us emphasize that these properties are present even at low densities (for which the composition is spatially uniform), in stark contrast to the familiar, analytic, Ornstein-Zernike form . Under the drive, composition variations are characterized by a preferred wavelength , where denotes the intensity of the drive. The second focus is a phase diagram in the density-drive plane, as we report the emergence of multiple stripes of A- and B-rich phases [of wavelength ], oriented perpendicular to the drive! By contrast, in the KLS model, the particle-rich (or hole-rich) regions form a single stripe parallel to the drive (especially in systems with aspect ratios). As in the KLS model, long-lived metastable states are observed as control parameters are changed across phase boundaries, a phenomenon reminiscent of hysteresis.

The remainder of this paper is organized as follows. In the next section we define the model. In Sec. III we examine the low-density, disordered phase, and discuss mean-field and field-theoretic approaches. Sec. IV reports our results for local and global ordering, correlations in the ordered phase, and interface roughness, followed by a summary and discussion of our findings in Sec. V .

## Ii Model

The WRLG is defined on a lattice of sites, each of which may be empty or occupied by a particle of type A or of type B; multiple occupancy is forbidden. It is convenient to introduce occupation variables , , or for site vacant, occupied by an A particle, or occupied by a B particle, respectively. To model the NN repulsion, if , none of its nearest neighbors may have . In equilibrium, the model may be studied in the canonical ensemble, i.e., with fixed numbers of A and B particles, and , respectively, or in the grand canonical ensemble, in which the associated chemical potentials, and are fixed. Here we limit our attention to the case of square lattice (). Thus, a site is labeled by integers , while is denoted by . Unless explicitly noted otherwise, our study focuses on systems with . For simplicity, we focus on systems with fixed , and choose the main control parameter to be the overall density: .

It is well to recall that the equilibrium thermodynamics of athermal systems is determined exclusively by maximization of entropy, subject to whatever constraints apply. Thus phase separation in the WR model occurs at densities high enough that the configurational entropy of the phase-separated, inhomogeneous system is higher than that of a homogenous one. As in the freezing of the hard-sphere fluid, a higher entropy is associated with a nominally more ordered phase because the disordered phase possesses relatively few configurations, due to excluded-volume constraints. To what extent such considerations apply to driven athermal systems is an open question: under a drive, the probability distribution is no longer uniform on configuration space, and the very definition of a thermodynamic entropy is in general unclear.

The driven WRLG is a stochastic process with a Markovian dynamics defined via a set of transition rates , from configuration (i.e., the set ) to configuration . Transitions occur via single-particle hopping, so that for to be nonzero, configurations and must differ by the exchange of a particle-hole pair. Specifically, each particle attempts to hop to a first or second neighbor, that is, a displacement of with , (excluding the no-hop case ). Any hopping move that results in a configuration satisfying the prohibition against NN A-B pairs is accepted. The attempt rates are parameterized in terms of (propensity to move along ) and (tendency to leave unchanged):

 w(Δnx,Δny)=⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩(1−a)(1±p)6,Δnx=±1, any Δnya2,Δnx=0,Δny=±10,Δnx=0=Δny. (1)

For and , detailed balance is satisfied, and the stationary distribution of the Markov process corresponds to the equilibrium distribution, i.e., all allowed configurations equally likely. The hopping rates are isotropic for and . (The special case of , in which particles are restricted a given column, is not considered here.) As is varied from zero to unity, the process interpolates from a case with symmetric hopping in (equilibrium) to one with maximal asymmetry, with no hopping in the direction. We consider a dynamics with next-nearest-neighbor as well as nearest-neighbor hopping to avoid the possibility of nongeneric patterns or dynamic arrest associated with a restriction to NN hopping szolnoki-szabo02 (), and because an extended set of hopping moves favors ergodicity. Our choice of equal hopping rates for jumps with the same value of is motivated by simplicity. In the simulation studies reported here, we set and only comment on some runs with .

### ii.1 Simulations

We perform Monte Carlo (MC) simulations the driven WRLG on the square lattice, using rectangular systems of sites with periodic boundaries in both directions. Unless otherwise noted, we use , with ranging from 30 to 400. We count one MC step as attempted particle moves. Studies of stationary properties use MC steps, preceded by an initial period of MC steps, which we found to be sufficient for relaxation. Averages and uncertainties (standard deviation of the mean) are calculated over 10-30 independent realizations.

## Iii Anomalies in the homogeneous phase: Singular structure factors and correlations

If the density is sufficiently low, the WRLG remains homogeneous whether driven or not. Simulations of the equilibrium case show this disordered phase to prevail for . As the drive parameter is increased from zero to unity, the change in is rather modest, attaining a value of about for . By contrast, the critical temperature of the driven KLS model is found to increase by as much as (in ) from equilibrium KLS (). In this section, we focus our attention on the low-density disordered state, in which we find anomalies in the correlation functions and structure factors, beyond those observed in the KLS system.

### iii.1 Simulation studies of the structure factor and correlation functions

One standard way to characterize collective behavior is through correlation functions () and their Fourier transforms, the structure factors (). For simplicity, we consider only equal-time, two-point correlations; since there are two species of particles, these are matrices. Instead of the densities of the two particle species, , symmetry guides us to consider their sum and difference, which will be referred to as “mass” and “charge” densities

 ρ(→x) ≡ ρA(→x)+ρB(→x), ψ(→x) ≡ ρA(→x)−ρB(→x),

respectively notation (). While this notation is unnecessary for reporting simulation data (as they are just, respectively, and at each ), it sets the stage for a field-theoretic approach, the simplest of which will be discussed below. The restrictions of our simulations correspond to the constraints and .

In terms of these densities, the matrix of correlations consists of and . The cross correlation, , is the difference which vanishes (statistically) in the symmetric systems we study. For simply, we write as and . In the steady state, translational invariance dictates that these are functions of , the distance between the two points in question. Thus,

 Gρ(→r) ≡ ⟨|σ(→x)||σ(→x+→r)|⟩−ρ2 Gψ(→r) ≡ ⟨σ(→x)σ(→x+→r)⟩

are independent of ; to improve statistics, we average over . When we form the covariance or connected correlation function by subtracting the associated expectations, , these ’s assume different values at the origin: vs. , since . Meanwhile, the sum vanishes for both functions. The structure factors are the Fourier transforms

 S(→q)=∑→rG(→r)ei→q⋅→r (2)

By averaging over to obtain , we also have, e.g., , which is times , the average over a run of the power spectrum beam (), of each configuration . Finally, the constraints on translate into , , and . These serve as useful checks on simulation data. In addition, the last constraint implies that, unlike the ordinary Ising model, cannot serve as a susceptibility-like quantity for detecting a second-order critical point. Instead, our assumes its maximum value at some non-vanishing . The maximum value behaves like the susceptibility of ferromagnetic systems, i.e., remaining finite (as ) in the disordered phase. For large enough , as we report in Sec. IV, this value diverges as and so is an ideal candidate for an order parameter. If the transition turns out to be critical (i.e., behaves much like second-order transitions in equilibrium statistical systems), then the large- behavior of can provide us with a critical exponent.

In general, the two-point equal-time correlation, , must be symmetric under parity () and statistically symmetric under the reflections ( and ). In equilibrium, there is an additional symmetry, the exchange (especially for samples with ). By contrast, the drive should break the last symmetry. Thus, in addition to studying and in the whole plane, we will pay special attention to the specific directions along and perpendicular to the drive. In particular, we will consider, e.g.,

 Gψ||(r)≡⟨σ(→x)σ(→x+r^x)⟩;  Gψ⊥(r)≡⟨σ(→x)σ(→x+r^y)⟩ (3)

and

 Sψ||(k)≡Sψ(q||=k,0);  Sψ⊥(k)≡Sψ(0,q⊥=k) (4)

Let us call attention to the different roles played by the quantities and . The former is a wave vector, being the conjugate to the vector , with components . The latter is a wavenumber, measuring the magnitude of along one of two axes, in the same spirit that is the magnitude of along one of the axes. Thus, and are not simply related to each other by Fourier transforms. Instead, is the transform associated with the average charge density over a column and its correlation function . Similarly, we can define the transverse counterparts: and . (Of course, translational invariance of the steady state provides the or independence of the ’s.) We will not discuss these functions explicitly below, but focus on the ’s, as the latter contain the same information about the system, displayed in a better form.

With this setup, let us turn to the remarkable behavior revealed by simulations. For simplicity, we focus mainly on and occasionally mention results for . Thus, we drop the subscripts on and as well. Considering first the ’s, we report results found mainly in a system. The reason for this choice is that the novel and most interesting features are best displayed here.

In Fig. 1, we show and immediately notice a feature quite distinct from the structure factor of the KLS model, namely, the presence of “twin peaks.” Despite being in the disordered, homogeneous phase, such prominent structures represent a major deviation from the simple Lorentzian form of the equilibrium WRLG, as well as the discontinuity singularity at found in the KLS model KLS (); schmittmann-zia (), . While the data show decreasing monotonically as increases from zero, rises to a maximum at wavenumber before dropping. Such a peak is indicative of a preferred wavelength induced by the drive: .

To be more confident of this unusual behavior in , we considered finite size effects, illustrated in Fig. 2. As the figure shows, the data points for and fall well within statistical errors of the sample, so that, for this density and drive, the thermodynamic limit is well represented by the system. Thus, we see that the small- region of is well fit by . Alternatively, we can write , which displays the presence of another length, , a candidate for a “correlation length.” Note that must be even in , as the equal-time two-point correlations in the steady state must be even in . Thus, the presence of peaks symmetric about is inevitable. Although, in principle, satisfies these symmetry conditions, Fig. 3 (diamonds) shows that there is, without doubt, a kink singularity () at , so that the quadratic expression proposed above is a much better fit to . In this regard, though a kink singularity is not necessary to generate period structures, the observed oscillations are best modeled by one. The inset of Fig. 2 shows that the data for are well fit by a quadratic expression, , with

 kmax≅0.56;  λ≅11.

This form implies that . Meanwhile, Fig. 3 (circles) shows that is consistent with the ordinary Lorenztian form. Of course, the drive induces anisotropy, so that in general. As in the KLS model, only one of the correlation lengths, , diverges as the critical point is approached, as discussed below. Our best phenomenological estimate for the structure factor (for small , for this sample) is

 S(→q)=ν||q2||+ν⊥q2⊥τ||q2||+τ⊥q2⊥+q2||(∣∣q||∣∣−kmax)2+2γ×q2||q2⊥+γ⊥q4⊥+... (5)

While such a form may seem unusual, it is based on a stochastic field theory which has proven successful in describing the KLS model JSLC () (the numerator and denominator being associated with the noise and deterministic part of a Langevin equation for the density field). The main novel feature here is the presence of . (We also note that at the critical point, it is that vanishes in KLS, whereas vanishes in the present case.) To retrieve the undriven case, we need only set , , and . Though this form most likely has limited validity (e.g., the presence of the singular expression for all ), there is no good theoretical basis for us to propose additional terms at present. Thus, we will discuss the salient features associated with this expression.

Recall that, in our notation, is given by Eq. (5) with and , while is obtained when we set and . As in the KLS model, the discontinuity singularity is induced by violation of detailed balance JSLC () and displayed through

 Δ≡ν||τ||+k2max−ν⊥τ⊥≠0

Considering the points nearest the origin (for this system size), we have and (Fig. 3), so that . This singularity is manifested in configuration space by decaying as at large distancesLRC (). Although such a power will dominate the ordinary exponential decay for sufficiently large , the crossover will depend on the ratio of the amplitudes of the two contributions. In this case, there appears to be only a small regime () where power laws finiteL (), are found in both and . While the data are consistent with being , the case for appears to support the power instead. If more extensive simulations bear out this scenario, then it will be a challenge to understand such anomalous, yet “generic” (in the sense of being present far from critical points) singularities.

Let us return to the more prominent feature, namely, the presence of a maximum at non-vanishing , accompanied by the singular at the origin. Not surprisingly, the main consequence is oscillations in , within an exponential envelope. In the case, , so that the oscillations implied by are not easily seen. These features are however prominent at larger densities, where is more comparable to . Meanwhile, as might be expected, no oscillations are found in . In Fig. 4 we show both and for systems with and (again, with , p=1 and ). The rapid increase in over that for , along with hardly any changes in , implies that the power-law decays in are mostly masked. To be more confident of the nature of such powers (expected or “anomalous”), we will need to perform longer simulations (for statistics) of larger systems (to open a larger window between and ).

As shown in Fig. 2, deep in the disordered phase (), finite-size effects appear to be minimal already for . It is clear that, for this density, there is no need to simulate larger systems to capture the essence of the prominent properties of the driven WRLG. By contrast, in the next section, we show that the maximum in increases with , so that, for sufficiently large densities, the value does not converge as . Instead, will diverge as the system size, so that can serve as an order parameter, in the same spirit as is a measure of the (square of the) spontaneous magnetization in the non-conserved Ising model in two or more dimensions. Along these lines, we mention for comparison, the divergence of or for the 2D Ising lattice gas and of alone for the KLS model. Different from these examples, here diverges at an -independent wave-vector over a sizable region of the - plane.

### iii.2 Discrete mean-field theory and its continuum limit

In this subsection, we present the simplest theoretical description, as the first, small step towards a more comprehensive theory. The idea here is to consider two continuous variables at each discrete site, each corresponding to the average occupation by the two species: . Next, we postulate the simplest gain/loss terms for the change in in a single time step. Instead of writing every term here, let us illustrate with just one example: the loss to due to particles hopping along (the argument is dropped for simplicity):

 (1−a)(1+p)6ρA(x,y)h(x+1,y)[1−ρB(x+1,y+1)][1−ρB(x+1,y−1)][1−ρB(x+2,y)]+ (6) +(1−a)(1−p)6ρA(x,y)h(x−1,y)[1−ρB(x−1,y+1)][1−ρB(x−1,y−1)][1−ρB(x−2,y)] (7)

Here, is the density of holes. The rest of the analysis is straightforward, though quite tedious. It is clear that such terms would be exact if we had written the averages of the products of the appropriate s. In a mean-field (MF) approximation, these are replaced by the appropriate products of the averages, e.g., . Beyond the site approximation pursued here, a (generally) more reliable (and complicated) approach is the pair approximation, which treats nearest-neighbor two-site joint probabilities as the basic elements marro (). Although we plan to implement the pair approximation in future work, it seems possible that only more sophisticated theories will be able to capture the periodic structure observed under a drive.

Since is a conserved density, the gain and loss terms are necessarily representable in terms of divergences of currents. Identifying the terms associated with the net hops from, say, to , we find the leading contribution by setting both densities to be the overall value: . This is the same approach used to find the current-density relationship in the KLS model. The result here is

 JMF(ρ)=(1−a)pρ2(1−ρ)(1−ρ2)3 (8)

where the various factors can be readily associated with those in, say, Eq. (6).

In Fig. 5, we show this (dashed line) and the simulation data (symbols) for a system of size with , as well as (dotted line) for comparison. It is remarkable how the properties of a complicated stochastic system are mostly captured (i.e., ) by such a simple minded approach. Near the limit, we note that the data jump to follow another “branch” of the - relationship, which is much closer to . This seemingly paradoxical result can also be understood, as follows. The jumps in the data are associated with “merging” transitions, where many strips merge into fewer. Eventually, the system settles into a single strip (for each species) configuration. In this limit, it is natural to regard the holes as “particles” and particles as “holes.” The idea is that only some of these “particles” are bound to interfaces which separate the - and -rich domains. Roughly speaking, of them remain mobile and are driven through the large domains of “holes,” interacting in a minimal way with the ones bound to the interfaces. Thus, we can apply to this system, with the density of “particles” being and (mobile) “holes” being . In the large limit, such an expression reduces precisely to . Without showing it explicitly, the agreement between data and , with , is excellent. The reason the current associated with single-stripe configurations lies much closer to the simple KLS prediction is that in this case the density of mobile vacancies is much nearer .

In Fig. 5, there is a range of densities (approximately ) in which both multistripe and single-stripe configurations exist. Although the former appear to be metastable, precise determination of the associated lifetimes is left for future study. On a lattice of sites, the maximum possible density is . Configurations with this density are frozen, since there are then no mobile vacancies. Outside this limit, however, we do not observe jamming in the sense found in the driven NNE model with NN hopping dynamics drnne ().

The next natural step in the attempt to formulate a theoretical approach is to write a continuum approximation for the discrete MF, while keeping the lowest few orders in the expansion of the displaced densities, e.g., . The final step is to limit our present study to the homogenous phase and consider the lowest orders in the expansion of the densities

 ρA,B(→x,t)=ρ2+ϕA,B(→x,t)+...

In this manner, we obtain expressions for . If we keep only the order linear in , we will find

 ∂t(ϕAϕB)=(D~D~DD)(ϕAϕB) (9)

where and are differential operators. Note that we have taken account of the symmetry under to write this form. Clearly, the sum and difference of the  fields diagonalize this matrix. There should be no confusion if we again use the notation for . In keeping with the discussions above, we will focus on the equation only, though the sum, , will play a role at the nonlinear level. The result of a tedious computation is, apart from the details of the coefficients, to be expected from the fact that is a conserved density. Thus, we find to be the divergence of a current density, obeying the appropriate space-time symmetries:

 ∂tψ=[(E1∂x+E3∂3x...)+(D||∂2x+D⊥∂2y+...)]ψ

For clarity, we have kept the parts odd in (which must be absent when ) separate from those even in (which represent the diffusive parts of this equation). Thus in Fourier space , the former set appears with and can be seen as the real part of a dispersion relation, , which plays the role of driven transport. The imaginary part plays the role of damping, providing the sign is appropriate. From these, we have also allowed for anisotropic diffusion (), as are not necessarily . Of course, the various coefficients are real functions of , and , the only control parameters here. At this level, all are expected to be analytic, so that we expect the ’s and ’s to be odd and even functions, respectively, of the drive . Naturally, the next step is to explore (linear) instability, by checking if any aspect of the damping vanishes. The only such occurrence is vanishing at , independent of , corresponding to the the unmixing transition in the equilibrium WRLG. Needless to say, this instability is of little relevance to our driven lattice gas. But, as a mean-field-like result, it may serve as a starting point for serious analysis.

To gain some insight into the structure factors and correlations, we must add noise terms to form a set of Langevin equations for the density fields, though these terms cannot be derived from the discrete mean-field equations above. Thus we write,

 ∂tψ=[(E1∂x+E3∂3x...)+(D||∂2x+D⊥∂2y+...)]ψ+η(→x,t). (10)

To account for the conservation law and possible anisotropy, we will assume that is delta-correlated, conserved, Gaussian noise, with zero mean and

 ⟨η(→x,t)η(→x′,t′)⟩=−[N||∂2x+N⊥∂2y]δ(→x−→x′)δ(t−t′).

Here again we have allowed for anisotropy, by writing two different coefficients. As a reminder, in the undriven limit, our dynamics satisfy detailed balance, which implies . Following standard routes pioneered by Martin-Siggia-Rose MSR (), Janssen HKJ (), and de Dominicis CdD (), we write the dynamic functional for both the original density fields and the associated response fields . At the quadratic level, corresponding to the linear Langevin example above, the and sectors decouple. The propagators and correlators (e.g., in the sector, the and terms) in Fourier space are, respectively, and . The most general terms consistent with symmetry, anisotropy, and analyticity are just , and . Thus, at this tree level, the equal-time correlation can be found easily; its transform reads,

 Stree(→q)=N||q2x+N⊥q2yD||q2x+D⊥q2y+O(q4)

While this approach may appear to be reasonable for describing the fluctuations of the homogeneous state, clearly, it fails to predict the most prominent feature: presence of a term and the associated maximum at . We are exploring non-linear terms in the Langevin equation (beyond quadratic in ) and whether a simple perturbative approach can generate such a kink singularity. We conjecture that at one-loop level, a term will emerge in the denominator, thereby justifying the form in Eq. (5). Work is in progress to examine a systematic approach, applying the Doi-Peliti formalism Doi (); Peliti (); Pert (); DPrev () to the defining microscopic dynamics to derive (instead of postulating) the appropriate functional. It is our hope that results from such an analysis will provide a better description of the phenomenological form above. The success of stochastic field theory in the past is, of course, the systematic study of universal properties in critical phenomena. Assuming that the transition from the homogeneous to the inhomogeneous state(s) is continuous, we believe this line of pursuit will lead to novel fixed point(s) and universality class(es).

## Iv Inhomogeneous states: Phase segregation, multiple stripes, correlations, and interfacial properties

In this section we report simulation results on the ordered phase of the driven WRLG. In equilibrium, for , the system exhibits separation into A- and B-rich phases above a certain critical density. We therefore begin by looking for signs of phase separation under a drive.

### iv.1 Emergence of striped phases under a drive

As shown in the preceding section, under a drive, the disordered system exhibits a preferred wavelength for charge density oscillations. We find that, for densities near the onset of phase separation, a striped pattern of wavelength appears, as is evident in configuration snapshots as well as quantitative measures such as the structure factor. With increasing density, first decreases, attaining a minimum for , and then grows, approaching the system size , as . (Recall that at least vacancies are required to satisfy the prohibition of NN A-B pairs.) We expect that in the limit , with fixed and , the wavelength converges to a size-independent function of the density and the drive parameters.

In equilibrium, of course, maximum entropy implies the minimum possible number of interfaces, so that is formally infinite for zero drive. (There is, of course, no striped pattern in equilibrium.) Consistent with this, grows systematically as the drive . For , 1/2, 1/3, 1/4 and 1/8, we find = 10, 15, 18, 22 and 36 respectively; these data follow where and are constants.

In contrast to the KLS model, stripes of distinct phases orient perpendicular to the drive, and have rough interfaces. For a wide range of densities, single-stripe configurations (i.e., ) are unstable, and the steady state exhibits multiple stripes; for the system sizes investigated, single-stripe configurations do appear to be stable for higher densities (specifically, for when and , and for , for and ). These observations are illustrated for and , in Figs. 6a and 7. For the densities shown, A- and B-rich regions are evident, organized into stripes oriented perpendicular to the drive. At the highest density shown, , each phase occupies a single stripe, a situation that persists for systems of size . The configurations shown in Figs. 6-8 were obtained using single-stripe initial conditions (perpendicular to the drive), thus demonstrating instability toward formation of multiple stripes for .

We note that, for and a nonzero drive, initial configurations with A- and B-regions separated by a boundary parallel to the drive rapidly reorganize so that resulting the stripe or stripes are perpendicular to it. After a short time ( MC steps) the initially flat interface develops waves which grow steadily in amplitude. When they reach a height , they reconnect via the periodic boundaries, forming stripes perpendicular to the drive. The latter are initially quite irregular, but (for sufficiently high densities) rapidly become ordered, so that the system is separated into clear bulk and interfacial regions.

While experience with the KLS model might lead one to expect interfaces parallel to the drive, it appears that without attractive interactions there is no mechanism for stabilizing and smoothing such an interface in the driven WRLG. On the other hand, sequences of the form A0B (where 0 represents a vacant site) along the drive represent barriers, behind which particles may form a queue (and analogously for B0A sequences). Formation of queues (together with the prohibition of A-B NN pairs) may stabilize stripes perpendicular to the drive. Why many stripes of the same phase form, instead of just one, is unclear. Figure 6c provides a hint of the mechanism: a fluctuation in the boundary has allowed B particles to invade an A stripe, leading to formation of a new stripe. Multi-stripe configurations may contain defects, as is evident in Figs. 6a and 8.

#### The preferred wavelength

As noted above, for a nonzero drive, the structure factor shows a distinct maximum for a certain wavenumber, , corresponding to the preferred wavelength . Power spectra of the columnwise charge density are particularly useful in determining (see Fig. 9). The latter is a decreasing function of drive at fixed density, as shown in Fig. 10. Since the density of vacancies is complementary to the particle density, and since each stripe requires at least vacancies for the interfaces, the number of stripes should decrease (and increase) as the density approaches unity, as verified in Fig. 11. This figure also shows that the dependence of on for different drive strengths is qualitatively similar, and that the dependence of on system size is weak.

Our results for derive from initial configurations (ICs) using both single and multiple stripes. For densities in the range , the different ICs yield consistent results for the stationary value of . For , however, the steady-state wavelength tends to remain the same as the initial value, , for the duration of the simulations, over a wide range of , hampering a precise determination of . For and , for example, we find for initial wavelengths in the range 22-26. For and the range of stable values broadens to 30-50. (While it remains possible that the uncertainty in would be smaller, using longer studies, it appears that the evolution of the number of stripes is very slow at densities .)

### iv.2 Local ordering

At intermediate densities, the driven system exhibits local ordering into A- and B-rich stripes, without global ordering of the pattern; an example is shown in Fig. 8. To quantify local ordering, we consider a function that essentially projects regions of size onto a patch of an ideal pattern. Consider a rectangle of sites, i.e., . Let be the excess number of A particles over B particles in the left half () and the excess number in the right half; define

 ψR=Q1−Q2λℓ. (11)

This quantity takes nonzero values when the region is centered on a stripe of width , but will be close to zero in random configurations, or in regions occupied by a single species. Thus a convenient measure of local ordering is , where the first average is over configurations in the stationary state at a given density, and the second is over configurations (of the same density) generated by inserting A and B particles, with equal probabilities, at random into , subject to the prohibition against A-B nearest-neighbor pairs. The probability density changes from unimodal to bimodal at a certain density, marking the growth of local order. In Fig. 12, for , the transition from unimodal to bimodal occurs near . This does not imply, of course, that there is no local ordering below this density. The behavior of is quite smooth over the range of densities of interest (see right inset of Fig. 12), suggesting that there is no phase transition associated with local ordering. For weaker drives, the transition from a unimodal to a bimodal distribution occurs at a similar density: for , and for . Thus, under a drive, short-range order appears at a density slightly above the equilibrium critical density, .

### iv.3 Order parameter and phase boundary

Given the evidence of phase separation discussed above, in the form of A- and B-rich stripes oriented perpendicular to the drive, we turn to a more quantitative discussion, which requires definition of an order parameter. The ordered phase is characterized by charge-density oscillations along the drive, i.e., a correlation function and an associated structure factor with a maximum (at ) proportional to . We therefore define the order parameter . The order parameter is plotted versus density in Fig. 13 for and ; the plot strongly suggests that there is a continuous phase transition at a density near 0.72.

To obtain a more precise estimate of the density marking the onset of global order, we perform a finite-size analysis of the order parameter . For fixed and (and ), we determine for a series of system sizes . In the disordered phase, we expect to decrease rapidly with , whereas in the ordered phase it must approach a nonzero limiting value as . A phase transition is marked by slow decay of the order parameter with , typically a power law, . Thus in a plot of versus for fixed density , we interpret upward (downward) curvature as a signal of the ordered (disordered) phase.

Estimating is complicated by the fact that as one varies the density, the preferred wavelength changes, as shown in Fig. 11. We obtain more orderly patterns (higher values of ) using system sizes that are integer multiples of the preferred wavelength. (Note that the latter need not be an integer. For and , for example, the preferred wavelength is 10.63.) A well ordered pattern consists of an integer number of wavelengths, which is only possible for specific system sizes. This leads to an irregular variation of with , as shown in Fig. 14. Small changes in (for example, from 140 to 144) can yield large (and reproducible) changes in . The various points on the graph of versus can nevertheless be bounded from above by a smooth “envelope.” The points falling on or near the envelope represent the most ordered cases; we therefore apply the curvature criterion described above to the envelope. For the data shown in Fig. 14, the values associated with points along the envelope fall in the range 10.50 - 10.77, with a mean of 10.63(5), providing an estimate of the preferred wavelength. We also note that in this case, the envelope is fairly well described by a power law, with an exponent .

We adopt the following procedure to determine as a function of :

1. Perform a survey of L values, to estimate the preferred wavelength .

2. Determine for a series of system sizes such that there is an integer with falling near .

3. If necessary, refine the estimate for the preferred wavelength and return to 2.

Varying the density from 0.71 to 0.78 (for fixed and ), we find the set of envelopes shown in Fig. 15. For the envelopes curve downward, while for there is no clear curvature; for the data appear to curve upward. Thus the onset of global order occurs at a density close to 0.77. For densities in the range 0.74-0.76, the order parameter appears to decay as a power law of the system size, with an exponent varying between about 0.173 and 0.143. For a small drive (, ), a similar analysis shows that the transition occurs at . Determining the precise value of as a function of the drive parameters, and verifying the existence of a “critical phase” in which the order parameter scales as an inverse power of , will require extensive studies of larger systems, a task we defer to future work.

### iv.4 Correlation functions

Along with the local and global order parameters discussed above, two-point correlation functions afford insight into the organization of the driven system. The charge-charge and density-density correlation functions are plotted in Figs. 16 and 17, respectively, for maximum drive. exhibits the oscillations expected given the stripe pattern. Interestingly, the decay of , which reflects the persistence of stripes in the direction perpendicular to the drive, decays in a manner qualitatively similar to the envelope of the oscillations in . The density-density correlation function along the drive, , exhibits weak oscillations at twice the spatial frequency as the charge-charge correlation, reflecting a reduction in density at the interfaces between stripes.

At the lowest density () shown in the figures, all four correlation functions decay in a manner best described by a stretched exponential, , with exponents ranging from 0.36 (for ) to 0.55 (for and . (For correlations along the drive, the envelope of the oscillations decays as a stretched exponential.) For density 0.75, stretched-exponential decay is again observed, with , except in the case of , which is better fit by a power law, , with . At density 0.8, approaches a nonzero limit of about 0.41 as , marking long-range coherence of the stripe pattern. The other correlation functions decay as power laws, with exponents (for ), 0.62 (for ) and 0.33 (). Thus, correlations decay more slowly at higher density. A nonzero limiting value of provides an alternative for detecting long-range order; we defer a systematic analysis using this criterion to future work.

### iv.5 Interface width

As shown above, the boundary between A- and B-rich phases defines an interface oriented, on average, perpendicular to the drive, but with significant fluctuations. Thus the scaling properties of the interface width with time and system size are of interest. We examine these properties in the single-stripe regime, in which identification of the interface is relatively simple. Given a configuration, we define for each in the direction perpendicular to the drive, and for along the drive, the function,

 s(ℓ,j)≡ℓ∑i=1σi,j. (12)

Since represents the excess of A particles over B particles in the first sites in row , the values of for which takes its maximum and minimum correspond to the interface positions in this row. (If the sites at which the maximum or minimum values occur are not unique, we use the smallest at which the extremum occurs.) Given the set of sites marking the maximum, we construct a random walk representation of the interface, , by setting and, for , taking , where the difference in the square brackets is calculated under periodic boundaries. A second interface is constructed using the positions at which takes its minimum. The width of each interface is given by,

 (13)

where .

We examine the scaling properties of the interface width as a function of time and system size, plotting, in Fig. 18 (main graph), the width versus time for system sizes , 100, 200 and 400, and (inset) the saturation (long-time) width versus system size. These data are for , , and density (single-stripe regime). The data for do not follow simple power laws, and cannot be collapsed onto a single curve by rescaling time and width. It is nevertheless possible to estimate the growth exponent (associated with the behavior ) using the most linear portions of the graphs (specifically, the second through fifth points in the present case), yielding , 0.205(4), 0.223(5) and 0.228(6), for system sizes 50 through 400, respectively. The expected scaling of the saturation width, is satisfied to reasonable approximation; for we find . We note that the exponent values are not very different from those of the Edwards-Wilkinson (EW) class, and ; it is conceivable that the observed deviations reflect finite-size effects. We defer a more precise determination of the interface width, including analyses of multi-stripe systems, to future work.

## V Summary and outlook

We study nonequilibrium stationary properties of a Widom-Rowlinson lattice gas subject to a drive favoring particle hopping along one axis and suppressing hopping to the opposite direction. The stationary properties are surprisingly different from those of the driven lattice gas with attractive interactions. As in the equilibrium WRLG wrlg (), there is a phase transition with segregation between particle species as the density is increased; the critical density for phase separation increases with drive. We find that even in the disordered phase, there is a preferred wavelength evident in the charge-charge correlation function along the drive direction, and that the associated structure factor does not take the usual Ornstein-Zernike form. With increasing density, local ordering into A- and B-rich regions occurs, with the appearance of stripes oriented perpendicular to the drive, but without global coherence. Increasing the density further, we observe a transition to a stripe pattern with long-range order. We characterize the system in terms of local and global order parameters, correlation functions and associated structure factors, and provide preliminary results on interface roughness and the particle current provoked by the drive.

Our study leaves many intriguing questions for future study. Among them, we highlight:

1. Why does phase separation result in stripes of a characteristic width, , instead of two regions, as happens in equilibrium? Simulations suggest that very broad stripes are unstable, but the underlying mechanism is unclear.

2. Why do the stripes form perpendicular to the drive? Again, we observe an instability in an interface oriented parallel to the drive - undulations on such an interface tend to grow - but a quantitative explanation is lacking.

3. Can one develop a hydrodynamic description or time-dependent Ginzburg-Landau theory capable of reproducing the phenomenology observed in simulations? Can one derive a quantitative prediction for the current? We find (Sec. III) that although simple mean-field analyses do not predict the phase diagram, they do yield quite good predictions for the current.

4. What is the nature of the phase transition to long-range order of the stripe pattern? Can it be connected to (equilibrium) transitions in smectic liquid crystals?

5. Do the scaling properties of the interfaces fall in the Edwards-Wilkinson class, or some other known class of growth processes?

6. To what extent does entropy maximization, which is the sole factor determining static properties at equilibrium, hold in the driven system, particularly for a weak drive? Related to this, is the zero-drive limit smooth or singular?

Although these and other questions remain open, our study demonstrates that the driven Widom-Rowlinson lattice gas exhibits surprising behaviors not previously observed in driven systems.

###### Acknowledgements.
This work was supported by CNPq and CAPES, Brazil.

### Footnotes

1. email: dickman@fisica.ufmg.br
2. email: rkpzia@vt.edu

### References

1. B. Schmittmann and R. K. P. Zia, Statistical Mechanics of Driven Diffusive Systems, Vol. 17 of Phase Transitions and Critical Phenomena, edited by C. Domb and J. L. Lebowitz (Academic Press, London, 1995).
2. J. Marro and R. Dickman, Nonequilibrium Phase Transitions in Lattice Models (Cambridge University Press, Cambridge, 1999).
3. S. Sasa and H. Tasaki, J. Stat. Phys. 125, 125 (2006).
4. S. Katz, J. L. Lebowitz, and H. Spohn, Phys. Rev. B28, 1655 (1983); J. Stat. Phys. 34, 497 (1984).
5. C. N. Yang and T. D. Lee, Phys. Rev. 87, 404 (1952).
6. Long range correlations and singular structure factors in connection with liquids in non-equilibrium conditions are known since 1979 Onuki79 (); Onuki81 (), and have been studied recently Kirkpatrick16 (). In the context of driven lattice gases, see e.g., Refs. schmittmann-zia () and KLS ().
7. A. Onuki and K. Kawasaki, Ann Phys, 121 456 (1979).
8. A. Onuki, K. Yamazaki and K. Kawasaki, Ann Phys, 131 217 (1981).
9. T. R. Kirkpatrick, J. M. Ortiz de Zárate, and J. V. Sengers, Phys. Rev. E93, 032117 (2016)
10. R. K. P. Zia, E. L. Praestgaard, and O. G. Mouritsen, Am. J. Phys. 70, 384 (2002)
11. H.K. Janssen, B. Schmittmann, Z. Phys. B: Cond. Mat. 64, 503 (1986) and K.-t. Leung and J. L. Cardy, J. Stat. Phys. 44, 567 (1986); Erratum J. Stat. Phys. 45, 1087 (1986)
12. K.-t. Leung, K. K. Mon, J. L. Vallés and R. K. P. Zia, “Suppression of interface roughness in driven nonequilibrium systems,” Phys. Rev. Lett. 61, 1744 (1988); Phys. Rev. B39, 9312 (1989).
13. R. K. P. Zia and K.-t. Leung, J. Phys. A24, L1399 (1991) and K.-t. Leung; R. K. P. Zia, J. Phys. A26, L737(1993).
14. P. Pradhan, R. Ramsperger, and U. Seifert, Phys. Rev. E84, 041104 (2011).
15. R. Dickman and R. Motai, Phys. Rev. E 89, 032134 (2014).
16. R. Dickman, New J. Phys. 18, 043034 (2016).
17. R. Dickman, Phys. Rev. E64, 016124 (2001).
18. R. Dickman, Phys. Rev. E 90, 062123 (2014).
19. R. Dickman and G. Stell, J. Chem. Phys. 102, 8674 (1995).
20. B. Widom and J. S. Rowlinson, J. Chem. Phys. 52, 1670 (1970).
21. J. S. Rowlinson and B. Widom, Molecular Theory of Capillarity (Clarendon, Oxford, 1982).
22. P. Nielaba and J. L. Lebowitz, Physica A244, 278 (1997).
23. G. Szabó, A. Szolnoki, and T. Antal, Phys. Rev. E 49, 299 (1994).
24. G. Szabó and A. Szolnoki, Phys. Rev. E 53, 2196 (1996).
25. A. Szolnoki and G. Szabó, Phys. Rev. E 65, 047101 (2002).
26. R. B. Potts, Proc. Camb. Philos. Soc. 48, 106 (1952).
27. M. Blume, V. J. Emery and R. B. Griffiths, Phys. Rev. A 4, 1071 (1971).
28. K. E. Bassler and R.K.P. Zia, Phys. Rev. E 49, 5871 (1994); G. Korniss, B. Schmittmann and R. K. P. Zia, J. Stat. Phys. 86, 721 (1997); B. Schmittmann, K. Hwang and R. K. P. Zia, Europhys. Lett. 19, 19 (1992); K. E. Bassler, B. Schmittmann and R. K. P. Zia, Europhys. Lett. 24, 115 (1993); G. Korniss, B. Schmittmann and R. K. P. Zia, Europhys. Lett. 45, 431 (1999); K. E. Bassler and D. A. Browne, Phys. Rev. Lett. 77, 4094 (1996), Phys. Rev. E 55, 5225 (1997).
29. K. Hwang, B. Schmittmann and R. K. P. Zia, Phys. Rev. E 48, 800 (1993).
30. R. K. P. Zia and B. Schmittmann, Z. Phys. B 97, 327 (1995).
31. G. Korniss, B. Schmittmann and R. K. P. Zia, Physica A 239, 111 (1997).
32. E. L. Praestgaard, B. Schmittmann and R. K. P. Zia, Eur. Phys. J. B 18, 675 (2000).
33. R. Dickman and R. Vidigal, J. Stat. Mech. Theory Exp. 2007, P05003 (2007).
34. P. Greulich and A. Schadschneider, Physica A 387, 1972 (2008).
35. Though we use to denote the local density, there should be no confusion with , the overall particle density. They are related by . Meanwhile, there should be no confusion that “charge” here is unrelated to the electric charge, but used as a way to distinguish from particles.
36. This average power spectrum is precisely the intensity of the scattered light when a beam is shined on a sample of our fluctuating lattice gas.
37. Due to both finite effects and the conservation law (which imposes non-vanishing values on at large ), extracting asymptotic behaviors is not straightforward. Thus, we cannot rely on data at distances comparable to
38. P. C. Martin, E. D. Siggia and H. A. Rose, Phys. Rev. A 8, 423 (1978).
39. H. K. Janssen, Z. Phys. B 23, 377 (1976).
40. C. De Dominicis, J. Physique 37, 247 (1976).
41. M. Doi, “Second quantization representation for classical many-particle system,” J. Phys. A 9, 1465 (1976).
42. L. Peliti, “Path integral approach to birth-death processes on a lattice,” J. Physique 46, 1469 (1985).
43. R. Dickman and R. Vidigal, “Path integrals and perturbation theory for stochastic processes,” Braz. J. Phys. 33, 73 (2003).
44. For a recent overview of this approach, see, e.g., U. C. Täuber, Critical dynamics: a field theory approach to equilibrium and non-equilibrium scaling behavior, (Cambridge University Press, Cambridge, 2014).