Exact fluctuating hydrodynamics of active lattice gases—typical fluctuations

While the hydrodynamic equations are exact in the L → ∞ limit, any finite value of L will lead to a fluctuating noisy dynamics. Our first result complements the hydrodynamic equations of [35] with Gaussian noise terms which account for the typical fluctuations δρ ∼ L^−1/2 around the most probable profile predicted by the deterministic equations (2) and (3). Adopting a representation in terms of the total particle density and polarization fields ρ = ρ₊ + ρ₋, m = ρ₊ − ρ₋ we find that the fluctuating hydrodynamical equations take the form:

$$\begin{equation}{\partial }_{t}\rho =D{\partial }_{x}^{2}\rho -\lambda {\partial }_{x}[m(1-\rho )]+\frac{1}{\sqrt{L}}\sqrt{D}\enspace {\partial }_{x}{\eta }_{\rho }\end{equation} \tag{ 5 }$$

$$\begin{equation}{\partial }_{t}m=D{\partial }_{x}^{2}m-\lambda {\partial }_{x}[\rho (1-\rho )]-2\gamma m+\frac{1}{\sqrt{L}}\left(\sqrt{D}\enspace {\partial }_{x}{\eta }_{m}+2\sqrt{\gamma }\enspace {\eta }_{K}\right).\end{equation} \tag{ 6 }$$

Here, η_ρ , η_m , and η_K are Gaussian white noises of zero mean and correlations given by

$$\begin{align}\langle {\eta }_{\rho }(x,t){\eta }_{\rho }({x}^{\prime },{t}^{\prime })\rangle & =2\rho (1-\rho )\enspace \delta (x-{x}^{\prime })\delta (t-{t}^{\prime }),\\ \langle {\eta }_{\rho }(x,t){\eta }_{m}({x}^{\prime },{t}^{\prime })\rangle & =2m(1-\rho )\enspace \delta (x-{x}^{\prime })\delta (t-{t}^{\prime }),\end{align} \tag{ 7 }$$

$$\begin{align}\langle {\eta }_{m}(x,t){\eta }_{m}({x}^{\prime },{t}^{\prime })\rangle & =2(\rho -{m}^{2})\delta (x-{x}^{\prime })\delta (t-{t}^{\prime }),\\ \langle {\eta }_{K}(x,t){\eta }_{K}({x}^{\prime },{t}^{\prime })\rangle & =\rho \enspace \delta (x-{x}^{\prime })\delta (t-{t}^{\prime }),\end{align} \tag{ 8 }$$

$$\begin{equation}\langle {\eta }_{m}(x,t){\eta }_{K}({x}^{\prime },{t}^{\prime })\rangle =\langle {\eta }_{\rho }(x,t){\eta }_{K}({x}^{\prime },{t}^{\prime })\rangle =0.\end{equation} \tag{ 9 }$$

The angular brackets denote an average over histories of the microscopic dynamics described in section 2. The noise terms are multiplicative, and as is usual in fluctuating hydrodynamics equations, it is straightforward to check that the scaling with L^−1/2 implies that the choice of time discretization, say Itō or Stratonovich, is not important. As we describe in section 4, the equations are obtained through a superposition of known results for the ABC lattice gas [56, 57] and a reaction process accounting for the tumbling. It is straightforward to check that by rescaling x by ℓ and t with the timescale 1/γ (after the diffusive rescaling) the noiseless hydrodynamic part is entirely controlled by the Péclet number $\text{Pe}=\lambda /\sqrt{\gamma D}$ [35]. The noise amplitudes in equations (5) and (6) then become proportional to $1/\sqrt{\ell L}$, indicating that the noise term keeps track of the lengthscale ℓ and the system size L.

At large L, the hydrodynamic equations augmented with the Gaussian noise terms equations (5)–(9) are enough to give the exact two-point functions of the model [24, 58] at leading order in L⁻¹. As detailed in section 5, we find to this order that when the system is homogeneous, namely outside the binodal, the static connected two-point correlation functions are given by:

$$\begin{align}{C}_{2}^{\rho }& \equiv \langle \delta \rho (x)\delta \rho ({x}^{\prime })\rangle =\frac{{\rho }_{0}(1-{\rho }_{0})}{L}\left[\delta (x-{x}^{\prime })-1\right]+\frac{{\text{Pe}}^{2}\enspace {\xi }^{2}{\rho }_{0}^{2}{(1-{\rho }_{0})}^{2}}{2\ell L(1-{\xi }^{2})}\\ & \quad {\times} \left[ \left({\mathrm{e}}^{-\frac{1}{\ell }\vert x-{x}^{\prime }\vert }-\ell \right)+\frac{1-2{\xi }^{2}}{\xi }\left({\mathrm{e}}^{-\frac{1}{\ell }\frac{\vert x-{x}^{\prime }\vert }{\xi }}-\ell \xi \right) \right]\end{align} \tag{ 10 }$$

$$\begin{align}{C}_{2}^{m}& \equiv \langle \delta m(x)\delta m({x}^{\prime })\rangle =\frac{{\rho }_{0}}{L}\left[\delta (x-{x}^{\prime })-1\right]+\frac{{\text{Pe}}^{2}\enspace {\xi }^{2}{\rho }_{0}^{2}(1-{\rho }_{0})(1-2{\rho }_{0})}{2\ell L(1-{\xi }^{2})}\\ & \quad {\times}\left({\mathrm{e}}^{-\frac{1}{\ell }\vert x-{x}^{\prime }\vert }-\frac{1}{\xi }\enspace {\mathrm{e}}^{-\frac{1}{\ell }\frac{\vert x-{x}^{\prime }\vert }{\xi }}\right),\end{align} \tag{ 11 }$$

$$\begin{equation}{C}_{2}^{\rho ,m}\equiv \left\langle \delta \rho (x)\delta m({x}^{\prime })\right\rangle =\mathrm{sign}(x-{x}^{\prime })\frac{\text{Pe}\enspace {\xi }^{2}{\rho }_{0}^{2}(1-{\rho }_{0})}{2\ell L(1-{\xi }^{2})}\left({\mathrm{e}}^{-\frac{1}{\ell }\vert x-{x}^{\prime }\vert }+\frac{1-2{\xi }^{2}}{{\xi }^{2}}\enspace {\mathrm{e}}^{-\frac{1}{\ell }\frac{\vert x-{x}^{\prime }\vert }{\xi }}\right),\end{equation} \tag{ 12 }$$

