Universal properties of repulsive self-propelled particles and attractive driven particles

Motility-induced phase separation (MIPS) is a nonequilibrium phase separation that has a different origin from equilibrium phase separation induced by attractive interactions. Similarities and differences in collective behaviors between these two types of phase separation have been intensely discussed. Here, to study another kind of similarity between MIPS and attraction-induced phase separation under a nonequilibrium condition, we perform simulations of active Brownian particles with uniaxially anisotropic self-propulsion (uniaxial ABPs) in two dimensions. We find that (i) long-range density correlation appears in the homogeneous state, (ii) anisotropic particle configuration appears in MIPS, where the anisotropy removes the possibility of microphase separation suggested for isotropic ABPs [X.-Q. Shi et al., Phys. Rev. Lett. 125, 168001 (2020)], and (iii) critical phenomena for the anisotropic MIPS presumably belong to the universality class for two-dimensional uniaxial ferromagnets with dipolar long-range interactions. Properties (i)-(iii) are common to the well-studied randomly driven lattice gas (RDLG), which is a particle model that undergoes phase separation by attractive interactions under external driving forces, suggesting that the origin of phase separation is not essential for macroscopic behaviors of uniaxial ABPs and RDLG. Based on the observations in uniaxial ABPs, we construct a coarse-grained Langevin model, which shows properties (i)-(iii) and corroborates the generality of the findings.


