Statistical Mechanics and Hydrodynamics of Self-Propelled Hard Spheres

TL;DR

This paper develops a statistical mechanics framework for self-propelled hard spheres, deriving key properties like pressure and phase separation thresholds, validated by simulations, advancing understanding of active matter systems.

Contribution

It introduces a novel kinetic theory approach to model interacting active particles with excluded volume, linking microscopic parameters to macroscopic behavior.

Findings

Derived steady state pressure of active fluid from microscopic parameters.

Identified critical density for motility-induced phase separation.

Validated theoretical results with overdamped simulations.

Abstract

Starting from a microscopic model of self-propelled hard spheres we use tools of non-equilibrium statistical mechanics and the kinetic theory of hard spheres to derive a Smoluchowski equation for interacting Active Brownian particles. We illustrate the utility of the statistical mechanics framework developed with two applications. First, we derive the steady state pressure of the hard sphere active fluid in terms of the microscopic parameters and second, we identify the critical density for the onset of motility-induced phase separation in this system. We show that both these quantities agree well with overdamped simulations of active Brownian particles with excluded volume interactions given by steeply repulsive potentials. The results presented here can be used to incorporate excluded volume effects in diverse models of self-propelled particles.

Equations

$$\frac{\partial {\mathbf{v}}_i}{\partial t}=\sum _{j\ne i}T(i,j){\mathbf{v}}_i+F{\hat{\mathbf{u}}}_i-\nabla V(\mathbf{r})-\xi {\mathbf{v}}_i+{\mathbf{\eta }}_i(t)$$

$$\frac{\partial {\omega }_i}{\partial t}=-{\xi }_R{\omega }_i+{\eta }_i^R(t).$$

$$T(i,j)=\sigma \int d\hat{\mathbf{\sigma }}\Theta (-{\mathbf{V}}_{ij}\cdot \hat{\mathbf{\sigma }})|{\mathbf{V}}_{ij}\cdot \hat{\mathbf{\sigma }}|\delta ({\mathbf{r}}_{ij}-\mathbf{\sigma })(b_{ij}-1)$$

$${\partial }_{t}f(\mathbf{r},\theta ,\mathbf{v},\omega ,t)+\mathcal{D}f(\mathbf{r},\theta ,\mathbf{v},\omega ,t)=0$$

$$\mathcal{D}=\mathbf{v}\cdot {\nabla }_{\mathbf{r}}+\omega {\partial }_{\theta }+F\hat{\mathbf{u}}\cdot {\nabla }_{\mathbf{v}}-{\nabla }_{\mathbf{r}}V(\mathbf{r})\cdot {\nabla }_{\mathbf{v}}-\xi {\nabla }_{\mathbf{v}}\cdot \mathbf{v}-{\xi }_R{\partial }_{\omega }\omega -k_{B}T\xi {\nabla }_{\mathbf{v}}^2-k_B{T}_R{\xi }_R{\partial }_{\omega }^2.$$

$${\partial }_{t}c+{\nabla }_{{\mathbf{r}}_1}\cdot {\mathbf{J}}^T+{\partial }_{\theta }{J}^R=0.$$

$$\displaystyle\lim_{t\gg 1/\xi}\mathbf{J}^{T}=-\frac{1}{\xi}\big{[}-F\hat{\mathbf{u}}c+c\nabla V(\mathbf{r})+\nabla\cdot\langle\mathbf{v}\otimes\mathbf{v}\rangle+\partial_{\theta}\langle\omega\mathbf{v}\rangle\big{]}$$

$$\displaystyle\lim_{t\gg 1/\xi_{R}}J^{R}=-\frac{1}{\xi_{R}}\big{[}\nabla\cdot\langle\omega\mathbf{v}\rangle+\partial_{\theta}\langle\omega^{2}\rangle\big{]}$$

$$f(\mathbf{r},\theta ,\mathbf{v},\omega ,t)=\frac{1}{2\pi k_{B}T}\sqrt{\frac{1}{2\pi k_B{T}_R}}c(\mathbf{r},\theta ,t)e^{-\frac{(\mathbf{v}-v_0\hat{\mathbf{u}}{)}^2}{2k_{B}T}}{e}^{-\frac{{\omega }^2}{2k_B{T}_R}}$$

$$\partial_{t}c+v_{0}\hat{\mathbf{u}}\cdot\nabla c-\xi^{-1}\nabla\cdot\big{[}c\nabla V(\mathbf{r})\big{]}-\nabla\cdot\tensor{\mathbf{D}}\cdot\nabla c-D_{r}\partial_{\theta}^{2}c=0$$

$${\partial }_{t}c+{\nabla }_{\mathbf{r}}\cdot {\mathbf{J}}^T-D_r{\partial }_{\theta }^{2}c=0$$

$$\mathbf{J}^{T}=v_{0}\hat{\mathbf{u}}c-\xi^{-1}c\nabla V(\mathbf{r})-\tensor{\mathbf{D}}\cdot\nabla c+\xi^{-1}\big{(}\mathbf{I}^{thermal}+\mathbf{I}^{sp}+\mathbf{I}^{cross}\big{)}.$$

$$I_{\alpha }^{thermal}=\frac{4{\pi }^2\sigma (k_{B}T{)}^3}{(2\pi k_{B}T+v_0^2{)}^2}{\int }_{\hat{\mathbf{\sigma }},{\theta }_2}{c}^{(2)}({\mathbf{r}}_1,{\theta }_1,{\mathbf{r}}_1-\mathbf{\sigma },{\theta }_2,t){\hat{\sigma }}_{\alpha }$$

$$I_{\alpha }^{sp}=\frac{\sigma v_0^6}{(2\pi k_{B}T+v_0^2{)}^2}{\int }_{\hat{\mathbf{\sigma }},{\theta }_2}\Theta (-{\hat{\mathbf{u}}}_{12}\cdot \hat{\mathbf{\sigma }})({\hat{\mathbf{u}}}_{12}\cdot \hat{\mathbf{\sigma }}{)}^2{c}^{(2)}({\mathbf{r}}_1,{\theta }_1,{\mathbf{r}}_1-\mathbf{\sigma },{\theta }_2,t){\hat{\sigma }}_{\alpha }$$

