Mode-coupling theory for the steady-state dynamics of active Brownian particles

TL;DR

This paper develops a mode-coupling theory to describe the steady-state dynamics of active Brownian particles, linking microscopic activity to macroscopic glassy behavior through a novel correlation function.

Contribution

It introduces a theoretical framework combining projection operators and mode-coupling approximations for active particles, incorporating nonequilibrium steady-state correlations.

Findings

Steady-state correlation functions quantify activity effects.

Activity influences short-time dynamics via spatial correlations.

Long-time glassy dynamics are affected by microscopic currents.

Abstract

We present a theory for the steady-state dynamics of a two-dimensional system of spherically symmetric active Brownian particles. The derivation of the theory consists of two steps. First, we integrate out the self-propulsions and obtain a many-particle evolution equation for the probability distribution of the particles' positions. Second, we use projection operator technique and a mode-coupling-like factorization approximation to derive an equation of motion for the density correlation function. The nonequilibrium character of the active system manifests itself through the presence of a steady-state correlation function that quantifies spatial correlations of microscopic steady-state currents of the particles. This function determines the dependence of the short-time dynamics on the activity. It also enters into the expression for the memory matrix and thus influences the long-time glassy dynamics.

Equations

$${\dot{\mathbf{r}}}_i$$

$${\dot{\theta }}_i$$

$${\mathbf{F}}_i=-\sum _{j\ne i}{\nabla }_{i}V(r_{ij}),$$

$$⟨{\left({\mathbf{r}}_i(t)-{\mathbf{r}}_i(0)\right)}^2⟩=4D_{t}t+2\frac{v_0^2}{D_r^2}\left(D_{r}t-1+e^{-D_{r}t}\right).$$

$$⟨{\left({\mathbf{r}}_i(t)-{\mathbf{r}}_i(0)\right)}^2⟩\approx 4D_{t}t+\frac{v_0^2}{D_r^2}{t}^{2}t\ll D_r^{-1}$$

$$D_{\text{eff}}=D_t+\frac{v_0^2}{2D_r}$$

$$T_{\text{eff}}=D_{\text{eff}}{\xi }_0=T+\frac{v_0^2}{2D_r{\xi }_0}.$$

$${\partial }_t{P}_N({\mathbf{r}}_1,{\theta }_1,...,{\mathbf{r}}_N,{\theta }_N;t)=\Omega P_N({\mathbf{r}}_1,{\theta }_1,...,{\mathbf{r}}_N,{\theta }_N;t),$$

$$\Omega =D_t\sum _{i=1}^N{\nabla }_i\cdot \left[{\nabla }_i-\beta {\mathbf{F}}_i\right]-\sum _{i=1}^N{v}_0{\nabla }_i\cdot \mathbf{e}({\theta }_i)+D_r\sum _i{\partial }_{{\theta }_i}^2.$$

$$\Omega P_N^{ss}({\mathbf{r}}_1,{\theta }_1,...,{\mathbf{r}}_N,{\theta }_N)=0.$$

$$P_N^{\text{ss}}({\mathbf{r}}_1,{\theta }_1,...,{\mathbf{r}}_N,{\theta }_N)\ne P_N^{\text{ss}}({\mathbf{r}}_1,...,{\mathbf{r}}_N)P_N^{\text{ss}}({\theta }_1,...,{\theta }_N),$$

$$P_N^{\text{ss}}({\mathbf{r}}_1,...,{\mathbf{r}}_N)$$

$$P_N^{\text{ss}}({\theta }_1,...,{\theta }_N)$$

$$P_N^{\text{ss}}({\theta }_1,...,{\theta }_N)={\left(2\pi \right)}^{-N}.$$

$$F(q;t)=\frac{1}{N}⟨\sum _i{e}^{-i\mathbf{q}\cdot {\mathbf{r}}_i(t)}\sum _j{e}^{i\mathbf{q}\cdot {\mathbf{r}}_j(0)}⟩.$$

$$F(q;t)=\frac{1}{N}⟨n(\mathbf{q})\exp \left(\Omega t\right)n(-\mathbf{q})⟩.$$

$$n(\mathbf{q})=\sum _l{e}^{-i\mathbf{q}\cdot {\mathbf{r}}_l},$$

$${\partial }_t{P}_N({\mathbf{r}}_1,{\theta }_1,...,{\mathbf{r}}_N,{\theta }_N;t)=$$

$${\mathbf{j}}_i({\mathbf{r}}_1,{\theta }_1,...,{\mathbf{r}}_N,{\theta }_N;t)=$$

$$j_i^{\theta }({\mathbf{r}}_1,{\theta }_1,...,{\mathbf{r}}_N,{\theta }_N;t)=$$

$${\mathbf{j}}_i^{\text{ss}}({\mathbf{r}}_1,...,{\mathbf{r}}_N;t)=\int d{\theta }_1...d{\theta }_N{\mathbf{j}}_i^{\text{ss}}({\mathbf{r}}_1,{\theta }_1,...,{\mathbf{r}}_N,{\theta }_N;t)$$

$$\equiv \int d{\theta }_1...d{\theta }_N\left(-D_t{\nabla }_i+{\xi }_0^{-1}{\mathbf{F}}_i+v_0{\mathbf{e}}_i\right)$$

$$\times P_N^{\text{ss}}({\mathbf{r}}_1,{\theta }_1,...,{\mathbf{r}}_N,{\theta }_N)=0.$$

$$v_0{⟨{\mathbf{e}}_i⟩}_{\text{lss}}=-{\xi }_0^{-1}{\mathbf{F}}_i+D_t{\nabla }_i\ln P_N^{\text{ss}}({\mathbf{r}}_1,...,{\mathbf{r}}_N)$$

$${⟨\dots ⟩}_{\text{lss}}=$$

$${\mathcal{P}}_{\text{lss}}{P}_N({\mathbf{r}}_1,{\theta }_1,...,{\mathbf{r}}_N,{\theta }_N;t)$$

