Universality class of the motility-induced critical point in large scale off-lattice simulations of active particles

We perform large-scale computer simulations of an off-lattice two-dimensional model of active particles undergoing a motility-induced phase separation (MIPS) to investigate the systems critical behaviour close to the critical point of the MIPS curve. By sampling steady-state configurations for large system sizes and performing finite size scaling analysis we provide exhaustive evidence that the critical behaviour of this active system belongs to the Ising universality class. In addition to the scaling observables that are also typical of passive systems, we study the critical behaviour of the kinetic temperature difference between the two active phases. This quantity, which is always zero in equilibrium, displays instead a critical behavior in the active system which is well described by the same exponent of the order parameter in agreement with mean-field theory.

Introduction. One of the pillars of statistical physics is the concept of universality in critical phenomena. In equilibrium systems, close to a second-order phase transition, universality can be ultimately attributed to the divergence of the correlation length of the order parameter. The behavior of this growing length-scale is found to be independent on the microscopic details of the systems but is determined only by few specific features, i.e. the spatial dimensionality and the symmetries of the order parameter as firstly hypothesized by Kadanoff [1]. Depending on such parameters it is possible to trace back the critical behaviour of disparate systems within few groups, called universality classes.
One of the biggest challenge of recent years is to transfer the vast knowledge acquired on universal behaviour of equilibrium systems into active matter physics. Active matter represents a peculiar class of non-equilibrium systems where the elementary units or agents are selfpropelled objects capable to convert energy in systematic movement [2,3]. The interacting agents are often complex biological objects that exhibit self-organized behavior at large scales giving rise to many and diverse living materials [4]. In particular, self-propulsion can trigger a feedback between motility and local density causing an effective attractive interaction between active particles [5]. This attractive force can bring to a phase separation in active fluids, remarkably similar to the gas-liquid coexistence in equilibrium systems, which is called Motility-Induced Phase Separation (MIPS) [6]. Since MIPS is a very general feature of active dynamics, observed independently on the details of the self- * claudio.maggi@roma1.infn.it † nicoletta.gnan@roma1.infn.it propulsion, it might play a role also in biological systems. For instance, it has been recently observed that multicellular aggregates of Myxococcus xanthus might take advantage of density-motility feedback for developing large scale collective behaviors that have been well described by MIPS [7]. In equilibrium physics the standard gas-liquid coexistence ends with a critical point that belongs to the Ising universality class [8,9]. A natural question is whether or not MIPS curve ends in a critical point and if there is a region close to phase separation in which the active critical behaviour that can be traced back to a specific universality class. Effective equilibrium approaches [10,11], and field-theoretic computations [12] have previously pointed to the Ising universality class. Concerning numerical simulations, despite numerous works have addressed the properties of phase-separated MIPS states [13][14][15][16][17][18][19][20], the study of the critical region and the determination of the critical properties still remain challenging and controversial. In particular it has been shown that active Brownian particles in two-dimensions (2d) display some critical exponents deviating considerably from the Ising ones [21]. However recent on-lattice simulations of an active model have shown that critical exponents are in good agreement with the Ising universality class [12] and suggested that simulations done off-lattice have been performed far from the scaling regime due to the small size of the systems. Since it is not a priori clear if both models lie in the same universality class, one must rely on large-scale computer simulations which are often a necessary tool to understand possible discrepancies between off-lattice and on-lattice critical exponents. This has been the the case, for example, for equilibrium Heisenberg fluids [22,23].
In this work, we report results of large-scale offlattice simulations of Active Ornstein-Uhlenbeck par- ticles (AOUPs) at criticality. Performing Finite-Size-Scaling (FSS) analysis, we show that the system's critical exponents agree with the Ising universality class. We also show the kinetic temperature difference between the two phases, emerging at criticality, display a critical behavior which is well described by the exponent of the order parameter in agreement with mean-field theory combined with a small-τ expansion of the AOUPs model. Model and Methods. We consider a system composed of N self-propelled AOUP disks in 2d [24,25]. This model is perhaps the simplest active particle model (due to the linearity of the process producing the "active noise") which has lead to numerous novel theoretical developments [26][27][28]. It has been shown [26] that AOUPs exhibit MIPS for large values of the persistence time of the activity, as also displayed in Fig. 1(a)-(d). The equations of motion of AOUPs reaḋ where we indicate with r i the i-th particle's position, with ψ i the self-propelling force, and with F i = j =i f ij the conservative force acting on the particle. We consider two-body interactions, i.e., f ij = −∇ ri φ(r ij ), with r ij = |r i − r j | and use a simple inverse power-law potential φ(r) = (r/σ) −12 /12 with a cut-off at r = 2.5 σ. Here σ represents the diameter of the particle and is set to 1. In Eq. (2) τ is the persistence time of the active force and η is a standard white noise source, i.e. η α i (t) = 0 and η α i (t)η β j (s) = 2D δ ij δ αβ δ(t − s), where the greek indices indicate Cartesian components. Here D is the diffusivity of the non-interacting particles which is related to the mean squared velocity v by D = v 2 τ . In all simulations we fix v = 1 and increase τ from small to large values to observe the transition as shown by the schematic phase diagram in Fig. 1(e). Eq.s (1),(2) have been integrated numerically using the Euler scheme with a time step ∆t = 10 −3 up to N t = 10 9 time steps for the largest systems which is adequate to observe full relaxation of the density auto-correlation function as shown in the Supplemental Material (SM). In addition we start each simulation from a random initial configuration and we perform up to N e = 10 8 equilibration steps which guarantees that all runs reach the steady-state, even close to the critical point. We perform averages over several independent runs and errors reported represent twice the standard error of the mean (see SM for details).
The active particles move in a rectangular box of size L x × L y with 1:3 ratio (L x = 3 L y ) and periodic boundary conditions. The simulated systems size are N = (7.5, 15, 30, 60) × 10 3 . To avoid spurious effects due to the presence of an interface [21,29] we compute the quantities of interest only in four sub-boxes of size L = L y /2 centered on the dense and diluted phases. These boxes are located at x = L x /2 ± L x /4 and at y = L y /2 ± L y /4. To ensure that these positions coincide with the locations of the dense and diluted phases, as in [21], for each configuration we first find the system center of mass along x (with periodic boundaries [30]), and we shift all particles so that the x-coordinate of the center of mass coincides with x = L x /2 + L x /4 as shown in Fig. 1(d).
We stress that the method also works at the critical point if the two phases are distinguishable enough. This is typically the case when the critical order parameter distribution shows two peaks as in many Ising [31] and Lennard-Jones systems [32]. This technique has been successfully applied to the 2d Ising model [21] and to the lattice active model [12].
This value approximately corresponds to the critical density estimated for the smallest investigated system size, which is found to be ρ c = 0.953(0.037) (errors from fit are given in brackets, see SM for details).
Results. Although finite systems can not develop any diverging correlation length, the finite-size scaling hypothesis allows us to systematically study the critical properties away from the thermodynamic limit [33]. Using the finite-size scaling ansatz, we assume that a generic observable O near the critical point behaves as , where ζ O is the critical exponent associated with the observable O, F O is a universal finite-size scaling function and ω is the power of the (subleading) correction-to-scaling exponent [33]. Here ν is the exponent associated with the divergence of the correlation length ξ as the control parameter is varied across the transition. In our active particle system the relaxation time of the noise τ is the control parameter, therefore we assume ξ ∼ (τ − τ c ) −ν . Using this and ignoring sub-leading corrections we get A particularly interesting observable is the fourth order cumulant of density fluctuations ∆ρ 2 2 / ∆ρ 4 (the Binder parameter [34][35][36]), where brackets indicate averages over configurations and over sub-boxes. The density fluctuations are computed in the four L×L sub-boxes described above, specifically ∆ρ where N b is the number of particles found in one single sub-box. For the Binder parameter we expect ζ O = 0 and thus it should be size-independent at τ = τ c . Exploiting this property we locate τ c = 16.361(0.058) and B = [ ∆ρ 2 2 / ∆ρ 4 ] τ =τc = 0.781(0.017) as the intersection of the data for N = 15 × 10 3 and N = 60 × 10 3 ( Fig. 2(a)). We choose to use these two sizes because N = 60 × 10 3 is the largest simulated size and N = 15 × 10 3 is two times smaller in linear size. The estimated value of B = 0.781(0.017) is lower than the corresponding value found in the triangular lattice gas (B = 0.8321(0.0023), see SM), but it is close to that found in the active lattice model [12] (B ≈ 0.75). Note that τ c = 16.36 is approximately the value at which cumulants of all sizes cross as shown in the inset of Fig. 2(a) where we report a magnification of the main panel in a small τ -interval around τ c (see SM for a systematic study of the crossing points). Fig. 2(b) shows a good data collapse of the cumulant data-points with the Ising exponent ν = 1. A direct way [21] to determine ν is to consider the size dependence of the slope of the cumulants at τ = τ c , this method yields ν = 1.03(0.10) as shown in the SM.
Next we test the scaling of the susceptibility Fig. 2(b)). Fig. 2(c) shows that scaling is very good using the Ising critical exponent γ = 7/4. In the SM we also show that, if we fit directly the size-dependent values of χ at τ c , we obtain γ = 1.84(0.20) which is compatible with the Ising γ. Note also that the χ in Fig. 2(c) does not show the typical peaked shape of the Ising model. This is due to the fact that the χ is obtained here by averaging together the values of N b both in the dense and diluted phase. In the SM we show that the χ (computed in the same way) for the lattice gas display a similar s-shaped curve as a function of the inverse temperature and that it scales with γ = 7/4. Furthermore, in Fig. 2(e) we consider the density difference between the boxes centered in the high and low-density phases (ρ h − ρ l ) (see Fig. 1(d)), which corresponds to the order parameter of the system. We find that this quantity displays a good scaling with an exponent Ising β = 1/8 ( Fig. 2(f)). It is worth to stress that, in the thermodynamic limit, the order parameter would be different from zero only for τ > τ c . It is however expected that, for finite systems, a smooth variation of the order parameter should be found also below τ c and that this should scale with the appropriate exponent. We have also checked that (ρ h − ρ l ), computed as in the active system, scales with β = 1/8 in the case of the 2d equilibrium lattice gas for temperatures above the critical temperature (see SM for details). A direct fit of the size-dependent critical (ρ h − ρ l ) gives β = 0.113(0.055), which is compatible with the Ising β = 0.125. To improve the accuracy of the β estimate we apply a finer technique finding the value of the exponent which optimizes the data collapse (see SM). This yields β = 0.133(0.022).
Next we check if our results are consistent with the Ising exponent η = 1/4 which controls the decay of the static structure factor S(k) near the critical point, i.e. S(k) ∼ k −2+η . To this aim we compute S(k) = A ρ * k ρ k where ρ k is the Fourier transform of the density fluctuations and the normalization factor A is chosen so that S(0) = ∆N 2 / N (N being the fluctuating number of particles in the sub-sytem considered). To avoid the interfaces and focus only onto the bulk phase we compute the S(k) for the particles in the diluted phase considering only those particles having |x − L x /4| < L/2, i.e. all particle in the left sub-boxes in Fig. 1(d). We show the resulting S(k) in Fig. 3: close to criticality S(k) becomes fairly linear in double log scale at low k. The data at low k are well fitted by the power law k −2+η with η = 1/4 (full line) which is appreciably different from the meanfield decay S(k) ∼ k −2 (dashed line). A direct powerlaw fit of these points gives 2 − η = 1.709(0.090) and η = 0.290(0.090) which are compatible with the Ising values 2 − η = 1.75 and η = 0.25.
Up to this point we have discussed quantities which display a critical behavior also in equilibrium fluids. We show now an observable that is zero in equilibrium while it exhibits a singular behavior in the active case. Since in active systems the instantaneous velocities are coupled to positions [26,37], whenever MIPS occurs dense regions of slow particles coexist with dilute regions of fast ones [20]. We then consider the average squared speed of particles in the dense and dilute sub-boxes that we indicate, respectively, with |ṙ| 2 h and |ṙ| 2 l . The quantity ∆T = 1 2 ( |ṙ| 2 l − |ṙ| 2 h ) can be seen as the (effective) kinetic temperature difference between the two phases and its behavior it is shown in Fig. 4(a). As expected ∆T decreases, as the two phases progressively mix upon low- theoretically derive the relation κ = β, within mean-field theory, by using a small-τ approximation of the AOUP model as shown in the SM. Discussion and Conclusions. In this article, we have studied the critical properties of an active system undergoing MIPS in two spatial dimensions. Performing largescale numerical simulations on GPU we have demonstrated that the critical behavior of the system agrees well with the Ising universality class. It is worth to stress the importance of simulating large system sizes: previous studies on off-lattice active models have reported different values of the critical exponents [21]. Although it has been speculated [12] that the limited sizes employed did not allow to observe the scaling regime, it is true that similar sizes have been exploited for the study of critical passive attractive liquids, finding numerical results compatible with the Ising universality class [29]. We instead suspect that for those sizes another correlation length, different from the critical one, may interfere with the scaling behaviour of the active system. A very recent work [38] has shown that the dense phase formed by active particles undergoing MIPS is made of a mosaic of hexatic micro-domains. We find that (see SM for discussion), already at the critical point, the hexatic correlation length is comparable with the size of the sub-boxes employed for the FSS analysis when the system size is small (N = 3750) justifying the choice of larger system sizes. Indeed for a size as small as N = 3750 we find that the crossing point of the Binder cumulant happens at quite low values, although a reasonable scaling is found also for this size (see SM).
Our large-scale simulation results are also consistent with recent works taking into account non-integrable active terms in a field-theoretical framework [39]. These results indicate that, when full MIPS is possible, these extra terms are irrelevant in a renormalization group sense and the system belongs to the Ising universality class. On the other hand, far from criticality, these non-equilibrium contributions could produce significant differences with respect to an equilibrium gas-liquid phase separation [20,[40][41][42][43]. Within this context it would be interesting to understand how one could make the active critical point unstable [39] by altering the microscopic interactions and/or the dynamics.

SUPPLEMENTAL MATERIAL Equilibration time and length of the simulation runs
We show here that the length of the simulation runs is long enough to allow the full relaxation of the density correlation function is the fluctuation of the number of particles in the sub-box. More specifically we analyze separately the relaxation in the dense and diluted phases considering the two sub-boxes on the left and two on the right as shown in Fig 1 of the main text. Figure S1(a) and (b) display C(t) for the largest system investigate (i.e N = 60 × 10 3 ) at different τ values both in the dense and diluted phase: in both cases C(t) decay to zero at all τ .

Averages and error estimation
All quantities appearing in the main text (i.e. the Binder cumulant, the susceptibility and the order parameter) are averaged over 16000 to 36000 configurations and over all sub-boxes. Individual configurations are taken at time intervals of duration ≈ 23 in reduced units, which is approximately equal the largest τ explored. Moreover for the values of τ close to τ c these quantities of interest are averaged over multiple initial random configurations (up to 12 for the largest system sizes).
To estimate the error on these quantities we proceed as follows: we first divide each run in time windows larger than the relaxation time of the density correlation function C(t) (Fig. S1) we than compute the observable in each time-window and compute the standard error of the mean over all time windows. In Fig. 2 of the main text we report the error as twice the standard error of the mean. We have also checked that by computing the average of the Binder cumulant over these time windows (instead that on all configurations) we get almost identical results than those reported in Fig. 2 of the main text.

Estimate of the critical ρ and τ
We roughly estimate the critical density by performing a density scan, at fixed τ , in the proximity of the critical point for the smallest system investigated (i.e. N = 3750). According to Ref. [36] the cumulant should exhibit a maximum at ρ = ρ c when plotted as a function of ρ and at fixed τ = τ c . Fig. S2(a) shows that the Binder parameter indeed displays a maximum if we fix τ = 16. This value has been chosen based on preliminary simulations and it is close to the critical value τ c = 16.36 estimated in the following. Note also that the cumulant varies much less upon changing density in this small interval than upon changing τ on a large interval as in Fig. 2(a) of the main text. This is shown in the inset of Fig. S2(a) plotting the same data of the main panel in the y-range (0.3, 1), i.e. the range of the binder cumulant as a function of τ . To extract ρ c we have fitted with a 2nd order polynomial the six closest points to the maximum in Fig. S2(a) finding ρ c = 0.953(0.037) where the fit error is reported in brackets. We set the value ρ c = 0.95 for all the sizes discussed in the main text neglecting the dependence of ρ c on the size of the system. The critical value τ , indicated by τ c , has been determined by finding the crossing of the cumulants for two sizes N = 15 × 10 3 and N = 60 × 10 3 . These have been chosen because N = 60 × 10 3 is the largest size simulated and N = 15 × 10 3 has a linear size which is two times smaller than the largest one. The procedure followed to find the intersection is illustrated in Fig. S2(b). We linearly interpolate the cumulant curves and find τ c as the x-intersection of the two lines, while the critical cumulant value B = [ ∆ρ 2 2 / ∆ρ 4 ] τ =τc is found from the intersection on the y-axis. We obtain τ c = 16.361(0.058) and B = 0.781(0.017). This B value is lower than the one found for the triangular lattice gas (B = 0.8321(0.0023), see section below) and for the square lattice gas (B ≈ 0.83 extracted from Ref. [21]). However it is close to B ≈ 0.75 which is the critical cumulant value of the active lattice model found in Ref. [12]. Note that previous studies on the Lennard-Jones fluid have also reported a lower value of the critical Binder parameter with respect to the Ising model [29].
To further check the size dependence of τ c and B, we have repeated this procedure also considering other system sizes, that are separated by a factor of 2 in linear size, i.  Fig. S2(c) and (d) we see that the values of both quantities at N = 7500 are closer to those of the largest system size (dashed lines) than to those corresponding to N = 3750. This suggests that, upon increasing the size, τ c (N ) and B(N ) progressively converge to the infinite-system critical values.

Direct estimates of the critical exponents
Here we show how we directly estimate the critical exponents. Following Ref. [21] we first focus on the dependence of the slope of the critical Binder cumulant on L whose scaling with size is controlled by the exponent ν: To use Eq. (S3) we evaluate the derivative of the cumulant by fitting the numerical data with a generalized logistic function of the form Where A 1 , A 2 , A 3 , A 4 , τ 0 and θ are fitting parameters. As shown in Fig. S3(a) this function fits well the data especially around τ c and allows us to estimate the derivative (S3). This derivative is reported in Fig. S3(b) as a function of the system size and it is indeed well fitted by ∼ L 1/ν with ν = 1. In fact a direct fit with a power law gives 1/ν = 0.968(0.096) (the fit error is indicated in brackets) and, by linear error propagation, ν = 1.03(0.10). Contrarily the best fit with the exponent ν = 1.5 proposed in Ref. [21] deviates considerably from the data.
Next we consider the size dependence of the susceptibility χ at τ c (i.e. χ τ =τc ), that should scale as χ τ =τc ∼ L γ/ν . To do this we simply linearly interpolate the χ at τ c at all sizes and report the results in Fig. S3(c). This quantity is also well fitted by the Ising exponent γ/ν = 7/4 = 1.75. A direct fit with a power law yields γ/ν = 1.787(0.090). Using ν found above and propagating also its error we get γ = 1.84(0.20). Also in this case by fixing the values of γ = 2.2 and τ = 1.5 (i.e. γ/ν ≈ 1.47) from Ref. [21] we obtain a worse fit of our data.  S4. (a) Order parameter as a function of size at τ = τc. Solid line is the best fit of data points which gives β/ν = 0.110(0.053) (i.e. β = 0.113(0.055)). Orange dashed line is the best fit fixing β/ν = 0.125 while the green dashed-dotted line is the fit with β/ν = 0.32 taken from Ref. [21]. (b) Collapsed order parameter data-points (colored dots) on the interpolating function (gray thick line) with exponent β/ν = 0.32 from [21], the distance of the points from the interpolating function is large resulting in a large E. (c) Same as (b) but with the fitted parameter β/ν = 0.130 minimizing the E function.
To directly estimate β we interpolate the points of the order parameter (ρ h − ρ l ) at τ c for all sizes and we plot them as a function of L in Fig. S4(a). It is evident that these points are quite noisy and also that β/ν is small. However we are still able to obtain an estimate of β compatible with the Ising value (β = 0.125) when we fit these points with a power law (ρ h − ρ l ) ∼ L −β/ν . This yields β/ν = 0.110(0.053) (full line in Fig. S4(a)) and β = 0.113(0.055) if the error on ν is propagated linearly. Note that this is close to the Ising value as shown by the orange line in Fig. S4(a) and appreciably smaller than the β = 0.45 of Ref. [21] shown by the green line Fig. S4(a).
To obtain a more accurate estimate of β we also implement a finer method which finds the exponent by minimizing the deviation between collapsed data. This type of technique has been applied in the past to extract the critical exponents of various spin models [44,45]. To practically apply this method we fix the values of τ c and ν to the values determined above (τ c = 16.36 and ν = 1.03). We then consider the following error function to be minimized: where L i andτ i = L i 1/ν (τ i − τ c ) are respectively the system size and the scaled control parameter of the i-th data-point, while m(τ i , L i ) = [ρ h − ρ l ] (τ =τi, L=Li) is the order parameter value at L i andτ i (the sum runs over all available data). The function G in Eq. (S5) is the scaling function describing the critical behavior of m whose analytic form is unknown. To circumvent this problem we evaluate the function G by interpolating the values of L β/ν m(τ , L) with a smooth function. We compute G by averaging over windows of fixed size ∆τ which we choose to be 10 times smaller than the overallτ range, i.e. ∆τ = max(|τ i |)/10. In this way a smoothed G can be evaluated at each desired valuẽ τ i . An example of the resulting G is plotted in Fig. S4(b) where we use the parameter β/ν = 0.32 of Ref. [21]. It is clear that, while the resulting G is smooth enough, the simulation data-points do not collapse well on the curve. The value of β/ν which minimizes the E function of Eq. (S5) is β/ν = 0.130(0.018) that results in a good data collapse as shown in Fig. S4(c) and is again compatible with the Ising value β/ν = 0.125. This gives β = 0.133(0.022), by linear error propagation, which is also compatible with the Ising value. Finally we estimate the exponent characterizing the critical behavior of the difference between the particle average squared speed of the dilute and dense phases, i.e. the quantity: ∆T = 1 2 ( |ṙ| 2 l − |ṙ| 2 h ). We use again the collapse optimization method introduced above for estimating β (see Eq. (S5)). In Fig. S5(a) we report the ∆T values if we use the exponent κ/ν = 0.32. It is evident that, with this exponent, the data-points do not collapse well on the interpolating function (gray curve in Fig. S5(a)). Minimizing the E-function we obtain κ/ν = 0.118(0.018) (i.e. κ = 0.122(0.022)) which gives a good data collapse as shown in Fig. S5 Derivation of the exponent identity κ = β To derive the exponent identity κ = β we consider an AOU particle in d = 1 and subjected to an external potential Φ(x). It is known [37,46] that for small τ the velocity distribution of this particle is is a zero-centered Gaussian with variance to first order in τ . Here v 2 = D/τ is the free particle mean squared velocity (which is assumed to be kept constant as in our simulations) and Φ (x) = ∂ x 2 Φ(x) is the potential curvature. We now assume that the total potential curvature Φ felt by the probe particle in x is generated by the interactions with other particles: is the second derivative of the pair interaction potential. This can be rewrit- where the integral extends over all space and we have introduced the density fieldρ(x) = i δ(x−x i ). By ignoring density fluctuations (mean-field approximation) we setρ(x) = ρ = const and we obtain Φ (x) = φ 2 ρ, where φ 2 = dx φ (x−x ) the mean potential curvature which is assumed to be positive. By using this in Eq. (S6) and neglecting higher order corrections we get: ẋ 2 = v 2 (1−τ φ 2 ρ). We now consider the difference between the averaged squared speed in the low and high density phases: If we now assume that (ρ h − ρ l ) ∼ (τ − τ c ) β near the critical point we have: with κ = β, which is the relation verified by the simulation data within errors. The derivation of this identity can be easily generalized to higher dimensions leading to the same result. To check whether the susceptibility and the order parameter behave in the same qualitative way in the active system and in equilibrium we consider a lattice gas on a triangular lattice in a rectangular geometry (similar to the one employed for the active system). We simulate systems of three different sizes composed by N = 192, 768 and 3072 sites. These sites are enclosed in a rectangular box of size (0, L x ) × (0, L y ) with L x = a N x and L y = a √ 3N y /2, where a = 1 is the lattice spacing. In all simulations we set N x = 3 N y and N = N x × N y is the total number of sites. Some near-critical configurations lattice gas simulated is shown in Fig. S6(a) and (b). By imposing periodic boundary conditions every site has 6 neighbours and the total lattice gas Hamiltonian (H lg ) is given by where J is the coupling constant set to 1 for convenience and n i is the occupancy of the i-th site which assume the values 0 or 1. The simulations conserve the total occupancy (i.e. i n i = const) by using a Kawasaki-type dynamics [47] in which a site can exchange its occupancy with any other site in the lattice in order to accelerate the approach to equilibrium. After an occupancy switch is proposed a standard Monte Carlo (MC) Metropolis rule is applied and the new configuration is accepted or rejected according to the energy change. All simulation results are obtained at fixed average occupancy i n i /N = 0.5 (i.e. at the critical occupancy), starting from a random configuration (i.e. at infinite temperature). It is possible to show, via the transformation n i = (1 + σ i )/2, that the model (S8) can be mapped onto the Ising model with spin σ i = ±1 on the triangular lattice having critical temperature T c = 4/ ln 3 ≈ 3.641 (for J = 1 and k B = 1) [48]. As a consequence, the T c of the lattice gas model turns out to be T c = 1/(ln 3) ≈ 0.91, i.e. an inverse critical temperature T c −1 ≈ 1.099 while the critical average occupancy is n c = 0.5.
The configurations of the lattice gas are analyzed as described in the main text for the active system. We start by shifting each configuration so that the its center of mass is positioned at x = 3L x /4 as also shown in Fig. S6(a) and (b). Subsequently the quantities of interest are averaged over all four L × L sub-boxes, where L = L y /2. The density ρ in one sub-box is computed as ρ = i n i /L 2 , where the prime indicates the sum runs only on those sites within the sub-box. Using this method we further check the correctness of the critical temperature by showing the Binder parameter and its good scaling with ν = 1 in Fig.s S6(c) and (d). By interpolating and averaging the values of the cumulants at the known value of T c for all sizes we get B = 0.8321(0.0023) which is close to the value of B of the square lattice gas found in Ref. [21]. In the main text we have mentioned that the χ, computed by averaging over all sub-boxes, does not show the typical peaked shape but rather forms a s-shaped curve when plotted as a function of the control parameter. This is the case also for the equilibrium lattice gas as shown in Fig. S6(e). In Fig. S6(f) we also show that this χ scales well with ν = 1 and γ = 7/4. In the main text we have also used the average difference of the density in the high-density phase ρ h and of the low-density phase ρ l as an order parameter. To check if this quantity behaves as expected at criticality also in the equilibrium case we compute ρ h and ρ l as the average density of the two sub-boxes on the right and on the left respectively. The resulting (ρ h − ρ l ) is shown in Fig. S6(b) as a function of T −1 . In Fig. S6(c) we show that we obtain a good data collapse by using the Ising exponents β = 1/8 and ν = 1. These data are clearly compatible with those presented for the off-lattice active system discussed in the main text, thus reinforcing the robustness of the analysis bringing to the Ising universality class in the case of the active system.

System sizes, hexatic order and velocity correlation length
We discuss here the data collapse for the smallest system simulated, i.e. N = 3750 (not included in the main text), which is comparable with the sizes used in a previous investigation on the critical behaviour of an off-lattice active system [21]. In Fig. S7 we report the data for this size for the cumulant, the susceptibility and the order parameter. It is evident that a reasonable data collapse with the exponents calculated above is found also for N = 3750 (see Fig. S7(b),(d) and (f)). However the crossing point of the Binder cumulant for this size seems significantly lower in height and τ than the larger sizes (see also Fig. S2(c) and (d)).
As mentioned in the main text we speculate that this could be due to the presence of another growing (but not diverging) correlation length. In the following we identify and compare two of them: the first related the hexatic order and the second associated to velocity correlations.
A very recent work [38] has shown that the dense phase formed by active particles undergoing MIPS is made of a mosaic of hexatic micro-domains whose size does not diverge. To compare the size of these regions with our smallest system size, near the critical point, we consider the state point τ = 16.5, ρ = 0.95 for N = 3750. In Fig. S8(a) we show a high-resolution density map of one configuration of this system (near τ c , already showing phase separation). This ρ-map is obtained by counting the number of particles in small squared bins of linear size s = 1. To characterize the hexatic order we calculate the parameter ψ 6j = N −1 j k e iθ jk for each particle. Here θ jk is orientation angle of the segment connecting the position of the j-th particle with its k-th (out of N j ) nearest neighbors found with a Voronoi tessellation. To visualize the regions with the same orientation we project ψ 6j onto the direction of the mean orientation N −1 i ψ 6i where the sum runs over all particles in the system. In Fig. S8(b) we show the ψ 6 -projection map obtained by averaging the ψ 6 -projection of the particles found in each small bin (white pixels correspond to empty bins). In Fig. S8(b) it is evident that, in the dense phase, hexatic domains (i.e. regions with the same color) have an extent comparable to the size L of the FSS analysis boxes (we have L ≈ 18 for N = 3750).
Recent works [49,50] have also shown that in active systems the colored noise induces an effective coupling between particles velocities. This effect gives rise to regions of densely packed particles with correlated speed and velocity orientation. We show here that, close to τ c , these regions have a size similar to the one of the hexatic regions. To visualize the extent of these velocity correlations we show in Fig. S8(c) the orientation map of particle velocities. This map is obtained by averaging the projected particle velocity vector on the x-axis, i.e. cos(ϑ j ) (where ϑ j is the orientation angle of the j-th particle velocity). Fig. S8(c) shows that the "islands" of velocity-correlated particles have a size comparable with the size of hexatic regions. Note however that when we consider a larger system (N = 60 × 10 3 and L ≈ 72) at the same τ and ρ that is phase separating (Fig. S8(d)) the extension of these correlated hexatic and velocity regions does not scale up but remains approximately of the same size (see Fig. S8(e) and (f)). To quantify this more precisely we compute the correlation function of the hexatic order parameter g 6 (r) = ψ * 6j ψ 6k |r k −rj |=r / |ψ 6j | 2 and the correlation function of the velocity orientation vector gv(r) = v j ·v k |r k −rj |=r . These functions are computed and reported in Fig. S9 considering only particles in the dense phase of the largest system. We find that both g 6 and gv decay to zero in an exponential-like fashion as shown in the double-log inset Fig. S9. We assume that both correlators are well described by a Ornstein-Zernike form in q-space (i.e. g(q) ∼ (ξ −2 + q 2 ) −1 ) and therefore we fit both data-sets with a function of the form g(r) = A K 0 (r/ξ) + B where K 0 is the modified Bessel function of the second kind ξ is the correlation length and A and B are amplitude and shift factors. The fit is quite good and reveals (in agreement with the qualitative map analysis discussed above) that the typical correlation lengths of the hexatic domains and velocity-oriented domains are, respectively, ξ 6 = 9.9(1.2) and ξv FIG. S9. Spatial correlation function of the hexatic order parameter ψ6 (blue points) and of the velocity orientation vector (orange points, see legend) for particles in the dense phase for the system with (N = 60 × 10 3 , τ = 16.5 and ρ = 0.95). The full lines are fits with the K0(r/ξ) Bessel function. The inset is the same of the main panel but on a double-log scale.