where ρ₀ is the mean particle density, δρ(x) = ρ(x) − ρ₀, δm(x) = m(x) with the mean-magnetization zero, and the (dimensionless) correlation length ξ is given by

$$\begin{equation}\xi =\frac{1}{\sqrt{2-{\text{Pe}}^{2}(1-{\rho }_{0})(2{\rho }_{0}-1)}}.\end{equation} \tag{ 13 }$$

The expressions are in excellent agreement with the Monte Carlo simulation results of the underlying microscopic dynamics of the lattice model, as illustrated in figures 3(b) and (c) [59]. The results (10)–(12) are valid for ℓ ≪ 1 and we give the expression for arbitrary ℓ in section 5.3. Several comments are in order. First, note that by setting the self-propulsion to zero by taking Pe = 0, we are left only with local correlations as expected—in this limit the model becomes identical to the simple symmetric exclusion model (SSEP) on a ring (see [25] for a review). Second, the integral of these expressions over each of the coordinates x, x' separately, vanishes to leading order in ℓ, as it should, due to total particle number conservation. Finally, when ℓ = 0, so that the diffusive length is microscopic, one is left with local contributions only, encoded in the delta function.

Zoom In Zoom Out Reset image size

When Pe ≠ 0, the model develops large-scale correlations of the order of the system size (recall that length scales are rescaled by the system size L). Indeed, in the hydrodynamic limit, any finite correlation length that does not scale with the system size will vanish and only manifest itself through the Dirac delta function term in the correlation function. It is well known that long-range correlations can appear in lattice gas systems when driven out of equilibrium, for example, by coupling the system to reservoirs of different densities at their boundaries [23]. The large-scale correlations that we observe here are, in contrast, a consequence of the active dynamics of the particles.

Interestingly there are two distinct length scales appearing in the correlation functions. One, ℓ, is independent of the density and stems from the diffusive length of the particles. Recall that here this length scale, which also controls the domain wall width in the phase-separated state, is assumed to be very small. Thus, the terms in the two-point functions associated with this scale are narrowly localized at x = x'. The second one, ξℓ, is related to the collective behavior leading to MIPS. In particular, ξ diverges upon approaching the critical MIPS point (ρ₀ = 3/4, Pe = 4) with a singular behavior of the form

$$\begin{equation}\xi ({\rho }_{0}=3/4,\text{Pe})=\frac{1}{\sqrt{4-\text{Pe}}}+\mathcal{O}\left[{(4-\text{Pe})}^{1/2}\right],\end{equation} \tag{ 14 }$$

corresponding to a mean-field like critical exponent 1/2. Thus, when approaching the critical point, the length scale ξℓ approaches the scale of the system, and the associated terms in the expressions of the correlation functions span the entire system. As detailed below, at this stage our small ℓ approximation fails, and one has to resort to the finite ℓ expressions presented in section 5.3. Note also that the amplitude of the two-point correlation function (10) diverges upon approaching the critical point as ∼ξ. Both these properties are consistent with a behavior expected from an equilibrium scalar field theory in one spatial dimension where only Gaussian fluctuations are taken into account. Indeed, here, the hydrodynamic scaling and the smallness of the noise ensure that, within our approach, only Gaussian fluctuations appear at leading order in 1/L, see also the discussion in section 3.3.

It is also interesting to compare the behavior to that of the ABC model which also presents a singular correlation function when it phase separates (as computed in reference [60]). In contrast to what we find, this singularity is not related to a divergent length scale. This is related to the fact that the quasi-potential of the ABC model is always long-ranged [60], as seen e.g. at equal densities [56, 57]. We return to the latter point below and a detailed comparison between the critical behavior in the present model and the phase transition in the ABC model is presented in section 5.3 for periodic systems.

While all of the correlation functions presented in equations (10)–(12) depend on the two length scales discussed above, their behavior is, as expected, distinct. Both ${C}_{2}^{\rho }$ and ${C}_{2}^{m}$ are even and continuous, but ${C}_{2}^{\rho }$ is monotonic and positive while ${C}_{2}^{m}$ is non-monotonic and changes sign as its argument is varied, see the inset of figure 3(b). Its extrema are located at |x| = x* ≡ 2ℓξ log ξ/(ξ − 1) ∼ ℓ log ξ. On the other hand, ${C}_{2}^{\rho ,m}$ is an odd function and is discontinuous at the origin (see figure 3(c)). These unusual properties can be best understood by rewriting equation (12) in terms of the +, − two-point functions

$$\begin{equation}{C}_{2}^{\rho ,m}=-\langle {\rho }_{+}(x){\rho }_{-}({x}^{\prime })\rangle +\langle {\rho }_{-}(x){\rho }_{+}({x}^{\prime })\rangle .\end{equation} \tag{ 15 }$$

With this decomposition, consider a fluctuation localized at position x with an enhanced + particle density. It is quite clear that it will be positively correlated with an enhanced − particle density fluctuation to the right of x. The reason being that the − particles to the right of x would be impeded by the surge in + particles to their left, causing them to accumulate. However, the effect on the particles to the left of x is exactly the opposite. There, the self-propulsion of the − particles is directed away from the local + particles surge. These features are manifested by the correlation functions between the densities of either positive and negative particles as shown in figure 3(a).

Finally, we note that the Gaussian fluctuating hydrodynamic equations also allow us to determine the dynamical two-point functions. The results, most conveniently written in Fourier space, are given by:

$$\begin{align}\begin{aligned}{~{C}}_{2}^{\rho }& =\langle \delta \rho (k,\omega )\delta \rho ({k}^{\prime },{\omega }^{\prime })\rangle =\frac{2{\rho }_{0}(1-{\rho }_{0})}{\ell LM(k,\omega )}\left[{k}^{2}{\omega }^{2}+\left\{{k}^{2}+2+\frac{2{\xi }^{2}-1}{{\xi }^{2}(2{\rho }_{0}-1)}\right\},{k}^{2}({k}^{2}+2)\right]\\ & \quad {\times}\delta (k+{k}^{\prime })\delta (\omega +{\omega }^{\prime }),\\ {~{C}}_{2}^{m}& =\langle \delta m(k,\omega )\delta m({k}^{\prime },{\omega }^{\prime })\rangle =\frac{2{\rho }_{0}}{\ell LM(k,\omega )}\left[({\xi }^{-2}-2)(1-2{\rho }_{0}){k}^{4}+({k}^{4}+{\omega }^{2})({k}^{2}+2)\right]\\ & \quad {\times}\delta (k+{k}^{\prime })\delta (\omega +{\omega }^{\prime }),\\ {~{C}}_{2}^{\rho ,m}& =\langle \delta \rho (k,\omega )\delta m({k}^{\prime },{\omega }^{\prime })\rangle =\frac{2ik\enspace \text{Pe}\enspace {\rho }_{0}(1-{\rho }_{0})}{\ell LM(k,\omega )}\left[2i{k}^{2}\omega (1-{\rho }_{0})-2{\rho }_{0}{k}^{2}({k}^{2}+2)\right]\\ & \quad {\times}\delta (k+{k}^{\prime })\delta (\omega +{\omega }^{\prime }),\end{aligned}\end{align} \tag{ 16 }$$