$${\mathcal{Q}}_{\text{lss}}=\mathcal{I}-{\mathcal{P}}_{\text{lss}},$$

$${\partial }_t{\mathcal{P}}_{\text{lss}}{P}_N({\mathbf{r}}_1,{\theta }_1,...,{\mathbf{r}}_N,{\theta }_N;t)=$$

$${\partial }_t{\mathcal{Q}}_{\text{lss}}{P}_N({\mathbf{r}}_1,{\theta }_1,...,{\mathbf{r}}_N,{\theta }_N;t)=$$

$${\mathcal{Q}}_{\text{lss}}{P}_N({\mathbf{r}}_1,{\theta }_1,...,{\mathbf{r}}_N,{\theta }_N;t=0)=0.$$

$$\mathcal{L}\mathcal{T}\left[{\partial }_t{\mathcal{P}}_{\text{lss}}{P}_N({\mathbf{r}}_1,{\theta }_1,...,{\mathbf{r}}_N,{\theta }_N;t)\right](z)=\left[{\mathcal{P}}_{\text{lss}}\Omega {\mathcal{P}}_{\text{lss}}+{\mathcal{P}}_{\text{lss}}\Omega {\mathcal{Q}}_{\text{lss}}\frac{1}{z-{\mathcal{Q}}_{\text{lss}}\Omega {\mathcal{Q}}_{\text{lss}}}{\mathcal{Q}}_{\text{lss}}\Omega {\mathcal{P}}_{\text{lss}}\right]{\mathcal{P}}_{\text{lss}}{P}_N({\mathbf{r}}_1,{\theta }_1,...,{\mathbf{r}}_N,{\theta }_N;z).$$

$${\mathcal{P}}_{\text{lss}}\Omega {\mathcal{P}}_{\text{lss}}{P}_N({\mathbf{r}}_1,{\theta }_1,...,{\mathbf{r}}_N,{\theta }_N;z)$$

Full text

I Introduction

Recently, there has been a lot of interest in the structure and dynamics of strongly interacting active matter systems Ramaswamyrev ; Catesrev ; Marchettirev ; Bechingerrev . This is motivated by a combination of experimental and simulational studies that uncovered fascinating phenomena with no analogs in equilibrium thermal (passive) systems. Very recent examples include transition from turbulent to coherent flows in confined three-dimensional active fluids Dogic2017 , spontaneously flowing two-dimensional crystals Dauchot2018 , and structure and dynamics of active systems on curved geometries Marchetti2017 ; Janssen2017 . The motivation comes also from the fact that activity can significantly and sometimes surprisingly modify phenomena occurring in thermal systems. For example, early works Henkes2011 ; BerthierKurchan ; Ni2013 ; Berthier2014 showed that introducing active forces changes the glassy behavior of an equilibrium fluid instead of wiping out the glassy phase completely, which could have been expected on the basis of a perhaps naive analogy with what happens when a shear flow is imposed on an equilibrium colloidal glass. The resulting non-equilibrium glassy behavior and the eventual glass transition exhibit features observed in equilibrium supercooled liquids Berthier2011 , including slowing down, transient localization of particles’ positions, non-exponential relaxations and dynamic heterogeneity. However, also present are unique non-equilibrium features such as non-trivial equal-time velocity correlations Szamel2015 and a variety of different effective temperatures Levis2015 ; Preisler2016 . It is this confluence of well-known but still not fully understood Berthier2011 glassy behavior and strongly non-equilibrium characteristics of active matter systems that makes non-equilibrium glassy behavior so fascinating.

Most of the experimental and many of the simulated active glassy systems are quite complex and are characterized by many independent parameters, which makes systematic studies quite involved. For this reason, in our initial investigations of non-equilibrium glassy behavior we focused on a a very simple active system, the so-called active Ornstein-Uhlenbeck particles (AOUPs) model Fodor2016 introduced independently in Refs. Szamel2014 and Maggi2015 . In this model, the dynamics is overdamped, there is no thermal noise, and the particles move under the combined influence of the inter-particle forces and the self-propulsion. The self-propulsion evolves according to the Ornstein-Uhlenbeck process, independently of the configuration of the particles. For a given interaction potential, the AOUP system is characterized by three parameters, the number density $\rho$, the single-particle effective temperature $T_{\text{eff}}$ characterizing the driving energy, and the persistence time of the self-propulsion $\tau_{p}$. In the limit of the vanishing persistence time, $\tau_{p}\to 0$, the AOUP system becomes equivalent to a thermal (passive) system at the temperature equal to the single-particle effective temperature. Thus, for the AOUP system the departure from equilibrium is characterized by a single parameter, the persistence time of the self-propulsion.

To analyze the dynamics of dense systems of AOUPs we used a combination of simulations and analytical theory. In collaboration with E. Flenner and L. Berthier, we uncovered the significance of the equal-time, steady-state correlation function of particles’ velocities Szamel2015 . This function exhibits strong wavevector dependence for long persistence times and becomes trivial, i.e. wavevector independent and related to the temperature, in the equilibrium limit $\tau_{p}\to 0$. Furthermore, we showed that increasing departure from equilibrium can result in both faster dynamics and fluidization of a glassy system Berthier2017 , which was also observed by other workers Ni2013 ; Berthier2014 ; Mandal2014 , and, unexpectedly, slower dynamics and glassification of a fluid system Flenner2016 ; Berthier2017 . We found that this highly non-trivial behavior of the effective glass transition line in the $T_{\text{eff}}-\rho$ plane correlates with the dependence of the steady-state structure on the departure from equilibrium Berthier2017 .

