An equation of state for active matter

We characterise the steady states of a suspension of two-dimensional active brownian particles (ABPs). We calculate the steady-state probability distribution to lowest order in Peclet number. We show that macroscopic quantities can be calculated in analogous way to equilibrium systems using this probability distribution. We then derive expressions for the macroscopic pressure and position-orientation correlation functions. We check our results by direct comparison with extensive numerical simulations. A key finding is the importance of many-body effective interactions even at very low densities.

Statistical mechanics tells us the probability that a system is in a certain state, as long as it is at equilibrium [1]. For instance, the degrees of freedom for any system in the canonical ensemble are sampled from the Gibbs-Boltzmann distribution, P eq ∝ e −U/(k B T ) (where U is the system's energy). P eq can then be used to determine macroscopic system properties, effectively reducing a large proportion of equilibrium statistical mechanics to a difficult exercise in integration. In contrast, the steady-state probability distribution of a non-equilibrium system in a non-equilibrium steady-state (NESS), P ss , is not known a-priori [2]. In fact, a NESS need not always exist. However, if a NESS exists and P ss can be determined, one should in principle be able to apply similar techniques to equilibrium statistical mechanics on nonequilibrium systems in a NESS. Herein lies the motivation of this letter: calculating macroscopic properties of non-equilibrium systems in a NESS using methods analogous to those of equilibrium statistical mechanics.
A popular subset of non-equilibrium phenomena, which we restrict our attention to here, are "active" matter systems [3,4], which never reach an equilibrium state due to the presence of internal (bulk) energy sources. Early studies of these active systems were motivated by biological processes on a wide array of length scales (from e.g. fluctuations of membranes [5] to flocking birds [6][7][8]), but has increasingly found application in synthetic man-made systems [9,10].
The active system we study here, known as Active Brownian Particles (ABPs) [11][12][13][14][15], is a suspension of spherical colloids in a solvent, which have a mechanism of self-propulsion driving them out of equilibrium [16]. ABPs are widely studied due to the relative simplicity of their equations of motion, while capturing the essentials of active matter systems.
Typically theoretical descriptions of active systems (such as ABPs) approximate the many-particle dynamics by an effective one-particle distribution function. For some macroscopic quantities, these single-particle distribution functions give qualitatively similar behaviour to explicit particle based simulations [12,17]. However quantitative comparison is only possible using (often many) phenomenological fitting parameters. Other macroscopic quantities (e.g. the pressure of an interacting active gas) which depend on particle correlations, are simply impossible to obtain using these approaches.
Theoretical frameworks which take account of the many-particle correlations, on the other hand, look formidably difficult. However recently, a generic approach to fluctuating dynamical systems highlighting the role of non-vanishing currents in non-equilibrium steadystates (NESS) was introduced, where the relaxation dynamics towards the NESS plays a key role [18]. In particular, [18] introduces the notion of "typical trajectories", which we use here to construct an ansatz for P ss amenable to analytical solution.
This letter has three main results. The first result is an explicit expression for the steady-state probability distribution of ABPs, P ss in 2D. The second result is an effective interaction potential for ABPs which is independent of particle orientation, but depends on many-body (as opposed to pair-wise) interactions. The third result uses P ss to obtain an equation of state, an expression for the active brownian swim pressure [19] at low Peclet number. The van der Waals equation is a landmark in statistical mechanics generalising the ideal gas law to an equation of state for realistic gases, expressing the pressure as an expansion in density [1]. Here we obtain an equivalent for active matter by obtaining an equation of state that is an expansion in both density and activity (here encoded in Peclet number). We also compute another macroscopic average, which probes local correlations between interparticle position and orientation. Comparison of our calculated swim pressure and the local correlation function to direct simulations yields good agreement between the two for a range of system densities.
We consider a collection of N active brownian particles on a plane of area A = L 2 with periodic boundary conditions (BCs). The 3N degrees of freedom satisfy the system of stochastic differential equations for their positions on the plane, r i , and their directions u i = (cos θ i , sin θ i ), of self-propulsion where W α are as usual Wiener processes. To be concrete, the interaction potential U is the sum of pair-wise (purely repulsive) Weeks-Chandler-Anderson (WCA) potentials V(r ij ) with energy scale . Length is measured in arXiv:2201.10813v2 [cond-mat.stat-mech] 8 Mar 2022 units of the interaction length σ, time in units of σ 2 /D t , and energy in units of D t γ, where D t and γ are the translational diffusion and drag, respectively. The remaining parameters in our system quantify the strength of the active force relative to fluctuations, f P ≡f P σ/D t γ, the relative magnitude of rotational diffusion, D r ≡D r σ 2 /D t , and the density of particles in our system, ρ ≡ N σ 2 /A. Our default parameters in this letter are /(D t γ) = 1 and D r = 3 (the latter corresponding to a no-slip boundary condition between the particles and the solvent in equilibrium systems). The eqns. (1) are equivalent to a Fokker-Planck equation for the probability density, P ss ( X, t) of degrees of freedom X = (r 1 , θ 1 , r 2 , θ 2 , · · · , r N , θ N ), subject to periodic boundary conditions. The steadystate probability density P ss ( X) is a solution of (∂ t P ss = 0). We define a new function ζ( X) via P ss ( X) ∝ e −U ( X)−f P ζ( X) , so that the equilibrium distribution is recovered when f P = 0. Inserting this into this Fokker-Planck equation maps the linear partial differential equation (PDE) for P ss ( X) to a non-linear PDE for ζ( X). This is helpful mathematically because ζ is less constrained than P and easier to approximate. To calculate ζ, we make some assumptions about its form. Symmetry under exchange of particles, translational and rotational invariance imply where r ij = |r ij |, r ij = r j − r i , u ij = u j − u i and w(r) is a scalar function (see Figure 1A). We neglect terms that are nonlinear in orientation, assuming local alignment interactions are weak. Finally, keeping terms to leading order in f P , one arrives at the ordinary differential equation (ODE) for w(r) where and ξ fluctuates around zero, hence we set ξ = 0 (see Methods). This ODE is solvable using standard techniques [20][21][22]. Thus, the problem of solving a 3N variable PDE for P ss ( X) in eqn 2 has been reduced to solving a single ODE. To solve for w(r) (and hence ζ( X)), we need two BC which we implement approximately as follows. First, we note that to be consistent with periodic BC's the range of w(r) (like that of V(r)) must be less than L/2, i.e. w(r) = 0, r > r c with r c < L/2 [23]. As we consider large systems, L 2 1/6 , we relax this and simply require that w(r) → 0 as r → ∞. Next, we note that the typical trajectories of a system in steadystate will generate (spatial) probability density currents J i = P ss V i [18] with These currents depend on P ss and specify a family of trajectories, i.e. V i . This implies that we can reverse this logic and use the typical trajectories to obtain P ss in principle. However since we need only one more BC, we only need to analyse a single trajectory. For this, we note that a trimer of three ABPs (labelled 1, 2, 3) in an equilateral triangle configuration (i.e. all separated by the zero-force radius r 0 = 2 1/6 + O(f P )) is meta-stable if their respective self-propulsion forces are directed to the centre of the trimer. This implies that the relative velocities of all three particles should be zero, i.e. V ij = V j − V i , then V 12 = V 13 = V 23 = 0 leading to the boundary condition 2 1/6 w (2 1/6 ) + w(2 1/6 ) = −1/3.
. While based on sound physical principles, the approximations we have applied are mathematically uncontrolled, and hence must be checked. We do this by empirical comparison with direct numerical simulations below.
We can now calculate macroscopic properties using steady-state distribution, P ss . The expectation value of any observable O is where h ij = V(r ij ) + f P w(r ij )r ij · u ij , and the normalisation constant (the "partition function") is This is our first main result. Furthermore, one can explicitly integrate out orientational degrees of freedom in P ss [24] to arrive at the marginal distribution Q ss (r 1 , . . . , r N ) ∼ e −U eff ({ri}) with an effective potential dependent only on particle positions [25], This effective potential constitutes the second main result of this letter. The self-propulsion f P = 0 introduces a minimum in the effective interaction [24]. It also demonstrates the importance of three-body and higher-body terms, which arise due to the coupling between position and orientation degrees of freedom. We note that the three-body effective interactions have also been observed in three dimensional simulations of ABPs [26] and an effective potential of the Active-Ornstein-Uhlenbeck model [14]. To obtain the average of a macroscopic observable where O is independent of orientational degrees of freedom, one may do so using e −U eff as the probability measure.
Next we use eqn. (8), to calculate the interacting swim pressure of ABPs [19,27]. We first remind the reader that the equation of state of a 2D fluid with temperature T at equilibrium whose particles interact via a pair potential, V(r) is p = k B T ρ + p v , with the virial interaction pressure, p v = − 1 4 ρ 2 ∞ 0 2πr 2 V (r)g (2) (r)dr [28] where g (2) (r) is the radial distribution function. We aim to obtain an equivalent for ABPs. Starting from the Kirkwood definition of the stress tensor [29], the total pressure for ABPs may split up as p = ρ 1 + f 2 P 2Dr + p v + p int s , with the first term corresponding to the active ideal gas pressure (since D t γ t = 1 in our units), the second is the virial interaction pressure p v [24] (present in both passive and active systems), and the last term is the interacting swim pressure p int s (which is non-zero only in the active case). The microscopic definition of p int s for the WCA potential V(r ij ) is [19] p int Using eqn. (8), we arrive at the expression for the steadystate swim pressure where a 2 (ρ) ≡ r 2 w(r)V (r)g (2) 0 (r)dr, (13a) are density dependent functions that are obtained from the passive (f P = 0) two-body: g (2) 0 (r) and three-body: G In Figure 1B, we compare eqn. (12), (black solid line) to p int s directly measured from simulation. We scale the theory by its value, p th 0 at ρ = 0.01 and the simulations by p sim 0 = p int s ρ 2 f 2 P (f P = 0.25, ρ = 0.01). With the scaling from eqn. (12), the simulation data does indeed collapse onto a a single curve in agreement with our calculation (the collapse is less noisy for higher densities due to more frequent collisions). We find p sim 0 /p th 0 ≈ 1.2. We present simulation results for four different active force values, f P = 0.25, 0.50, 1.0, 2.0, as a function of density. We expect the theory to begin to deviate from simulations for f P > 1. The theoretical prediction of cubic dependence in density is surprising as both a 2 (ρ) and ρa 3 (ρ) are density dependent and so any higher dependence on ρ must cancel. In fact, keeping only pair correlations gives the red dashed line which shows completely different behaviour. Clearly, three-body correlations control the dependence of p int s on density. Finally at large ρ the theory curve diverges. We expect this is due to higherbody correlations which are in principle calculable in our framework.
Finally, we compute the local correlation function defined as Physically, this quantity represents angle-averaged correlations between inter-particle orientation and interparticle displacement, for a given spacing r between pairs of particles. It is interesting to measure this function as it probes local structure of the ABPs, and so gives a more detailed picture of whether our calculation captures the mesoscale structure as well as the thermodynamic behaviour of the system. We find C 1 also depends on twobody and three-body terms, 0 (r, s)ds + . . . (15) In Figure 2, we measure C 1 (r) in simulation and compare our results to those predicted by eqn. (15). Again we find that with two-body terms only that the theoretical calculation does not agree with simulations and only compares well (for all r) on inclusion of the three-body terms (without fitting parameters).
In conclusion, we have presented the first ab-initio cal-culation of the many-body steady-state probability density of states of Active Brownian Particles. This is based on a sequence of approximations that we check by numerical simulations. It is promising that despite the approximations the calculation has managed to capture the behaviour of a number of local and global observables of the system without any free parameters. This indicates that these approximations are founded on sound physical principles and capture the essential features of the non-equilibrium steady-state of ABPs at low f P . It also suggests that they can form the starting point for a more detailed theory of these kinds of system that can be systematically improved by tightening the approximations and the assumptions behind them. The