where $M(k,\omega )={[{k}^{2}({k}^{2}+{\xi }^{-2})-{\omega }^{2}]}^{2}+4{({k}^{2}+1)}^{2}{\omega }^{2}$. Most notably, by examining scaling properties of the correlation function we find that the dynamical scaling exponent for our model is given by z = 4 on scales smaller than the correlation length and z = 2 on larger scales, in agreement with a mean-field model B and consistent with the mean-field divergence of the static correlation function. Numerical results supporting this are presented below.

As detailed below in section 5 the correlation functions are obtained by linearizing the noisy hydrodynamic equations about a homogeneous profile. The divergence of the amplitude of the fluctuations exhibited by the two-point function (10) as criticality is approached indicates that the assumption of small fluctuations may not be self-consistent near the critical point. To check this we employ a Ginzburg criterion and check if:

$$\begin{equation}\langle \delta \phi {(\xi )}^{2}\rangle \ll 1.\end{equation} \tag{ 17 }$$

Here the average density fluctuation inside a correlation length is δϕ(ξ) ≡ (ξℓ)⁻¹ ∫_ξℓ dx  δρ(x), and we assume that the average density is of the order of 1. We note that this analysis is equivalent to checking whether ⟨[δρ(x) − δρ(0)]²⟩ ≪ 1 holds for any value of x. In the small ℓ limit, and using equation (16) the condition reads

$$\begin{equation}\frac{{\mathrm{P}\mathrm{e}}^{2}\enspace \xi }{\ell }{\rho }_{0}^{2}{(1-{\rho }_{0})}^{2}\ll L.\end{equation} \tag{ 18 }$$

Using the asymptotics (14) and omitting $\mathcal{O}\left(1\right)$ pre-factors, we find that the above condition can be rewritten as

$$\begin{equation}\frac{1}{L\ell }\ll \sqrt{{\text{Pe}}_{\mathrm{c}}-\text{Pe}}.\end{equation} \tag{ 19 }$$

The Ginzburg criterion implies that for any finite ξ, or Pe_c − Pe, a large enough L can be chosen so that it is satisfied. This specifies how the thermodynamic limit has to be taken so that the calculations of the correlation functions, around the homogeneous profiles, are valid. Note that exactly at the critical point the calculation is not self-consistent. While beyond the scope of this manuscript, it is interesting to ask if this is due to a breakdown of the macroscopic hydrodynamic description used here to describe the system or due to the contribution of non-linearities to the fluctuations.

The steady-state probability distribution of observing density and polarization profiles ρ(x) and m(x) takes a large-deviation form at large L as ${\mathrm{e}}^{-L\mathcal{F}[\rho (x),m(x)]}$ where $\mathcal{F}[\rho (x),m(x)]$ is an analog of a free-energy density. Our linearized theory detailed in section 5 allows us to probe the large-deviation function $\mathcal{F}$ up to Gaussian order. In particular, the deviations from the most probable profiles, outside the binodal curve, take the form

$$\begin{equation}\mathcal{F}[\rho (x),m(x)]=\frac{1}{2}\int \mathrm{d}k\left[\begin{matrix}\hfill \delta ~{\rho }(k)\hfill \\ \hfill \delta ~{m}(k)\hfill \end{matrix}\right]{\mathbb{C}}_{k}^{-1}\left[\begin{matrix}\hfill \delta ~{\rho }(-k)\hfill \\ \hfill \delta ~{m}(-k)\hfill \end{matrix}\right]+o\left[{(\delta \rho ,\delta m)}^{2}\right],\end{equation} \tag{ 20 }$$

in Fourier space, where $~{f}$ denotes the Fourier transform of f, ${\mathbb{C}}_{k}$ is a Hermitian matrix with diagonal elements given by (38) and (39), and the off-diagonal elements given by (40) and its complex conjugate [61]. Similarly, the free energy $\mathcal{F}[\rho (x)]$ corresponding to total density fluctuations is obtained by integrating over the magnetization fluctuations:

$$\begin{equation}\mathcal{F}[\rho (x)]=\frac{1}{2}\int \mathrm{d}k\frac{\vert \delta ~{\rho }(k){\vert }^{2}}{{C}_{2}^{\rho }(k)}+\mathcal{O}(\delta {\rho }^{4})\end{equation} \tag{ 21 }$$

with ${C}_{2}^{\rho }(k)$ given by (38).

The long-wavelength physics is captured by the small-k expansion, where we find that $\mathcal{F}[\rho (x)]$ takes the form of the standard Landau theory near the critical point

$$\begin{equation}\mathcal{F}[\rho (x)]\simeq {\left(2\pi \right)}^{2}{\int }_{0}^{1}\mathrm{d}x\left\{\frac{\alpha }{{\xi }^{2}}\delta {\rho }^{2}+\beta {\ell }^{2}{\left(\frac{\mathrm{d}\delta \rho }{\mathrm{d}x}\right)}^{2}\right\},\end{equation} \tag{ 22 }$$

with

$$\begin{equation}\alpha =\frac{1}{2{\rho }_{0}\left(1-{\rho }_{0}\right)\left[2+{\text{Pe}}^{2}\left(1-{\rho }_{0}\right)\right]};\qquad \beta =\frac{{\text{Pe}}^{2}\left(2+{\xi }^{-2}\right)}{2{\left[2+{\text{Pe}}^{2}\left(1-{\rho }_{0}\right)\right]}^{2}}.\end{equation} \tag{ 23 }$$

