Driven WidomRowlinson lattice gas
Abstract
In the WidomRowlinson lattice gas, two particle species (A, B) diffuse freely via particlehole exchange, subject to both onsite exclusion and prohibition of AB nearestneighbor 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 Brich domains separated by holerich 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 KatzLebowitzSpohn (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 nonvanishing 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 (holerich) 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 densitydrive 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 farfromequilibrium thermodynamics and statistical mechanics, in particular in the study of phase transitions schmittmannzia (); marro () and steadystate thermodynamics (SST) ST (). In this context, most studies have focussed on the lattice gas with attractive nearestneighbor (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 schmittmannzia (); marro (), such as generic longrange correlations LRC (); Onuki79 (); Onuki81 (); Kirkpatrick16 (), ‘negative response’ ZiaPM (), nonIsing 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 WidomRowlinson (WR) model or WidomRowlinson 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 discretespace version of the original (continuousspace) WidomRowlinson model WR70 (); RW82 (). Each site may be vacant, or occupied by a particle belonging to species A or species B. Nearestneighbor AB pairs are prohibited, as is multiple occupancy of a site. (The latter restriction departs from the original, continuousspace 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, lowdensity phase and coexisting A and Brich 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 Isinglike critical point, as expected on grounds of symmetry. Subsequent studies of the WRLG have focussed on multispecies versions NielabaLebowitz (). 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 sublatticeordered phases are equivalent under a unit translation in the or direction, whereas the symmetry connecting the A and Brich phases in the WRLG is global exchange of A and B particles. In the driven NNE lattice gas with first and secondneighbor hopping (the model studied by Szolnoki and Szabo, szolnokiszabo02 ()), 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 singlecomponent 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 nearestneighbor 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 lowdensity regions, there is a continuous, Isinglike transition to a phase with sublattice ordering, as in equilibrium szolnokiszabo02 (). In the present work, the dynamics includes both NN and secondneighbor 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 spin1 or 3state system, a class that includes well known systems such as the Potts Potts () and BlumeEmeryGriffiths BEG () models. Some driven 3state 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 meanfield 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 nonvanishing 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, OrnsteinZernike 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 densitydrive plane, as we report the emergence of multiple stripes of A and Brich phases [of wavelength ], oriented perpendicular to the drive! By contrast, in the KLS model, the particlerich (or holerich) regions form a single stripe parallel to the drive (especially in systems with aspect ratios). As in the KLS model, longlived 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 lowdensity, disordered phase, and discuss meanfield and fieldtheoretic 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 phaseseparated, inhomogeneous system is higher than that of a homogenous one. As in the freezing of the hardsphere fluid, a higher entropy is associated with a nominally more ordered phase because the disordered phase possesses relatively few configurations, due to excludedvolume 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 singleparticle hopping, so that for to be nonzero, configurations and must differ by the exchange of a particlehole pair. Specifically, each particle attempts to hop to a first or second neighbor, that is, a displacement of with , (excluding the nohop case ). Any hopping move that results in a configuration satisfying the prohibition against NN AB pairs is accepted. The attempt rates are parameterized in terms of (propensity to move along ) and (tendency to leave unchanged):
(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 nextnearestneighbor as well as nearestneighbor hopping to avoid the possibility of nongeneric patterns or dynamic arrest associated with a restriction to NN hopping szolnokiszabo02 (), 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 1030 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 lowdensity 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 equaltime, twopoint 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
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 fieldtheoretic 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,
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
(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 susceptibilitylike quantity for detecting a secondorder critical point. Instead, our assumes its maximum value at some nonvanishing . 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 secondorder transitions in equilibrium statistical systems), then the large behavior of can provide us with a critical exponent.
In general, the twopoint equaltime 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.,
(3) 
and
(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 (); schmittmannzia (), . 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 equaltime twopoint 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
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
(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
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 nonvanishing , 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 powerlaw 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 (), finitesize 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 nonconserved 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 wavevector over a sizable region of the  plane.
iii.2 Discrete meanfield 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):
(6)  
(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 meanfield (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 nearestneighbor twosite 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 currentdensity relationship in the KLS model. The result here is
(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 singlestripe 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 singlestripe 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
In this manner, we obtain expressions for . If we keep only the order linear in , we will find
(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 spacetime symmetries:
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 meanfieldlike 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 meanfield equations above. Thus we write,
(10) 
To account for the conservation law and possible anisotropy, we will assume that is deltacorrelated, conserved, Gaussian noise, with zero mean and
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 MartinSiggiaRose 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 equaltime correlation can be found easily; its transform reads,
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 nonlinear terms in the Langevin equation (beyond quadratic in ) and whether a simple perturbative approach can generate such a kink singularity. We conjecture that at oneloop 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 DoiPeliti 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 Brich 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 AB pairs.) We expect that in the limit , with fixed and , the wavelength converges to a sizeindependent 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, singlestripe configurations (i.e., ) are unstable, and the steady state exhibits multiple stripes; for the system sizes investigated, singlestripe 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 Brich 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. 68 were obtained using singlestripe 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 Bregions 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 AB 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. Multistripe 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 steadystate 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 2226. For and the range of stable values broadens to 3050. (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 Brich 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
(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 AB nearestneighbor 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, shortrange 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 Brich 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 chargedensity 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 finitesize 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 :

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

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

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.740.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, twopoint correlation functions afford insight into the organization of the driven system. The chargecharge and densitydensity 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 densitydensity correlation function along the drive, , exhibits weak oscillations at twice the spatial frequency as the chargecharge 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, stretchedexponential 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 longrange 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 longrange order; we defer a systematic analysis using this criterion to future work.
iv.5 Interface width
As shown above, the boundary between A and Brich 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 singlestripe 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,
(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 (longtime) width versus system size. These data are for , , and density (singlestripe 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 EdwardsWilkinson (EW) class, and ; it is conceivable that the observed deviations reflect finitesize effects. We defer a more precise determination of the interface width, including analyses of multistripe systems, to future work.
V Summary and outlook
We study nonequilibrium stationary properties of a WidomRowlinson 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 chargecharge correlation function along the drive direction, and that the associated structure factor does not take the usual OrnsteinZernike form. With increasing density, local ordering into A and Brich 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 longrange 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 timedependent GinzburgLandau 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 meanfield 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 longrange 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 EdwardsWilkinson 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 zerodrive limit smooth or singular?
Although these and other questions remain open, our study demonstrates that the driven WidomRowlinson lattice gas exhibits surprising behaviors not previously observed in driven systems.
Acknowledgements.
This work was supported by CNPq and CAPES, Brazil.Footnotes
 email: dickman@fisica.ufmg.br
 email: rkpzia@vt.edu
References
 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).
 J. Marro and R. Dickman, Nonequilibrium Phase Transitions in Lattice Models (Cambridge University Press, Cambridge, 1999).
 S. Sasa and H. Tasaki, J. Stat. Phys. 125, 125 (2006).
 S. Katz, J. L. Lebowitz, and H. Spohn, Phys. Rev. B28, 1655 (1983); J. Stat. Phys. 34, 497 (1984).
 C. N. Yang and T. D. Lee, Phys. Rev. 87, 404 (1952).
 Long range correlations and singular structure factors in connection with liquids in nonequilibrium conditions are known since 1979 Onuki79 (); Onuki81 (), and have been studied recently Kirkpatrick16 (). In the context of driven lattice gases, see e.g., Refs. schmittmannzia () and KLS ().
 A. Onuki and K. Kawasaki, Ann Phys, 121 456 (1979).
 A. Onuki, K. Yamazaki and K. Kawasaki, Ann Phys, 131 217 (1981).
 T. R. Kirkpatrick, J. M. Ortiz de Zárate, and J. V. Sengers, Phys. Rev. E93, 032117 (2016)
 R. K. P. Zia, E. L. Praestgaard, and O. G. Mouritsen, Am. J. Phys. 70, 384 (2002)
 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)
 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).
 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).
 P. Pradhan, R. Ramsperger, and U. Seifert, Phys. Rev. E84, 041104 (2011).
 R. Dickman and R. Motai, Phys. Rev. E 89, 032134 (2014).
 R. Dickman, New J. Phys. 18, 043034 (2016).
 R. Dickman, Phys. Rev. E64, 016124 (2001).
 R. Dickman, Phys. Rev. E 90, 062123 (2014).
 R. Dickman and G. Stell, J. Chem. Phys. 102, 8674 (1995).
 B. Widom and J. S. Rowlinson, J. Chem. Phys. 52, 1670 (1970).
 J. S. Rowlinson and B. Widom, Molecular Theory of Capillarity (Clarendon, Oxford, 1982).
 P. Nielaba and J. L. Lebowitz, Physica A244, 278 (1997).
 G. Szabó, A. Szolnoki, and T. Antal, Phys. Rev. E 49, 299 (1994).
 G. Szabó and A. Szolnoki, Phys. Rev. E 53, 2196 (1996).
 A. Szolnoki and G. Szabó, Phys. Rev. E 65, 047101 (2002).
 R. B. Potts, Proc. Camb. Philos. Soc. 48, 106 (1952).
 M. Blume, V. J. Emery and R. B. Griffiths, Phys. Rev. A 4, 1071 (1971).
 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).
 K. Hwang, B. Schmittmann and R. K. P. Zia, Phys. Rev. E 48, 800 (1993).
 R. K. P. Zia and B. Schmittmann, Z. Phys. B 97, 327 (1995).
 G. Korniss, B. Schmittmann and R. K. P. Zia, Physica A 239, 111 (1997).
 E. L. Praestgaard, B. Schmittmann and R. K. P. Zia, Eur. Phys. J. B 18, 675 (2000).
 R. Dickman and R. Vidigal, J. Stat. Mech. Theory Exp. 2007, P05003 (2007).
 P. Greulich and A. Schadschneider, Physica A 387, 1972 (2008).
 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.
 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.
 Due to both finite effects and the conservation law (which imposes nonvanishing values on at large ), extracting asymptotic behaviors is not straightforward. Thus, we cannot rely on data at distances comparable to
 P. C. Martin, E. D. Siggia and H. A. Rose, Phys. Rev. A 8, 423 (1978).
 H. K. Janssen, Z. Phys. B 23, 377 (1976).
 C. De Dominicis, J. Physique 37, 247 (1976).
 M. Doi, “Second quantization representation for classical manyparticle system,” J. Phys. A 9, 1465 (1976).
 L. Peliti, “Path integral approach to birthdeath processes on a lattice,” J. Physique 46, 1469 (1985).
 R. Dickman and R. Vidigal, “Path integrals and perturbation theory for stochastic processes,” Braz. J. Phys. 33, 73 (2003).
 For a recent overview of this approach, see, e.g., U. C. Täuber, Critical dynamics: a field theory approach to equilibrium and nonequilibrium scaling behavior, (Cambridge University Press, Cambridge, 2014).