On the theory side, we derived an approximate theory for the steady-state dynamics of AOUPs Szamel2016 . The basic assumption of this theory was the absence of the steady state currents, or more precisely, that the velocities of individual particles vanish after averaging over the self-propulsions. The derivation was done in two steps. First, we approximately integrated out the self-propulsions and obtained an effective equation of motion for the many-particle distribution of particles’ positions. This equation featured a time-dependent diffusivity matrix. The time-dependence was crucial for retaining proper (ballistic) short-time dynamics of the AOUPs. Next, using the effective many-particle equation of motion we derived an approximate equation of motion for the intermediate scattering function. In this last step we used a factorization approximation analogous to that employed in the mode-coupling theory of glassy dynamics and the glass transition Goetzebook . The resulting self-consistent equation of motion for the intermediate scattering function resembled an equation of motion for an under-damped colloidal system (without hydrodynamic interactions). This could have been expected since the AOUP dynamics is ballistic at short times and diffusive at long times, as in an under-damped colloidal system. However, we found that the equation of motion for the intermediate scattering function depends, in a highly non-trivial way, on the activity of the system. This dependence manifests itself through the correlation function of particles’ velocities, which entered into an analogue of the frequency matrix term (and thus determined the short-time dynamics) and into the vertices of the approximate expression for the memory function.

In the present contribution we extend the derivation presented in Ref. Szamel2016 to the most often studied active system, the so-called active Brownian particles (ABPs) model tenHagen ; FilyMarchetti . In this model, the dynamics is overdamped but there is also thermal noise. Thus, the particles move under the combined influence of the inter-particle forces, thermal noise originating from fluctuations of the solvent, and the self-propulsion. The magnitude of the self-propulsion is fixed and its direction changes via rotational diffusion. For a given interaction potential, the ABP system is characterized by four parameters, the number density $\rho$, the translational diffusion coefficient $D_{t}$, which depends in the temperature $T$ characterizing the thermal noise, the magnitude of the self-propulsion $v_{0}$ and the rotational diffusion constant specifying the evolution of the direction of the self-propulsion $D_{r}$. For real active colloidal particles both $D_{t}$ and $D_{r}$ are determined by the hydrodynamics of the solvent and are, therefore, related but we will treat them as independent model parameters, as was done in some simulational Ni2013 and theoretical Voigtmann2017 studies.

Compared to the AOUP system, the ABP model introduces two complications. First, we need to include thermal noise (translational diffusion) in addition to the self-propulsion. Second, in the ABP model the relation between the variable describing the state of the self-propulsion and the particle motion is non-linear and thus even the solution of the single-particle motion is highly non-trivial Franosch2016 .

To put our work in the context of earlier microscopic theoretical investigations we would like to mention two other theories for the dynamics of the ABP system Voigtmann2017 ; Farage and a theory for a generalized AOUP model that includes thermal noise Feng2017 . All these theoretical approaches share one important feature with our theory for the AOUP system Szamel2016 and with the theory presented here: all these theories rely upon a factorization approximation and thus are of mode-coupling flavor. However, they differ in their focus and in the way they include the self-propulsion. In the next two paragraphs we will delineate these differences. We shall also mention a theory of Nandi and Gov NandiGov2017 , which falls in between microscopic theories of Refs. Voigtmann2017 ; Farage ; Feng2017 and the theoretical analysis of an effective $p$-spin-like active model of Ref. BerthierKurchan .

The difference between the theories of Farage and Brader Farage and of Liluashvili et al. Voigtmann2017 and both our theories and the approach of Feng and Hou Feng2017 is that the former theories follow the philosophy of the “integration-through-transients” approach to the dynamics of colloidal systems under shear developed by Fuchs and Cates Fuchsrev . In this approach one assumes that the system was in an equilibrium state in the infinitely distant past and then the drive, in this case the activity, was turned on. The advantage of this approach is that, in principle, it allows one to calculate both equal time and dynamic properties of an active system within a single theory. Within this approach one easily derives approximate equations of motion for transient correlation functions. However, it is much more difficult to obtain equations of motion for time-dependent correlation functions in a steady state Krueger2009 . In addition, it is not clear how one could use this approach to describe a fully athermal system that does not have any dynamics without the activity. In contrast, our theories and the theory of Feng and Hou focus on the steady-state dynamics of the active system. Their disadvantage is that a separate investigation, either analytical or simulational, is needed to supply these theories with necessary equal-time steady-state correlation functions.

On the other hand, the difference between theories of Farage and Brader Farage and Feng and Hou Feng2017 and both our theories and the theory of Liluashvili et al. Voigtmann2017 is in the way the self-propulsion is treated. In the former theories the active system is approximately replaced by a thermal (passive) system using either a very simple time coarse-graining procedure Farage or a procedure introduced by Fox Fox in the context of approximating a system with a colored noise by an equivalent system with white noise. In contrast, both our theories and the theory of Liluashvili et al. include, in an approximate way, the time evolution of the self-propulsion and the resulting time-delayed response. In the theory of Ref. Szamel2016 and in the theory presented here this is achieved by retaining the time dependence in the many-particle diffusivity matrix. The theory of Liluashvili et al. treats the self-propulsion explicitly and thus constitutes the most accurate description of the self-propulsion.

The paper is organized as follows. In the next section, Sec. II, we introduce and briefly discuss the ABP system. In Sec. III we discuss the basic assumption of our theory and its physical meaning. In Sec. IV we present the derivation of the approximate theory for the long-time dynamics of the intermediate scattering function of the ABP model. We end with discussion in Sec. V.

II Active Brownian particles

We consider a two-dimensional comment2d system of $N$ interacting, self-propelled particles in a volume (area) $A$. The average density $\rho=N/A$. The particles interact via a spherically symmetric potential $V(r)$. They move in a viscous medium that is characterized by the friction coefficient of a single particle, which we denote by $\xi_{0}$. We assume that the friction felt by a particle is independent of the particle density and configuration, and thus we neglect hydrodynamic interactions Dhont . Each particle moves under the combined influence of interparticle force derived from the potential $V(r)$, random, thermal forces due to the fluctuations of the solvent and self-propulsion commentSP . We assume that the velocity due to the self-propulsion has a constant magnitude $v_{0}$ and that it evolves in time via rotational diffusion. The corresponding equations of motion read,

[TABLE]