As expected, the 'mass' term in equation (22) vanishes as the critical point is approached and ξ → ∞. Namely, the form of the large deviation function is that of a Gaussian field theory. Note, however, that if a standard Landau–Ginzburg free energy is rescaled so that x → x/L the coefficient of ${(\frac{\mathrm{d}\delta \rho }{\mathrm{d}x})}^{2}$ would scale as 1/L and vanishes in the large L limit since the cost of a domain wall is finite. In the model we study, in contrast, the cost of a domain wall scales with L for any non-vanishing ℓ.

It is interesting to compare the result to those obtained previously for the ABC model which also exhibits phase separation in one dimension. There the large deviation function has been computed at equal densities [62], close to the equal density point [62], and for arbitrary densities but at Gaussian order [60]. In striking contrast, the large deviation function of the ABC model does not display similar criticality associated with a diverging correlation length (see also section 5.3 for a comparison of the divergent features of the correlation functions between the ABC model and the model studied here).

Next, we detail how the results described above are obtained.

To obtain the conservative noise, we set the tumbling rate of the model to zero. In this case, the diffusive part of the model maps to a driven ABC model [62]. In fact, to compute the covariance of the flux noise terms, it is sufficient to consider the non-driven case where λ = 0 (as generic in MFT, in our weak-drive limit, the noise is independent of the drive [26]). In this limit, the model coincides with an ABC model with symmetric dynamics, first presented in [56] (corresponding to q = 1 in their notations). Within the model, each site can be occupied by either an A, B, or C particle. The dynamics are such that chiral pairs AB, BC, CA exchange with rate q (in arbitrary units) into BA, CB, AC while the reverse exchanges occur with rate q. The model is known to exhibit phase separation for any finite q < 1 [56, 57] and when the rate q is scaled with the system size as q = e^−β/L the phase transition occurs at a finite value of β = β_c [62]. The mapping follows by identifying +, − and vacant states with A, B, and C particle types respectively. The fluctuating hydrodynamics for this model was studied in [60, 64]. The derivation there was based on yet another mapping, with each of the three species A, B, and C behaving as an SSEP, a gas composed of particles interacting only through excluded volume interactions [25]. The corresponding flux noise terms have a covariance:

$$\begin{align}\langle {\eta }_{{\pm}}(x,t){\eta }_{{\pm}}({x}^{\prime },{t}^{\prime })\rangle & =2{\rho }_{{\pm}}(1-{\rho }_{{\pm}})\enspace \delta (x-{x}^{\prime })\delta (t-{t}^{\prime }) ,\\ \langle {\eta }_{+}(x,t){\eta }_{-}({x}^{\prime },{t}^{\prime })\rangle & =-2{\rho }_{+}{\rho }_{-}\delta (x-{x}^{\prime })\delta (t-{t}^{\prime }).\end{align} \tag{ 27 }$$

Transforming to the density and polarization fields ρ and m, this yields the conservative part of the fluctuating hydrodynamics presented in equations (5)–(8). Note that in (27) the same-species diagonal variance terms coincide with those of the SSEP, as one might anticipate. At zero bias λ = 0, the different species are neutral to one another, i.e. the dynamics restricted to that of only one of the species coincides with that of freely evolving gas of SSEP particles. Nevertheless, the cross-correlation between the fluxes of different species is non-vanishing and negative. This stems from the swapping of adjacent particles of different signs. Indeed, such moves contribute to oppositely oriented fluxes of the two particle types which leads to negative correlations. With this, we now turn to discuss the non-conserving part of the dynamics.

The reaction dynamics corresponding to tumbles can be accounted for within the fluctuating hydrodynamics of lattice gas models using several equivalent methods such as the Doi–Peliti formalism [65–67], a large deviation formalism for reactive lattice gas [27, 68–71], or for typical fluctuations using a van Kampen expansion [72]. Here, where the focus is on typical fluctuations, we use a simpler approach.

Specifically, we note that the coarse grained scaling allows us to treat the total tumble term K(x, t) ≡ K₋(x, t) − K₊(x, t) as the difference of two uncorrelated Poisson processes K₊(x, t) and K₋(x, t) which specify the tumbling rates of + and − particles, respectively. According to equation (24), $\mathcal{K}(x,t)\equiv {\int }_{x}^{x+{\Delta}x}L\enspace \mathrm{d}{x}^{\prime }{\int }_{t}^{t+{\Delta}t}\mathrm{d}{t}^{\prime }\enspace K\left({x}^{\prime },{t}^{\prime }\right)$ specifies the total number of particles that tumbled into the + state subtracted by the number of particles that tumbled out of this state, inside a mesoscopic compartment of length Δx during a time interval Δt. Since the densities evolve over macroscopic time and length scales, they can be approximated as constant inside Δx during the time interval Δt if we take Δx such that LΔx ∼ L^δ ≫ 1 with 1 < δ < 0. Hence, the random variable $\mathcal{K}(x,t)$ is simply the difference of two independent Poisson processes with rates γρ₋(x, t)LΔx and γρ₊(x, t)LΔx evolving over a time Δt.

Next, using the fact that the Poisson distribution can be well-approximated by a Gaussian distribution in the large L limit, we can specify the fluctuations by calculating the first two cumulants of $\mathcal{K}(x,t)$. These are easily found to be

$$\begin{equation}\langle \mathcal{K}(x,t)\rangle =\gamma \left[{\rho }_{-}(x,t)-{\rho }_{+}(x,t)\right]L{\Delta}x{\Delta}t\end{equation} \tag{ 28 }$$

$$\begin{equation}\left\langle {\mathcal{K}}^{2}(x,t)\right\rangle ={\langle \mathcal{K}(x,t)\rangle }^{2}+\gamma \left[{\rho }_{-}(x,t)+{\rho }_{+}(x,t)\right]L{\Delta}x{\Delta}t,\end{equation} \tag{ 29 }$$

which then implies that

$$\begin{equation}K\left(x,t\right)=\gamma \left[{\rho }_{-}(x,t)-{\rho }_{+}(x,t)\right]+\sqrt{\frac{\gamma }{L}\left[{\rho }_{-}(x,t)+{\rho }_{+}(x,t)\right]}\enspace \eta (x,t),\end{equation} \tag{ 30 }$$

where $\eta \left(x,t\right)$ is a zero mean and unit variance white Gaussian field. Defining ${\eta }_{K}(x,t)\equiv \sqrt{\rho }\eta (x,t)$ with the variance