$$I_{\alpha }^{cross}=\frac{\sigma }{(2\pi k_{B}T+v_0^2{)}^2}{\int }_{\hat{\mathbf{\sigma }},{\theta }_2}{c}^{(2)}({\mathbf{r}}_1,{\theta }_1,{\mathbf{r}}_1-\mathbf{\sigma },{\theta }_2,t)$$

$$\times [2\pi (k_{B}T{)}^2{v}_0^2+\pi k_{B}T{v}_0^4(({\hat{\mathbf{u}}}_1\cdot \hat{\mathbf{\sigma }}{)}^2+({\hat{\mathbf{u}}}_2\cdot \hat{\mathbf{\sigma }}{)}^2)-\sqrt{8\pi (k_{B}T{)}^3}{v}_0^3{\hat{\mathbf{u}}}_{12}\cdot \hat{\mathbf{\sigma }}]{\hat{\sigma }}_{\alpha }$$

$$P={\int }_0^{\infty }\rho (x){\partial }_{x}V(x)dx.$$

$$P=P_{0}+\int d\theta\bigg{[}\frac{v_{0}}{D_{r}}\hat{u}_{x}(I_{x}^{thermal}+I_{x}^{sp}+I_{x}^{cross})\bigg{|}_{x=0}+\int_{0}^{\infty}(I_{x}^{thermal}+I_{x}^{sp}+I_{x}^{cross})dx\bigg{]}.$$

$$c^{(2)}(\mathbf{r}_{1},\theta_{1},\mathbf{r}_{2},\theta_{2},t)=g(\mathbf{r}_{1},\mathbf{r}_{2}\big{|}\rho(\mathbf{r},t))c(\mathbf{r}_{1},\theta_{1},t)c(\mathbf{r}_{2},\theta_{2},t),$$

$$g\left(\mathbf{r}_{1},\mathbf{r}_{2}\big{|}\rho(\mathbf{r},t)\right)=g(\sigma|\rho)=\frac{1-\frac{7}{16}\phi}{(1-\phi)^{2}},$$

$$P=P_0+{\lambda }_{D}g(\sigma ){\rho }_0^2-g(\sigma )\frac{\xi v_0}{D_r}{\lambda }_I{\rho }_0^2.$$

$$P=k_{B}T\rho_{0}\bigg{(}1+\frac{\pi\sigma^{2}g(\sigma)}{2}\rho_{0}\bigg{)}$$

$$P=\frac{v_0^2}{2D_r{\xi }^{-1}}{\rho }_0+\frac{v_0^2}{2}{\rho }_0+\frac{3\pi }{16}{v}_o^2{\sigma }^{2}g(\sigma ){\rho }_0^2-\frac{2}{3}\frac{v_0^3\sigma }{D_r}g(\sigma ){\rho }_0^2$$

$$\rho (\mathbf{r},t)=\int d\theta c(\mathbf{r},\theta ,t).$$

$$P_{\alpha }(\mathbf{r},t)=\int d\theta {\hat{u}}_{\alpha }c(\mathbf{r},\theta ,t).$$

$${\partial }_t\rho =-\nabla \cdot \mathbf{J}$$

$${\partial }_t{P}_{\beta }=-{\partial }_{\alpha }{J}_{\alpha \beta }-D_r{P}_{\beta }$$

$$\left(\begin{array}{c}{J}_{\alpha }\\ J_{\alpha \beta }\end{array}\right)=\int d\theta \left(\begin{array}{c}1\\ {\hat{u}}_{\beta }\end{array}\right)J_{\alpha }^T.$$

$$c(\mathbf{r}-\mathbf{\sigma },{\theta }_2,t)=c(\mathbf{r},{\theta }_2,t)-{\sigma }_{\alpha }{\partial }_{\alpha }c(\mathbf{r},{\theta }_2,t)+\dots$$

$${\partial }_t\rho +v_0\nabla \cdot \mathbf{P}-\nabla \cdot [\mathcal{D}(\rho )\nabla \rho ]+D_p\nabla \nabla :\mathbf{PP}=0$$

$$\partial_{t}P_{\gamma}+D_{r}P_{\gamma}+\nabla_{\gamma}\mathcal{P}(\rho)-D_{\perp}\nabla^{2}P_{\gamma}+2\lambda_{I}\big{[}P_{\gamma}(\nabla\cdot\mathbf{P})+(\mathbf{P}\cdot\nabla)P_{\gamma}\big{]}$$

$$=$$

Full text

I Introduction

In recent years, the study of active materials has been on the forefront of research in soft condensed matter physics. A particular class of active materials is one composed of self propelled particles that convert energy from a local bath into persistent motion. Self propelled particles describe systems on many length scales, ranging from bacteria Cates (2012); Berg (2001), synthetic colloids F. et al. (2015); Hong et al. (2007); Palacci et al. (2013), vibrated granular rods Narayan et al. (2007), to schools of fish Marchetti et al. (2013); Ballerini et al. (2008). Despite the fact that there exists continual energy consumption and dissipation at the level of individual particles, collections of self-propelled particles form large-scale, stable, coherent structures Ramaswamy (2010). Phenomena exhibited include athermal phase separation among purely repulsive particles Buttinoni et al. (2013); Cates and Tailleur (2015); Redner et al. (2013); Putzig and Baskaran (2014), anomalous mechanical properties Yang et al. (2014); Takatori et al. (2014); Solon et al. (2015a); Mallory et al. (2014), emergent structures and pattern formation Hagan and Baskaran (2016); Gopinath et al. (2012); Farrell et al. (2012).

Theoretical progress in understanding active materials has been propelled by the study of minimal models for self-propelled particles. The most widely studied of these minimal models is known as an active Brownian particle (ABP). ABPs travel with constant velocity along their body axis and the orientation of their body axis changes through ordinary rotational diffusion. While numerical investigations of ABPs abound in the literature and have led to a significant understanding of the dynamics of this system, analytical progress for interacting particles has been difficult. The reason that analytical descriptions of the collective behavior remain elusive is that interactions have a finite collision time, which in turn would result in intractable many body effects. To overcome this challenge, authors have taken phenomenological approaches such as incorporating the effect of interactions in a density dependence of the self propulsion speed Cates and Tailleur (2015, 2013); Cates (2012); Tailleur and Cates (2009); Solon et al. (2015b) and a mean field version of dynamical density functional theory Bialke et al. (2013); Speck et al. (2014, 2015).