In Eq. (1), $\mathbf{r}_{i}$ is the position of particle $i$ and $\mathbf{F}_{i}$ is the force acting on particle $i$ originating from the interactions,

[TABLE]

where $\boldsymbol{\nabla}_{i}$ denotes a partial derivative with respect to $\mathbf{r}_{i}$, $\boldsymbol{\nabla}_{i}=\frac{\partial}{\partial\mathbf{r}_{i}}$. Furthermore, in Eq. (1), $\boldsymbol{\eta}_{i}$ is the Gaussian thermal noise with zero mean and variance $\left<\boldsymbol{\eta}_{i}(t)\boldsymbol{\eta}_{j}(t^{\prime})\right>=2\xi_{0}T\mathbf{I}\delta_{ij}\delta(t-t^{\prime})$ with $T$ being the temperature (we use units such that the Boltzmann constant $k_{B}=1$) and $\mathbf{I}$ being the unit tensor. Finally, in Eq. (1), $\mathbf{e}(\theta_{i})$ is the unit vector specifying the direction of the self-propulsion, $\mathbf{e}(\theta_{i})=(\cos\theta_{i},\sin\theta_{i})$. The single particle translational diffusion coefficient is proportional to the temperature, $D_{t}=T/\xi_{0}$. In Eq. (2), $\eta_{\theta i}$ is an internal Gaussian noise with zero mean and variance $\left<\eta_{\theta i}(t)\eta_{\theta j}(t^{\prime})\right>=2D_{r}\delta_{ij}\delta(t-t^{\prime})$, with $D_{r}$ being the rotational diffusion coefficient. As mentioned earlier, we treat $D_{t}$ and $D_{r}$ as independent parameters.

We emphasize that the direction of each self-propulsion evolves independently of positions of the particles and of (direct) forces acting on the particles.

Without interactions, particles evolving according to Eqs. (1-2) perform a somewhat complicated random walk, with non-trivial higher-order cumulants, but with a relatively simple mean-square displacement tenHagen

[TABLE]

One should note that the expression for the mean-square displacement is the same as for a simpler system of AOUPs with additional thermal noise. According to Eq. (4), the self-propulsion contributes to the mean-square displacement at the level of $t^{2}$,

[TABLE]

and the long-time motion, $t\gg D_{r}^{-1}$, is diffusive with diffusion coefficient $D_{\text{eff}}$,

[TABLE]

Comparing expression (6) with the well-known formula for the diffusion coefficient of a Brownian particle moving in a viscous medium with friction constant $\xi_{0}$, $D_{\text{Brownian}}=T/\xi_{0}$, we can define the single-particle effective temperature,

[TABLE]

We note that in the limit of fast evolution of the self-propulsions, $D_{r}\to\infty$, $v_{0}\to\infty$, $v_{0}^{2}/D_{r}=\text{ const.}$, the mean-square displacement as given by Eq. (4) becomes linear in time. In fact, in this limit the system becomes equivalent to a thermal system with the temperature given by Eq. (7).

For further theoretical analysis it is convenient to replace the description of the system’s dynamics in terms of the equations of motion (1-2) by the equivalent description in terms of an $N$-particle joint probability distribution of positions and self-propulsions $P_{N}(\mathbf{r}_{1},\mathbf{e}_{1},...,\mathbf{r}_{N},\mathbf{e}_{N};t)\equiv P_{N}(\mathbf{r}_{1},\theta_{1},...,\mathbf{r}_{N},\theta_{N};t)$. This distribution evolves in time with evolution operator $\Omega$,

[TABLE]

which can be derived from equations of motion (1-2),

[TABLE]

We assume that there exists a steady state. In other words, we assume there exists a probability distribution $P_{N}^{ss}(\mathbf{r}_{1},\theta_{1},...,\mathbf{r}_{N},\theta_{N})$ such that

[TABLE]

We expect that there are correlations between positions and self-propulsion velocities BialkeEPL and thus

[TABLE]

where $P_{N}^{\text{ss}}(\mathbf{r}_{1},...,\mathbf{r}_{N})$ and $P_{N}^{\text{ss}}(\theta_{1},...,\theta_{N})$ are the steady-state distributions of positions and self-propulsions,

[TABLE]

In general, neither the joint steady-state distribution $P_{N}^{\text{ss}}(\mathbf{r}_{1},\theta_{1},...,\mathbf{r}_{N},\theta_{N})$ nor the steady-state distributions of positions $P_{N}^{\text{ss}}(\mathbf{r}_{1},...,\mathbf{r}_{N})$ are known exactly (for approximate theories for the latter distribution see Refs. Maggi2015 ; FarageKB ). However, the steady-state distribution of self-propulsions has a trivial form,

[TABLE]

The main object of our theory is the intermediate scattering function, $F(q;t)$, which describes the time dependence of the collective density fluctuations,

[TABLE]

Here and in the following the brackets $\left<\dots\right>$ denote averaging over a steady-state distribution of positions and self-propulsions.

Evolution operator (9) allows us to rewrite the definition of the intermediate scattering function (15),

[TABLE]

In Eq. (16) $n(\mathbf{q})$ is the Fourier transform of the microscopic density,

[TABLE]

We emphasize that in Eq. (16) and in all similar formulas the steady-state distribution stands to the right of the quantity being averaged, and all operators act on it too.

III The main assumption: absence of average currents

We follow our earlier theory for AOUPs Szamel2016 and assume that in the steady state, the currents vanish after integrating out the self-propulsions. This assumption will allow us to approximate our system by a passive system with a time-dependent diffusivity matrix.

To make our assumption explicit we first rewrite the equation of motion for the joint probability distribution of positions and self-propulsions, Eq. (8), in the form of a continuity equation,

[TABLE]

where current densities are defined as

[TABLE]

[TABLE]

Current densities (III-III) are microscopic quantities (i.e in principle they depend on positions and self-propulsions of all the particles), which may be non-zero in a system without detailed balance. We assume that in the steady state, the current density in the position space integrated over self-propulsions vanishes,

[TABLE]