I. INTRODUCTION
Liquid-gas or liquid-liquid phase separation is a typical collective phenomenon that has been observed in a wide range of systems from polymer solution [1] to biological materials [2,3].Basically, equilibrium phase separation is caused by attractive interactions between molecules or particles [1], and the corresponding critical phenomena have been considered to belong to the Ising universality class [4][5][6].In contrast, in nonequilibrium systems, depending on how the detailed balance is broken, the critical exponents for phase separation can deviate from the Ising model values [7], and phase separation can emerge from different mechanisms such as chemical reactions [8] and coupling to multiple heat baths [9].A comprehensive understanding of the seemingly broad spectrum of nonequilibrium phase separation requires theoretical studies from a unified viewpoint.
Self-propulsion is another way to break the detailed balance [37][38][39].In a crowd of self-propelled particles, or active matter, collective phenomena ranging from giant number fluctuations [40,41] to active turbulence [42] have been found using biological [43][44][45][46][47] and artificial [48][49][50][51][52][53][54][55][56] systems.In particular, as shown in simulations [57][58][59][60][61] and experiments [62], self-propelled particles with repulsive interactions can undergo phase separation, which is called motilityinduced phase separation (MIPS) [63].No attractive interactions are necessary for MIPS, which is distinct from equilibrium phase separation or nonequilibrium phase separation under external driving.MIPS has been studied in comparison with equilibrium phase separation, and similarities and differences between them have been reported [see Figs.1(a) and (c)].For example, the global phase diagrams for MIPS [64,65] and equilibrium phase separation are similar if we exchange the axis of self-propulsion strength for MIPS with that of attractive interaction strength for equilibrium phase separation.In addition, the lever rule [1], which is common to equilibrium phase separation, holds for MIPS in particle models [66], and consistently, effective free energy has been proposed based on coarse-grained models [67,68].In contrast, it is still unclear whether the critical phenomena for MIPS belong to the Ising universality class [69][70][71][72].Furthermore, as a unique feature of MIPS, the nucleation of persistent gas bubbles that can lead to microphase separation has been found [66,[73][74][75].
In the previous work [76], one of the authors has proposed another kind of similarity between the anisotropic version of MIPS and attraction-induced phase separation under external driving.Briefly, it has been found that a lattice gas model with spatially anisotropic self-propulsion exhibits a variety of col- lective behaviors: long-range density correlation, anisotropic phase separation, and critical phenomena with the universality class expected to be the same as that for uniaxial dipolar ferromagnets.All these behaviors have also been seen in RDLG, which indicates a connection between repulsively interacting particles with anisotropic self-propulsion and attractively interacting particles under external driving.However, the generality of such observations is still unclear beyond the considered lattice gas model.In particular, though persistent gas bubbles have been observed in active Brownian particles (ABPs) [66], a prototypical model of MIPS [59], the fate of gas bubbles under spatial anisotropy has not been investigated.More broadly, systematic studies of the effect of spatial anisotropy on active matter are still scarce [77][78][79][80][81].
In this paper, toward a comprehensive understanding of the relation between the anisotropic MIPS and attractioninduced phase separation under external driving, we consider ABPs with anisotropic self-propulsion.In Fig. 1, we show typical particle configurations obtained from model simulations for the above-mentioned four types of phase separation: attraction/motility-induced phase separation with isotropic/anisotropic dynamics.In each panel of Fig. 1, we also schematically show the single particle motion and typical configuration of small clusters, which can grow up to a macroscopic scale and lead to phase separation.Our present focus is on the relation between the two types of anisotropic phase separation in the right panels of Fig. 1.
We find that, as expected from the previous study [76], uniaxial anisotropy dramatically changes the collective behaviors and causes long-range correlation, anisotropic phase separation, and critical phenomena that are presumably in the same universality class as that for uniaxial dipolar ferromagnets.Furthermore, uniaxial anisotropy suppresses the growth of gas bubbles in MIPS [66] and stabilizes macroscopic phase separation.Developing a coarse-grained model for particles with anisotropic self-propulsion, we corroborate the generality of the observed phenomena.

II. MICROSCOPIC MODELS
In this section, we explain the numerical implementation of uniaxial ABPs and RDLG, which are anisotropic extensions of the isotropic ABPs and equilibrium lattice gas, respectively.We also present phase diagrams for the two models, which provide preliminary insights into collective behaviors.

A. Active Brownian particles with uniaxial anisotropy
For uniaxial ABPs, N particles are confined in [0, L x ] × [0, L y ] with periodic boundary conditions.The state of the ith particle is specified by position r i and polarity angle θ i .The time evolution of (r i , θ i ) is governed by where a ∈ {x, y}, , and n i := (cos θ i , sin θ i ).Also, ξ i , ξ ⊥ i , and ξ θ i are Gaussian white noises with zero mean and unit variance.The translational noise η a i , representing thermal noise, is added to satisfy the detailed-balance condition when the self-propulsion force F 0 n b i is absent.We assume the two-body interaction as V(r) = (k/2)(σ − r) 2 for r < σ and V(r) = 0 otherwise.The potential for the polarity angle, U(θ), is added to model the effect of spatial anisotropy on self-propulsion, and ǫ (≥ 0) represents the strength of anisotropy.In this work, we use a simple potential function, U(θ) = − cos(2θ), which enhances the alignment of polarity along the x-axis (i.e., θ = 0 or π).Note that the polarity angle of each particle can take any value between 0 and 2π, in contrast to the previous model [76], in which the polarity angle is restricted to 0 or π.We also stress that we consider anisotropy of the self-propulsion direction, not of the particle shape.
In the case of ǫ = 0, the properties of the model [Eq.( 1)] have been studied in Ref. [66].In particular, the anisotropic mobility tensor µ ab i has been used to enhance the nucleation of gas bubbles in the phase-separated state.In our numerical simulations, we follow Ref. [66] and use anisotropic µ ab i for the case of ǫ > 0. While the anisotropy of µ ab i results in the polarity-dependent response of particle motion to the force, it does not induce spatially anisotropic particle motion along a fixed axis, in contrast to the spatial anisotropy caused by ǫ.Thus, in the following, we use the term "anisotropy" to refer to the effects of ǫ.
Figure 2(a) displays snapshots with two sets of parameters, which show that this model undergoes anisotropic phase separation.We stress that there is no attractive interaction in uniaxial ABPs, just like isotropic ABPs.As suggested in Fig. 1(d), this phase separation originates from the self-propulsion of each particle.We also present the phase diagram in Fig. 2(c); phase separation emerges for large Pe, which is also the same as in isotropic ABPs.Thus, this phase separation is regarded as the anisotropic extension of isotropic MIPS [see Fig. 1(c)].

B. Randomly driven lattice gas
For RDLG, we consider N particles on a square lattice with system size (L x , L y ) in units of the lattice constant.The state of the ith site is specified by occupation number n i , and the set of n i represents the configuration of the whole system.We assume exclusion between particles so that each site can be occupied by at most one particle, i.e., n i ∈ {0, 1}.We also consider attractive interaction between neighboring particles, which is represented by the following Hamiltonian: (2) The state of the system is updated in three steps: 1. We randomly choose two adjacent sites, (i, j), and calculate the energy difference (∆H) between the original configuration and the new configuration obtained by exchanging the state of the ith site with the state of the jth site.
2. If sites (i, j) are located along the x-axis, the new configuration is accepted with probability min(1, e −β∆H ).
3. If sites (i, j) are located along the y-axis, the new configuration is accepted with probability min(1, e −β(∆H+Eη) ), where E is the strength of the driving force, and η is a random number drawn from a Gaussian distribution with zero mean and unit variance.
For step 3, the random driving force is applied along the y-axis.We basically set the parameters to J = 4 and E = 100 and control β and ρ := N/L x L y .Note that our numerical implementation leads to the same type of macroscopic behaviors as those in the previous studies [27,84].
Since the driving force along the y-axis (i.e., Eη) competes with and effectively weakens the attractive interaction (i.e., ∆H), the motion of particles that interact with the neighbors is enhanced along the y-axis.In particular, E = 100 is practically equivalent to the limiting case with E = ∞, where the configuration is updated regardless of the value of ∆H.This limiting case has been commonly used in simulations of DLG and RDLG [7].
Figure 2(b) displays snapshots with two sets of parameters.This model undergoes phase separation induced by the attractive interaction [Fig.1(b)] though the motion of each particle is affected by the random driving force.We present the phase diagram in Fig. 2(d); phase separation is controlled by inverse temperature β, just like in equilibrium particle systems with attractive interactions [see Fig. 1(a)].

C. Orientation of phase separation
Self-propulsion is favored along the x-axis in uniaxial ABPs, while the driving force is applied along the y-axis in RDLG.Despite this difference in the direction of the enhanced particle motion, the dense and dilute regions are segregated along the x-axis in both uniaxial ABPs and RDLG [see Figs.2(a) and (b)].Such a coincidence of the collective behavior can be interpreted from a microscopic viewpoint as follows.For uniaxial ABPs, self-propulsion induces persistent collision of particles along the x-axis, leading to effective adhesion between particles along the x-axis.Since this type of collision is less probable along the y-axis, particles can move more freely along the y-axis.Thus, particle clusters that are caused by the effective adhesion should be elongated along the y-axis, which results in the segregation along the x-axis [see Fig. 1(d)].Note that similar cluster patterns have been recently found in simulations of ABPs with anisotropic selfpropulsion [81].For RDLG, the driving force enhances the free motion of particles along the y-axis.Thus, particle clusters caused by the attractive interaction should be elongated in the y direction, leading to the segregation along the x-axis [see Fig. 1(b)].See Appendix A for further comparisons between uniaxial ABPs and RDLG.

III. PROPERTIES OF HOMOGENEOUS STATE
Hydrodynamic descriptions are helpful in understanding the collective behavior of particles.For RDLG, homogeneous state properties have been studied using a linear coarsegrained model [27,85]: Here, φ(r, t) is the density fluctuation field, ξ(r, t) is a Gaussian noise with ξ a (r, t) = 0 and ξ a (r, t)ξ b (r ′ , t ′ ) = δ ab δ(r − r ′ )δ(t − t ′ ).In the isotropic limit (K xx = K xy = K yx = K yy = K, a x = a y = a, and D x = D y = D), Eq. ( 3) is reduced to the so-called model B [86], where H is a coarse-grained Hamiltonian: Thus, Eq. ( 3) is regarded as an extension of model B to an anisotropic system that respects the symmetry of particle dynamics in RDLG.
In the following, we demonstrate that the homogeneous states of uniaxial ABPs and RDLG exhibit the same type of long-range correlation as a generic feature of the nonequilibrium collective dynamics, which can be explained by Eq. ( 3).In Appendix D, using the well-known correspondence between RDLG and uniaxial dipolar ferromagnets [27,32], we further establish the connection between uniaxial ABPs and dipolar ferromagnets.

A. Long-range density correlation
The steady-state long-range correlation of a conserved quantity has been recognized as a general feature of nonequilibrium systems with anisotropic dynamics [11,13].Specifically, the fluctuation of a conserved quantity, which we denote as δA(r) here, decays as where • is an ensemble average in the steady state, and c eq and c neq are constants.The first term represents an exponential decay that also appears in equilibrium systems, while the second term is a nonequilibrium correction that leads to the long-range correlation with a power-law decay.The presence of long-range correlation (i.e., c neq 0) is ubiquitous in nonequilibrium systems with spatial anisotropy.
In uniaxial ABPs and RDLG, the self-propulsion and driving force violate the detailed balance in a spatially anisotropic way, respectively.Thus, the long-range correlation of the density field, which is a locally conserved field, is expected to appear in both systems.Though RDLG has been known to show the long-range correlation [11,27,85], for completeness, we explain the results for uniaxial ABPs and RDLG in parallel.Assuming small self-propulsion Pe in uniaxial ABPs and low inverse temperature β in RDLG [see the plus sign (+) in Figs. and respectively.Here, ρ(r) := N i=1 δ(r−r i ), δρ(r) := ρ(r)− ρ(r) , and ρ(k) is the Fourier transformation of ρ(r).
We show the heatmaps of S (k) for uniaxial ABPs and RDLG in Figs.3(a) and (b), respectively, both of which exhibit owl-like or butterfly-like patterns [85].Analytically, the observed pattern of S (k) can be characterized by the discontinuity at the origin in the Fourier space, i.e., lim This discontinuity of S (k) reflects the power-law decay of C(r) in the real space [85].As shown in Figs.

B. Linear coarse-grained model
According to the previous studies, the owl-like pattern of the structure factor observed in RDLG [Fig.3(b)] can be reproduced by the linear coarse-grained model [Eq.( 3)] [85].The similar pattern observed in uniaxial ABPs [Fig.3(a)] suggests that uniaxial ABPs and RDLG share the same macroscopic dynamics described by Eq. (3).To confirm the validity of Eq. (3) for both uniaxial ABPs and RDLG, we examine the structure factor for the coarse-grained density fluctuation, S lin (k) := | φ(k)| 2 /(L x L y ), and φ(k) is the Fourier transformation of φ(r).From Eq. (3), we can obtain [7,85] (10) For uniaxial ABPs, we fit the simulation data of 10), using D x , D y , a x , a y , K xy , and K yy as fitting paramters with K xx = 1.The fitting results are as follows: D x = 0.0287, D y = 0.00600, a x = 0.0990, a y = 0.0778, K xy = 0.525, K yy = 0.145.(11) In Figs.4(a) and (c), we plot the observed S (k) (with dots) and the fitted S lin (k) (with lines).The results show that Eq. ( 10) quantitatively reproduces the observed behavior of the structure factor for small |k|, which reflects the long-wavelength density fluctuation.We also fit the simulation data of RDLG in the same way as used for uniaxial ABPs.The fitting results are as follows: In Figs.4(b) and (d), we compare the observed S (k) and the fitted S lin (k), which show quantitative agreement as expected.
As discussed in previous studies of DLG and RDLG [7], we can derive the asymptotic behavior of the long-range part of the correlation function, C lin (r), which is the inverse Fourier transformation of S lin (k).From Eq. (10) we can obtain which is also consistent with the power-law decay of C(r) observed in uniaxial ABPs [Fig.3(c)] and RDLG [Fig.3(d)].

IV. PHASE SEPARATION PROPERTIES
As briefly explained in Sec.II, uniaxial ABPs and RDLG undergo anisotropic phase separation (Fig. 2).In this section, we investigate the properties of phase separation of uniaxial ABPs in more detail.We focus on the nucleation of persistent gas bubbles and the possibility of microphase separation, which have been found in recent studies [66].

A. Anisotropy-induced removal of gas bubbles
In Fig. 5 procedure for drawing this figure is given in Appendix E. From this figure, we find that for ǫ = 0, numerous gas bubbles are nucleated within the liquid phase.Throughout this paper, we use a "gas bubble" to refer to a connected region of the gas phase surrounded by the largest liquid phase.Note that we regard the largest gas phase as the gas reservoir and not as the gas bubble (see Appendix E 2 for the method to detect gas bubbles).As ǫ increases, the number of gas bubbles decreases.For sufficiently large values of ǫ (e.g., ǫ = 0.02), the presence of gas bubbles becomes less evident.To quantitatively characterize this observation, we define the bubble fraction as where S := L x L y and S bubble is the total area occupied by gas bubbles.We plot f b as a function of ǫ in Fig. 5(b), which shows that the fraction of gas bubbles monotonically decreases as ǫ increases.For sufficient large ǫ, f b reaches zero, indicating the absence of gas bubbles.This observation demonstrates that the uniaxial self-propulsion prevents the nucleations of gas bubbles.
In isotropic ABPs (i.e., ǫ = 0), the nucleation of gas bubbles has been examined in Ref. [66], which has revealed a connection between the existence of gas bubbles and a novel type of phase separation called microphase separation [73].To briefly explain the previous results in Ref. [66], we focus on the size distribution of gas bubbles divided by the total liquid area, n(a)/S liq , where a is the area of a single bubble.In Fig. 5(d), we plot n(a)/S liq for isotropic ABPs.We find that n(a)/S liq for large a fits well with the power-law decay observed in the reduced bubble model [66]: Considering that the bubble fraction, f b , and the size distribution, n(a), are related as [87] we can derive the system size dependence of f b as Here, χ liq := S liq /S and χ gas := 1−χ liq represent the area fractions of the liquid and gas phases, respectively, and are nearly independent of the system size, S .Thus, as S increases, f b is expected to increase until it reaches the area fraction of the gas phase, χ gas .This implies that the whole gas phase exists as persistent gas bubbles surrounded by the liquid phase.This state has been defined as the microphase-separated state [66].
As seen in Fig. 5(a), we find that gas bubbles are still observed for small but finite ǫ.We consider whether the size distribution of such gas bubbles can show the power-law decay as observed in isotropic ABPs (i.e., ǫ = 0).In Fig. 5(d), we plot n(a)/S liq for ǫ = 0.002.In contrast to the isotropic case, the bubble size distribution does not show the power-law behavior.Note that this result is not attributed to the finite-size effect since n(a)/S liq for different system sizes fall on a universal curve.More specifically, n(a)/S liq for ǫ = 0.002 decays faster than a −2 .From Eq. ( 17), f b is expected to converge to zero in the large system size limit, implying that uniaxial ABPs undergo macroscopic phase separation rather than microphase separation.Thus, we confirm that the type of phase separation significantly changes by the anisotropic self-propulsion.We also plot the ǫ dependence of n(a)/S liq for (L x , L y ) = (2880, 1440) in Fig. 5(e), which shows that the functional form of n(a)/S liq is changed by a small amount of ǫ.This suggests that microphase separation can be prohibited even for extremely small ǫ (e.g., ǫ = 0.0005), though we need a more detailed finite-size scaling analysis to draw a conclusion.
We comment on possible gas bubbles in RDLG.Note that previous studies on RDLG have not reported any possibility of microphase separation.As shown in Fig. 6(a), the nucleation of gas bubbles is hardly observed in typical snapshots for large systems, and macroscopic phase separation is expected to appear regardless of the strength of anisotropy.The bubble fraction, f b , plotted in Fig. 6(b) suggests that the nucleation of gas bubbles is suppressed by anisotropic driving force E in a similar way to uniaxial ABPs.

B. Nonlinear coarse-grained model
Though the linear coarse-grained model [Eq.(3)] succeeds in explaining the homogeneous state far from the critical point as discussed in Sec.III, it cannot describe phase separation since nonlinear terms are not included.In previous studies on isotropic ABPs [66], the qualitative features of microphase separation and the mechanism behind the observed persistent gas bubbles have been demonstrated using a coarse-grained model called Active Model B+ (AMB+) [73,88].To discuss the observed suppression of gas bubbles by the anisotropic self-propulsion from a general perspective, we consider an anisotropic extension of AMB+: which is also regarded as a nonlinear extension (i.e., adding the b, λ, and ζ terms) of Eq. ( 3).The b term can be derived from a coarse-grained Hamiltonian, and the λ and ζ terms reflect the violation of the time-reversal symmetry [73].To improve numerical stability, the higher-order gradient term with a small K ′ is also introduced.This term is irrelevant in the RG sense (see Appendix G 2 for the detail) and is not expected to affect the qualitative phase behavior.For simplicity, the effect of anisotropy is minimally retained in the difference between a x and a y .
We explain the isotropic limit (a x = a y ) with the present parameter set.In the low-density case, we observe phase separation with persistent gas bubbles [Fig.7(a), left], which is similar to the behavior of uniaxial ABPs [Fig.5(a), left].In the high-density case, we observe microphase separation, where gas bubbles are present throughout the system [Fig.7(b), left].Such phase behaviors are consistent with the previous observations in the isotropic AMB+ [73].
We consider the effect of anisotropy on phase separation with gas bubbles [Fig.7(a)].Similarly to the observation in uniaxial ABPs [Fig.5(b)], we find the suppression of bubble fraction f b as shown in Fig. 7(c).This suggests that the minimal extension of AMB+ (i.e., a x a y ) is sufficient to explain the qualitative behavior of uniaxial ABPs.We next examine the effect of anisotropy on microphase separation [Fig.7 Let us focus on the case with a x < 0 < a y [see the right panels of Figs.7(a) and (b)] to consider why strong anisotropy suppresses gas bubbles and stabilizes macroscopic phase separation.We neglect the noise term in Eq. ( 18) by the meanfield approximation, which has been used in the previous studies [67,68,73].Then, the linearized equation for φ − φ 0 is obtained in the Fourier space as From a x < 0 < a y , K > 0, and K ′ > 0, we see that the most unstable wavevector is along the k x -axis.Thus, we approximately neglect the modulation in the y direction and replace Eq. ( 18) by Here, chemical potential µ is a local quantity, in contrast to the isotropic limit (a x = a y ), where nonlocality of chemical potential can lead to phase separation with gas bubbles and microphase separation [73].Thus, macroscopic phase separation is expected to appear for a x < 0 < a y .20)] based on the RG analysis.