In this paper, we aim to capture these many body effects by applying the well developed theoretical framework of hard-sphere liquids van Beijeren and Ernst (1979); Hansen and McDonald (1990) to a fluid consisting of ABPs. Starting with underdamped Langevin equations to describe the dynamics of self-propelled particles that interact with hard core instantaneous collisions, we systematically coarse grain the microdynamics to obtain a description applicable on longer length and time scales after which we take the limit of large friction to obtain an effective description for an overdamped systems of ABPs. The primary theoretical result in this work is a first principles derivation of the statistical mechanics of self-propelled hard spheres. To illustrate the utility of this result, we compute two well studied emergent properties in fluids of ABPs: the mechanical pressure and the phase boundary that determines the onset of athermal phase separation into a dense liquid and a dilute gas. We then compare our results to those in the literature obtained from numerical studies and more phenomenological theoretical approaches and use this comparison to place the context of our work within the existing body of work.

II Theoretical Framework

II.1 Microdynamics

We consider self propelled hard disks in two dimensions with particle diameter $\sigma$, unit mass, and moment of inertia $I$. The microstate of $N$ such hard disks is given by $\Gamma=\{\mathbf{r}_{1}(t),\dots,\mathbf{r}_{n}(t),\mathbf{v}_{1}(t),\dots,\mathbf{v}_{n}(t),\omega_{1}(t),\dots,\omega_{n}(t),\theta_{1}(t),\dots,\theta_{n}(t)\}$. Here the variables ($\mathbf{r}_{i}$, $\mathbf{v}_{i}$, $\theta_{i}$, $\omega_{i}$) are the position, velocity, orientation, and angular velocity respectively of the $i$ th particle. The equations of motion in this case are given by $\partial_{t}\mathbf{r}_{i}=\mathbf{v}_{i}$ and $\partial_{t}\theta_{i}=\omega_{i}$ where the linear and angular velocities evolve according to

[TABLE]

[TABLE]

Here, the unit vector $\hat{\mathbf{u}}=(\cos(\theta),\sin(\theta))$ is the orientation of the particle’s body axis along which a propulsive force of strength $F$ acts and $V(\mathbf{r})$ is some external potential. The binary collision operator $T(i,j)$ generates the instantaneous linear momentum transfer between disks at contact and is given by van Beijeren and Ernst (1979)

[TABLE]

where $\hat{\bm{\sigma}}$ is the unit normal at the point of contact directed from disk $j$ to disk $i$ and $\mathbf{V}_{ij}=\mathbf{v}_{i}-\mathbf{v}_{j}$. The operator $b_{ij}$ replaces pre-collisional velocities with post-collisional velocities, e.g., $b_{12}\mathbf{v}_{1}=\mathbf{v}_{1}-(\mathbf{V}_{12}\cdot\hat{\bm{\sigma}})\hat{\bm{\sigma}}$   for collisions which conserve energy and momentum. The spatial delta function ensures that particles are in contact. The prefactors that depend on the relative velocities $\mathbf{V}_{ij}$ ensures that the incoming flux of colliding particles is taken into account correctly. The random forces $\bm{\eta}_{i}(t)$ and $\eta_{i}^{R}(t)$ are Gaussian white noise variables with correlations given by $\langle\eta_{i}^{\alpha}(t)\eta_{j}^{\beta}(t^{\prime})\rangle=2k_{B}T\xi\delta_{\alpha\beta}\delta_{ij}\delta(t-t^{\prime})$ and $\langle\eta_{i}^{R}(t)\eta_{j}^{R}(t^{\prime})\rangle=2k_{B}T_{R}\xi_{R}\delta_{ij}\delta(t-t^{\prime})$ respectively. Latin indices label the particle number, while Greek indices label the vector components of the noise. The noise amplitudes depend on parameters $T,T_{R}$ that have the units of temperature. These parameters need not be the same as $T_{R}$ is an intrinsic quantity describing the reorientation of the active drive.

Thus, Eqs.(1-2) are the Langevin equations for $N$ interacting self-propelled hard disks. For a single particle in the high friction limit these equations would reduce to the well studied Redner et al. (2013); Fily and Marchetti (2012); Cates (2012); Romanczuk et al. (2012); Cates and Tailleur (2015); Marchetti et al. (2013); Yang et al. (2014); Takatori et al. (2014); Solon et al. (2015a); Mallory et al. (2014); Farage et al. (2015); Marini Bettolo Marconi and Maggi (2015); Cates and Tailleur (2013); Solon et al. (2015b); Bialke et al. (2013); Tailleur and Cates (2009) overdamped Langevin equations $\partial_{t}\mathbf{r}=v_{0}\hat{\mathbf{u}}+\bar{\bm{\eta}}$ and $\partial_{t}\theta=\bar{\eta}^{R}$ where the self-propulsion velocity is given by $v_{0}=F/\xi$ and the noise terms are given by $\bar{\bm{\eta}}=\bm{\eta}/\xi$ and $\bar{\eta}^{R}=\eta^{R}/\xi_{R}$.

II.2 Statistical Mechanics

Illustration of Coarse-graining Technique: We seek to derive the statistical mechanics of this system of self-propelled hard disks in the large friction limit. One could derive the statistical mechanics from the overdamped Langevin equation mentioned at the end of the previous subsection rather straightforwardly. In this work, we are using hard disks to be able to tractably incorporate the physics of excluded volume interactions into the statistical mechanics. In this case, the large friction limit needs to be taken at the level of the Fokker-Planck equation. In order to illustrate the technique involved, let us begin by considering a collection of noninteracting particles (described by Eqs.(1-2) without the binary collision operator). In this case the statistical mechanics is given by the one particle probability distribution function $f(\mathbf{r},\theta,\mathbf{v},\omega,t)$ of finding a particle with some position $\mathbf{r}$, orientation $\theta$, velocity $\mathbf{v}$, and angular velocity $\omega$ at time $t$. This PDF obeys the Fokker-Planck equation given by Baskaran and Marchetti (2010)

[TABLE]

and the Fokker-Planck operator $\mathcal{D}$ is given by

[TABLE]