Assumption (III) implies the following expression for the local steady-state average of the self-propulsion

[TABLE]

where the local steady-state average is defined as

[TABLE]

Eq. (22) expresses a balance of the self-propulsion acting on particle $i$ and the sum of the total potential force acting on this particle and the averaged force due to the solvent fluctuations, for a given configuration, i.e. for a given set of the positions of the particles.

We are not aware of any study that specifically focused on the existence of non-trivial steady-state currents in high density active systems without aligning interactions commentcurrents . We note that the assumption (III) is made at the level of $N$-particle quantities. Thus, its direct simulational verification seems rather difficult. However, it might be possible to define and measure reduced (few-particle) current densities. Work in this direction is planned and the results will be reported in the future.

IV Mode-coupling theory

IV.1 Effective evolution operator for particles’ positions

We again follow our earlier theory Szamel2016 and start the analysis of the ABP system dynamics by deriving an approximate equation of motion for the $N$-particle distribution of particles’ positions. We note that since we retain the time-dependence that originates from the evolution of particles’ self-propulsions, this step is exact. The reason for introducing approximations is to make the formal expressions explicit. Thus, in principle, we are not restricted to fast evolution of self-propulsions. Only if we neglect the time delay caused by evolution of self-propulsions on finite time scale and introduce an effective passive system with Markovian dynamics we are assuming that the evolution of self-propulsions is fast compared to the evolution of particles’ positions.

We start by introducing a projection operator that acts on an $N$-particle probability distribution of self-propulsions and positions and projects it on a local steady-state distribution, i.e. on a distribution in which self-propulsions have a steady-state distribution for a given sample of positions,

[TABLE]

We note that by integrating $\mathcal{P}_{\text{lss}}P_{N}(\mathbf{r}_{1},\theta_{1},...,\mathbf{r}_{N},\theta_{N};t)$ over self-propulsions we get the probability distribution of particles’ positions, $P_{N}(\mathbf{r}_{1},...,\mathbf{r}_{N};t)$.

Next, we define the orthogonal projection,

[TABLE]

and write down equations of motion for $\mathcal{P}_{\text{lss}}P_{N}(\mathbf{r}_{1},\theta_{1},...,\mathbf{r}_{N},\theta_{N};t)$ and $\mathcal{Q}_{\text{lss}}P_{N}(\mathbf{r}_{1},\theta_{1},...,\mathbf{r}_{N},\theta_{N};t)$,

[TABLE]

[TABLE]

Since our goal is to calculate the intermediate scattering function, Eq. (16), which is a function of positions only, we can restrict ourselves to initial states that satisfy the following condition,

[TABLE]

Then we solve Eqs. (26-27) for the Laplace transform, $\mathcal{LT}$, of $\partial_{t}\mathcal{P}_{\text{lss}}P_{N}(\mathbf{r}_{1},\theta_{1},...,\mathbf{r}_{N},\theta_{N};t)$ and we obtain

[TABLE]

The first term inside the brackets on right-hand-side of Eq. (29) reads

[TABLE]

We see that if current densities vanish in the steady state, Eq. (III), the first term at the right-hand side of Eq. (30) vanishes and the second term can be re-written in terms of the effective steady state force,

[TABLE]

Furthermore, combining the assumption that currents vanish with Eq. (22) one can show that

[TABLE]

Similarly, one can show that

[TABLE]

To proceed, we will need to deal with the projected evolution operator $\mathcal{Q}_{\text{lss}}\Omega\mathcal{Q}_{\text{lss}}$ in Eq. (29). This operator describes evolution in the space orthogonal to the local steady-state space. The simplest possible approximation is to assume that this evolution is entirely due to the free relaxation of the self-propulsions. In this case $\mathcal{Q}_{\text{lss}}\Omega\mathcal{Q}_{\text{lss}}$ is approximated as follows

[TABLE]

Approximation (34) is equivalent to assuming that the relaxation in the space orthogonal to the local steady-state space is the same in non-interacting and interacting systems. In particular, the approximation (34) neglects the influence of the correlations between self-propulsions and positions on the evolution in the space orthogonal to the local steady-state space.

Combining approximation (34) with Eqs. (32-33) we get the following approximate equality

[TABLE]

Now, we expand $\left[z-D_{r}\sum_{j=1}^{N}\partial_{\theta_{i}}^{2}\right]^{-1}$ and integrate by parts. In this way we get the following expression for the second term inside the brackets on right-hand-side of Eq. (29),

[TABLE]

Combining the right-hand-sides of Eqs. (31) and (36) allows us to identify the effective evolution operator $\Omega^{\text{eff}}(z)$ governing the time-dependence of the distribution of particles’ positions,

[TABLE]

The dependence of the effective evolution operator $\Omega^{\text{eff}}(z)$ on $z$ allows us to retain, albeit in an approximate way, the finite relaxation time of the self-propulsions. This is in contrast with approximate approaches of Farage and Brader Farage and of Feng and Hou Feng2017 . In particular, the evolution operator corresponding to the former approach can be recovered from Eq. (37) by taking the $z\to 0$ limit, i.e. by neglecting the finite relaxation time of the self-propulsions, and at the same time neglecting local steady-state correlations between the self-propulsions, $\left<\mathbf{e}_{i}\mathbf{e}_{j}\right>_{\text{lss}}-\left<\mathbf{e}_{i}\right>_{\text{lss}}\left<\mathbf{e}_{j}\right>_{\text{lss}}\rightarrow\frac{1}{2}\delta_{ij}\mathbf{I}$, where $\mathbf{I}$ is the unit tensor.

The above described approximations become exact in the limit of fast evolution of the self-propulsions, $D_{r}\to\infty$, $v_{0}\to\infty$, $v_{0}^{2}/D_{r}=\text{ const}$. In this case, the ABP system becomes equivalent to a thermal system at temperature $T+v_{0}^{2}/(2D_{r}\xi_{0})$.