V. CRITICAL PROPERTIES
Since uniaxial ABPs and RDLG share the common properties in the homogeneous and phase-separated states (see Secs.III and IV), we expect that the critical point for anisotropic phase separation in each model belongs to the same universality class.In the following, we support this expectation using the RG analysis of the coarse-grained model [Eq.(18)] and the finite-size scaling analysis of simulation data for uniaxial ABPs.

A. Renormalization group analysis of coarse-grained model
We consider the critical phase transition between the homogeneous and phase-separated states in the coarse-grained model [Eq.( 18)] under sufficiently large anisotropy with a x < a y .We first review the previous RG analyses of Eq. ( 18) for K ′ = λ = ζ = 0 [26,27,31,32].Retaining only the relevant variables in the RG sense, we can obtain a model that is equivalent to a coarse-grained model of uniaxial dipolar ferromagnets, which have dipolar long-range interactions [26,27,31,32] (see Appendix D for the detail).At the two-loop level, the critical exponents for the coarsegrained model of uniaxial dipolar ferromagnets have been obtained [26,27,32] as Here, β is the exponent for the onset of the order parameter, and ν x and ν y (≃ 2ν x ) are the exponents for the divergent correlation lengths along the x-and y-axes, respectively.For RDLG, the finite-size scaling analysis of simulation data has been performed to obtain the critical exponents [27] as These values coincide with the RG results [Eq.(20)] within the numerical error, suggesting that the critical point for anisotropic phase separation in RDLG belongs to the universality class of uniaxial dipolar ferromagnets.
Considering nonzero λ and ζ to discuss the phase behavior of uniaxial ABPs (see Sec. IV), we can show that λ and ζ are irrelevant variables in the RG sense (see Appendix G 2 for the detail).This suggests that the introduction of small λ or ζ does not affect the critical properties of anisotropic phase separation, and the critical exponents remain the same as those given in Eq. (20).Thus, like RDLG, the critical point for anisotropic phase separation in uniaxial ABPs is expected to belong to the universality class of uniaxial dipolar ferromagnets.Note that the irrelevance of λ or ζ is further supported by the suppression of gas bubbles under strong anisotropy (see Fig. 7).