Let us define the particle concentration $c(\mathbf{r},\theta,t)=\int d\mathbf{v}d\omega f(\mathbf{r},\theta,\mathbf{v},\omega,t)$, a translational current $\mathbf{J}^{T}=\int d\mathbf{v}d\omega\mathbf{v}f(\mathbf{r},\theta,\mathbf{v},\omega,t)$, and a rotational current $J^{R}=\int d\mathbf{v}d\omega\omega f(\mathbf{r},\theta,\mathbf{v},\omega,t)$. By taking appropriate velocity moments of Eq.(4), one finds that the concentration field obeys a conservation law of the form

[TABLE]

In the large friction limit, the currents are given by

[TABLE]

[TABLE]

with $v_{0}=F/\xi$. This assumption says that in this large friction regime, there exists separation of time scales for spatial relaxations and velocity relaxations, the latter being fast and so at late times, the velocity distribution has relaxed to its local equilibrium form Baskaran and Marchetti (2010). The velocity averages can then be evaluated and the dynamical equation reads

[TABLE]

In the above $D_{r}=\frac{k_{B}T_{R}}{\xi_{r}}$ and the diffusion tensor is of the following form $D_{\alpha\beta}=D_{\parallel}\hat{u}_{\alpha}\hat{u}_{\beta}+D_{\perp}(\delta_{\alpha\beta}-\hat{u}_{\alpha}\hat{u}_{\beta})$. With $D_{\parallel}=v_{0}^{2}\xi^{-1}+k_{B}T\xi^{-1}$ and $D_{\perp}=k_{B}T\xi^{-1}$. Note that in the completely thermal limit we recover ordinary isotropic diffusion $D_{\alpha\beta}=\xi^{-1}k_{B}T\delta_{\alpha\beta}=D_{t}\delta_{\alpha\beta}$. In the purely self-propelled limit we obtain $D_{\alpha\beta}=\xi^{-1}v_{0}^{2}\hat{u}_{1\alpha}\hat{u}_{1\beta}$ as was reported in Baskaran and Marchetti (2010) for self-propelled hard rods. This procedure yields the well studied Smoluchowski equation for ABPs that others have obtained Romanczuk et al. (2012); Solon et al. (2015b); Bialke et al. (2013).

Result: Now, repeating the calculation but retaining the hard core interactions (see Appendix A), the Smoluchowski equation is of the following form

[TABLE]

Note that the rotational part of the Smoluchowksi equation remains unchanged because no torques are exerted in a collision of smooth disks. However, the translational current has the form

[TABLE]

The first three terms of Eq.(11) are equivalent to the terms in Eq.(9) and the additional terms in Eq.(11) are collisional contributions to the translational current. They are given by

[TABLE]

[TABLE]

and where $\hat{\mathbf{u}}_{12}=\hat{\mathbf{u}}_{1}-\hat{\mathbf{u}}_{2}$.

Discussion:

The above collisional contributions in Eq.(10) all have the form of a mean field force, $I_{\alpha}=\int_{\hat{\bm{\sigma}},\theta_{2}}F_{\alpha}(\theta_{1},\theta_{2},\hat{\bm{\sigma}})c^{(2)}(\mathbf{r}_{1},\theta_{1},\mathbf{r}_{1}-\bm{\sigma},\theta_{2},t)$. Where $F_{\alpha}(\theta_{1},\theta_{2},\hat{\bm{\sigma}})$ is some orientation dependent force density and $c^{(2)}(\mathbf{r}_{1},\theta_{1},\mathbf{r}_{1}-\bm{\sigma},\theta_{2},t)$ is the two body probability distribution function of finding particle 1 with position $\mathbf{r}_{1}$ and orientation $\theta_{1}$ and particle 2 with position $\mathbf{r}_{1}-\bm{\sigma}$ and orientation $\theta_{2}$ at a time $t$ (see fig.1).

In the completely thermal limit ($v_{0}=0$) the only term that would contribute to the collision integral is Eq.(12a) and would reproduce the statistical mechanics of thermal hard spheres in the overdamped limit. This is precisely the contribution one would find starting from the revised Enskog theory Balescu (1975); Lachowicz et al. (1991); Chapman and Cowling (1970). 3. 3.

Eq.(12b) is a contribution from collisions arising from self-propulsion alone. Within this equation is a theta function which constrains the range of allowed orientations of the self-propulsion direction of the second particle. That is, the theta function is nonzero only for orientations that would result in a collision. 4. 4.

Eq.(12c) is a collisional contribution that arises from the coupling of thermal noise and self-propulsion and contains the leading order contribution to the translational current from the collisions among particles (see Appendix A).

In this section we have outlined the systematic derivation of the statistical mechanics of self-propelled hard spheres. The results described above should be useful to describe any collection of active particles that interact through strongly repulsive short range interactions, as will be shown in later sections. We now illustrate the applicability of the derived statistical mechanics by investigating some macroscopic properties of a fluid of active particles.

III Steady State Mechanical Properties

As a first illustration, we seek to use the statistical mechanics developed above to calculate a steady state property of this system. Let us consider the pressure of an active fluid. As has been shown in Solon et al. (2015a) this is indeed a state variable for self-propelled particles in the absence of any torques as is the case for smooth hard disks considered here. Mechanically one can then define pressure by considering the system as being confined by a wall at some position $x_{w}\gg 0$ (see fig.2).

represented by some confining potential $V(x)$ and where $x=0$ is taken to be deep in the bulk of the active particle fluid. We assume that deep in the bulk of the fluid the density has a constant bulk value $\rho_{0}$, far beyond $x_{w}$ the density vanishes, and that the system stays uniform in the $\hat{y}$-direction. By Newton’s 3rd Law, the pressure can be computed from the total force acting on the wall,

[TABLE]

Using the procedure developed in Solon et al. (2015a) (see Appendix B for details), we arrive at the following expression for the pressure

[TABLE]

In the above, P_{0}=\big{(}\frac{v_{0}^{2}\xi}{2D_{r}}+\frac{v_{0}^{2}}{2}+k_{B}T\big{)}\rho_{0} is the pressure of noninteracting self-propelled particles and $I$’s are the collision kernels given in Eqs. (12a-12c).

To make any further analytic progress we must evaluate integrals over the collisional contributions which involve the two body distribution $c^{\left(2\right)}$. We now make the ansatz that this two body distribution can be written as a functional of the one body distributions in the following way,

[TABLE]