$$\begin{equation}\langle {\eta }_{K}(x,t){\eta }_{K}({x}^{\prime },{t}^{\prime })\rangle =\rho \enspace \delta (x-{x}^{\prime })\delta (t-{t}^{\prime }),\end{equation} \tag{ 31 }$$

leads to the announced result (8).

Finally, we remark that this derivation also elucidates the validity of the Gaussian approximation. To see this, we evaluate the third cumulant of $\delta \mathcal{K}(x,t)$ which reads

$$\begin{equation}\langle \delta {\mathcal{K}}^{3}(x,t)\rangle =\gamma \left[{\rho }_{+}(x,t)-{\rho }_{-}(x,t)\right]L{\Delta}x{\Delta}t,\end{equation} \tag{ 32 }$$

and is not accounted for by Gaussian random noise. To incorporate this into K(x, t) of equation (25), we need to add a term of order L^−2/3, which is sub-leading when compared to the Gaussian contribution which is of order L^−1/2. Note, that to account for large but rare fluctuations of K(x, t), which go beyond the Gaussian approximation we present here and are described by large deviations, one needs to account for the non-Gaussian contributions to the statistics of K [73].

To combine both the results presented above we consider a coarse-grained description of the biased diffusion dynamics given by (27) in the presence of the additional tumbling reactions. The former relies on the dynamics obeying 'local equilibrium' described by a local product measure parameterized by the local particle densities [22]. This local equilibrium property is in fact necessary for establishing the more basic deterministic description (2) and (3) and was proven to hold in [34]. That is, as one might expect, introducing a slow tumbling reaction to the ABC dynamics does not violate the local equilibrium property that holds on mesoscopic scales.

Thus, we see that the two dynamics can be superimposed. Since tumbles occur independently of the particle diffusion, the noise terms of the two dynamics are uncorrelated. Under the appropriate change of variables to ρ and m variables, one arrives at equations (5)–(8).

In this section, we derive the expressions for the different correlation functions shown above. Our main focus is the case ℓ ≪ 1 where two distinct phases coexist in the system within the binodal. However, in subsection 5.3 we also consider the case when ℓ is not small. While in this case one cannot identify two distinct phases, the stability of the homogeneous phase can be studied and the correlation functions in it can be evaluated.

In the hydrodynamic limit L ≫ 1, the two-point functions are found by solving for small fluctuations $\rho ={\rho }_{0}+\delta \rho /\sqrt{L}+\cdots \enspace$ and $m=\delta m/\sqrt{L}+\cdots \enspace$, around the steady-state homogeneous profiles. For convenience we rescale the variables so that x → ℓx, t → γt. Then −1/ℓ < x < 1/ℓ and for small ℓ we can use the continuum Fourier transforms

$$\begin{equation}\delta ~{m}\left(k,t\right)=\frac{1}{2\pi }\int \mathrm{d}x\enspace {\text{e}}^{-\text{i}kx}\delta m\left(x,t\right);\qquad \delta ~{\rho }\left(k,t\right)=\frac{1}{2\pi }\int \mathrm{d}x\enspace {\mathrm{e}}^{-\text{i}kx}\delta \rho \left(x,t\right)\end{equation} \tag{ 33 }$$

in equations (5) and (6) with k continuous. We then arrive at the linear system of equations

$$\begin{equation}{\partial }_{t}\left(\begin{matrix}\hfill \delta ~{\rho }(k,t)\hfill \\ \hfill \delta ~{m}(k,t)\hfill \end{matrix}\right)=-{\mathbb{M}}_{k}\left(\begin{matrix}\hfill \delta ~{\rho }(k,t)\hfill \\ \hfill \delta ~{m}(k,t)\hfill \end{matrix}\right)+\mathbb{R}(k,t).\end{equation} \tag{ 34 }$$

Here, the matrix ${\mathbb{M}}_{k}$ and the random noise $\mathbb{R}(k,t)$ are

$$\begin{equation}{\mathbb{M}}_{k}=\left(\begin{matrix}\hfill {k}^{2}\hfill & \hfill iCk\hfill \\ \hfill iEk\hfill & \hfill {k}^{2}+2\hfill \end{matrix}\right);\qquad \mathbb{R}(k,t)=\left(\begin{matrix}\hfill Aik\enspace {~{\eta }}_{1}(k,t)\hfill \\ \hfill \sqrt{{k}^{2}+2}B\enspace {~{\eta }}_{2}(k,t)\hfill \end{matrix}\right)\end{equation} \tag{ 35 }$$

with $A=\sqrt{2{\rho }_{0}\left(1-{\rho }_{0}\right)/(\ell L)}$, $B=\sqrt{2{\rho }_{0}/(\ell L)}$, $C=\text{Pe}\left(1-{\rho }_{0}\right)$ and $E=\text{Pe}\left(1-2{\rho }_{0}\right)$, and where ${~{\eta }}_{1,2}$ are Gaussian white noises with variance

$$\begin{equation}\left\langle {~{\eta }}_{i}(k,{t}_{1}){~{\eta }}_{j}({k}^{\prime },{t}_{2})\right\rangle ={\delta }_{i,j}\delta (k+{k}^{\prime })\delta ({t}_{1}-{t}_{2}).\end{equation} \tag{ 36 }$$

The steady-state correlations can now be readily calculated using the solution of equation (34)

$$\begin{equation}\left(\begin{matrix}\hfill \delta ~{\rho }(k,t)\hfill \\ \hfill \delta ~{m}(k,t)\hfill \end{matrix}\right)={\int }_{-\infty }^{t}\mathrm{d}{t}^{\prime }\enspace {\mathrm{e}}^{-{\mathbb{M}}_{k}(t-{t}^{\prime })}\mathbb{R}(k,{t}^{\prime }).\end{equation} \tag{ 37 }$$

Note that the zeroth mode k = 0 remains identically zero at all times due to the conservation of the total number of particles in the system. This implies that the equal time correlation function vanishes at k = k' = 0. Accounting for this using a delta function whose magnitude is proportional to the finite inter-modes spacing l gives

$$\begin{equation}{C}_{2}^{\rho }(k,{k}^{\prime })=\left\{\frac{{A}^{2}\left[4+(6+CE){k}^{2}+2{k}^{4}\right]+{B}^{2}{C}^{2}(2+{k}^{2})}{4(1+{k}^{2})(2+CE+{k}^{2})}-\frac{2{A}^{2}+{B}^{2}{C}^{2}}{2\left(2+CE\right)}l\delta \left(k\right)\right\}\delta (k+{k}^{\prime }),\end{equation} \tag{ 38 }$$