B. Connection to uniaxial dipolar ferromagnets
To study the critical point for anisotropic phase separation in uniaxial ABPs, we perform simulations with a fixed strength of anisotropy, ǫ = 0.01.Here, we assume that the critical exponents are not affected by the specific value of ǫ.First, assuming the law of rectilinear diameter [6,89], we estimate the critical density as ρ c = 0.71 (see Appendix F 1 for the detail).Next, we perform simulations with ρ = ρ c = 0.71 to identify the universality class of the critical point using the anisotropic finite-size scaling analysis, which has been widely applied to critical phenomena in externally driven sys-tems [7,22,90,91].Since the liquid and gas phases are separated along the x-axis for large Pe [Fig.2(a)], the degree of phase separation can be measured by an order parameter, The finite-size scaling hypotheses for m and the Binder ratio, U := m2 2 / m4 , are given as and respectively.Here, τ := Pe − Pe c is the distance from the critical point, and M and U are scaling functions.Equations ( 23) and ( 24) are extensions of the scaling hypotheses for isotropic systems with ν x = ν y [22], and the values of ν x and ν y can be different in anisotropic systems such as uniaxial ABPs and RDLG.For ν x ν y , to perform the finite-size scaling analysis, we need to vary the system size with L y /L x ν y /ν x fixed.Though ν y /ν x should be determined in principle by the finite-size scaling analysis, we choose ν y /ν x = 2, which has been commonly used for RDLG based on the RG analysis [26,27].Following this choice, we perform simulations with five different system sizes satisfying L y /L x 2 = 1/24 2 : (L x , L y ) = (180, 56.25), (210, 76.5625), (240, 100), (300, 156.25), and (360, 225).
The results of the finite-size scaling analysis are summarized in Fig. 8 (See Appendix F 2 for the detailed procedure).Varying Pe from 11.5 to 13.0, we find that U as a function of Pe for different system sizes approximately crosses at a unique point [Fig.8(a)], which suggests the presence of the critical point, Pe c .By fitting U(τ, L x ) and m (τ, L x ) with secondorder polynomials, we obtain Pe c as Pe c = 12.408(5) (25) and the critical exponents as Using these obtained values, we find that the rescaled plots of U and m collapse onto universal curves [Figs.8(b) and (c)], which validates the anisotropic finite-size scaling hypotheses given by Eqs. ( 23) and (24).The obtained β and ν x [Eq.( 26)] agree with the RG result for the coarse-grained model [Eq.(20)] and the simulation result of RDLG [Eq.(21)] within the error margin.This indicates that the critical phenomena in uniaxial ABPs belong to the universality class of uniaxial dipolar ferromagnets, as expected from the RG analysis (see Sec. V A).To check the consistency of the obtained values of β and ν x , we plot the L x dependence of ∂U/∂Pe and m at Pe = 12.415 (≃ Pe c ) in Figs.8(d) and (e).According to Eqs. ( 23) and ( 24), the slopes of ∂U/∂Pe and m on the logarithmic scale are 1/ν x and −β/ν x , respectively.Indeed, Figs. 8(d) and (e) show that the slopes are comparable to the counterparts for the two-loop RG result [Eq.( 20)].