where g(\mathbf{r}_{1},\mathbf{r}_{2}\big{|}\rho(\mathbf{r},t)) is a functional of the density field. For thermal hard spheres, this function $g$ is the equilibrium pair correlation function and it provides an exact functional relationship between the one and two body distributions. In the case of self-propelled particles, this ansatz cannot be exact as we expect orientational correlations to play some role. Such orientational correlations were in fact characterized in Bialke et al. (2013) through the use of density functional theory and simulations. A systematic estimation of $g$, even in the limited form we have chosen is a hard problem that is beyond the scope of the present work. In the following, we use the form of $g$ associated with thermal hard spheres at contact. That is, we assume orientational correlations can be neglected and that positional correlations are accounted for in the same way as for thermal hard spheres. The test of the validity of this assumption will be the comparison to numerical simulations considered later in this presentation. For the rest of the paper we use the well known estimate of the contact pair-correlation function known as the Carnahan-Starling pair-correlation function in 2 dimensions Balescu (1975).

[TABLE]

where $\phi=\frac{\pi}{4}\rho\sigma^{2}$ is the packing fraction. This estimate is known to give an accurate description of the fluid phase of hard spheres and consequently, this approximation will not capture crystallization effects. We also note that one can choose other estimates of the pair correlation function, such as the Hypernetted Chain Balescu (1975), to approximate density correlations in different parameter regimes. With the ansatz above, the computation of the pressure in Eq.(14) reduces to integrals over the one particle distribution function $c\left(r,\theta\right)$ which is in turn the steady state solution to the Smoluchowski equation Eq.(10). Using the standard procedure Solon et al. (2015b); Baskaran and Marchetti (2008) of representing the distribution as a harmonic expansion in terms of the angular moments, $c=\frac{\rho}{2\pi}+\frac{1}{\pi}\mathbf{P}\cdot\hat{\mathbf{u}}+\tensor{\mathbf{Q}}:(\hat{\mathbf{u}}\hat{\mathbf{u}}-\frac{1}{2}\mathbf{I})+\dots$ where $\mathbf{P}=\int d\theta\hat{\mathbf{u}}c$, is the first moment, $\tensor{\mathbf{Q}}=\int d\theta(\hat{\mathbf{u}}\hat{\mathbf{u}}-\frac{1}{2}\mathbf{I})c$ is the second moment, and assuming a low-moment closure (i.e., truncating the moment expansion at some order, see appendix C for details) we can now evaluate the integrals in Eq.(14) with the result

[TABLE]

In the above, the constants $\lambda_{D}$ and $\lambda_{I}$ control the strength of the collisional contributions to the pressure and depend on the microscopic parameters (see appendix D for the explicit forms). In order to understand the structure of this result, it is useful to consider some limits of Eq.(17). First, in the absence of the self-propulsion (i.e. $v_{0}=0$) we have.

[TABLE]

This is precisely what one would find when calculating the pressure for a hard sphere gas when using Revised-Enskog theory Baskaran et al. (2007); Garzó and Dufty (1999). It consists of the ideal gas term plus an additional correction due to the hard core interactions. In the athermal limit (i.e. $k_{B}T=0$) we obtain

[TABLE]

This limit ($k_{B}T=0$) has been simulated in Solon et al. (2015a) using the Weeks-Chandler-Andersen (WCA) potential. Even though this interaction has a finite collision time we find that our expression captures the computed pressure well for low to moderate densities (see fig.3), thus illustrating the validity of our results for systems with short-range strongly repulsive interactions.

Finally, we compare our result for the pressure with those already in the literature for an overdamped system of ABPs interacting with repulsive potentials. In Solon et al. (2015c), the authors use fluctuating hydrodynamics to derive expressions for the pressure of interacting ABPs in terms of correlation functions between moments of $c(\mathbf{r},\theta,t)$. They find that the pressure can be written as the sum of three terms $P=P_{0}+P_{I}+P_{D}$ where $P_{0}$ is the ideal active pressure and the ”indirect pressure” $P_{I}$ and ”direct pressure” $P_{D}$ are given respectively by $P_{I}=\frac{v_{0}}{D_{r}}\int d^{2}r^{\prime}F_{x}(\mathbf{r}^{\prime})\langle\rho(\mathbf{r}^{\prime})P_{x}(0)\rangle$ and $P_{D}=\int_{0}^{\infty}dx\int d^{2}r^{\prime}F_{x}(\mathbf{r^{\prime}}-\mathbf{r})\langle\rho(\mathbf{r}^{\prime})\rho(\mathbf{r})\rangle$ where $F_{x}$ is the component of the interaction force between particles and the angular brackets here represent an average over the noise. In our theory, we start with a noise averaged dynamical equation. The analogous contributions for the pressure in our case are as follows. The second term in Eq. (14) is the indirect pressure $P_{I}$ and is indeed the evaluation of the corresponding correlation function over the solution of the Smoluchowski equation for the case of hard spheres. The third term in Eq. (14) is the direct pressure $P_{D}$ evaluated for our hard core interactions. We also note that in Solon et al. (2015c) the authors prove that the sum $P_{0}+P_{I}$ is equivalent to the ”swim pressure” investigated by Yang et al. (2014); Takatori et al. (2014) and the same equivalence holds between our theory and the swim pressure as well.

Summarizing, in this section, we have used the statistical mechanics of self propelled hard spheres to tractably evaluate the interaction contributions to the pressure of an active fluid. We have found good agreement with numerical simulations of the pressure of ABPs interacting via a steeply repulsive potential with a finite collision time, thus illustrating the usefulness of the derived statistical mechanics for a variety of model active fluids.

IV Hydrodynamics and Phase Separation

As a second illustration of the utility of the nonequilibrium statistical mechanics constructed here, we derive and characterize a dynamical description of an active fluid on length scales long compared to the particle diameter and time scale long compared to the mean free time. In this regime, the relevant dynamical variables are the conserved quantities and quantities associated with any possible broken symmetries that couple to them. The only conserved quantity is the density of particles given by the zeroth moment of the probability distribution

[TABLE]

The other relevant dynamical quantities are higher angular moments of concentration field. Of these moments, the only relevant quantity that couples to the hydrodynamic field is the polarization described by the first moment of the probability distribution.

[TABLE]

We seek to identify the dynamical equations obeyed by these two quantities.

To derive the continuum equations for the relevant macroscopic fields we must take the corresponding moments of the Smoluchowski equation and they are of the form

[TABLE]

[TABLE]

where

[TABLE]

