Critical behavior in active lattice models of motility-induced phase separation

Abstract Lattice models allow for a computationally efficient investigation of motility-induced phase separation (MIPS) compared to off-lattice systems. Simulations are less demanding, and thus, bigger systems can be accessed with higher accuracy and better statistics. In equilibrium, lattice and off-lattice models with comparable interactions belong to the same universality class. Whether concepts of universality also hold for active particles is still a controversial and open question. Here, we examine two recently proposed active lattice systems that undergo MIPS and investigate numerically their critical behavior. In particular, we examine the claim that these systems and MIPS in general belong to the Ising universality class. We also take a more detailed look on the influence and role of rotational diffusion and active velocity in these systems. Graphic Abstract


Introduction
Non-equilibrium active systems composed of selfpropelled particles offer a wide range of interesting behavior and applications [1][2][3]. A fundamental phenomenon is the so-called motility-induced phase separation (MIPS) [4]: At large propulsion speeds and low rotational diffusion, self-propelled particles block each other due to excluded volume and form initial clusters. If the timescale for the rotational diffusion of the particle orientations at the border of such a cluster is larger than the time it takes to enrich the cluster with additional particles, a dynamical instability leading to non-equilibrium phase separation is induced. Although the phase-separated state resembles passive liquid-gas separation with dense domains surrounded by an active gas, no explicit attractive interactions are required. Still, the phase behavior is very similar, with the binodal curve of coexisting densities terminated by a critical point below which the system remains homogeneous for all densities. The question whether or not the behavior close to such a non-equilibrium critical point is universal, and whether it can be attributed to one of the standard universality classes, is not only of fundamental interest but has also stirred up an ongoing controversy [5][6][7][8].
Numerical studies of active Brownian particles (ABPs) [9][10][11][12][13][14][15][16] are a common approach to investigate MIPS in a simple continuous system. These particles are modeled as disks interacting with each other via a purely repula e-mail: thomas.speck@uni-mainz.de (corresponding author) b e-mail: virnau@uni-mainz.de sive Weeks-Chandler-Anderson potential in the framework of an overdamped Langevin equation. In addition, they are propelled with constant speed along their orientation, which is subject to rotational diffusion. For this system, we have determined the location of the critical point in two dimensions and reported critical exponents, which are incompatible with any of the known universality classes [5]. To gain access to the critical point and the critical exponents, we have proposed a novel method to sample subboxes that minimizes the influence of interfaces on density fluctuations. In contrast, subsequent numerical investigations of related but different models have come to a different conclusion, supporting Ising universality in two dimensions for off-lattice active Ornstein-Uhlenbeck particles [8] and a lattice variant of ABPs [7] (and for ABPs in three dimensions [17,18]). Following generic arguments of renormalization [19], however, all these models (in two dimensions) should fall into the same universality class and thus exhibit the same critical exponents. Indeed, in a first renormalization study of an active field theory ("active model B+" [20,21]) it was found that the critical behavior is controlled by the Wilson-Fisher fixed point [6]. What, then, is the reason for the reported differences? Regarding geometry and subboxes, all three numerical works have employed the same method (for details, see Sec. 2.2). One reason could be insufficient statistics, or insufficient range of system sizes leading to a biased estimate of critical exponents. Or, more intriguingly, are there additional features that characterize universality classes in active matter? We stress that MIPS of repulsive particles is a genuine non-equilibrium phenomenon. The effect of self-propulsion on the critical behavior in models that exhibit phase separation already under equilibrium conditions has been studied for Lennard-Jones interactions [22] and a three-dimensional Asakura-Oosawa model [23] driven by a Vicsek-type force [24][25][26] and found to be compatible with the 3d Ising universality class [27].
In this manuscript, we take another step toward a comprehensive understanding of critical behavior in active matter. To this end, we numerically investigate different variants of two-dimensional active lattice gases with excluded volume and dynamics that break detailed balance. While a range of active lattice gas models has been investigated [28][29][30][31][32][33], we focus on two variants that mimic active Brownian particles. In particular, we study two lattice geometries (the square and hexagonal lattice) and two implementations of the dynamics, either treating rotation and translation serial [7] or concurrently [34]. Our numerical results indicate that details of the dynamics have an influence on the critical behavior and question the proposition that MIPS falls into Ising universality.