VI. DISCUSSION
In this paper, to investigate the relation between MIPS and nonequilibrium phase separation caused by attractive interactions, we have studied the collective properties of 2D uniaxial ABPs, in which self-propulsion along the x-axis is favored.Performing simulations, we have found three distinctive features of uniaxial ABPs: (i) generic long-range density correlation in the homogeneous state, (ii) anisotropic phase separation with suppressed nucleation of gas bubbles in contrast to isotropic ABPs, and (iii) critical phenomena that presumably belong to the universality class of 2D uniaxial ferromagnets with dipolar long-range interactions.Since properties (i)-(iii) are common to RDLG, in which phase separation is induced by attractive interactions under external driving, we have established the connection between collective behaviors of uniaxial ABPs and RDLG.Additionally, we have constructed a nonlinear coarse-grained model [Eq.(18)] and substantiated the generality of properties (i)-(iii).
The critical exponents for the models related to this study are summarized in Table I, which points out that the critical behaviors of 2D uniaxial ABPs are close to those of the 3D Ising model rather than the 2D Ising model.This property is consistent with the previous study concerning 2D uniaxial ferromagnets with dipolar long-range interactions [35,36].For 2D uniaxial dipolar ferromagnets, the effective increase in dimensionality has been attributed to the consequence of the long-range correlation caused by the dipolar interactions.For 2D uniaxial ABPs, the long-range density correlation arising from the anisotropic nonequilibrium dynamics (see Sec. III) effectively increases the dimensionality from two to three, according to the analogy with uniaxial dipolar ferromagnets (see Appendix D for the detail).
Our results suggest that the origin of phase separation (i.e., self-propulsion or attractive interaction) is not essential for the collective behaviors of particles with anisotropic dynamics [Figs.1(b) and (d)].In contrast, for isotropic systems [Figs.1(a) and (c)], the collective phenomena of selfpropelled particles can be distinct from those of attractively interacting particles.Specifically, in 2D isotropic ABPs, persistent gas bubbles or microphase separation can appear (see Sec. IV) [66,73], and the universality class for critical phenomena can be different from the 2D Ising class [88].Further studies are required to elucidate the condition for such differences in isotropic systems.
Recently, a wide range of active matter phases has been realized using biological [42][43][44][45][46][47] and artificial [48][49][50][51][52][53][54][55][56]62] systems, especially under anisotropic conditions [92].The connection between uniaxial ABPs and RDLG suggests that active matter can serve as a platform for materializing the properties predicted for externally driven systems.Though we have focused on uniaxial anisotropy in this study, it will be interesting to examine whether the collective behaviors of the standard DLG can be observed in ABPs with unidirectional anisotropy, which can be relevant to biological systems with chemical gradients.We then observe the density ρ l and ρ h in the high-and lowdensity regions.In the phase-separated state, the values of ρ l and ρ h give the coexisting (binodal) curve, which is drawn in Fig. 2(c) and (d).
Appendix C: Parameter details of Figs. 3 and 4 We set the simulation box to L x = L y = 360.The particle number is set to N = 92016 for uniaxial ABPs and N = 64800 for RDLG, which respectively correspond to the density of 0.710 and 0.50.We start from the initial state in which the particles are randomly located with zero overlaps.We perform the relaxation run for 10 8 time steps (i.e., time = 10 8 dt = 2.0× 10 6 ) for uniaxial ABPs and for 4.0×10 6 Monte Carlo steps for RDLG.After that, we observe the structure factor S (k).The real-space density correlation ρ(r)ρ(0) is calculated by the inverse Fourier transformation of the structure factor S (k).
We take the time average in the steady state and the ensemble average over different noise realizations.For uniaxial ABPs, the ensemble average is performed over 28 different noise realizations, and the time average is performed over 400 samples obtained every 10 6 time steps (i.e., time = 10 6 dt = 20000).For RDLG, the ensemble average is performed over 96 different noise realizations, and the time average is performed over 400 samples obtained every 20000 Monte Carlo steps.