The fluxes in Eq.(24) are again moments of the two particle distribution as in the case of the pressure calculation above. In order to evaluate these fluxes, we proceed as in the preceding section by assuming the simplest phenomenological closure of the Smoluchowksi equation, that the two particle distribution function can be written as the product of one particle distributions $c^{(2)}(\mathbf{r}_{1},\theta_{1},\mathbf{r}_{1}-\bm{\sigma},\theta_{2},t)=g(\sigma)c^{(1)}(\mathbf{r}_{1},\theta_{1},t)c^{(1)}(\mathbf{r}_{1}-\bm{\sigma},\theta_{2},t)$ and give the one particle distribution function a series representation in terms of its angular moments (Appendix B). As before, we are using the Carnahan-Starling estimate for $g(\sigma)$. Since we are interested in a long wavelength description of the system the nonlocal dependence of the concentration field is expanded in gradients

[TABLE]

Using the above gradient expansion coupled with a low moment closure we arrive the following equations

[TABLE]

[TABLE]

[TABLE]

where the explicit expressions for all the macroscopic parameters and the functions $\mathcal{D}(\rho)$, $K\left(\rho\right)$ and $L_{\alpha\beta}\left(P_{\gamma}\right)$ in terms of the microscopic parameters of the model are given in Appendix D.

These macroscopic equations are complex and nonlinear, with the effect of the repulsive interactions showing up in the coefficients and in the detailed form of the nonlinearities. While careful study of the phase behavior predicted by these equations is warranted, we defer this to future work and make only a few remarks about the structure of these equations. The density equation above has a similar form to those written down for non-interacting ABPs but with a density dependent diffusion coefficient $\mathcal{D}(\rho)$ and a term analogous to the curvature induced flux term found in hydrodynamic theories of orientable active particles Aditi Simha et al. (2003), signifying the fact that the orientation comes with a physical velocity and hence its fluctuations can result in a diffusive flux. The polarization equation does not have a homogeneous nonlinearity reflecting the fact that the interactions among smooth particles are non-aligning, but still has the complex nonlinearities one expects in Toner-Tu type hydrodynamic theories Toner and Tu (1998). Note that $v_{0}\mathbf{P}$ is a measure of the collective self-propulsion velocity of the system and thus encompasses a compressible flow. This is reflected by the presence of the hydrostatic pressure through $\mathcal{P}=\frac{D_{r}}{v_{0}\xi}P$, together with additional Euler order terms $K(\rho)\nabla\cdot\mathbf{P}-\lambda_{2}\mathbf{P}\cdot\nabla\rho$ in the dynamical equation for this flow.

Note that the relaxation of time of Eq.(27) is given by $t=1/D_{R}$. In the rest of this section, we focus on the behavior of this system on times much longer than this characteristic relaxation time. For such times it is reasonable to assume that the polarization has relaxed to its steady state value, i.e., $\partial_{t}\mathbf{P}=0$. Solving for $\mathbf{P}$ in Eq.(27) and substituting into Eq.(26) we find, neglecting higher order gradient terms, the following diffusion equation for the density,

[TABLE]

where the effective density dependent diffusion coefficient is given by

[TABLE]

In the limit $k_{B}T\rightarrow 0$, this effective diffusion coefficient becomes (putting in the explicit forms for $\mathcal{D}(\rho)$ and $\lambda_{I}$ from Appendix D)

[TABLE]

The last term in Eqs.(29-30) is negative and there exists a critical density $\rho_{c}$ or equivalently a critical packing fraction $\phi_{c}=\rho_{c}\frac{\pi}{4}\sigma^{2}$ above which this diffusion coefficient becomes negative. This signals the onset of clustering in the system that has been referred to in the literature as Motility Induced Phase Separation (MIPS). Using the Carnahan-Starling estimate (Eq.(24)) for $g(\sigma)$ one can identify this critical density as a function of system parameters and this is shown by the black line in fig.4.

To better understand this critical density, it is useful to take some limits. First consider the limit of high Peclet number, in this limit the critical density predicted by Eq.(29) is of the form to $\phi_{c}=\frac{1}{\text{Pe}}\frac{\pi(37+21\alpha)}{224\alpha}$, where $\alpha=\xi/D_{r}$. In the limit $k_{B}T\rightarrow 0$, (i.e., no translational diffusion in the overdamped limit) we obtain the following expression for the critical density,

[TABLE]

which again for high enough Peclet number goes as $\phi_{c}\sim\text{Pe}^{-1}$. This $1/\text{Pe}$ behavior has also been seen through phenomenological theories Cates and Tailleur (2015); Marchetti et al. (2016) and in kinetic estimates of the critical density for the onset of MIPS Redner et al. (2013). Finally in fig.4, we compare the estimate provided by this theory against the data from simulations of a system of ABPs interacting through a WCA potential treating $\alpha=\frac{\xi}{D_{r}}$ as a fit parameter and we find good agreement with the data, again illustrating the validity of the result to systems with short range strongly repulsive interactions.

V Summary & Discussion

In this paper we have provided a derivation of the statistical mechanics for self-propelled hard disks starting from first principles. We then considered two applications of the derived statistical mechanics. First is computing the steady state mechanical pressure of a hard sphere active fluid. In the athermal limit ($k_{B}T=0$) we find good agreement with existing numerical simulations Solon et al. (2015a) of ABPs for low to moderate densities. In the absence of activity ($v_{0}=0$) we reproduce the pressure one finds for thermal hard spheres by using Revised Enskog theory Baskaran et al. (2007). We then derived the hydrodynamic equations describing self-propelled hard spheres and identified the critical density for the onset of MIPS in this system. We find that our prediction fits well the numerical data Redner et al. (2013) from ABP simulations and agrees with earlier phenomenological estimates Redner et al. (2013); Marchetti et al. (2016); Cates and Tailleur (2015) in the large Peclet number regime.

This work uses the theory of hard spheres to circumvent the finite collision time problem and hence many-body effects present in arbitrary repulsive potentials. While an idealization, the results obtained in this paper should be useful for any system of self-propelled particles interacting through strongly repulsive potentials and can be useful for incorporating excluded volume effects into diverse models for active systems. This can be seen in the two applications presented, where we have obtained good agreement with numerical simulations of the overdamped dynamics of ABPs interacting the WCA potential. The statistical mechanics presented here transcends to the two applications we have used to illustrate its utility and is potentially useful in diverse active materials modeling.