Critical behavior
Before embarking on the computational study, let us recall some of the properties close to a critical point. We consider systems that undergo phase separation with two coexisting phases having different densities ρ. The two phases are identified with gas (ρ gas ) and liquid (ρ liq ). The (average) order parameter is the difference, m = ρ liq − ρ gas . As we approach the critical point, the gap in m closes and follows a path through the critical point. Hence, we observe whereby τ measures the distance to the critical point (typically the reduced temperature) and β is the corresponding critical exponent. The transition is continuous, with m > 0 for τ > 0 and m = 0 in the homogeneous phase for τ < 0. In addition, both the susceptibility χ and the correlation length ξ diverge at the critical point, defining two more exponents. Of particular importance is Ising universality in equilibrium systems with shortrange interactions and scalar order parameter, for which in two dimensions the exponents can be obtained analytically [35][36][37] Note that these three critical exponents are not independent but obey the hyperscaling relation Arguments why this relation might still be valid for driven active systems are sketched in "Appendix" A. The diverging correlation length ξ implies that the critical behavior is modified in finite systems, where the correlation length is bounded by the system size l. One of the most remarkable successes of computational statistical physics is that the critical behavior can still be extracted from simulations of finite systems [38,39]. To locate the critical point, we turn to Binder's cumulant ratio which becomes independent of l exactly at the critical point. Note that m l = (N l − N l )/l 2 in the lattice gas formulation with N l being the number of particles.
Plotting the ratio Q l as a function of some parameter for different l thus allows-notwithstanding systematic effects as discussed below-to locate the critical point from the intersection of curves. Moreover, the derivative dQ l /dτ | τ =0 = 1/ν yields the inverse of the critical exponent ν. Once we have located the critical point, we can extract β from plotting m l as a function of τ . Finally, we exploit the scaling form χ l = l γ/νχ (l/ξ) for the susceptibility with scaling functionχ that depends on system size through the ratio l/ξ. Plotting the susceptibility (obtained from the fluctuations of the order parameter) as a function of l allows to extract the ratio γ/ν. We thus have access to the three critical exponents ν, γ and β.

Simulations
While the ensemble of choice for simulations of critical behavior in equilibrium is the grand canonical ensemble, for driven active systems breaking detailed balance this route is not available due to the absence of a comprehensive framework in which a chemical potential is defined (although attempts have been made [40,41]). Therefore, we closely follow the method and analysis proposed in Ref. [5]. All simulations were performed in a periodic box with 1:3 geometry. In such elongated boxes, the dense phase nucleates to a slablike struc- In order to measure dense and dilute phase, as well as to avoid interface regions between the two phases, we place two quadratic subboxes of size l above each other right in the center of mass to sample the dense phase (see Fig. 1). Another set of two subboxes is placed shifted away by one half of the simulation box's width from the center of mass to sample the dilute phase. In total, we sample the number of particles N l within the four subboxes of size l in a simulation box of total size 2l × 6l. Note that N l is counted for each of the four subboxes separately, resulting in four measurements for each snapshot. By adjusting the size of the simulation box through the subbox size, we couple the maximum correlation length to l and achieve a clear crossing point of Q l for different l. The susceptibility is evaluated as χ l = (N l − N l ) 2 / N l . Coexistence densities of dense and dilute phase (ρ liq and ρ gas ) are obtained as plateau values of density profiles generated from a simulation box of size 252×84 (corresponding to l = 42) at activities slightly above the tentative critical point.
As reference, in Fig. 2 we show our analysis applied to the 2d Ising model. Since the 2d Ising model is subject to critical slowing down, we get somewhat poorer statistics (see "Appendix" B). Nevertheless, the analysis yields the following critical exponents β 0.113(1), γ/ν 1.751(2), 1/ν 0.97 (3) (6) in good agreement with the analytical values. Errors in this and later sections refer to statistical errors obtained by splitting respective data sets into three parts and calculating the standard error of the mean. While γ/ν agrees exactly, there is a slight underestimation for 1/ν and a noticeable underestimation for β. While the error of 1/ν might be within statistical uncertainties, the deviation observed for β is more systematic: The power-law scaling is clearly only valid very close to the critical point. However, measuring points closer to the critical point than shown in Fig. 2 is challenging. Below we show that the active lattice models follow a power-law behavior over a wider range.

Model description
We first turn to the model studied in Ref. [7] employing a hexagonal lattice, which we will refer to as model I. On a hexagonal lattice each particle has six neighboring sites and six discrete directions it can be orientated toward (Fig. 3). Specifically, each Monte Carlo (MC) step works as follows: 1. A particle is picked at random. 2. A Gaussian distributed random number (with standard deviation σ and zero mean) is drawn and rounded to the nearest integer n. The current orientation of the particle is adjusted by that integer (n = 1 means one step clockwise, n = −1 means one step counterclockwise, n = 2 means two steps clockwise and so on), cf.
Sketch to illustrate serial model I on a hexagonal lattice. a First, the orientation (arrow) of the particle is updated drawing a random number from a Gaussian distribution. b Then, a move along the particle's orientation is attempted with rate w+, diffusion along another direction is attempted with rate wt each. Note that the adjustment of orientation in step 2 is always accepted and a translation does not change the orientation of the particle. Since w + and w t are fixed, the "activity" of the system is solely adjusted via the rotational diffusion, which is defined by the width σ of the Gaussian distribution. A low value for σ corresponds to low rotational diffusion and therefore highly persistent motion [see Fig. 3c]. Note that the probability to keep the current orientation is not linear in σ, especially not around the estimated values for the critical points. It is also important to note that in contrast to model II discussed below, rotation (step 2) and translation (step 4) are always performed in series.