Appendix D: Relation to equilibrium uniaxial dipolar ferromagnet
For RDLG, it is known that the specific patterns of structure factor S (k) involving the long-range correlations are analogous to the long-range nature of the uniaxial dipolar system.Here, we give the definition of uniaxial dipolar ferromagnet [35] and briefly discuss the analogy between the density correlation of uniaxial ABPs and the spin correlation of uniaxial dipolar ferromagnet.
We start with the Heisenberg model with the short-range exchange interaction and long-range dipolar interaction.The Heisenberg spin S R is defined on the two-dimensional square lattice {R = (n x , n y ) | n x , n y = 0, ±1, ±2, • • • }, where the lattice constant is set to 1.The Hamiltonian H of this model consists of the short-range exchange interaction and long-range dipolar interaction, which is expressed as where R δ runs over all nearest-neighbor pairs.Let us impose the uniaxial condition where the Heisenberg spin S R is restricted to pointing in the direction of the y-axis: S R = (0, S R , 0).The model reduces to the Ising model with anisotropic interaction: This model is called the uniaxial dipolar ferromagnet.
In the Fourier space, the dipolar part of the Hamiltonian is expanded near the k = 0 as where {a i } i=1,••• ,4 is a set of numerical constants depending on the lattice structure.By expanding the short-range part of Hamiltonian in the same way, we rewrite the Hamiltonian as where we ignore the higher-order terms in S k .The values of the numerical factor are given in Ref. [35].
The equilibrium state of this system is described by the canonical ensemble.In the disordered state, the linear approximation leads to the static spin-spin correlation: with This form is the special case of Eq. ( 10), indicating that uniaxial ABPs acquire dipolar-like long-range natures.As discussed in Appendix G 2, this feature determines the universality class of critical phenomena.this observation, we can infer that the critical Péclet number, Pe c , lies between 12.0 < Pe c < 13.0 and the critical density, ρ c , is estimated as ≈ 0.708.