Acknowledgments

We thank Yaouen Fily and Gabe Redner for providing the numerical data. We acknowledge support from the NSF DMR-1149266, the Brandeis MRSEC DMR-1420382, and IGERT DGE-1068620.

Appendix A Derivation of the Smoluchowski equation

We provide a complete derivation of the overdamped dynamics of self-propelled hard disks. Our starting point will be the coupled Langevin equations

[TABLE]

[TABLE]

where $T(i,j)$ is the binary collision operator given in Eq.(3). The noise function $\bm{\eta}_{i}$ and $\eta_{i}^{R}$ are both Gaussian white noise variables with zero mean and the same correlations defined in the main body of text . The phase space variables associated with the above Langevin dynamics are the position $\mathbf{r}_{n}(t)$, the velocity $\mathbf{v}_{n}(t)$, the angular velocity $\bm{\omega}_{n}(t)$ and the orientation $\theta_{n}(t)$ of the n particles. We begin by considering an arbitrary phase space function $A(\Gamma)$ evaluated at some phase point $\Gamma_{t}=\left\{\mathbf{r}_{1}(t),\dots,\mathbf{r}_{n}(t),\mathbf{v}_{1}(t),\dots,\mathbf{v}_{n}(t),\bm{\omega}_{1}(t),\dots,\bm{\omega}_{n}(t),\theta_{1}(t),\dots,\theta_{n}(t)\right\}$ which has evolved from an initial configuration $\Gamma=\left\{\mathbf{r}_{1},\dots,\mathbf{r}_{n},\mathbf{v}_{1},\dots,\mathbf{v}_{n},\bm{\omega}_{1},\dots,\bm{\omega}_{n},\theta_{1},\dots,\theta_{n}\right\}$. A phase space function at some time t can be expressed in terms of a generator for the dynamics in the following way

[TABLE]

In the above equation, $L$ is the sum of generators for translations and rotations for each particle coordinate. For a system with pairwise additive, conservative, non-singular interactions the operator L can be written as

[TABLE]

where the single particle component for self propelled particles with a propulsion force $F^{sp}$ along its body axis is given by

[TABLE]

The inclusion of the last two terms in the above generator results from using Ito calculus to correctly incorporate the stochastic component of the dynamics. Eq.(35) is the generalization of the completely deterministic case, in which the dynamics are governed by Hamilton’s equations. One can verify that, for Hamiltonian systems, the generator L changes the positions and velocities according to Hamilton’s equations. In this case, L is the linear operator representing the Poisson bracket of its operand with the Hamiltonian for the system Dufty and Baskaran (2005). However, we are considering the stochastic dynamics of elastic hard spheres, thus the interactions must be taking into account via the binary collision operator $T(i,j)$. In this case the operator L becomes

[TABLE]

The above formalism completely determines the stochastic dynamics of any observable $A(\Gamma)$ for given initial conditions in phase space. In this study, we would like to write down not the stochastic observable, but its ensemble averaged value for a given ensemble of initial conditions $\hat{\rho}(\Gamma)$. This is represented by the following phase space average

[TABLE]

One can equivalently treat the phase space density as the dynamical variable Zwanzig (2001)

[TABLE]

and taking the time derivative of both equation yields

[TABLE]

[TABLE]

[TABLE]

Where $\mathcal{L}$ is the adjoint operator to L. The result of this simple manipulation is a Liouville-like equation for the phase space probability density

[TABLE]

with the adjoint operator is given by

[TABLE]

In Eq.(40), the single particle component has been identified by an integration by parts, while the collision operator $\bar{T}(i,j)$ has been explicitly constructed using the restituting collisions Baskaran and Marchetti (2010)

[TABLE]

In the above, $b_{ij}^{-1}$ is the generator of restituting collisions $(b_{ij}^{-1}A(x_{i}^{\prime},x_{j}^{\prime})=A(x_{i},x_{j}))$ which replaces post collisional velocities with its pre-collisional values. Finally, we average over the noise $\rho=\langle\hat{\rho}\rangle$. The resulting equation is given by

[TABLE]

with the operator $\mathcal{L}$ given by

[TABLE]

The above equation governs the time evolution of the phase space density of N self-propelled hard spheres. To proceed, we now introduce reduced distribution functions

[TABLE]

and consider the first equation of the resulting hierarchy

[TABLE]

Where the one particle operator $\mathcal{D}$ in Eq.(45) is given by

[TABLE]

The goal of this section is to obtain a dynamical equation for the local concentration field by

[TABLE]

For convenience, we introduce a translational and rotational current defined as the velocity moments of the 1-particle distribution function

[TABLE]

[TABLE]

Taking velocity moments of Eq.(45), we arrive at the following equations

[TABLE]

[TABLE]

[TABLE]

In the above,

[TABLE]

[TABLE]

[TABLE]

are the second velocity moments of the probability distribution function. Also present are the following terms

[TABLE]

[TABLE]

Eq.(56) arises from the linear momentum transfer due to collisions. The integral in Eq.(57) equates to zero because these are smooth hard discs and therefore, no angular momentum is transferred in a collision. The translational and rotational currents in the above equations are subject to frictional damping and as such relax on time scales of order $m/\xi$. On time scales $t\gg m/\xi$ the flux can be approximated as

[TABLE]

[TABLE]

Velocity Integration

To complete the derivation of the Smoluchowski equation, we must evaluate the velocity integrals in Eqs.(58-59). As outlined in section II of the main text, we assume that on these timescales the velocity distributions have relaxed to their local equilibrium form. Explicitly we have that

[TABLE]

where N is a normalization factor. With this distribution, it is readily shown that

[TABLE]

[TABLE]

[TABLE]

The above averages are exact for the one-body terms, but for the collision integral (which depends on the orientation of the colliding particles and an integration over the two body distribution), the exact evaluation of the velocity integrals is not possible. To continue further, we make the following asymptotic approximation that accurately captures the physics in the limits that $v_{0}^{2}\gg k_{B}T/m$ (or $k_{B}T/m\gg v_{0}^{2}$)

[TABLE]

The delta function enforces that the particle will have a self-propulsion velocity $v_{0}$ along its body axis $\hat{\mathbf{u}}$ and reproduces the the correct distribution in the athermal limit. The Maxwellian accounts for thermal noise and reproduces the statistical mechanics of thermal hard spheres in the limit of zero self-propulsion.