$$\begin{equation}{C}_{2}^{m}(k,{k}^{\prime })=\left[\frac{{A}^{2}{E}^{2}{k}^{2}+{B}^{2}(2+{k}^{2})(2+CE+2{k}^{2})}{4(1+{k}^{2})(2+CE+{k}^{2})}-\frac{{B}^{2}}{2}l\delta \left(k\right)\right]\delta (k+{k}^{\prime }),\end{equation} \tag{ 39 }$$

$$\begin{equation}{C}_{2}^{\rho ,m}(k,{k}^{\prime })=\frac{ik(2+{k}^{2})({A}^{2}E-{B}^{2}C)}{4(1+{k}^{2})(2+CD+{k}^{2})}\delta (k+{k}^{\prime }).\end{equation} \tag{ 40 }$$

Performing the inverse Fourier transform on equations (38)–(40), one arrives at equations (10)–(12). Notice that the finite ℓ correction introduced at the origin k = k' = 0 via the additional delta functions, results in a constant offset of the real space two point function equations (10) and (11). This $\mathcal{O}\left(1\right)$ offset ensures that the integral over each of the real space coordinates vanishes, to leading order in ℓ, as it should, due to total particle number conservation.

For dynamic correlation functions, it is useful to use the Fourier transforms

$$\begin{align}\begin{aligned}\delta ~{m}(k,\omega )& =\frac{1}{{(2\pi )}^{2}}\int \mathrm{d}x\enspace \mathrm{d}t\enspace {\text{e}}^{-\text{i}(kx-\omega t)}\delta m(x,t);\\ \delta ~{\rho }(k,\omega )& =\frac{1}{{(2\pi )}^{2}}\int \mathrm{d}x\enspace \mathrm{d}t\enspace {\text{e}}^{-\text{i}(kx-\omega t)}\delta \rho (x,t).\end{aligned}\end{align} \tag{ 41 }$$

Inserting these into equations (5) and (6), one obtains at leading order in L:

$$\begin{equation}\left(\begin{matrix}\hfill \delta ~{\rho }\hfill \\ \hfill \delta ~{m}\hfill \end{matrix}\right)={\mathbb{M}}_{k,\omega }^{-1}\mathbb{R}(k,\omega ).\end{equation} \tag{ 42 }$$

with

$$\begin{equation}{\mathbb{M}}_{k,\omega }=\left(\begin{matrix}\hfill -i\omega +{k}^{2}\hfill & \hfill iCk\hfill \\ \hfill iEk\hfill & \hfill -i\omega +{k}^{2}+2\hfill \end{matrix}\right);\qquad \mathbb{R}(k,\omega )=\left(\begin{matrix}\hfill Aik\enspace {~{\eta }}_{1}(k,\omega )\hfill \\ \hfill \sqrt{{k}^{2}+2}B\enspace {~{\eta }}_{2}(k,\omega )\hfill \end{matrix}\right).\end{equation} \tag{ 43 }$$

Here, the random noises satisfy

$$\begin{equation}\left\langle {~{\eta }}_{i}(k,\omega ){~{\eta }}_{j}({k}^{\prime },{\omega }^{\prime })\right\rangle ={\delta }_{ij}\delta (k+{k}^{\prime })\delta (\omega +{\omega }^{\prime }).\end{equation} \tag{ 44 }$$

Solving for $\delta ~{\rho }$ and $\delta ~{m}$ and averaging over noise realizations one finds the two-point correlation functions equation (16). As expected these expressions agree with equations (10)–(12) after performing the inverse Fourier transform and taking the equal-time limit. By rescaling x → bx and t → b^z t and observing the scaling behavior of the dynamic correlation given in equation (16) we identify the dynamic scaling exponent z as follows. As b becomes large, the correlation function converges to the following asymptotic form given as

$$\begin{equation}{~{C}}_{2}^{\rho }\simeq \frac{2{\rho }_{0}(1-{\rho }_{0})}{\ell L}\frac{{k}^{2}{\omega }^{2}+2{b}^{2z}{k}^{2}\left\{2+\frac{1-2{\xi }^{2}}{(1-2{\rho }_{0}){\xi }^{2}}\right\}}{4{\omega }^{2}+{\left[{b}^{z-4}{k}^{2}({k}^{2}+{b}^{2}{\xi }^{-2})-{b}^{-z}{\omega }^{2}\right]}^{2}}{b}^{z-1}\delta (k+{k}^{\prime })\delta (\omega +{\omega }^{\prime }),\end{equation} \tag{ 45 }$$

where the rescaling leads to k → b⁻¹ k and ω → b^−z ω. The correlation function shows different scaling behaviors depending on whether the length scale of interest is larger than or smaller than the correlation length ξ. For the length scales much larger than the correlation length (k ≪ ξ⁻¹), the terms in the square brackets behave as ${[{b}^{z-2}{k}^{2}{\xi }^{-2}]}^{2}$. Taking large b, the correlation function becomes a scaling function of ${~{C}}_{2}^{\rho }\propto {k}^{-2}[4{(\omega /{k}^{2})}^{2}+{\xi }^{-4}]$ only if z = 2, which implies the diffusive scaling. In contrast, if the length scale of interest is much smaller than the correlation length (k ≫ ξ⁻¹), due to divergence of the correlation length nearby the critical point for instance, the terms in the square brackets converges to ${[{b}^{z-4}{k}^{4}]}^{2}$. In this case, the correlation function becomes ${~{C}}_{2}^{\rho }\propto {k}^{-6}[4{(\omega /{k}^{4})}^{2}+1]$ when z = 4.

To check this numerically, it is more convenient to directly study the eigenvalues of ${\mathbb{M}}_{k}$. As equation (34) indicates, the eigenvalues determine the time scale for the relaxation of different modes. According to equation (35), the eigenvalues are given as

$$\begin{equation}{\lambda }_{{\pm}}=({k}^{2}+1)\left[1{\pm}\sqrt{1-\frac{{k}^{2}({k}^{2}+{\xi }^{-2})}{{({k}^{2}+1)}^{2}}}\right].\end{equation} \tag{ 46 }$$

Here, we are interested in the long-wavelength dynamics governed by small-k modes where

