{'fraction_non_alphanumeric': 0.03785448479404594, 'fraction_numerical': 0.012275974164848796, 'mean_word_length': 4.541597534960891, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Microscopic theory for the phase separation of self-propelled repulsive disks

Motivated by recent experiments on colloidal suspensions, we study analytically and numerically a microscopic model for self-propelled particles lacking alignment interactions. In this model, even for purely repulsive interactions, a dynamical instability leading to phase separation has been reported. Starting from the many-body Smoluchowski equation, we develop a mean-field description based on a novel closure scheme and derive the effective hydrodynamic equations. We demonstrate that the microscopic origin of the instability is a force imbalance due to an anisotropic pair distribution leading to self-trapping. The phase diagram can be understood in terms of two quantities: a minimal drive and the force imbalance. At sufficiently high propulsion speeds there is a reentrance into the disordered fluid.

Introduction. -Living "active matter" [1] ranging from bacterial suspensions [2] to flocks of birds [3] is an emerging paradigm on the interface of physics, chemistry and biology. These systems are composed of identical subunits, which are able to show collective dynamical behavior such as swarming [4], active clustering [5], and even micro-bacterial turbulence [6]. The fact that every individual by itself is far from equilibrium allows a virtually unlimited number of propagation rules to be devised and, moreover, the emergence of purely dynamical phases. This is in contrast to the phenomenon of equilibrium phase separation [7] such as the liquid-vapor transition, which can be understood on thermodynamic grounds going back to the seminal work of van der Waals [8]. For one-component systems as considered here, equilibrium phase separation requires attractive interactions.
Theoretical progress in statistical physics is often based on insights gained from the study of models that are simple but still capture an essential feature of real, complex systems. Self-propelled particles constitute one such class of models that has been very successful in describing nonequilibrium collective behavior. For example, a large class of models [9,10]-the most famous of which is the Vicsek [4] model-describe particles that move with constant velocity and align their orientations with the average orientation of neighboring particles. Such an alignment can also arise due to volume exclusion as has been observed for granular rods [11]. In the Vicsek model, a dynamical phase transi-tion from an ordered into an unordered phase takes place as a function of the orientational noise strength.
Following earlier work [12,13], suspensions of selfpropelled spherical colloidal particles have recently been realized experimentally. The propulsion mechanism is diffusiophoresis: either the catalytic decomposition of hydrogen peroxide on a platinum hemisphere [14] or hematite in conjunction with light [15], or the reversible local demixing of a near-critical water-lutidine mixture [16,17]. The non-equilibrium behavior of these driven suspensions is unexpectedly rich and includes a stable fluid of "living" clusters that form, break, and merge but with a steady mean size [15]; and a phase separation, i.e., the largest cluster grows until it is composed of a finite fraction of the particles [17]. In a model of run-and-tumble bacteria, Tailleur and Cates have shown theoretically that such a dynamical instability resembling liquid-vapor phase separation can be attributed to a density-dependent mobility (or propulsion speed) [18][19][20]. Non-equilibrium phase separation has been shown to also occur in computer simulations of a minimal model for self-propelled repulsive disks, where particle orientations are independent and can be modelled as free rotational diffusion [17,21,22]. Although recently a link of run-and-tumble motion to active Brownian motion has been made [23], a connection between the microscopic details and the large-scale dynamical instability leading to phase separation has been missing.
In this Letter, we close this gap between numerical simulations and experiments of colloidal suspensions on one side and phenomenological models positing a densitydependent mobility on the other side by deriving an explicit expression for the effective swimming speed from first principles. To this end, we start from the full manybody Smoluchowski equation of the minimal model and obtain the effective evolution equation for a single tagged particle. We employ a novel scheme to close the ensuing hierarchy of evolution equations and discuss the approximations involved. The effective swimming speed of the tagged particle is reduced compared to its free swimming speed due to a force imbalance that is quantified by the pair distribution function. We finally test our results employing Brownian dynamics computer simulations for different repulsive pair potentials.
Model. -We consider a suspension of N colloidal particles with number densityρ moving in two dimensions. The coupled equations of motion arė where the Gaussian noise ξ k models the interactions with the solvent. The noise has zero mean and temporal shortranged correlations Particles interact through a pair potential u(r) with total potential energy U = k<k u(|r k − r k |). In addition to the conservative force due to the potential energy, every particle is propelled with a constant speed v 0 in the direction characterized by the angle ϕ k the particle orientation encloses with the x-axis. The speed v 0 corresponds to the free swimming speed in a dilute suspension. It characterizes the propelling force without specifying details of the actual propulsion mechanism. Throughout, we employ reduced units with a length scale a related to the particle size, unit of energy given by the thermal energy, and unit of time τ 0 = D 0 /a 2 , where D 0 is the bare diffusion coefficient (as measured in a dilute passive suspension). We assume that particle orientations undergo free diffusion with rotational diffusion coefficient D r , Real colloidal particles are of course three-dimensional objects (i.e., spheres) and thus rotate in three dimensions even if their translational motion is confined to (quasi) two dimensions. However, for the sake of simplicity, here we assume that particles rotate only about the fixed z-axis. Moreover, for particles obeying the no-slip boundary condition rotational and translational diffusion coefficient are coupled. Nevertheless, we will treat D r as a free parameter for most of our considerations.
Derivation. -The time evolution of the normalized joint probability distribution ψ N ({r k , ϕ k }; t) is governed by the Smoluchowski equation The first step is to derive an approximate equation of motion for the projected density ψ 1 (r 1 , ϕ 1 ; t) = dr 2 · · · dr N dϕ 2 · · · dϕ N N ψ N (6) of a single particle. Since particles are identical, in the following we simply designate particle 1 as the tagged particle and drop the subscript. Performing the integration, Eq. (5) becomes with mean force exerted by the surrounding particles onto the tagged particle. Here, ψ 2 (r, ϕ, r ; t) is the two-body density distribution function to find another particle at r (with arbitrary orientation) together with the tagged particle at r with orientation ϕ. Following this scheme leads to the wellknown BBGKY hierarchy of coupled equations [24]. In order to proceed, we need to find an approximate closure. The two-body density ψ 2 can be decomposed into the product of the conditional probability g(|r − r |, θ|r, ϕ) to find another particle at position r given there is a selfpropelled particle at r with orientation e times the probability to find a particle at r, ψ 2 (r, ϕ, r ; t) = ψ 1 (r, ϕ; t)ρg(|r − r |, θ|r, ϕ; t). (9) Here, θ is the angle enclosed by the displacement vector r − r (pointing away from the tagged particle) and the orientation e.
Inserting the decomposition Eq. (9) into Eq. (8), the projection of the force onto the orientation reads e · F = −ρζψ 1 with coefficient where u(r) is the pair potential. So far, we have made no approximations. To proceed, we assume that the coefficient ζ is independent of the position of the tagged particle, which is equivalent to assuming that the system is homogeneous. Moreover, we neglect the time-dependence of the pair distribution function, g(r, θ|r, ϕ; t) ≈ g(r, θ).
These are sufficient conditions to study the onset of a p-2 possible dynamical instability of an initially homogeneous suspension. We decompose the mean force in the (nonorthogonal) basis spanned by the orientation e and the density gradient ∇ψ 1 , and expand to linear order of |∇ψ 1 | with result Hence, Eq. (7) now reads which is the desired closed equation for the temporal evolution of the tagged particle density. Here, the effective enters. In the following, we assume that D is a constant coefficient that does not depend on the position of the tagged particle. It thus corresponds to the long-time diffusion coefficient of the passive suspension (v 0 = 0). Eq. (12) together with Eq. (13) constitutes our first main result.
The physical picture we have in mind is the following: As shown in Fig. 1(a), the propulsion of particles leads to a pair distribution g(r, θ) that is anisotropic, i.e., there is a depletion zone behind the tagged particle and an excess zone in front of the particle. The effective swimming speed v is reduced due to this anisotropy. This is intuitively clear since particles block each other in a dense suspension due to volume exclusion or repulsive interactions. The projection coefficient ζ quantifies how strongly the tagged particle is slowed down by its neighbors as a function of the particle densityρ. Note that for passive particles with v 0 = 0 the pair distribution g(r) becomes isotropic and, therefore, ζ = 0 vanishes.
Dynamical instability. -To make the connection to previous work on mobility-induced phase separation [23], we now cast Eq. (12) into the form of the more familiar effective hydrodynamic equations (see also Refs. 25,26). To this end we introduce the particle density and the orientational field as the first moment of the tagged particle orientation. Plugging in Eq. (12), we find which couples the temporal evolution of the density to the orientational field. The orientational field in turn couples to an integral involving  Fig. 1: (a) Anisotropic pair distribution function g(r, θ) for a tagged particle (white circle) at the origin with its orientation pointing along the y-axis (arrow). There is a larger probability to find another particle in front of the tagged particle compared to finding it behind. For a closure, we drop the second term corresponding to the second moment. Again, this approximation is only justified at the onset of a large-scale instability close to the homogenous state. We thus arrive at The stationary solution of the two effective hydrodynamic equations (16) and (18) is ρ(r, t) =ρ and p(r, t) = 0. Clearly, for a constant effective speed v this solution is always stable. Now suppose that ρ(r, t) corresponds to a slowly varying density profile such that within the short range of the pair potential the density is approximately constant. We can then replace the global densityρ in the expression for the effective swimming speed Eq. (13) by the local density with v(ρ) = v 0 − ρζ. Following Ref. 23, we consider the limit of large length scales and time scales much longer than 1/D r . Eq. (18) then implies the quasi-stationary solution i.e., the orientational field is adiabatically enslaved to the density and points along the density gradient. Plugging this result back into Eq. (16) leads to the diffusion equation with collective diffusion coefficient D(ρ) governing the density relaxation. Instability of the homogenous suspension is signaled by D(ρ) < 0. From this criterion we find the instability region ζ − ζ ζ + bounded bȳ p-3  Fig. 1(b). The solid line encloses the instability region. The vertical dashed lines indicate the corresponding critical speeds vc.
Eq. (21) constitutes our second main result. Fig. 1(b) shows the instability region assuming that both the propulsion speed v 0 and the coefficient ζ can be controlled independently. Of course, keeping all other parameters fixed, ζ is a function of v 0 and thus describes a curve in this plot. Our theory is only able to predict the onset of this dynamical instability. Below we will demonstrate using computer simulations that phase separation ensues from this instability. Hence, phase separation is predicted to occur within the shaded region if ζ(v 0 ) crosses the boundary. A curve ζ(v 0 ) can enter the instability region for v 0 > v * and, as we will demonstrate below, leave it for larger speeds corresponding to a reentrance into the disordered fluid phase. From Eq. (12) we can also determine the self-diffusion coefficient i.e., the long-time diffusion coefficient of the interacting suspension can be calculated as the diffusion coefficient of a free particle but employing the reduced swimming speed v instead of v 0 . The swimming speed v can be calculated via Eq. (13) and is not a fit parameter as in Ref. 21.
Brownian dynamics simulations. -We test our theoretical predictions with Brownian dynamics simulations. We consider up to N = 4900 particles in a quadratic box and employ periodic boundary conditions. The equations of motion Eq. (1) are integrated with time step 5 · 10 −6 -10 −5 depending on the speed v 0 . Motivated by recent experiments [17], we first study a suspension of (almost) hard disks employing the WCA potential [27] u(r) = 4ε δ r 12 − δ r 6 + ε (r 2 1/6 δ) where we set the potential strength ε = 100. Particles have diameter unity in reduced units. We set 2 1/6 δ = 1 so that interactions are only present for overlapping particles. We study two area fractions φ = 0.4 and φ = 0.5, where φ =ρπ/4. Moreover, we assume that the no-slip boundary condition holds, which determines the rotational diffusion coefficient as D r = 3 in reduced units.
We perform a cluster analysis based on a simple overlap criterion. Two particles are considered as "bonded" if their distance is smaller than unity (they overlap), and a cluster is the set of all mutually bonded particles. As a geometrical order parameter, in Fig. 2(a) we plot the mean size P of the largest cluster as a function of the propulsion speed v 0 . For both densities there is a transition from the unordered, homogeneous phase to an ordered phase, where the largest cluster occupies a finite fraction of the system. This drive-induced transition occurs as we increase the speed v 0 and resembles a second order transition, i.e., the order parameter P does not jump but increases continuously.
In order to gain a first insight into the universality class of the transition, we have performed finite-size scaling [28]. To this end, we have simulated smaller systems with particles N = L 2 for L = 40, 50, 70. We employ L as the relevant system size and not the actual length of the simulation box since the order parameter P is based on the fixed particle size. Defining ≡ v 0 /v c − 1, finitesize scaling predicts the order parameter to behave like P L (v 0 ) = L −β/νP (L 1/ν ) with universal scaling functioñ P (x). Fig. 2(b) shows that the data for different system sizes indeed collapses onto a single curve using the equilibrium Ising exponents ν = 1 and β = 1/8 even though their use lacks a firm theoretical basis. The quality of the collapse is somewhat better for v 0 v c . From the collapse we have estimated the critical speeds v c 50 (φ = 0.4) and v c 38 (φ = 0.5).
From the numerical data, we calculate the projection coefficient ζ by evaluating the integral Eq. (10). In Fig. 2(c), the long-time self-diffusion coefficient D lt is plotted as a function of the effective speed v = v 0 −ρζ. Below the transition v 0 < v c , the diffusion coefficient is indeed very well described by D lt = D + v 2 /(2D r ) ≈ v 2 /(2D r ) as predicted and confirms the validity of Eq. (12). For v 0 > v c we observe that the effective speed v stays approximately constant while the diffusion coefficient D lt grows further. Note that in the phase-separated suspension the global speed v and diffusion coefficient D lt are still calculated from all particles. Although particles are slower while part of a large cluster, there is a constant particle exchange between clusters and dilute phase leading to a monotonously increasing D lt as a function of v 0 . Finally, in Fig. 2(d) we plot the force coefficient as a function of reduced speed v 0 /v * , where v * follows from Eq. (22) with the long-time diffusion coefficient D ( 0.46, 0.36 for φ = 0.4, 0.5) measured at equilibrium (v 0 = 0). In agreement with Fig. 2(a), for low speeds the value of ζ corresponds to the homogeneous suspension. The speeds where ζ crosses the stability boundary agree well with the finite-size scaling estimates of the critical speeds v c .
The observed phase separation is a robust phenomenon that does not depend on the employed pair potential. For a demonstration we have gathered numerical data for three further potentials: (H) soft spheres with a harmonic repulsion u(r) = ε(r − 1) 2 for r 1 and u(r) = 0 otherwise [21], (GCM) the Gaussian core model u(r) = εe −r 2 , and (Y) the screened Coulomb or Yukawa potential u(r) = εe −κ(r−1)−1 /r with κ = 5. For the last two potentials we assume the effective particle diameter to be unity. In this section, we relax the coupling of rotational and translational diffusion by setting D r = 3 · 10 −5 and neglecting the translational noise altogether, assuming that the effective noise induced by the propulsion is dominant. To be as general as possible, we absorb the strength of the potential in the modified time unit τ 0 /ε.
We study the three potentials for different speeds v 0 at densities: (H) φ = 0.7, (GCM) φ = 0.091, and (Y) φ = 0.2. Fig. 3(a) shows the structure factor S(q) ≡ kk e iq·(r k −r k ) /N for different speeds v 0 . For all three potentials, we see an abrupt change for small q from approaching a finite value S(0) to a power law, whereby the latter signals phase separation. In Fig. 3(b), the force coefficient ζ is plotted as a function of the speed v 0 , where open symbols correspond to the homogeneous, and closed symbols to the phase separated suspension as determined from the behavior of the structure factor at small q. Since translational noise is neglected, v * 0 is very small and the suspension clusters almost immediately. Going to larger speeds v 0 , we observe a reentrance into the fluid phase. Finally, Fig. 3(c) shows a snapshot for the harmonic potential (H) at v 0 = 0.2 in the instability region showing a single large cluster surrounded by a dilute fluid phase. In contrast, in Fig. 3(d) at higher speeds, the cluster has dissolved and the suspension is homogeneous again. Hence, beside the ordering transition, also a reentrance into the . Every particle is colored according to Eq. (25) with ∆t = 25 quantifying the persistence of particle motion with respect to the initial particle orientation.
disordered fluid at higher propulsion speeds can be observed. The color code corresponds to the quantity measuring the persistence of particle motion over the time interval ∆t with respect to the initial orientation. For ∆t 1/D r we expect for a free particle α k ∼ 1 while for α k < 0 a particle has actually moved against the propelling force. We find that in dense regions particles are indeed less likely to move along their orientations.
Discussion. -What is the origin of the dynamical instability? The emerging microscopic picture is the following: particles collide and, due to the persistence of their swimming motion, they block each other forming a small cluster [14,15,17]. Collisions with other particles lead to a growth of the cluster. For a particle situated in the rim of the cluster to become free, it has to wait a time ∼ 1/D r for its orientation to point outward. Depending on the ratio between this waiting time and the collision time controlled by the propulsion speed, there is either a stable steady state corresponding to a homogeneous suspension with dynamically evolving clusters, or clusters grow until the largest cluster occupies a finite fraction of the system. This microscopic picture is confirmed by the effective large scale evolution equations. The relaxation of a density fluctuation is dominated by the flux p along the local mean orientation, which is given by p ∼ −∇(vρ). Note that the "potential" is the product of swimming speed v and density ρ. Now image a small region in which the density is increased but at the same time the swimming speed is decreased. The density of the surrounding region is smaller but the speed of particles is higher such that the product vρ might be larger than that of the dense region. In this case there is a flux towards the denser region and phase separation sets in.
Conclusions. -We have studied a dynamical instability in a minimal model for self-propelled repulsive disks for densities below the freezing transition [29]. The instability results in phase separation and intrinsically requires the system to be driven away from thermal equilibrium. It cannot be modelled as an effective isotropic attraction that shifts an equilibrium coexistence line [30]. We have shown that instead the relevant mechanism is a force imbalance due to the propulsion, which implies an effective swimming speed that depends on the density.
Starting from the Smoluchowski equation, we have derived a closed equation of motion for the tagged particle density by decomposing the mean force due to the surrounding particles into two independent contributions: along the particle orientation and along the density gradient. We have employed Brownian dynamics simulations to calculate the two free parameters entering this equation: the force coefficient ζ and the passive long-time diffusion coefficient D. Using these values as input, our theory is able to reproduce the single-particle long-time diffusion coefficient in the homogeneous suspension and to predict the onset of the instability. For different potentials we have confirmed numerically that there is also a reentrance into the disordered fluid at sufficiently high propulsion.
Here we have studied a two dimensional system but our approach also extends to three dimensions. For the sake of clarity and brevity, we have neglected hydrodynamic interactions. At least pair-wise hydrodynamic interactions are easily incorporated in our framework and will change the coefficients D, D r , and ζ. However, the qualitative predictions, and in particular the instability diagram Fig. 1(b), will remain unchanged. Models that, e.g., imply a reorientation of particles during collisions instead of a persistent motion (see, e.g., Ref. 31) might move out of the instability region either because v * is increased or the force coefficient ζ is reduced. While we have used computer simulations in order to calculate ζ, as a next step it will be desirable to go to the next level of the BBGKY hierarchy in order to calculate the anisotropic pair distribution, and therefore the force coefficient ζ, from first principles. * * * This work has been supported financially by the DFG within SFB TR6 (project C3) and by the ERC (advanced grant INTERCOCOS under project number 267499).