Estimation of critical exponents
Based on the estimation of the critical density in the previous section, we set the density to ρ = 0.710 and change the Péclet number, Pe, from Pe = 11.5 to Pe = 13.0.As explained in main text, we set the system sizes to (L x , L y ) = (180, 56.25), (210, 76.5625), (240, 100), (300, 156.25), and (360, 225).We show the typical time evolution of the ensemble average of the order parameter for (L x , L y ) = (240, 100), (300, 156.25), and (360, 225) in Fig. 10.This figure confirms that our simulation achieves the steady state after a sufficiently long relaxation run.Using the data within the red region, we take the time and ensemble averages for the order parameter m and the Binder Parameter U := m2 2 / m4 .The ensemble average is taken over 800 different noise realizations for (L x , L y ) = (300, 156.To estimate the critical exponents, we use the anisotropic finite-size scaling hypothesis Eqs. ( 23) and (23).We refer to Ref. [22] for a more detailed discussion of the anisotropic finite-size scaling.Since the scaling functions M and U are analytic, we can expand m (τ, L x ) and U(τ, L x ) around τ = 0 as For simulations of the coarse-grained model [Eq.( 18)], ), we discretize time as t = n∆t and spatial coordinates as x = i∆x and y = j∆y with periodic boundary conditions.Accordingly, we replace φ(x, y, t) by φ n i, j and ξ a (x, y, t) by (∆x∆y∆t) −1/2 ξ n a,i, j , where ξ n a,i, j is a Gaussian noise with ξ n a,i, j = 0 and ξ n a,i, j ξ n ′ b,i ′ , j ′ = δ ab δ ii ′ δ j j ′ δ nn ′ .Using the explicit Euler method, we replace Eq. (G1) by where [F(φ)] n i, j is the discretized form of the right-hand side of Eq. (G1).To determine [F(φ)] n i, j , we use the second-order central finite difference for the differential operators that appear in Eq. (G1) (i.e., ∂ x , ∂ y , ∂ x 2 , and ∂ y 2 ), such as , where φ 0 is the spatial average of φ(r, t).As the initial state for all the simulations, we use a phase-separated state, φ init (r) := −2sgn(φ 0 ) exp[−(x − L x /2) 4 /(L x /4) 4 ] − C, where C is a constant to set the spatial average of φ init (r) to φ 0 (Fig. 11).
We define the liquid and gas phases as the spatial regions satisfying φ(r) > 0 and φ(r) < 0, respectively.In the same way as applied to uniaxial ABPs (see Appendix E 2), a Julia package (JuliaImages.jl)is used to detect the connected regions of the gas phase.The size of each gas phase, a, is defined as the area of the regions that satisfy φ(r) < 0 and are connected to each other along the x or y-axis.The bubble fraction, f b , which is plotted in Figs.7(c) and (d), is calculated as f b := S gas − a max /(L x L y ), where S gas and a max are the total and maximum areas of the gas phase, respectively, and • • • means the average over samples.
To characterize the steady state using bubble fraction f b and order parameter m, independent samples are taken with different noise realizations.For the low-density condition with φ 0 = −0.1, which is used for Figs.7(a 192,192).To obtain the expectation values, we take the average over independent samples as well as the time average over 51 points in the last half of the total time.
We show the typical time evolution of φ in the liquid and gas phases (φ liq and φ gas , respectively), averaged over space and independent samples [Figs.12(a-c) and (e-g)].The points in the red region in Figs.12(a-c) and (e-g) are used in time averaging to obtain the a y dependence of φ liq and φ gas , which is plotted in Figs.12(d

Renormalization group analysis
Assuming anisotropic systems with a y > 0, we consider the critical phase transition between the homogeneous state and anisotropic phase separation that occurs as a x is changed.Applying the approach by Martin, Siggia, Rose, Janssen, and de Dominicis (MSRJD) [93][94][95][96] to Eq. (G1), we can obtain the probability density for a dynamical path of configurations {φ(r, t)} t∈[0,T ] as  which coincides with the effective action for the randomly driven lattice gas [26,27,31,32].

FIG. 1 .
FIG. 1. Four types of phase separation.The row and column correspond to the type of phase separation (attraction-or motility-induced) and the type of dynamics (isotropic or anisotropic), respectively.In each panel, a typical particle configuration obtained from model simulations is shown with schematic figures of the single particle motion and small cluster formation.(a) Brownian particles follow overdamped dynamics with attractive interactions (wavy lines) and random forces.(b) In RDLG, particles stochastically move with attractive interactions (wavy lines) and external driving force (red arrow) along an axis (i.e., y-axis in the figure).(c) ABPs show self-propelled motion (red arrow) with repulsive interactions and random forces.(d) Uniaxial ABPs show anisotropic self-propelled motion favored along an axis (i.e., x-axis in the figure) with repulsive interactions and random forces similar to (c) [see Eq. (1) for the detail].

FIG. 3 .
FIG. 3. Singular structure factors in the homogeneous states.(a, b) Heatmap of structure factor S (k) for (a) uniaxial ABPs and (b) RDLG.(c, d) Density correlation functions C(x, 0) (yellow) and C(0, y) (purple) for (c) uniaxial ABPs and (d) RDLG.In the insets, the absolute value is plotted on the log-log scale.The parameters used for (a, c) and (b, d) are the same as those for the left panels of Figs.2(a) and (b), respectively.The system size is set to L x = L y = 360 for both models.

FIG. 4 .
FIG. 4.  Quantitative comparison between the simulated structure factor and the theoretical expression [Eq.(10)],where the same data as plotted in Fig.3is used.(a, b) Structure factor S (k x , k y ) with k y = 4π/L y , 8π/L y , and 12π/L y for (a) uniaxial ABPs and (b) RDLG.(c, d) Structure factor S (k x , k y ) with k x = 4π/L x , 8π/L x , and 12π/L x for (c) uniaxial ABPs and (d) RDLG.In all figures, the colored dots represent the simulation results, and the black lines represent the theoretical expression with the best-fit parameter.
2(c) and (d)], we focus on typical homogeneous states [Figs.2(a) and (b)].We calculate the structure factor and the two-point correlation function, which are defined as 3(c) and (d), the correlation function [C(x, y = 0) (yellow) and C(x = 0, y) (purple)] indeed shows a power-law decay as ∼ r −2 , which implies the long-range density correlation.The negative correlation observed in C(x, y = 0) suggests the formation of transient clusters elongated along the y-axis.This orientation of clusters is consistent with the configurations in phase separation shown in Figs.2(a) and (b) (see Sec. II C).
FIG. 6. Absence of gas bubbles in RDLG.(a) Typical snapshots in the steady state for three values of E. The colors represent the particle density from 0 (blue) to 1 (red).(b) Bubble fraction f b as a function of E. In all figures, the parameters are chosen as ρ = 0.5, β = 0.556, and (L x , L y ) = (720, 360).
(b)].We find that microphase separation discontinuously changes into macroscopic phase separation, indicated by the abrupt change in f b [Fig.7(d)].In addition, we define an order parameter for macroscopic phase separation along the x-axis as m := S (k x = 2π/L x , 0), where the structure factor is defined as S (k) := | φ(k)| 2 /(L x L y ) with φ(k) := d 2 r e −ik•r φ(r).As shown in the inset of Fig. 7(d), the discontinuous change in m also suggests the discontinuous transition between microphase separation and macroscopic phase separation.

