Velocity distribution in active particles systems

We derive an analytic expression for the distribution of velocities of multiple interacting active particles which we test by numerical simulations. In clear contrast with equilibrium we find that the velocities are coupled to positions. Our model shows that, even for two particles only, the individual velocities display a variance depending on the interparticle separation and the emergence of correlations between the velocities of the particles. When considering systems composed of many particles we find an analytic expression connecting the overall velocity variance to density, at the mean-field level, and to the pair distribution function valid in the limit of small noise correlation times. Finally we discuss the intriguing analogies and main differences between our effective free energy functional and the theoretical scenario proposed so far for phase-separating active particles.

A simple derivation of the stationary velocity distribution of a collection of active particles subject to non equilibrium colored noise is provided in this note. Moreover simulation results for the 2dimensional case are reported.

I. DERIVATION OF THE STATIONARY VELOCITY DISTRIBUTION
The eective dynamics for space coordinates of an assembly of active spheres 1,2 iṡ where the velocities u i evolve according to the law: The force F i = −∇ i U acting on the i-th particle is conservative and associated to the potential U(r 1 , . . . , r N ), γ is the drag coecient, whereas the stochastic vectors η i (t) are Gaussian and Markovian processes distributed with zero mean and moments η i (t)η j (t ) = 2δ ij δ(t − t ). where d is the spatial dimensionality. The coecient D a due to the activity is related to the correlation of the Ornstein-Uhlenbeck process u i (t) via where d is the spatial dimension. To simplify the notation we switch from r i to an array x i as done in previous publications:ẋ Dierentiate againẍ Let us introduce the variable v i and recast eq.(5) as: A. One particle Before presenting the multidimensional result, we digress to illustrate the kinetic method of solution in a simple one-dimensional case. We begin with a single particle in one dimension and drop the index i. We dierentiate eq. (1) with respect to time and introduce the velocity variable v =ẋ so that instead of the original system (1) and (2) we have:v We obtain the Kramers equation for the phase-space distribution f (x,ẋ; t): with Γ(x) = 1 + τ γ ∇ 2 U(x)). By multiplying and integrating over v eq. (9) and considering only time independent solutions f 0 (x, v) one obtains: Such an integro-dierential equation can be solved by the following Γ(x)v 2 whose width depends on the particle position and the average velocity v vanishes. After substituting the factorization f 0 into (9) and evaluating the velocity varianceˆd we arrive at the following dierential equation determining the steady state coordinate distribution P(x) : which is identical to the dierential equation determining the stationary coordinate distribution in the unied color noise approximation (UCNA): with and T s = D a γ and Z 1 a normalization constant dened as: The method can be easily extended to the multidimensional case.

MANY PARTICLE SYSTEMS
We now generalize the kinetic argument above illustrated to a many particle system and write the following multidimensional Kramers equation 3 describing the evolution of the phase-space distribution with the non dimensional friction matrix Γ ik dened as The probability conservation and the momentum balance equations are straightforwardly obtained by projection, i.e. by integrating eq.(16) over the dN dimensional velocity space after multiplying by 1 and v i , respectively: where P N ({x i }, t) is the marginalized distribution giving the distribution of positions of the particles regardless their velocities: Equations (18)-(19) are an exact consequence of (16), but require the knowledge of p ij which can be obtained by continuing the projection procedure in velocity space to higher order in v i and closing the hierarchy by a suitable truncation ansatz, such as the introduction of phenomenological transport coecients. On the other hand, if we limit ourselves to study the stationary regime the approach is fruitful and leads to a dierent interpretation of the MUCNA equations. Let us rst, notice that the operator featuring in the right hand side of (16) posseses a null eigenvalue whose associated eigenfunction (the ground state) is the multivariate Gaussian velocity distribution : In other words Π({v i }|{x i }) is the conditional probability distribution of velocities given the positions take on the values {x i }. We now construct a time-independent trial phase-space distribution having a factorized form: and corresponding to zero current J i (x 1 , . . . , x N ). For generic potentials the velocity moments of G are non constant in space, but depend on coordinates x i and dierent velocity components may be correlated and are given by the formula: so that the "pressure " tensor reads Using eq. (19) in conjunction with the condition of vanishing currents, J i = 0, we have the following balance equations: Eliminating the pressure tensor using (25) we arrive at the dierential equation determining P st ({x i }) : Finally, by multiplying by the matrix Γ equation (27) we arrive at an equation identical to the MUCNA equation (see in the case D t = 0, so that we may identify the congurational part of the trial solution P st with P N derived by the MUCNA method 2 : where Z N is a normalization constant.
We turn now to the formula for the velocity covariance matrix for many particles. Using the Gaussianity of the distribution we can immediately write the averages of the velocity as where · stands for the double average over velocities and over space since the velocity distribution depends on the positions, unlike the equilibrium case.
In the general case we must evaluate the matrix elements Γ ij . As shown in the appendix this can be done to rst order in τ /γ with the result for the velocity self-correlations: Finally, averaging over the whole system and introducing the positional pair correlation function g(r − r ) one obtains In analogy with granular gases we can dene the second moment of the velocity distribution to be the average kinetic temperature, T k of the system, via

III. APPROXIMATION FOR THE DETERMINANT AND VELOCITY CORRELATIONS
The exact evaluation of the determinant Γ associated with the Hessian matrix is beyond the authors capabilities and we look for approximations in order to evaluate the eective forces. We consider the associated determinant in the case of two spatial dimensions and vanishing external potential: It is interesting to remark that the o-diagonal elements contain only one term, while the diagonal elements and their neighbors contain N terms. Thus in the limit of N → ∞ we expect that the matrix becomes eectively diagonal. However, even with such a limit, we see that to order τ /γ the inverse matrix, which is directly related to velocity correlator v αi v βj (see (25)), reads: We, now, assume translational invariance and perform the average over the spatial distribution P N (r 1 , . . . , r N ), and introduce the pair correlation as P 2 (r 1 , r 2 ) = ρ 2 N (N −1) g(r 1 − r 2 ) . The sought average is: . One can see that in the limit of large N the o-diagonal terms become neglegible because they carry factors 1/N and one obtains a block diagonal matrix whose blocks are 2 × 2 matrices: Evaluation of the ensemble average of the determinant to linear order in τ /γ and velocity variance Notice that to linear order τ /γ the ensemble average of the determinant is: To order τ /γ: We simulate a 2d system driven by GCN and composed by N = 1000 particles with periodic boundary conditions. The particles interact via the purely repulsive pair potential ϕ(r 1 − r 2 ) = |r 1 − r 2 | −12 . We x τ = 1 and we explore several values of density ρ = N/L 2 and diusivity D. The results are reported in Fig. 1. Fig. 1(a) shows the normalized velocity variance as a function of density for dierent values of D. It is seen that, as in the 1d case, the quantity ẋ 2 /(D/τ ) is a decreasing function of ρ. However in the 1d case ẋ 2 /(D/τ ) is also a decreasing function of D (at xed ρ), while in 2d a non-monotonic behavior in D can be observed at high ρ. Fig. 1(b) shows the amplitude of the Fourier-transformed density uctuations at low q (q ≈ (20σ) −1 ) which shows a very similar behaviour to the 1d case. Fig. 1(c) and (d) show two snapshots at high D and intermediate densities where evident clustering of the particles is observed but not a full phase separation. This is the case also upon further increasing τ to very high values as shown in Fig. 2.