$$\begin{align}\begin{aligned}{\lambda }_{+}& =2+\left(2-\frac{1}{2{\xi }^{2}}\right){k}^{2}-\frac{1}{8}{\left(2-\frac{1}{{\xi }^{2}}\right)}^{2}{k}^{4}+\mathcal{O}({k}^{6}),\\ \quad {\lambda }_{-}& =\frac{1}{2{\xi }^{2}}{k}^{2}+\frac{1}{8}{\left(2-\frac{1}{{\xi }^{2}}\right)}^{2}{k}^{4}+\mathcal{O}({k}^{6}).\end{aligned}\end{align} \tag{ 47 }$$

Note that λ₊ → 2 and λ₋ → 0 as k → 0, so that the long-time dynamics of the fields are governed by λ₋ [74]. In figure 4, we plot the relaxation dynamics of the principal mode when the system is initially prepared in a square-wave form with ρ(x) = 1 for 0 ⩽ x < ρ₀ and ρ(x) = 0 for ρ₀ < x < 1, and m(x) = 0. The simulation results, shown in symbols, display good agreement with the analytical relaxation ${\mathrm{e}}^{-{\lambda }_{-}t}$ with λ₋ from equation (47).

Zoom In Zoom Out Reset image size

As the critical point is approached, the correlation length diverges and λ₋ becomes proportional to k⁴ as shown by equation (47), indicating the critical slowing down addressed with the correlation function above. Unfortunately, a direct numerical confirmation of the critical slowing down turned out to be beyond the reach of our numerics. Nearby the critical point, the Ginzburg criterion equation (18) requires L ≫ ξ for the results to hold while the hydrodynamic time scale, required for the simulations, diverges as L². This makes it difficult to capture the long-time behavior of large systems in numerics using our numerical capabilities.

In problems presenting a phase coexistence or a discontinuous phase transition, finite systems can present different features than those in the thermodynamic limit. It is thus of interest to study how finite-size effects manifest themselves in the system studied here. As we have shown, the finite system results also allow us to compare the critical behavior that we report here to that of the ABC model, which was studied in a finite system in [60].

To do this, it is instructive to first discuss the deterministic steady state of the system at finite ℓ, as described by the noiseless hydrodynamics (2) and (3). At finite ℓ the system no longer displays MIPS. Nevertheless, at large enough Péclet numbers, the homogeneous state becomes linearly unstable along a (modified) 'spinodal' line [35] which is found by analyzing the eigenvalues of the linearized dynamics equation (46).

Noting that at finite ℓ the discrete nature of the Fourier modes, k = 2πℓn with n a non-zero integer, has to be accounted for (the zeroth mode k = 0 is stationary due to total particle number conservation) the first mode to become unstable is k₁ = 2πℓ. This happens when λ₋ changes sign at ${k}_{1}^{2}+{\xi }^{-2}{< }0$ (where ξ in equation (13) is purely imaginary), which sets the instability condition

$$\begin{equation}{\text{Pe}}^{2}{ >}\frac{2+{\left(2\pi \ell \right)}^{2}}{(1-{\rho }_{0})(2{\rho }_{0}-1)}.\end{equation} \tag{ 48 }$$

Note that at finite ℓ, the instability sets in at larger Péclet values when compared to the small ℓ limit. In particular, there is no instability at Pe = 4 and density ρ = 3/4. However, as in the small ℓ limit, the homogeneous solution can become metastable before the 'spinodal' curve (48) is reached. While at vanishing ℓ this happens along the binodal curve depicted in figure 2 [35], at finite ℓ at the moment there is no prescription for finding this line in a systematic manner. Therefore, in what follows, we are able to calculate the correlations about the flat state which may be, depending on the region in the phase diagram and outside the binodal, stable or meta-stable. Interestingly, even when the homogeneous state is metastable, for large L, it is long-lived. Since the fluctuations which drive the system out of the metastable state scale as L^−1/2, the mean time for the system to relax to the globally stable state is expected to grow exponentially with L as exp(αL) with α determined by the most probable path connecting the two states. The scaling with L arises from the Gaussian cost of fluctuations, illustrated, for example, in the large-deviation function equation (20).

The derivation proceeds along the lines of section 5 where the continuous Fourier transform is replaced by a Fourier series for the fields f(x) = δρ(x), δm(x) in a system of finite length 1/ℓ:

$$\begin{equation}{~{f}}_{n}=\ell {\int }_{0}^{1/\ell }\mathrm{d}x\enspace f(x){\mathrm{e}}^{2\pi \mathrm{i}n\ell x};\qquad f(x)=\sum\limits _{n=-\infty }^{\infty }{~{f}}_{n}\enspace {\mathrm{e}}^{-2\pi \mathrm{i}n\ell x}.\end{equation} \tag{ 49 }$$

With equations (5) and (6), these give a discrete version of equation (34)

$$\begin{equation}{\partial }_{t}\left(\begin{matrix}\hfill \delta {~{\rho }}_{n}(t)\hfill \\ \hfill \delta {~{m}}_{n}(t)\hfill \end{matrix}\right)=-{\mathbb{M}}_{n}\left(\begin{matrix}\hfill \delta {~{\rho }}_{n}(t)\hfill \\ \hfill \delta {~{m}}_{n}(t)\hfill \end{matrix}\right)+{\mathbb{R}}_{n}(t),\end{equation} \tag{ 50 }$$

where the matrix ${\mathbb{M}}_{n}$ and the noise ${\mathbb{R}}_{n}(t)$ are given by

$$\begin{equation}{\mathbb{M}}_{n}=\left(\begin{matrix}\hfill {k}_{n}^{2}\hfill & \hfill iC{k}_{n}\hfill \\ \hfill iE{k}_{n}\hfill & \hfill {k}_{n}^{2}+2\hfill \end{matrix}\right);\qquad {\mathbb{R}}_{n}(t)=\left(\begin{matrix}\hfill Ai{k}_{n}\enspace {~{\eta }}_{1,n}(t)\hfill \\ \hfill \sqrt{{k}_{n}^{2}+2}B\enspace {~{\eta }}_{2,n}(t)\hfill \end{matrix}\right).\end{equation} \tag{ 51 }$$