The finite time scale relaxation of the self-propulsions is retained in Eq. (37) in the simplest possible way. One could try using a more sophisticated approximation in place of Eq. (34), for example by following the procedure described in Appendix A of Ref. Szamel2016 . However, the resulting expressions would rapidly become rather difficult to handle. We note that, in principle, within the theory of Liluashvili et al. Voigtmann2017 the relaxation of self-propulsions is handled much more accurately commentlowd . However, additional technical approximations used in Ref. Voigtmann2017 make their description of the relaxation of self-propulsions similar to ours. Specifically, both in the theory of Ref. Voigtmann2017 and in our approach the matrix representing the self-propulsion part of the evolution operator is effectively approximated by its lowest order non-trivial elements.

The effective evolution operator $\Omega^{\text{eff}}(z)$ allows us to re-write the Laplace transform of the intermediate scattering function,

[TABLE]

Here $\left<...\right>_{\mathbf{r}}$ denotes averaging over the steady-state distribution of particles’ positions. Eq. (38) will be the starting point for the standard projection operator derivation of the memory function representation for $F(q;z)$ in subsection IV.3.

IV.2 Short-time dynamics and the importance of the correlations

of particles’ currents

Before turning to the memory function representation, which focuses on the long-time dynamics, we briefly examine the short-time dynamics of the intermediate scattering function and identify the contribution of the correlations of the particles’ currents.

To evaluate the short-time behavior of the intermediate scattering function we expand expression (16) in powers of $t$,

[TABLE]

The first term at the right-hand-side of Eq. (39) is the steady state static structure factor,

[TABLE]

the second term, using the assumption of vanishing currents, gives

[TABLE]

and the third term, also under the assumption of vanishing currents, reads

[TABLE]

where $\hat{\mathbf{q}}$ is a unit vector, $\hat{\mathbf{q}}=\mathbf{q}/q$. The first term in Eq. (IV.2) originates from thermal fluctuations while the second is the lowest order in time contribution of the self-propulsion.

Eq. (IV.2) shows that at short times the self-propulsion always speeds up the relaxation. The contribution of the self-propulsion can be expressed in terms of the correlation function of particles currents,

[TABLE]

The interpretation of function $\omega_{\parallel}(q)$ comes from the fact that

[TABLE]

can be interpreted as a current of particle $i$ averaged over thermal noise.

The same result for the short-time dynamics are obtained if one starts from expression (38) for the intermediate scattering function. In this case it is easiest to proceed in the Laplace space,

[TABLE]

Using the effective evolution operator (37) in Eq. (38) we obtain the following expression for the first derivative,

[TABLE]

where $\Omega_{t}$ is the translational diffusion part of the effective evolution operator,

[TABLE]

Furthermore, for the second derivative we get

[TABLE]

which can be shown to coincide with the second derivative obtained from the third term of the Taylor expansion, Eq. (IV.2).

IV.3 Memory function representation

In this subsection we rewrite the formal expression (38) for the intermediate scattering function in terms of the so-called frequency matrix and irreducible memory matrix. The latter quantity contains all the unknown non-trivial dynamic information about the system. The resulting expression for the density correlation function in terms of the frequency matrix and the memory matrix is known as the memory function representation. Some formal manipulations in this subsection are the same as in Sec. VI of Ref. Szamel2016 . They are included here for completeness.

To derive the memory function representation of $F(q;z)$ we use the projection operator approach Goetzebook ; CHess ; SL . We define a projection operator on the microscopic density

[TABLE]

We emphasize that projection operator $\mathcal{P}_{n}$ is defined in terms of the steady-state distribution, unlike in the approaches of Farage and Brader Farage and of Liluashvili et al. Voigtmann2017 . Next, we use the identity

[TABLE]

where $\mathcal{Q}_{n}$ is the projection on the space orthogonal to that spanned by the microscopic density,

[TABLE]

to rewrite the Laplace transform of the time derivative of $NF(q;t)$ in the following way

[TABLE]

The important part of the first term on the right-hand-side of the last equality sign in Eq. (IV.3) is the matrix element of the effective evolution operator, $\left<n(\mathbf{q})\Omega^{\text{eff}}(z)n(-\mathbf{q})\right>_{\mathbf{r}}$, which can be expressed in terms of the frequency matrix $\mathcal{H}(q;z)$,

[TABLE]

The frequency matrix is given by the following expression

[TABLE]

The important part of the second term at the right-hand-side of the last equality sign in Eq. (IV.3) can be expressed in terms of reducible CHess ; Kawasaki memory matrix $\mathcal{M}(q;z)$,

[TABLE]

Explicitly, the memory matrix is given by the following expression

[TABLE]

We can now rewrite the Laplace transform of the intermediate scattering function in terms of the frequency and memory matrix,

[TABLE]

where $S(q)$ is the steady-state structure factor,

[TABLE]

The second equality sign in Eq. (58) follows from the fact that for self-propulsion-independent quantities averaging over particles’ positions is equivalent to averaging over the full steady-state distribution of positions and self-propulsions.

Next, following Cichocki and Hess CHess and Kawasaki Kawasaki we introduce an irreducible memory matrix. First, we define the irreducible evolution operator $\Omega^{\mathrm{irr}}(z)$,

[TABLE]

where the subtraction term $\delta\Omega^{\text{irr}}(z)$ reads

[TABLE]

Next, we define the irreducible memory matrix $\mathcal{M}^{\text{irr}}(q;z)$, which is given by the expression analogous to Eq. (IV.3) but with the projected evolution operator $\mathcal{Q}_{n}\Omega^{\text{eff}}(z)\mathcal{Q}_{n}$ replaced by irreducible evolution operator $\Omega^{\text{irr}}(z)$,

[TABLE]

Finally, we use an identity similar to Eq. (IV.3),

[TABLE]

and we derive the following relation between $\mathcal{M}(q;z)$ and $\mathcal{M}^{\text{irr}}(q;z)$,

[TABLE]

Combining Eqs. (57) and (63) we arrive at the following representation of the intermediate scattering function in terms of the irreducible memory matrix,