Analysis and results
By closely following the analysis described in Sec. 2.2, we determine the critical point σ cr,I 0.3048 as the average of the cumulant ratio crossings (Fig. 4) for the four largest system sizes under consideration (l = 24,30,36,42). This value is in agreement with the results published in Ref. [7], which has analyzed systems of comparable system sizes. Excluding the two smallest boxes (l = 12, 18) 1 , the cumulant ratios for the bigger boxes cross within a small interval as expected for critical scaling. Figure 5a-c shows results for the order parameter, the susceptibility and the derivative of the cumulant ratio. Fitting power laws yields the following exponents and thus ν 0.97 and γ 1.63. While the agreement with the corresponding 2d Ising values is reasonable for ν (ν Ising = 1.0) and γ (γ Ising = 1.75), the exponent β differs by more than 25% from β Ising = 0.125. This disagreement is also clearly visible in Fig. 5a.
and thus ν 0.98 and γ 1.66. These critical exponents are very similar to the hexagonal case and within numerical uncertainties, indicating that the influence of the underlying lattice is negligible as one would expect.

Model description
The second model is based on a square lattice and has been proposed in Ref. [34]. As illustrated in Fig. 6, there are now six possible moves: either rotation of the particle orientation clockwise or counterclockwise by 90 • with weight w 1 [Fig. 6a], or translation along the orientation with weight w + or any of the three other directions with weight w t [Fig. 6b]. In contrast to model I, the weight w 1 = 0.1 for rotation is now kept constant and we vary w + with w t = 1. Moreover, in each MC step one of the moves is selected according to its weight. Hence, the particle can either rotate or move in one time step, which we term concurrent. Sketch to illustrate model II on a square lattice. a The particle orientation is turned clockwise or counterclockwise by 90 • with rate w1. b A move along the orientation is attempted with rate w+, diffusion into any other direction with rate wt each. Note that only one of these moves is attempted in each time step. The probability for each move is given by the respective rate divided by the sum of all rates Figure 7 shows the crossings of the cumulant ratios Q l for different box lengths l. The crossings start to converge for the bigger boxes with l ≥ 24. Hence, we only take these system sizes into account and determine the critical point to be at w +,cr 4.76. Corresponding results for the critical exponents are displayed in Fig. 8, for which we find

Alternative determination of β
The accuracy of measuring the exponent β from the density difference is limited by the fact that it becomes more and more difficult to reliably estimate this difference as we approach the critical point. For the Ising model (cf. Fig. 2), this has led to a noticeable deviation from the known analytical value. There is an alternative method using the density fluctuations (ρ l − ρ l ) 2 ∼ l −2β /ν at the critical point which is equivalent to (and therefore not independent from) how γ/ν is determined from χ l , see also Ref. [8]. The densities for the subboxes are given by ρ l = N l /l 2 . This analysis is shown in Fig. 9.
For the 2d Ising model we now obtain 2β /ν 0.249(1), which is in excellent agreement with the analytical value. The values for the other models are included in Table 1, but it is clear that they substantially deviate from 2d Ising universality.