The last step is to perform the velocity averaging in Eq.(56). Explicitly we must evaluate the following,

[TABLE]

This evaluation of this integral leads to Eqs.(12a-12c) for $I^{thermal}_{\alpha}$, $I^{active}_{\alpha}$, and $I^{cross}_{\alpha}$. To evaluate $I^{cross}_{\alpha}$, which comes from the cross terms in the multiplication of the velocity distributions, one must make an additional approximation. Contained within the integral for the cross terms, one finds terms like

[TABLE]

Since the integral is over all possible velocities, to good approximation we have

[TABLE]

Using this approximation, the integrals can be computed. The integrals leading to $I^{thermal}_{\alpha}$ and $I^{active}_{\alpha}$ can be computed exactly. These are the most important contributions since, in the limits of purely self-propelled or purely thermal, these are the only contributing factors. This completes the derivation of the Smoluchowski equation. Combining Eq.(50), Eqs.(58-59) and the velocity averages yields Eqs.(10-11) in the main text.

Appendix B Derivation of Pressure

To evaluate this expression we start with the Smoluchowski equation (Eq.(10)). In steady state, the Smoluchowksi equation takes the following form

[TABLE]

Because of the translational invariance in $\hat{y}$, in steady state, the only spatial dependence in the above is through the x-coordinate. Integrating Eq.(68) over the orientations $\theta$, on can see clearly that the resulting equation is of the form $\partial_{x}J^{T}=0$ with $J^{T}$ being the particle current. Since this system has impermeable boundary conditions ($J^{T}=0$ at the wall), the only admissible solution is that $J^{T}=0$ everywhere. Written explicitly we have,

[TABLE]

Rewriting the left hand side of Eq.(69) in terms of the density by noting that $\rho=\int cd\theta$ and by integrating over x, we obtain the pressure (Eq.(13)) written in terms of the concentration field

[TABLE]

Now let us multiply Eq.(68) by $\hat{u}_{x}$ and integrate over $\theta$

[TABLE]

and now integrate Eq.(71) over x. Note that the right hand side of Eq.(71) is given as a total derivative and therefore trivially integrated. Aside from the interaction terms in Eq.(71), the only surviving term is the angular integral over $\hat{u}_{x}\hat{u}_{x}c$ which is proportional $\rho/2$. The remaining terms vanish because there can be no orientational order in the bulk of the fluid $(x=0)$ or at $x=\infty$. With these facts, we have that

[TABLE]

Using Eq.(72) to eliminate the first term on the right hand side of Eq.(70) we recover Eq.(14)

[TABLE]

Appendix C Evaluation of Mean-Field Force and Low Moments Closure

In this section we give the explicit evaluation of the mean-field force (Eqs.(12a-12c)) and an outline of the low moment closure procedure used in the text. To construct the low moment closure we first represent the concentration field as the following harmonic expansion Ahmadi et al. (2006)

[TABLE]

where the irreducible tensors $T^{m}_{i_{1}\dots i_{m}}(\theta)$ are equivalent to the spherical harmonics but expressed here in Cartesian coordinates. For illustrative purposes, the first four irreducible tensors are given by

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The $m$ th order moment $a^{m}_{i_{1}\dots i_{m}}(\mathbf{r},t)$ of the concentration is given by

[TABLE]

The dynamical equation for the $m$ th moment can then be obtained by taking moments of the Smoluchowksi equation (Eq.(10)), with the result

[TABLE]

We see that when written in this form, the density and polarization are just given as the zeroth and first moment of the concentration field. The low moment closure is the approximation that $c(\mathbf{r}_{1},\theta_{1},t)$ can be expressed as a functional of the first two moments $c^{(1)}(\theta_{1};\{\rho(\mathbf{r},t),P_{\alpha}(\mathbf{r},t)\})$. The motivation for this closure is the following: One can see from Eq.(80) that all moments greater than the zeroth moment have a finite relaxation rate given by $m^{2}Dr$. Thus, these moments will relax to values that depend on gradients of the local concentration. When the values of the higher moments are substituted into the density equation they result in terms that are irrelevant compared to the terms resulting from the polarization equation, i.e., contain more powers of gradients and fields. Therefore, we close the expansion by assuming that the second and higher moments can be neglected. This implies that $\int d\theta\hat{u}_{i}\hat{u}_{j}c=\frac{1}{2}\delta_{ij}\rho$ and $\int d\theta\hat{u}_{i}\hat{u}_{j}\hat{u}_{k}c=\frac{1}{4}(\delta_{ij}P_{k}+\delta_{ik}P_{j}+\delta_{jk}P_{i})$. We note that this is a valid closure in the absence of aligning interactions because the density and polarization represent the relevant macroscopic fields. To accurately capture the effects of aligning interactions one must include the second moment as is done in Farrell et al. (2012); Baskaran and Marchetti (2008); Hancock and Baskaran (2015). With this we can now evaluate the mean field force. The mean field force consists of three parts $I_{\alpha}^{thermal}$, $I_{\alpha}^{sp}$, and $I_{\alpha}^{cross}$. To evaluate these quantities we first make the functional ansatz used in the main body of the text, namely that $c^{(2)}(\mathbf{r}_{1},\theta_{1},\mathbf{r}_{1}-\bm{\sigma},\theta_{2},t)=g(\sigma)c(\mathbf{r}_{1},\theta_{1},t)c(\mathbf{r}_{1}-\bm{\sigma},\theta_{2},t)$. We then gradient expand the nonlocal dependence of the concentration field $c(\mathbf{r}-\bm{\sigma},\theta_{2},t)=c(\mathbf{r},\theta_{2},t)-\sigma_{\alpha}\partial_{\alpha}c(\mathbf{r},\theta_{2},t)+\dots$ , which is valid in a long wavelength (hydrodynamic) description of the system. To first order in gradients the integrals can be evaluated. The result of this integration is

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Appendix D Constants and functions used in the hydrodynamics and pressure

For convenience, let us call $A=(2\pi k_{B}T+v_{0}^{2})^{-1}$ and $a=(1+\frac{v_{0}^{2}}{2k_{B}T})$. The constants in the hydrodynamic equations and pressure are then given by

[TABLE]

[TABLE]

[TABLE]

[TABLE]

and the functions are given by

[TABLE]

[TABLE]

[TABLE]