[TABLE]

Eq. (64) constitutes the memory function representation for the intermediate scattering function of the ABP system. It has the same structure as the memory function representation for the intermediate scattering function of the AOUP model. The difference between the two models is buried in the expressions for the frequency and irreducible memory matrices.

IV.4 Mode-coupling approximation for the irreducible memory matrix

Here, we derive an explicit approximate expression for the irreducible memory matrix using a factorization approximation of the type used in the mode-coupling theory of the glass transition. To this end we follow the steps of the derivation of the mode-coupling theory for systems evolving with Brownian dynamics SL . The derivation consists of three steps.

First, we project onto the subspace of density pairs,

[TABLE]

Here $n_{2}(\mathbf{q}_{1},\mathbf{q}_{2})$ is the Fourier transform of the microscopic two-particle density,

[TABLE]

and $\left[\left<\mathcal{Q}_{n}n_{2}(\mathbf{q}_{1},\mathbf{q}_{2})\mathcal{Q}_{n}n_{2}(-\mathbf{q}_{3},-\mathbf{q}_{4})\right>_{\mathbf{r}}\right]^{-1}$ is the inverse of the correlation matrix of microscopic pair densities.

Second, we factorize averages resulting from substituting projection (65) into the expression for the memory function and at the same time replace the irreducible operator $\Omega^{\text{irr}}(z)$ by effective evolution operator $\Omega^{\text{eff}}(z)$. We should emphasize that this factorization has to be done in the time domain,

[TABLE]

Here $\mathcal{LT}^{-1}$ denotes the inverse Laplace transform and $\left\{3\leftrightarrow 4\right\}$ means the preceding expression with labels $3$ and $4$ interchanged. Consistently with Eq. (IV.4) we also factorize the steady-state correlation matrix of microscopic pair densities and for its inverse we get

[TABLE]

Third, we approximate the vertex functions. Due to the presence of the current correlations in one part of the vertex, this last step is a bit more complex than the approximation used in the derivation of the standard mode-coupling theory SL . We will explain it on the example of the left vertex, $\mathcal{V}_{l}$.

The left vertex is given by the following formula

[TABLE]

According to Eq. (IV.4) the left vertex consists of two terms. The first term originates from the translational diffusion part of the effective evolution operator. It has the same form as the vertex of the standard mode-coupling theory SL ,

[TABLE]

where $c(q)$ is the direct correlation function HansenMcDonald

[TABLE]

We note that in the derivation of formula (IV.4) a convolution approximation is used, which is equivalent to neglecting the three-particle direct correlation function. The same approximation is used in the derivation of the standard mode-coupling theory SL .

The second term originates from the activity-related part of the effective evolution operator. It is proportional to function $\mathcal{V}_{l}^{a}$, which is defined as follows,

[TABLE]

$\mathcal{V}_{l}^{a}(\mathbf{q};\mathbf{q}_{1},\mathbf{q}_{2})$ has the form very similar to the form of the vertex of the mode-coupling theory for the AOUP model Szamel2016 and can be analyzed in the same way. The result is

[TABLE]

where a new function $\mathcal{C}(q)$ reads

[TABLE]

Again, we note that the derivation of formula (IV.4) involves a generalization of the convolution approximation introduced in Ref. Szamel2016 .

Combining Eqs. (70) and (IV.4) we get the following approximate expression for the left vertex,

[TABLE]

The right vertex can be analyzed in the same way resulting in

[TABLE]

We note that both $\mathcal{V}_{l}$ and $\mathcal{V}_{r}$ are $z-$ or time-dependent. The time dependence originates from the finite relaxation rate of the self-propulsion.

To simplify the formulae below, we factor out $\pm iN\rho S(q_{1})S(q_{2})$ from the left and right vertex and we define the vertex function $\mathcal{V}(\mathbf{q};\mathbf{q}_{1},\mathbf{q}_{2})$, which is analogous to the standard vertex functions of mode-coupling theories,

[TABLE]

Now, combining the three steps constituting the mode-coupling approximation, expressing the vertices in terms of $\mathcal{V}^{t}$ and $\mathcal{V}^{a}$ and taking the thermodynamic limit we arrive at the following expression for the irreducible memory matrix,

[TABLE]

where $\ast$ denotes convolution in the time domain.

Formula (IV.4) describes the standard mode-coupling dynamic feedback mechanism: the time-delayed internal friction arising due to interparticle interactions decays due to the relaxation of the two-particle density, which is included at the level of factorization approximation (IV.4). Through formula (IV.4), slow decay of the density fluctuations feeds back into slow decay of the irreducible memory matrix, which quantifies the internal friction.

The non-equilibrium nature of the active system is included at two levels. First, since the microscopic dynamics consists of random thermal motion and motion due to the independently evolving self-propulsion, the vertices exhibit two contributions corresponding to these two microscopic modes of motion. Second, the contribution due to the self-propulsion $\mathcal{V}^{a}$, involves equal-time correlation function of particles’ currents, $\omega_{\parallel}(q)$. This function sets the overall scale of the self-propulsion contribution (together with $D_{r}$). It also enters into the expression for the new function $\mathcal{C}(q)$, Eq. (74), which plays the role of the direct correlation function in $\mathcal{V}^{a}$.

IV.5 Long-time dynamics and ergodicity-breaking transition

To obtain the time dependence of the intermediate scattering function one needs to solve the combined set of Eqs. (64) and (IV.4). In the case of AOUPs, it is possible to rewrite the analogous set of equations and to get a self-consistent equation of motion for the intermediate scattering function that has the form very similar to that of the mode-coupling equation for the intermediate scattering function of an under-damped thermal colloidal system Szamel2016 . In the present case, with microscopic dynamics due to two independent mechanisms, this does not seem possible. The exception is the asymptotic long-time dynamics close to the ergodicity-breaking transition, which is discussed below.