Discussion and conclusions
Our results for the critical exponents are summarized in Table 1. We have also added the corresponding values for active Brownian particles as determined in Ref. [5]. For all lattice models studied here, we find values for ν that are in very good agreement with Ising universality (< 3% smaller) and values for γ/ν that are in good agreement (< 5% smaller). These values are in agreement with plots shown in Ref. [7], which concludes that Ising universality holds. This conclusion seems   (7) 2.10 0.391 (6) The first three columns are the estimated values. The fourth column is the hyperscaling relation Eq. (4), which is approximately obeyed by all models. The last column shows results for the alternative determination of β (cf. Sec. 5). Errors refer to statistical errors obtained by splitting respective data sets into three parts and calculating the standard error of the mean questionable when taking the exponent β into account, which deviates substantially. Indeed, the determination of β is technically the most challenging. However, note that the hyperscaling relation Eq. (4) places a strong constraint on the exponents. From γ/ν 1.68 and ν 0.98, we can obtain an estimate for β 0.157 that is in excellent agreement with the numerically estimated values for model I on both lattice geometries, supporting that reduction of γ (and ν) is not a statistical effect but systematic. In our analysis, we cannot completely exclude the possibility that the system sizes under consideration are still too small and larger sizes would lead to a slightly shifted critical point. This could in principle also shift the critical exponents, particularly β/ν, while the influence on β and 1/ν is less pronounced.
The value for β estimated for model II is even larger. However, in this case the hyperscaling relation is only fulfilled approximately, which might indicate that β is too large. We have observed that obtaining "good" crossings of the cumulant ratio in this model is more challenging, which might be because the speed is changed in contrast to the rotational diffusion in model I. Moreover, the distance to the critical point is larger since the determination of ρ liq − ρ gas requires stable liquid slabs with a well defined plateau of the density profile. The alternative method of Sec. 5 yields a smaller β 0. 19. We notice that the ratio γ/ν has become even smaller, moving away from the Ising value. While this seems to indicate an influence of the different dynamic rules on the critical behavior, we cannot rule out that the scaling closer to the critical point changes (but note that the smaller value of γ/ν accommodates a larger β). Even further from Ising universality are off-lattice active Brownian particles, where also the exponent ν now changes substantially from ν = 1 to ν 1.5. Still, the hyperscaling relation is again approximately fulfilled, indicating that the exponents are consistent. We cannot exclude the possibility that corrections to scaling are relevant and modify these exponents in a way that is compatible with scaling relations [43].
Based on our numerical results, we find the general conclusion from Ref. [7] that MIPS belongs to the 2d Ising universality class to be somewhat premature. Our results even cast some serious doubts on the weaker claim that model I exhibits 2d Ising behavior. At this point, we would like to emphasize that the numerical evidence presented in Ref. [7] is based on figures similar to our Figs. 5(b) and (c) in which the slopes for the 2d Ising values were drawn on top of the simulation values suggesting excellent agreement. However, the authors neither provide values for γ or ν, nor did they mention the discrepancy for the exponent β.
Instead, we see mounting evidence that the critical behavior for models exhibiting MIPS is at least to some degree model-dependent. Whether or not there is an underlying Ising universality or any universality at all, and to which extent deviations occur and why still remains an interesting and challenging question for simulations and theory alike.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this arti-cle are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

A Hyperscaling relation
The scaling relation Eq. (4) is typically derived from a free energy following arguments originally developed by Widom [44]. Since active matter is steadily driven away from thermal equilibrium, its behavior is not governed by such a free energy. However, the property of the free energy that is mostly exploited in deriving scaling laws is that of a generating function, and some of the results can be transferred to non-equilibrium systems.
To this end, consider the cumulant generating function for the order parameterm =m(C) summing over all possible configurations C with probability p(C; τ ) depending on the control parameter τ . The auxiliary field h allows to obtain the average m and susceptibility χ as A system with linear extend l and correlation length ξ in d dimensions can be viewed as n (l/ξ) d independent systems, for which the joint probability p(C) = n i=1 p ξ (Ci) becomes a product of probabilities p ξ (C ; τ ) for the configuration in a smaller system with linear extend ξ. Witĥ In the second step, we invoke the usual scaling hypothesis positing a scaling functionφ(x) of the combined argument (4)]. Hence, this scaling relation only requires extensivity and homogeneity of the generating function, two properties that are not restricted to equilibrium.

B Critical slowing down
As mentioned above, the statistical quality of the data obtained for the active lattice models is way better than that obtained for the 2d Ising model, despite comparable simulation run time. We determine the autocorrelation function for N − N in order to estimate effects of critical slowing  A comparable amount of data is used for each of the other simulation points. The same run time is applied to the smaller simulation boxes, therefore more data is available for these systems. The difference in MC steps performed is a result of different computational effort for the particular models. The equilibration time is not included down. The strongest effects can be found for the biggest systems; therefore, Fig. 10 only shows the results for l = 42 (simulation box of size 84 × 252). Each model is evaluated at its respective critical point. The autocorrelation function is calculated for each individual subbox separately and then averaged over all subboxes and all simulation data obtained. The time lag Δt is given in units of 1,000 Monte Carlo (MC) steps. For each MC step, each lattice site is on average picked once. So in total Δt = 1 is in this case equivalent to performing 21,168,000 simulation steps. Table. 2 provides an overview of the total amount of simulation data evaluated for this work. Figure 10 shows that the correlation for the 2d Ising model is quiet persistent, indicating critical slowing down. On the other hand, the correlations decay rather quickly for the active lattice models and are more or less gone after Δt = 100. Consequently, critical slowing down is not an issue for the active lattice models. The main reason for these differences is the acceptance schemes in the models. For the Ising model, the metropolis criterion has to be fulfilled for a spin swap. In comparison, each move to a free lattice site is accepted for the active lattice models and there is no interaction except for hard repulsion.