with k_n = 2πℓn and $\left\langle {~{\eta }}_{i,{n}_{1}}\left({t}_{1}\right){~{\eta }}_{j,{n}_{2}}\left({t}_{2}\right)\right\rangle =\ell {\delta }_{i,j}{\delta }_{{n}_{1},-{n}_{2}}\delta \left({t}_{1}-{t}_{2}\right)$. Solving the coupled equations and performing the inverse Fourier transform one finds,

$$\begin{equation} {C}_{2}^{\rho }(x,{x}^{\prime })=\frac{{\rho }_{0}(1-{\rho }_{0})}{L}\left[\delta (x-{x}^{\prime })-1\right]+\frac{1}{L}\frac{{\text{Pe}}^{2}\enspace {\xi }^{2}{\rho }_{0}^{2}{(1-{\rho }_{0})}^{2}}{2\ell ({\xi }^{2}-1)}\left[f(\vert x-{x}^{\prime }\vert )+{g}_{\xi }(\vert x-{x}^{\prime }\vert )\right]\end{equation} \tag{ 52 }$$

where we defined

$$\begin{equation}f(\vert x-{x}^{\prime }\vert )\equiv -\left[\frac{\mathrm{cosh}\left(\frac{\vert x-{x}^{\prime }\vert }{\ell }-\frac{1}{2\ell }\right)}{2\enspace \mathrm{sinh}\left(\frac{1}{2\ell }\right)}-\ell \right],\quad \text{and}\end{equation} \tag{ 53 }$$

Here |ξ| is the modulus of ξ, and k₁ = 2πℓ. Note that equation (56) is the analytic continuation of equation (54) obtained from the 0 < ξ⁻² regime to the ξ⁻² < 0 one where ξ becomes purely imaginary. Also, notice that the integral of equation (52) over each of the coordinates x, x' separately, vanishes due to total particle number conservation. We compare the expressions for finite ℓ correlation to simulation results of the microscopic dynamics, all of which show excellent agreement as presented in figure 5.

Zoom In Zoom Out Reset image size

Equation (52) can now be used to describe the finite ℓ rounding of the critical behavior predicted by equation (10). Consider first the limit where we take ℓ → 0 before approaching the critical point $\left({\rho }_{\mathrm{c}}=3/4,\text{Pe}=4\right)$. In this limit, we recover the infinite system result since equation (52) converges to equation (10) and the region for equation (56) shrinks to zero since k₁ → 0. Here once ξ becomes of the order of 1/ℓ, the small ℓ approximation breaks down and one must resort to the finite ℓ expression equation (52) where no singularity is present at Pe = 4 as discussed above. For Pe < 4 where ξ⁻² > 0 expression equation (54) holds. At Pe = 4 and ρ₀ = 3/4 the correlation function takes the form of the parabola equation (55) in contrast to the ℓ → 0 where there is no valid analytical expression. Then for Pe > 4 such that $-{k}_{1}^{2}{< }{\xi }^{-2}{< }0$ expression equation (56) holds. These results are illustrated in figure 5.

Note that beyond Pe = 4 for finite ℓ the correlation length spans the whole system size. To see this it is useful to evaluate the span of the two-point function using the first moment

$$\begin{equation}\langle x\rangle \equiv {\int }_{0}^{1}\mathrm{d}x\enspace {~{C}}_{2}^{\rho }\left(x\right)x,\end{equation} \tag{ 57 }$$

with

$$\begin{equation}{~{C}}_{2}^{\rho }\left(x\right)\equiv \frac{{C}_{2}^{\rho }\left(x\right)-{\mathrm{min}}_{x}\enspace {C}_{2}^{\rho }\left(x\right)}{{\int }_{0}^{1}\mathrm{d}x\left[{C}_{2}^{\rho }\left(x\right)-{\mathrm{min}}_{x}\enspace {C}_{2}^{\rho }\left(x\right)\right]}\end{equation} \tag{ 58 }$$

a normalized and shifted two point function so that it serves as a proper distribution function. This measure grows from a small value which is of order ℓ for Pe < 4 to the order of the system size at Pe = 4 and remains at the order of the system size up to the linear instability point. This is illustrated in figure 6 and compared to the ℓ → 0 case where the span ⟨x⟩ diverges at Pe = 4.

Zoom In Zoom Out Reset image size

Finally, the results above also allow one to compare the phase transition to that found in the version of the ABC model studied in [62]. In the ABC model, each site can be occupied by either an A, B, or C particle and the exchange dynamics lead to phase separation. For average densities r_A,B,C of the three species, and denoting ${\Delta}=1-2({r}_{A}^{2}+{r}_{B}^{2}+{r}_{C}^{2})$, it was shown that a continuous transition into a phase separated state occurs at $\beta ={\beta }_{\mathrm{c}}\equiv 2\pi /\sqrt{{\Delta}}$ [62]. Moreover, in the homogeneous phase the connected two-point correlation function is [60]:

$$\begin{equation}{\langle {\rho }_{A}(x){\rho }_{A}({x}^{\prime })\rangle }_{c}=\frac{1}{L}{r}_{A}(1-{r}_{A})\left(\delta (x-{x}^{\prime })-1\right)-\frac{3}{L}{r}_{A}^{2}{r}_{B}{r}_{C}\left(\frac{\mathrm{cos}\frac{\beta \sqrt{{\Delta}}(1-2x-2{x}^{\prime })}{2}}{\sqrt{{\Delta}}\enspace \mathrm{sin}\frac{\beta \sqrt{{\Delta}}}{2}}-\frac{2}{\beta {\Delta}}\right).\end{equation} \tag{ 59 }$$

Note that at the critical point the amplitude of the non-local contribution (59) to the correlation function diverges due to its denominator, in the same way as is observed in equation (56). However, this phase transition is not associated with a divergent correlation length, in contrast to the critical features of the active model studied in this paper. Indeed, in the ABC model, equation (59) shows that the correlation length spans the entire system also away from criticality, and does not grow as β → β_c, while in our case ξℓ does become formally divergent as Pe → Pe_c (i.e. tends to the full system size). In comparison to the model we study here, in the ABC model the domain wall width, which is known to behave as 1/|ln q| ∝ L/β is always of the order of the system size. In our case, the domain wall width is controlled by ℓ which can be held at any value as the critical point is approached. In particular, in the ℓ → 0 limit the domain wall vanishes. Then our system falls into the more conventional phase separation picture. In contrast, in the ABC model the domain wall is not an independent parameter and takes a large value at the critical point. In sum, in comparison to the ABC model, the model presents more standard features of critical phenomena.