FIG. 8 .
FIG. 8. Finite-size scaling analysis for uniaxial ABPs.The parameters are chosen as ρ = 0.71, (µ , µ ⊥ , µ θ ) = (1, 0.25, 1.5), and ǫ = 0.01.(a) The Binder ratio U as a function of Pe for different system sizes.(b) U and (c) the rescaled order parameter m as functions of the rescaled Pe with the best-fitted critical exponents (β/ν x = 0.540, 1/ν x = 1.54).(d) ∂U/∂Pe and (e) m against L x near the critical point in the log-log plot.In (d) and (e), the red dashed lines represent ∂U/∂Pe ∝ L x 1/νx and m ∝ L x −β/νx , respectively, with the critical exponents used in (b) and (c), and the blue dashed lines are counterparts for the expected universality class [Eq.(20)] based on the RG analysis.

P 0 dt d 2 r 3 − b y ∂ y 2 φ 3 +
FIG.12.Typical time and a y dependence of φ in the liquid and gas phases (φ liq and φ gas , respectively), averaged over space and independent samples.We show the time evolution with (a) a y = −0.25,(b) a y = 0, and (c) a y = 0.2 for φ 0 = −0.1 and (L x , L y ) = (256, 128), as well as (e-g) the counterparts for φ 0 = 0.4 and (L x , L y ) =(192, 192).The values at equally spaced 51 time points within the red region in (a-c) and (e-g) are used in time averaging to obtain the a y dependence of φ liq and φ gas , which is plotted in (d) and (h).

FIG. 13 . 3 +
FIG.13.Typical time dependence of bubble fraction f b , averaged over independent samples.The same parameters as used in Figs.12(a-c) and (e-g) are used, and the error bar represents the standard error.The values at equally spaced 51 time points within the red region are used in time averaging to obtain the a y dependence of f b , which is plotted in Figs.7(c) and (d).

TABLE II .
Basic features in microscopic implementation of uniaxial ABPs, RDLG, isotropic ABPs, and equilibrium LG.