Generically, equations similar to (64) and (IV.4) predict that as equal-time correlations grow (which can happen by lowering the temperature, increasing the number density, or, in the present case, by manipulating the self-propulsion), the memory matrix grows and leads to the slowing down of the time evolution of the intermediate scattering function. Eventually, the the equations predict an ergodicity breaking transition, at which the intermediate scattering function ceases decaying to zero. Close to the transition and at long times (which corresponds to small Laplace variable $z$), Eqs. (64) and (IV.4) can be approximated by a simpler set of equations. These asymptotic equations are usually expressed in terms of the so-called normalized correlator, $\phi(q;t)=F(q;t)/S(q)$, and normalized irreducible memory function, $m(q;t)$. They have the following form,

[TABLE]

[TABLE]

Eqs. (79-IV.5) have the structure very similar to that of the self-consistent equations of motion for the long-time dynamics near the ergodicity breaking transition described by the standard mode-coupling theory. This implies that all analytical results based of the standard mode-coupling theory can be used for the present theory for the dynamics of the ABP system. In particular, the only quantity that one needs to calculate in order to predict mode-coupling exponents is the so-called exponent parameter $\lambda$ Goetzebook . This parameter can be calculated from the solution of the self-consistent equations for the order parameter at the ergodicity breaking transition.

The reduced memory function $m(q;z)$ differs from the corresponding quantity of the standard mode-coupling theory by the fact that the vertex involves a weighted average of the direct correlation function and the new function $\mathcal{C}(q)$, which involves $S(q)$ and $\omega_{\parallel}(q)$. We emphasize that it is the non-equilibrium steady-state structure factor that enters into Eq. (IV.5) and determines $c(q)$ and $\mathcal{C}(q)$. However, unlike in the standard mode-coupling theory, the structure factor itself does not completely determine the system’s dynamics, since the memory function involves also the correlation function of particles’ currents $\omega_{\parallel}(q)$.

Finally, let us assume that the system, as described by Eqs. (79-IV.5), undergoes an ergodicity breaking transition. At such a transition, the normalized correlator ceases decaying to zero. The long-time limit of the correlator, $\lim_{t\to\infty}\phi(q;t)=f(q)$, is the order parameter of the non-ergodic state. Eqs. (79-IV.5) lead to the following set of self-consistent equations for the order parameter, $f(q)$,

[TABLE]

where $m(q)$ is given by the following equation

[TABLE]

Again, self-consistent equations (81-IV.5) for the order parameter are very similar to the equations derived in the standard mode-coupling theory. The non-equilibrium character of the ABP system manifests itself in the second equation.

We note that, unlike in the integration-through-transients theory of Liluashvili et al. Voigtmann2017 , the ergodicity-breaking transition is determined only by equal-time quantities characterizing the steady state of the ABP system.

Finally, we recall that in the limit of rapidly relaxing self-propulsion, $D_{r}\to\infty$, $v_{0}\to\infty$, $v_{0}^{2}/D_{r}=\text{ const.}$, the ABP system becomes equivalent to a thermal system at temperature equal to the the single-particle effective temperature, $T_{\text{eff}}=T+v_{0}^{2}/\left(2D_{r}\xi_{0}\right)$. We note that in this limit the static structure factor becomes equal to the equilibrium structure factor at $T_{\text{eff}}$. Moreover, $\mathcal{C}(q)$ becomes equal to the direct correlation function. With these two changes, Eqs. (79-IV.5) coincide with standard mode-coupling equations for a thermal system at temperature $T_{\text{eff}}$.

V Discussion

We presented here a theory for the steady state dynamics of dense systems of active Brownian particles. The derivation identified a function that quantifies correlations of steady state particles’ currents, which influences both short- and long-time dynamics of the ABP system.

The present theory for the dynamics of ABP systems relies upon a factorization approximation. Thus, like our previous theory for the dynamics of AOUP systems, it belongs to the class of mode-coupling-like theories. Our theory predicts that the dynamics upon approaching an ergodicity breaking transition is qualitatively similar to that predicted by the standard mode-coupling theory close to the corresponding mode-coupling transition. Quantitative details, including the location of the transition, the exponent parameter $\lambda$ and the mode-coupling exponents, depend on both the steady-state structure factor and the correlation function of steady-state currents.

In the limit of vanishing self-propulsion the theory reduces itself to the standard mode-coupling theory for colloidal glassy dynamics and the colloidal glass transition. In the limit of vanishing thermal fluctuations the theory becomes equivalent to our earlier theory for the glassy dynamics of the AOUP model (at the level of our description of the time evolution of the self-propulsion these models are equivalent). Finally, in the limit of rapidly varying self-propulsion the theory becomes equivalent to the standard mode-coupling theory for the thermal system at a temperature equal to the effective temperature. On the other hand, our effective evolution operator (37) becomes ill-defined in the limit of very slowly evolving self propulsion, $D_{r}\to 0$. We note that the last limit is rather difficult since none of the present theories for the dynamics of dense active systems seem to work well in this limit, with the exception of the theory of Liluashvili et al.

According to the present theory, the location of the active ergodicity-breaking transition is determined by equal-time steady-state correlation functions only. This is in contrast with the approach of Liluashvili et al., in which the equation determining the ergodicity-breaking transition involves also a time integral over the decaying elements of the matrix correlator, which is the basic fundamental quantity of that theory.

The main fundamental assumption behind the present theory is the vanishing of stead-state currents after integration over the self-propulsions. Preliminary analytical results suggest that the presence of the currents may lead to wiping out the ergodicity-breaking transition. A combined analytical and simulational study in this direction is planned for the near future.

Finally, we note that, somewhat unexpectedly, the theory of Liluashvili et al. predicts that in addition to density correlations also some cross-correlations between the density and the self-propulsion slow down and get arrested at the ergodicity-breaking transition. While our theory cannot address such cross-correlation functions directly, we feel that this is another area where a combined analytical and simulational study would be very interesting.

Acknowledgments

I thank Elijah Flenner for comments on the manuscript. I gratefully acknowledge the support of NSF Grant No. CHE 1800282.