Mode Coupling Theory for Nonequilibrium Glassy Dynamics of Thermal Self-Propelled Particles

TL;DR

This paper develops a mode coupling theory to analyze the nonequilibrium glassy dynamics of thermal self-propelled particles, revealing how propulsion and temperature influence the glass transition and dynamics.

Contribution

It introduces a theoretical framework for nonequilibrium glassy dynamics of self-propelled particles, extending mode coupling theory to active matter with thermal noise.

Findings

Critical density increases with propulsion force and effective temperature.

Critical effective temperature rises with persistence time of self-propulsion.

Theory recovers passive system results and can extend to complex active-passive mixtures.

Abstract

We present a promising mode coupling theory study for the relaxation and glassy dynamics of a system of strongly interacting self-propelled particles, wherein the self-propulsion force is described by Ornstein-Uhlenbeck colored noise and thermal noises are included. Our starting point is an effective Smoluchowski equation governing the distribution function of particle's positions, from which we derive a memory function equation for the time dependence of density fluctuations in nonequilibrium steady states. With the basic assumption of absence of macroscopic currents and standard mode coupling approximation, we can obtain expressions for the irreducible memory function and other relevant dynamic terms. With these equations obtained, we study the glassy dynamics of this thermal self-propelled particles system by investigating the Debye-Waller factor f_{q} and relaxation time \tau_{\alpha} as functions of the persistence time \tau_{p} of self-propulsion, the single particle effective temperature T_{\text{eff}} as well as the number density \rho. Consequently, we find the critical density \rho_{c} for given \tau_{p} shifts to larger values with increasing magnitude of propulsion force or effective temperature, in good accordance with previous reported simulation works. In addition, the theory facilitates us to study the critical effective temperature T_{\text{eff}}^{c} for fixed \rho as well as its dependence on \tau_{p}. We find that T_{\text{eff}}^{c} increases with \tau_{p} and in the limit \tau_{p}\to0, it approaches the value for a simple passive Brownian system as expected. Our theory also well recovers the results for passive systems and can be easily extended to more complex systems such as active-passive mixtures.

Equations

$${\dot{\mathbf{r}}}_i\left(t\right)$$

$$⟨{\mathbf{\xi }}_i\left(t\right){\mathbf{\xi }}_j\left(t^{\prime }\right)⟩=2D_{t}1{\delta }_{ij}\delta \left(t-t^{\prime }\right),$$

$$\dot{{\mathbf{f}}_i}\left(t\right)=-{\tau }_p^{-1}{\mathbf{f}}_i\left(t\right)+\mathbf{\eta }{}_i(t)$$

$$⟨{\mathbf{\eta }}_i\left(t\right){\mathbf{\eta }}_j\left(t^{\prime }\right)⟩=2D_{f}1{\delta }_{ij}\delta \left(t-t^{\prime }\right)$$

$$⟨{\mathbf{f}}_i(t){\mathbf{f}}_j(t^{\prime })⟩=D_f{\tau }_p{e}^{-|t^{\prime }-t|/{\tau }_p}1{\delta }_{ij}$$

$$⟨\delta r^2\left(t\right)⟩=\frac{6D_f{\tau }_p^2}{{\gamma }^2}\left[t+{\tau }_p\left(e^{-t/{\tau }_p}-1\right)\right]+6D_{t}t$$

$$D_0=D_t+D_f{\tau }_p^2/{\gamma }^2$$

$$T_{\text{eff}}=T+k_B{D}_f{\tau }_p^2/\gamma =\left(\frac{D_0}{D_t}\right)T.$$

$$⟨{\mathbf{p}}_i\left(t\right){\mathbf{p}}_j\left(t^{\prime }\right)⟩\simeq \frac{1}{3}{e}^{-2D_r|t-t^{\prime }|}1{\delta }_{ij}$$

$$u\left(r\right)=\{\begin{array}{ll}4\epsilon \left[{\left(\frac{\sigma }{r}\right)}^{12}-{\left(\frac{\sigma }{r}\right)}^6\right]+\epsilon & r\le 2^{1/6}\sigma \\ 0& r>2^{1/6}\sigma \end{array}$$

$$\frac{\partial }{\partial t}\Psi \left({\mathbf{r}}^N,t\right)=\hat{\Omega }\Psi \left({\mathbf{r}}^N,t\right)$$

$$\hat{\Omega }=\sum _{j=1}^N{\nabla }_j\cdot D_j({\mathbf{r}}^N)\left[{\nabla }_j-\beta {\mathbf{F}}_j^{\mathrm{eff}}\left({\mathbf{r}}^N\right)\right].$$

$$D_j\left({\mathbf{r}}^N\right)=D_t+\frac{D_f{\tau }_p^2/{\gamma }^2}{1-{\tau }_p\beta D_t{\nabla }_j\cdot {\mathbf{F}}_j^{\text{mbox}}}.$$

$${\mathbf{F}}_j^{\text{mbox}eff}\left({\mathbf{r}}^N\right)=\frac{D_t}{D_j\left({\mathbf{r}}^N\right)}\left[{\mathbf{F}}_j\left({\mathbf{r}}^N\right)-\frac{1}{\beta D_t}{\nabla }_j{D}_j\left({\mathbf{r}}^N\right)\right]$$

$${\hat{\Omega }}_{{\tau }_{p\to 0}}=D_0\sum _{j=1}^N{\nabla }_j\left({\nabla }_j-\beta \frac{D_t}{D_0}{\mathbf{F}}_j\right)=D_0\sum _{j=1}^N{\nabla }_j\left({\nabla }_j-{\beta }_{\text{eff}}{\mathbf{F}}_j\right)$$

$$\hat{\Omega }{P}_s\left({\mathbf{r}}^N\right)=-\sum _i{\nabla }_i\cdot {\mathbf{J}}_i^s=0$$

$${\mathbf{J}}_i^s=-D_j\left({\mathbf{r}}^N\right)\left[{\nabla }_j-\beta {\mathbf{F}}_j^{\mathrm{eff}}\left({\mathbf{r}}^N\right)\right]P_s\left({\mathbf{r}}^N\right)$$

$${\hat{\Omega }}^{†}=\sum _{j=1}^N\left({\nabla }_j+\beta {\mathbf{F}}_j^{\mathrm{eff}}\right)D_j\left({\mathbf{r}}^N\right)\cdot {\nabla }_j$$

$$F_q\left(t\right)=\frac{1}{N}⟨{\text{ρ}}_{\mathbf{q}}^{\ast }\left(e^{{\hat{\Omega }}^{†}t}{\rho }_{\mathbf{q}}\right)⟩=\frac{1}{N}⟨{\text{ρ}}_{-\mathbf{q}}\left(e^{{\hat{\Omega }}^{†}t}{\rho }_{\mathbf{q}}\right)⟩$$

$${\rho }_{\mathbf{q}}=\sum _{j=1}^N{e}^{-i\mathbf{q}\cdot {\mathbf{r}}_j}$$

$$F_q\left(0\right)=\frac{1}{N}⟨{\text{ρ}}_{-\mathbf{q}}{\rho }_{\mathbf{q}}⟩=S\left(q\right)$$

$$F_q\left(t\right)=\frac{1}{N}⟨{\text{ρ}}_{-\mathbf{q}}{e}^{\hat{\Omega }t}{\rho }_{\mathbf{q}}⟩$$

$$F_q^s\left(t\right)$$

$$=$$

$$\frac{\partial }{\partial t}{F}_q(t)+{\omega }_q{F}_q(t)+{\int }_0^{t}du{\tilde{M}}^{\mathrm{irr}}\left(q,t-u\right)\frac{\partial }{\partial u}{F}_q(u)=0$$

$$\omega_{q}=-\left\langle\rho_{\mathbf{q}}^{*}\left(\hat{\mbox{$\Omega$}}^{\dagger}\rho_{\mathbf{q}}\right)\right\rangle\left\langle\rho_{\mathbf{q}}^{*}\rho_{\mathbf{q}}\right\rangle^{-1}=\frac{q^{2}\sum_{j}\left\langle D_{j}\left(\mathbf{r}^{N}\right)\right\rangle}{NS(q)}=\frac{q^{2}\bar{D}}{S\left(q\right)}$$

$$\bar{D}=N^{-1}\sum _j⟨D_j\left({\mathbf{r}}^N\right)⟩.$$

$${\tilde{M}}^{\mathrm{irr}}\left(q,t\right)=\frac{\rho \bar{D}}{16{\pi }^3}\int d\mathbf{k}{\left[\left(\hat{\mathbf{q}}\cdot \mathbf{k}\right)C_2\left(\mathbf{q};\mathbf{k}\right)+\left(\hat{\mathbf{q}}\cdot \mathbf{p}\right)C_2\left(\mathbf{q};\mathbf{k}\right)\right]}^2{F}_k\left(t\right)F_p\left(t\right).$$

$$C_2\left(\mathbf{q};\mathbf{k}\right)={\rho }^{-1}\left[1-\frac{D_0}{\bar{D}}\frac{S_2\left(p\right)}{S\left(p\right)}{S}^{-1}\left(k\right)\right]$$

$$S_2(k)=\frac{1}{D_{0}N}⟨\sum _{i,j}{D}_j\left({\mathbf{r}}^N\right)e^{-i\mathbf{k}\cdot {\mathbf{r}}_j+i\mathbf{k}\cdot {\mathbf{r}}_i}⟩$$

$${\tilde{M}}^{irr}\left(q,t\right)=\frac{\rho D_t}{16{\pi }^3}\int d\mathbf{k}{\left[\left(\hat{\mathbf{q}}\cdot \mathbf{k}\right)c\left(k\right)+\left(\hat{\mathbf{q}}\cdot \mathbf{p}\right)c\left(k\right)\right]}^2{F}_k\left(t\right)F_p\left(t\right)$$

Full text

I Introduction

The collective behaviors of systems containing active (self-propelled) particles have gained extensive attention in recent years due to its great importance both from a fundamental physics perspective and for understanding many biological systems Vicsek and Zafeiris (2012); Marchetti and Joanny (2013); Bechinger et al. (2016). A wealth of new nonequilibrium phenomena have been reported, such as active swarming, large scale vortex formationSchaller et al. (2010); Sumino1 et al. (2012), phase separationSchwarz-Linek et al. (2012); Fily and Marchetti1 (2012); Buttinoni et al. (2013); Wittkowski et al. (2014); Das et al. (2014); Bechinger et al. (2016), etc, both experimentally and theoretically. Recently, a new trend in this field has been the glassy dynamics and the glass transition in dense assemblies of self-propelled particles and their comparisons to their corresponding phenomena in equilibrium systemsBerthier and Kurchan (2013); Ni et al. (2013); Berthier (2013); Mandal et al. (2014); Berthier (2014); Sharma et al. (2017). Experiments demonstrated that active fluids may show dynamic features such as jamming and dynamic arrest that are very similar to those observed in glassy materials. For instance, migrating cells exhibited glassy dynamics, such as dynamic heterogeneity as the cell density increasesAngelini et al. (2011), amorphous solidification process was found in collective motion of a cellular monolayerKrakoviack (2005), and glassy behaviors could even be found for ant aggregates in large scalesWhitesides and Grzybowski (2002), to list just a few.

Computer simulations also demonstrated that nonequilibrium glass transition or dynamic arrest behavior does occur in dense suspensions of self-propelled particles. So far, mainly two types of self-propulsion particle systems have been studied. One is the Rotation diffusional Active Brownian (RAB) particles system, where each particle is subjected to a self-propulsion force with constant amplitude $v_{0}$ but a randomly changing direction evolving $via$ rotational diffusion with diffusion coefficient $D_{r}$. Ni et al. ten Hagen et al. (2011); Ni et al. (2013, 2014)studied the glassy behavior of this RAB system of hard-sphere particles, finding that the critical density for glass transition shifts to larger density as the active force increases, thus pushing the glass transition point to the limit of random packing. The other one is the so-called Active Ornstein-Uhlenbeck(AOU) particles system, wherein the self-propulsion force is realized by a colored noise described by the OU process. In contrast to the RAB model, thermal noise is ignored in the AOU system, such that the system can never reach the equilibrium state determined by the canonical distribution. For this athermal system, an effective temperature $T_{\text{eff}}$ can be introduced to quantify the strength of the self-propulsion force, and a persistence time $\tau_{p}$ controls the duration of the persistent self-propelled motion. Berthier and co-workersBerthier (2014); Levis and Berthier (2014); Flenner et al. (2016) had performed detailed studies about the structural and glassy dynamics of this AOU model in two and three dimensions. Similarly, the glass transition shifts to larger densities compared to the equilibrium one when the magnitude or persistence time of Cthe self-propulsion force increases. Besides the studies of self-propelled particles, using molecular dynamics simulations of a model glass former, Mandal et al.Mandal et al. (2014) showed that the incorporation of activity or self-propulsion can induce cage breaking and fluidization, resulting in a disappearance of the glassy phase beyond a critical force. And related to the glassy dynamics, it was shown that particle activity can shift the freezing density to larger valuesJulian Bialke (2012) and particularly, hydrodynamic interactions can further enhance this effectLi et al. (2015).

Besides experimental and simulation studies mentioned above, on the theoretical side, important progresses have also been achieved in recent yearsDing et al. (2015). Starting from a generalized Langevin equation with colored non-thermal noise, Berthier and KurchanBerthier (2013) predicted that dynamic arrest can occur in systems that are far from equilibrium, showing that non-equilibrium glass transition moves to lower temperature with increasing activity and to higher temperature with increasing dissipation in spin glasses. Farage and Brader attempted to extend the mode coupling theory to the limiting case of the RAB modelFarage and Brader (2014), wherein activity leads to a higher effective particle diffusivity. They then started from the effective Smoluchowski equation governing the many-particle distribution function and obtained a memory function equation for the equilibrium intermediate scattering function, showing that self-propulsion could shift the glass transition to larger densityFarage et al. (2015). In a recent work, Nandi proposed a phenomenological extension of random first order transition theory to study glass transition of the RAB model, showing that more active systems are stronger glass formers Nandi (2016). Very recently, Szamel et al.Berthier (2014); Szamel et al. (2015) presented an elegant theoretical modeling of the structure and glassy dynamics of the athermal AOU system. In their approach, they first integrated out the self-propulsion and then used the projection operator method and a mode-coupling-like approximation to derive an approximate equation of motion for the collective intermediate scattering functions, defined upon the nonequilibrium steady state distribution. In particular, this work highlighted the importance of the steady state correlations of particle velocities, which played a crucial role to understand the relaxation dynamics of the system. Nevertheless, extension of this framework to more general cases with thermal noise included is not present yet.

In the present work, we have developed an alternative mode coupling theory for the nonequilibrium glass dynamics of a general system of active particles, wherein the self-propulsion force is described by an OU process, and besides, thermal noise is included. To be specific, we call this active OU particles system with thermal noise presented as the AOU-T system, where ’T’ apparently stands for thermal noise. Our starting point is an effective Smoluchowski equation(SE) obtained via Fox approximation method, which was recently adopted by Farage et al. to study the effective interactions among RAB particlesFarage et al. (2015). This effective SE allows us to derive a memory function equation for the nonequilibrium steady state collective intermediate scattering function $F_{q}\left(t\right)$, as well as that for the self-intermediate scattering function $F_{q}^{s}\left(t\right)$. With the basic assumptions of absence of macroscopic currents and standard mode coupling approximation, we are able to get the expressions for the irreducible memory functions and other relevant variables. Particularly, we find that the dynamics is governed by an averaged diffusion coefficient $\bar{D}$ and a nontrivial steady state structure function $S_{2}\left(q\right)$, both depending on the effective temperature $T_{\text{eff}}$, the persistence time $\tau_{p}$ as well as the number density $\rho$. With $\bar{D}$, $S_{2}\left(q\right)$ and the nonequilibrium structure factor $S\left(q\right)$ as inputs, we can calculate $F_{q}\left(t\right)$ and $F_{q}^{s}\left(t\right)$ for different parameter settings and investigate the glass transition behaviors. We calculate the critical density $\rho_{c}$ for glass transition as a function of the effective temperature $T_{\text{eff}}$, the magnitude of propulsion force $v_{0}$ and the persistence time $\tau_{p}$, by investigating the Debye-Waller factor $f_{q}$ as well as the relaxation time $\tau_{\alpha}$. Consequently, $\rho_{c}$ shifts to larger values with increasing propulsion force $v_{0}$ or effective temperature $T_{\text{eff}}$, in good accordance with previous simulation works. The theory also facilitates us to study the critical temperature $T_{\text{eff}}^{c}$ for fixed number density $\rho$, as well as its dependence on $\tau_{p}$.

The remainder of paper is organized as follows. In Section II, we present descriptions of the AOU-T model and our theory, the latter being the main result of the present work. Section III includes the numerical results predicted by our theory and conclusions in Section IV.

II Model and Theory

II.1 AOU-T Model for self-propelled particles

We consider a system of $N$ interacting, self-propelled particles in a volume $V$. The particles move in a viscous medium with single particle friction coefficient $\gamma$ and hydrodynamic interactions are neglected. As mentioned in the introduction, the self-propulsion force is given by a colored noise described by OU process and thermal noises are included. The equations of motion for these AOU-T particles are thus given by

[TABLE]

where $\mathbf{r}_{i}$ denotes the position vector of particle $i$, the force $\mathbf{F}_{i}=-\sum_{j\neq i}\nabla_{i}u\left(r_{ij}\right)$ originates from the interactions where is the pair-potential $u\left(r_{ij}\right)$, $\mathbf{f}_{i}$ is the self-propulsion force, and $\bm{\xi}_{i}\left(t\right)$ is the thermal noise with zero mean and variance

[TABLE]

where $D_{t}=k_{B}T/\gamma$ with $k_{B}$ the Boltzmann constant and $T$ the ambient temperature, and $\bm{1}$ denotes the unit tensor. The equations of motion for the self-propulsion force $\mathbf{f}_{i}$ are given by

[TABLE]

where $\tau_{p}$ is the persistence time of self-propulsion and $\bm{\eta}{}_{i}(t)$ is a Gaussian white noise with zero mean and variance $D_{f}$

[TABLE]

Accordingly, the correlation function of the force $\mathbf{f}_{i}$ reads

[TABLE]

For an isolated particle, the mean square displacement can be obtained as

[TABLE]

In the short time limit $t\ll\tau_{p}$, the particle’s motion is diffusive with $\left\langle\delta r^{2}\left(t\right)\right\rangle=6D_{t}t$, which is at variance with the AOU system, wherein the term $6D_{t}t$ is absent and the motion is ballistic $\left\langle\delta r^{2}\left(t\right)\right\rangle=3D_{f}\tau_{p}\gamma^{-2}t^{2}$ for $t\ll\tau_{p}$. For long time, the motion is also diffusive but with $\left\langle\delta r^{2}\left(t\right)\right\rangle=6\left(D_{t}+D_{f}\tau_{p}^{2}/\gamma^{2}\right)t$, implying that the long diffusion coefficient is given by

[TABLE]

This allows us to introduce a single-particle effective temperature

[TABLE]

In the limit of vanishing $\tau_{p}\rightarrow 0$, $\left\langle\mathbf{f}_{i}(t)\mathbf{f}_{j}(t^{\prime})\right\rangle\rightarrow 2D_{f}\tau_{p}^{2}\bm{1}\delta\left(t-t^{\prime}\right)\delta_{ij}$ and the system becomes equivalent to a Brownian system at the effective temperature $T_{\text{eff}}$ with diffusion coefficient $D_{0}$. For the AOU model where the thermal noise is absent, the effective temperature is simply given by $T_{\text{{eff}}}=D_{f}\tau_{p}^{2}/\gamma$Szamel (2014) .

We note here that the AOU-T model can be mapped onto the RAB model at a coarse-grained levelFarage et al. (2015). The equation of motion for the RAB particles is $\dot{{\bf r}}_{i}=v_{0}{\bf p}_{i}+\gamma^{-1}{\bf F}_{i}+\bm{\xi}_{i}$, where $v_{0}$ denotes the magnitude of the propulsion force and $\mathbf{p}_{i}$ is the unit vector of the direction of particle $i$. $\mathbf{p}_{i}$ changes randomly with time via rotational diffusion, $\dot{{\bf p}}_{i}=\bm{\zeta}_{i}\times{\bf p}_{i}$, where $\bm{\zeta}_{i}$ is also a Gaussian white noise with zero mean and correlation $\left\langle\bm{\zeta}_{i}\left(t\right)\bm{\zeta}_{j}\left(t^{\prime}\right)\right\rangle=2D_{r}\bm{1}\delta_{ij}\delta\left(t-t^{\prime}\right)$ where $D_{r}$ denotes the rotational diffusion coefficient. Recently, it was shown that $\mathbf{p}_{i}$ can be approximated by a colored noise with persistence time $\tau_{p}=\left(2D_{r}\right)^{-1}$ if one average over the angular degree of freedom, i.e.,

[TABLE]

Comparing Eq.(9) with (5), we see that $v_{0}\mathbf{p}_{i}$ has same correlation property as $\mathbf{f}_{i}$ by the mapping $\tau_{p}\to\left(2D_{r}\right)^{-1}$ and $D_{f}\tau_{p}\to v_{0}^{2}/3$.

In the studies of the AOU model, the authors usually used $T_{\text{eff}}$ as one of the independent parameters together with the persistence time $\tau_{p}$ and number density $\rho=N/V$. In the simulation works of the RAB model, however, the authors often used the magnitude of the propulsion force $v_{0}$ as an independent parameter. Note that in the simulation work performed by Ran Ni et al.Ni et al. (2013), they have adopted Stokes-Einstein relation to set the rotational diffusion coefficient $D_{r}=3D_{t}/\sigma^{2}$ where $\sigma$ is the particle diameter. In the dimensionless version by setting $D_{t}=1$, $\gamma=1$ and $\sigma=1$, this means that $D_{r}$ is fixed and the persistence time is $\tau_{p}=\left(2D_{r}\right)^{-1}=1/6$. In more general cases, the coupling of $D_{t}$ and $D_{r}$ may not hold and one can thus set $\tau_{p}$ as a free independent parameter. In the present work, we will set $v_{0}=\sqrt{3D_{f}\tau_{p}}$ or $T_{\text{eff}}$, $\tau_{p}$, and $\rho$ as independent parameters if not otherwise specified.

For simplicity, we consider a one-component pure-repulsive Lenard-Jones(LJ) system of self-propelled particles. The pair potential is given by

[TABLE]

where $\varepsilon$ is the strength of the potential. Here we set $\epsilon=1k_{B}T$, where $k_{B}T=D_{t}\gamma$ is the unit of energy. Moreover, $\sigma\gamma/k_{B}T$ is the unit of time. The number density $\rho$ is set to be large enough such that the phase separation dynamics is not relevant and we mainly focus on the glassy dynamics. Simulations are performed in a cubic box with $L=10$ and periodic boundary conditions.

II.2 Effective Smoluchowski Equation

The AOU-T model described in Eq.(1) is non-Markovian due to the colored noise term $\mathbf{f}_{i}$. Consequently, it is not possible to derive an exact Fokker-Planck equation (FPE) for the time evolution of the probability distribution function $\Psi\left(\mathbf{r}^{N},t\right)$, which gives the probability that the system is at a specific configuration $\mathbf{r}^{N}=\left(\mathbf{r}_{1},\mathbf{r}_{2},\dots,\mathbf{r}_{N}\right)$ at time $t$. Nevertheless, one may obtain an approximate FPE for such a colored noise system by applying the method introduced by FoxFox (1986), where a perturbative expansion in powers of correlation time is partially resumed using functional calculus. The resulting FPE thus defines implicitly a Markovian process, and is shown to be rather accurate for short correlation time of the colored noise. Very recently, T. Farage and co-workers had applied such a method to the RAB model and obtained an effective FPE for $\Psi\left(\mathbf{r}^{N},t\right)$, and analyzed the effective interaction among active particles in the low density limit. Since the FPE only involves the distribution in the configuration space $\mathbf{r}^{N}$, it is also known as Smoluchowski equation (SE).

Here we use the same method to obtain the approximate SE of the AOU-T model. The procedure is similar to that in Ref.Farage et al. (2015) noting that the AOU-T model can be mapped to the RAB model as discussed in the last subsection. For self-consistency, the main steps with necessary illustrations are given in Appendix A. Finally, we can obtain the effective SE as

[TABLE]

where $\hat{\Omega}$ denotes the effective Smoluchowski operator given by

[TABLE]

Herein, $D_{j}\left(\mathbf{r}^{N}\right)$ is a configuration-dependent instantaneous diffusion coefficient of particle $j$ which is given by

[TABLE]

with $\beta=\left(k_{B}T\right)^{-1}$. $\mathbf{F}_{j}^{\text{eff}}\left(\mathbf{r}^{N}\right)$ defines an instantaneous effective force subject to particle $j$, which also depends on the configuration, given by

[TABLE]

Note that for passive particles, one has $D_{f}\tau_{p}=0$ such that $D_{j}\left(\mathbf{r}^{N}\right)=D_{t}$ and $\mathbf{F}_{j}^{\text{eff}}\left(\mathbf{r}^{N}\right)=\mathbf{F}_{j}\left(\mathbf{r}^{N}\right)$ as expected. In the limit $\tau_{p}\to 0$,corresponding to a white noise $\mathbf{f}_{i}$, we have $D_{j}\left(\mathbf{r}^{N}\right)=D_{t}+D_{f}\tau_{p}^{2}/\gamma^{2}=D_{0}$ and $\mathbf{F}_{j}^{\text{eff}}=\left(D_{t}/D_{0}\right)\mathbf{F}_{j}$. In this latter case, the effective Smoluchowski operator is given byFarage and Brader (2014)

[TABLE]

where $\beta_{\text{eff}}=\left(k_{B}T_{\text{eff}}\right)^{-1}$, and the system reduces to $N$ interacting Brownian particles at the effective temperature $T_{\text{eff}}$.

We assume that the system will reach a nonequilibrium steady state (NESS) $P_{s}\left(\mathbf{r}^{N}\right)$ in the long time limit, which satisfies

[TABLE]

where the steady state current $\mathbf{J}_{i}^{s}$ is given by

[TABLE]

For a passive system, $P_{s}\left(\mathbf{r}^{N}\right)$ will be given by the canonical equilibrium distribution $P_{s}^{\text{eq}}\left(\mathbf{r}^{N}\right)=\exp\left(-\beta U\left(\mathbf{r}^{N}\right)\right)/Z$ where $U\left(\mathbf{r}^{N}\right)=\frac{1}{2}\sum_{j\neq i}u\left(r_{ij}\right)$ is the system potential and $Z$ is the partition function. But for the active system studied in the present work, the explicit form of $P_{s}\left(\mathbf{r}^{N}\right)$ is hard to obtain. Nevertheless, in the case $\tau_{p}\to 0$, $P_{s}^{\tau_{p}\to 0}\left(\mathbf{r}^{N}\right)\sim\exp\left(-\beta_{\text{eff}}U\left(\mathbf{r}^{N}\right)\right)$ satisfies $\hat{\Omega}_{\tau_{p\to 0}}P_{s}^{\tau_{p}\to 0}\left(\mathbf{r}^{N}\right)=0$ indicating that the system can be described by an effective equilibrium distribution at an effective temperature $T_{\text{eff}}$.

For latter purposes, it is convenient to introduce an adjoint operator of the Smoluchowski operator as

[TABLE]

which satisfies $\int d\mathbf{r}^{N}f^{*}\left(\hat{\Omega}g\right)=\int d\mathbf{r}^{N}\left(\hat{\Omega}^{\dagger}f\right)^{*}g$ for any functions $f\left(\mathbf{r}^{N}\right)$ and $g\left(\mathbf{r}^{N}\right)$. For the collective dynamic behaviors of the system, one can then define the collective intermediate scattering function asFarage and Brader (2014)

[TABLE]

where

[TABLE]

is the Fourier transform with wave vector $\mathbf{q}$ of the density variable $\rho\left(\mathbf{r}\right)=\sum_{j=1}^{N}\delta\left(\mathbf{r}-\mathbf{r}_{j}\right)$ and $q=\left|\mathbf{q}\right|$. In particular, one must emphasize that the brackets $\left\langle\right\rangle$ in Eq.(18) denotes the ensemble average over the NESS distribution $P_{s}\left(\mathbf{r}^{N}\right)$, rather than over the equilibrium one $P_{s}^{\text{eq}}\left(\mathbf{r}^{N}\right)$. For $t=0$, $F_{q}\left(t\right)$ is related to the non-equilibrium static structure factor

[TABLE]

where again $\left\langle\cdots\right\rangle$ denotes averaging over the NESS. Nevertheless, for the non-equilibrium system studied here, $S\left(q\right)$ can not be calculated by analytical methods like the Ornstein-Zernike (OZ) equations and must be obtained by direct simulations.

Note that $F_{q}\left(t\right)$ can also be written as

[TABLE]

wherein the operator $\hat{\Omega}$ acts on all the functions on its right side including $P_{s}\left(\mathbf{r}^{N}\right)$, while in Eq.(18) the adjoint operator $\hat{\Omega}^{\dagger}$ only acts on $\rho_{\mathbf{q}}$. We also consider a closely related function, $F_{q}^{s}\left(t\right)$, called self-intermediate scattering function

[TABLE]

where $\rho_{q}^{s}$ is the Fourier transform of microscopic tagged particle (tracer) density $\rho_{q}^{s}=e^{-i\mathbf{q}\cdot\mathbf{r}_{s}}$.

II.3 Memory Function Equations

In this subsection, we derive a formal expression for the collective (and self-) intermediate scattering functions Eqs. (21) and (II.2) in terms of the so-called irreducible memory function. This can be done most easily in the Laplace domain, and the details are given in the Appendix B. Consequently, the equation for the time evolution of $F_{q}\left(t\right)$ is given by

[TABLE]

where

[TABLE]

is the frequency term where $D_{j}\left(\mathbf{r}^{N}\right)$ is given by Eq.(13) and

[TABLE]

denotes an averaged single-particle diffusion coefficient in the NESS.

The irreducible memory function $\tilde{M}^{{\rm irr}}\left(q,t\right)$ is given by

[TABLE]

with $\mathbf{p}=\mathbf{q}-\mathbf{k}$, $p=\left|\mathbf{p}\right|$. Herein, a pseudo-correlation function $C_{2}\left(\mathbf{q},\mathbf{k}\right)$ is introduced which is defined as

[TABLE]

where

[TABLE]

denotes a static structure function involving the coupling of the instantaneous diffusion coefficient $D_{j}\left(\mathbf{r}^{N}\right)$ and density fluctuation $e^{-i\mathbf{k}\cdot\left(\mathbf{r}_{j}-\mathbf{r}_{i}\right)}$.

Eqs. (24) to (29) contribute to the main theoretical results of the present paper. The equation for $F_{q}\left(t\right)$, (24), has the same form as that for an equilibrium colloid systemNägele (1996). However, Eqs.(25) to (29) contain important new features that are specific to the AOU-T system. The frequency $\omega_{q}$ depends on the parameter $\bar{D}$, which denotes an averaged effective diffusion coefficient of a particle in the NESS. The irreducible memory function, Eq.(27), has similar form as that for passive colloid systems, except that a new pseudo-direct correlation function $C_{2}\left(\mathbf{q};\mathbf{k}\right)$ is introduced in replace of the usual direct correlation function $c(k)=\rho^{-1}\left[1-S^{-1}\left(k\right)\right]$. The definition of $C_{2}\left(\mathbf{q};\mathbf{k}\right)$ now involves another function $S_{2}\left(k\right)$, defined by Eq.(29), which resembles the structure factor $S\left(k\right)$ but with $D_{j}\left(\mathbf{r}^{N}\right)$ involved. Since $D_{j}$ is a configuration-dependent function, it cannot be drawn out of the summation $\sum_{i,j}$ in Eq.(29). Interestingly, for a homogeneous passive system, $D_{j}=D_{t}=D_{0}$, hence $\omega_{q}=q^{2}D_{t}S^{-1}\left(q\right)$, $S_{2}\left(k\right)=N^{-1}\left\langle\sum_{i,j}e^{-i\mathbf{k}\cdot\mathbf{r}_{j}+i\mathbf{k}\cdot\mathbf{r}_{i}}\right\rangle\equiv S\left(k\right)$ and $C_{2}\left(\mathbf{q};\mathbf{k}\right)$ simply reduces to $c\left(k\right)$. In this case, Eq. (27) becomes

[TABLE]

which reduces exactly to that of a passive colloid systemCichocki and Hess (1987); Kawasaki (1994). Note that in the limit $\tau_{p}\to 0$, we have $D_{j}=\overline{D}=D_{0}$ such that $S_{2}\left(k\right)=S\left(k\right)$ and $C_{2}\left(\mathbf{q};\mathbf{k}\right)=c\left(k\right)$ also hold. In this case, $\omega_{q}=q^{2}D_{0}S^{-1}\left(q\right)$ and the equations describe the dynamics of an equivalent Brownian system with effective diffusion coefficient $D_{0}$ as described above.

In general cases, $D_{j}$ is dependent on the particle positions, thus it cannot be drawn out from the summation of $S_{2}\left(k\right)$ in Eq.(29). Interestingly, if we approximately replaces $D_{j}$ by its ensemble average value $\left\langle D_{j}\right\rangle$ in the summation, we can obtain

[TABLE]

in the second equality, we use the fact that $\left\langle D_{j}\right\rangle=\bar{D}$ for a homogenous system. We find then

[TABLE]

and the memory function Eq.(27) reduces to that for a passive system with effective diffusion coefficient $\bar{D}$. Note that this approximation holds if the coupling of $D_{j}$ and density fluctuation $e^{-i{\bf k}\cdot\left({\bf r}_{j}-\mathbf{r}_{l}\right)}$ is weak or the fluctuation of $D_{j}$ is very small. In the latter sections of the present paper, we will show by simulations that $S_{2}\left(k\right)\simeq\left(\bar{D}/D_{0}\right)S\left(k\right)$ is a very good approximation when $k$ is large, whereas for small $k$ $S_{2}\left(k\right)D_{0}/\bar{D}S\left(k\right)$ does show apparent structures.

In the next section, we will adopt our above theoretical results to study the glassy behaviors of the one-component LJ active system described by the AOU-T model. In the dimensionless unit, $\gamma=D_{t}=k_{B}T=1$ and we choose $D_{f}$, $\tau_{p}$ together with the number density $\rho$ as adjustable parameters. As already discussed in the model description part, now the effective temperature is given by $T_{\text{eff}}=1+D_{f}\tau_{p}^{2}$, while the amplitude of active force is quantified by $v_{0}=\sqrt{3D_{f}\tau_{p}}$. To compare our results with those simulation works of Ni and others, we will choose $v_{0}$ and $\tau_{p}$ as independent variables together with $\rho$. Nevertheless, we will also study the behavior of the system by choosing $T_{\text{eff}}$ and $\tau_{p}$ as independent free parameters since it has been a regular choice in recent studiesSzamel et al. (2015); Flenner et al. (2016). To begin, we will run the system until it reaches the steady state from which we can get the parameter $\bar{D}$ and the function $S_{2}\left(k\right)$, with which the memory function Eq. (24) can be numerically calculated. We can then investigate the time dependence of $F_{q}\left(t\right)$ to address the glass transition issue.

For the self-intermediate scattering function $F_{q}^{s}\left(t\right)$, the memory function equation reads(see the Appendix B)

[TABLE]

where $\omega_{q}^{s}=q^{2}\bar{D}$ and

[TABLE]

If $S_{2}\left(k\right)\simeq\left(\bar{D}/D_{0}\right)S\left(k\right)$, the second term in the bracket can be neglected, and the equation reduces to the equilibrium version.

It would be instructive here to compare our theoretical results with those in the literature. As mentioned in the introduction, Farage and BraderFarage and Brader (2014) had tried to develop a MCT for the RAB model in the limit $\tau_{p}\to 0$. In this circumstance, the effective Smoluchowski operator is actually given by Eq.(15). Starting from this effective operator, they obtained a memory function for the collective scattering function $F_{q}^{\text{eq}}\left(t\right)$, but defined for the equilibrium distribution. In our work, the effective Smoluchowski operator is now extended to finite (small) $\tau_{p}$, and importantly, the scattering function is now defined for the NESS which is more relevant for the active system as pointed out by SzamelSzamel et al. (2015). The extension to finite $\tau_{p}$ and using a nonequilibrium function makes it feasible to compare with simulation results. Surely, for a nonequilibrium MCT theory, some static functions such as $\bar{D}$, $S_{2}\left(k\right)$ in the present work must be obtained from simulations, which is currently not avoidable.

On the other hand, Szamel et al.Szamel (2016) had made important progress in the theoretical modeling of active particle systems very recently. In particular, they mainly focused on athermal system, the AOU model, which is applicable for large colloidal systems wherein thermal noise might be ignored compared to the self-propulsion. Their treatment followed a quite different way as in the present work, where they performed a projection onto the local steady state defined by the self-propulsion force $\mathbf{f}_{i}$. With the assumption of vanishment of system currents in the local steady state and mode coupling approximation, they were able to obtain an effective Smoluchowski operator, which is time dependent, and the memory function for the nonequilibrium scattering function $F_{q}\left(t\right)$. Importantly, their theory involved a function $\omega_{||}\left(q\right)$ which highlights the role of the velocity correlations. In particular, this theory reproduced a nontrivial non-monotonic dependence of the relaxation time $\tau_{\alpha}$ with $\tau_{p}$ if $T_{\text{eff}}$ is fixed which was observed in their simulations for a standard LJ system, although the theory apparently overestimated $\tau_{\alpha}$ in the $\tau_{p}\to 0$ limit. In our present work, we have considered the AOU-T model where thermal noise is taken into account. We have not tried to extend Szamel’s method to this thermal situation, which might be hard to realize, rather we have adopted a different scheme. Given that the Fox’s method is applicable, the effective Smoluchowski operator given by Eq.(12) provides the starting point for the derivation. This approach actually involves a type of coarse-gaining over time, wherein the effects of colored noise has been replaced by an effective white one but with configuration-dependent correlation functions. As shown in our theory above, the dynamics is then mainly determined by the effective diffusion coefficient $\bar{D}$ and a static structure function $S_{2}\left(k\right)$ wherein both involves the instantaneous diffusion coefficient $D_{j}\left(\mathbf{r}^{N}\right)$. Interestingly, although our method are quite different with that of Szamel, we note that $\omega_{||}\left(q\right)\tau_{p}$ in their work plays the same role as $\bar{D}$ in ours. We also note that in a recent paper, Marconi et al.Maggi et al. (2015) had studied the velocity correlations in the AOU-model, finding an expression very similar to $D_{j}\left(\mathbf{r}^{N}\right)$ under mean-field approximation.

III Numerical Results and Discussions

III.1 Static Properties

As discussed above, to solve $F_{q}\left(t\right)$, we must obtain the parameter $\bar{D}$ and the pseudo-structure factor $S_{2}\left(k\right)$ in the NESS $via$ direct numerical simulations. In Fig.1(a), the dependence of $\bar{D}$ on the effective temperature $T_{\text{eff}}$ is presented, for different fixed values of $\tau_{p}$ and $\rho$. As can be seen, $\bar{D}$ increases monotonically with $T_{\text{eff}}$, which is reasonable since $\bar{D}$ denotes a kind of averaged diffusion coefficient that should be larger for a higher temperature. If $T_{\text{eff}}$ is fixed, $\bar{D}$ decreases with increasing $\tau_{p}$ and the variation of $\bar{D}$ with $T_{\text{eff}}$ becomes less sharp, i.e., $\left(\partial\bar{D}/\partial T_{\text{eff}}\right)$ decreases. Such qualitative behaviors are robust with the change of number density $\rho$, despite that the value of $\bar{D}$ decreases slightly with increasing $\rho$ for given values of $\tau_{p}$ and $T_{\text{eff}}$. In Fig.1(b), we have also plotted $\bar{D}$ as a function of persistent time $\tau_{p}$ for different values of $T_{{\rm eff}}$ and $\rho$. In this case, we can see that $\bar{D}$ decreases with $\tau_{p}$, tending to approach $D_{0}=k_{B}T_{{\rm eff}}/\gamma$ when $\tau_{p}\rightarrow 0$ and close to $D_{t}$ at a large $\tau_{p}$ value. Besides, at the lower density $\bar{D}$ has a slightly larger value as shown in (a).

In Fig.2(a), the non-equilibrium static structure factors $S\left(k\right)$ obtained from direct simulations are drawn for different particle activities $v_{0}$ for fixed $\rho=1.12$ and $\tau_{p}=0.167$. The value of $\rho$ is chosen such that the system is close to the glass transition and that of $\tau_{p}$ is consistent with the setting in the simulation work of NiNi et al. (2013). It seems that $S\left(k\right)$ does not change much with the variation of $v_{0}$, except that the main peak decreases slightly and shifts a little to right with increasing $v_{0}$. This decreasing of the main peak indicates that the structure becomes looser with increasing active force. The other peaks at larger values of $k$ show little discrepancy for different $v_{0}$. Such observations are in qualitative agreements with the simulation results obtained by Ni. Since we have fixed $\tau_{p}$, the effective temperature $T_{\text{eff}}\sim 1+D_{f}\tau_{p}^{2}$ changes in the same tendency as $v_{0}$, such that Fig.2(a) also shows the change of $S\left(k\right)$ with $T_{\text{eff}}$. In Fig.2(b), $S\left(k\right)$ for different $\tau_{p}$, but with fixed $T_{\text{eff}}$ have been presented. In this case, one can see that the main peak increases apparently with increasing $\tau_{p}$ and also shifts a little bit to smaller values of $k$. Since $T_{\text{eff}}\sim 1+v_{0}^{2}\tau_{p}/3$, increasing $\tau_{p}$ with fixed $T_{\text{eff}}$ corresponds to decreasing $v_{0}$, this observation is consistent with Fig.2(a). The second and third peak also show observable differences with the variation of $\tau_{p}$ in that the peak gets higher and moves to smaller values of k with increasing $\tau_{p}$.

As discussed in the last section, an important new feature of our theory is the introduction of the function $S_{2}\left(k\right)$, which couples the instantaneous diffusion coefficient $D_{j}$ and the density fluctuations. It is therefore instructive for us to investigate how $S_{2}\left(k\right)$ looks like. In Fig.3, we have plotted $S_{2}\left(k\right)$ with the same parameter settings as in Fig.2. As shown in Fig.3(a), the particle activity (or effective temperature) drastically influences $S_{2}\left(k\right)$, with the main peak decreasing considerably with increasing $v_{0}$ or $T_{\text{eff}}$. Compared to Fig.2(a), the value of $S_{2}\left(k\right)$ is much smaller than $S\left(k\right)$, which reflects the fact that $D_{j}$ is generally less than $D_{0}$. The behaviors of $S_{2}\left(k\right)$ for fixed $T_{\text{eff}}$ but with variant $\tau_{p}$ are shown in Fig.3(b). In this latter case, we see that the main peak now slightly reduces with increasing $\tau_{p}$ and it seems to saturate for large $\tau_{p}$, which are at variance with the observations in Fig.2(b). The apparent discrepancies between $S_{2}\left(k\right)$ and $S\left(k\right)$ indicate that our theory may show interesting new features.

Another new feature of our theory is the pseudo-correlation function $C_{2}\left(\mathbf{q},\mathbf{k}\right)$, which plays a similar role to $c\left(k\right)$ in the irreducible memory function $M^{irr}\left(q,t\right)$. As discussed above, $C_{2}\left(\mathbf{q},\mathbf{k}\right)$ reduces to $c\left(k\right)$ if $S_{2}\left(p\right)D_{0}/\bar{D}S\left(p\right)\simeq 1$. In Fig.4, the dependence of $S_{2}\left(k\right)D_{0}/\bar{D}S\left(k\right)$ on $k$ has been presented, for fixed $\rho$ with varying $T_{{\rm eff}}$ and $\tau_{p}$. Interestingly, we find that it is approximately one if $k$ is larger than $2\pi/\sigma$ which is approximately the peak position for $S\left(k\right)$. Nevertheless, for small values of $k$, $S_{2}\left(k\right)D_{0}/\bar{D}S\left(k\right)$ shows some structures. Specifically, $S_{2}\left(k\right)D_{0}/\bar{D}S\left(k\right)$ becomes much less than one for small $\tau_{p}$ if $T_{\text{eff}}$ is fixed. Such a feature may lead to enhancement of the irreducible memory function $M_{s}^{irr}\left(q,t\right)$ shown in Eq.(34) with decreasing $\tau_{p}\to 0$ if $T_{\text{eff}}$ is fixed. This would lead to the increment of $\tau_{\alpha}$,if other effects are not accounted for. Note that, however, $\bar{D}$ increases with decreasing $\tau_{p}$ with constant $T_{\text{eff}}$ as shown in Fig.1, such that $\tau_{\alpha}$ would decrease with decreasing $\tau_{p}$ in terms of this effect. Therefore, it might be possible that the relaxation time $\tau_{\alpha}$ shows some re-entrance behaviors in the small $\tau_{p}$ region, similar to that reported for the AOU modelBerthier (2014); Szamel et al. (2015).

III.2 Intermediate Scattering Function

With the static properties obtained above, particularly $\bar{D}$ and $S_{2}\left(k\right)$, we are ready to investigate the behavior of the intermediate scattering function $F_{q}\left(t\right)$ by numerically solving the memory functions Eqs. (24) and (33). In Fig.5(a), the normalized scattering functions $\phi_{q}\left(t\right)=F_{q}\left(t\right)/S\left(q\right)$ are shown for different values of $v_{0}$ (or $T_{\text{eff}}$) and number density $\rho$, wherein we have chosen $q=7.5$ which is around the first peak of $S\left(q\right)$. The results for two densities $\rho=1.05$ and 1.10 are plotted, and the value of $\tau_{p}$ is fixed to be 0.167. For the higher density $\rho=1.07,$ one can see that $F_{q}\left(t\right)$ finally reaches a plateau in the long time limit for $v_{0}=0$ (or $T_{\text{eff}}=1$), indicating that the system reaches the glassy state. For a nonzero value of $v_{0}$ as shown in the figure, $F_{q}\left(t\right)$ will finally relax to zero for large $t$ indicating that the system is in a liquid state, and the relaxation time decreases apparently with increasing $v_{0}$. Therefore, activity will push the glass transition to higher number density, in consistent with the simulation results of the RAB model and other related models. For a smaller $\rho=1.05$, the system is in the liquid state for $v_{0}=0$ and the relaxation of $F_{q}\left(t\right)$ also becomes faster with increasing $v_{0}$ or $T_{\text{eff}}$. The behaviors of the self-scattering function $F_{q}^{s}\left(t\right)$ are similar as shown in Fig.5(b). While for $v_{0}=0$ the tracer particle is trapped and $F_{q}^{s}\left(t\right)$ reaches a non-zero value for $t\to\infty$, it relaxes to zero for $v_{0}=10$ and $20$ with the relaxation time $\tau_{\alpha}$ decreases apparently with increasing $v_{0}$.

The limiting value $f_{q}=\lim_{t\to\infty}\phi_{q}\left(t\right)$ at the plateau defines the so-called Debye-Waller factor. A non-zero value of $f_{q}$ indicates that the system is in the glassy state. With the increase of $\rho$, $f_{q}$ may change from zero to an apparent nonzero value, and the value $\rho_{c}$ thus corresponds to the glass transition point. One may also fix $\rho$ but vary $T_{\text{eff}}$, then $f_{q}$ may become nonzero for $T_{\text{eff}}$ less than some critical value $T_{\text{eff}}^{c}$ , which defines a critical temperature for glass transition. According to the MCT Eq. (24), the Debye-Waller factor $f_{q}$ follows

[TABLE]

where

[TABLE]

This equation can be solved numerically and self-consistently to get $f_{q}$ for given control parameters $v_{0}$ (or $T_{\text{eff}}$), $\tau_{p}$ as well as $\rho$.

In Fig.6(a), the dependence of $f_{q}$ on the number density $\rho$ is presented for different $v_{0}$ (or $T_{\text{eff}}$) with given $\tau_{p}=0.001$. Clearly, $f_{q}$ changes abruptly from zero to a large nonzero value at a critical density $\rho_{c}$, indicated for example by the vertical dashed line for $v_{0}=0$ at about $\rho_{c}\simeq 1.064$. For given $\tau_{p}$, the curve shifts to larger values of $\rho$ with increasing $v_{0}$, indicating that that glass transition is pushed to higher values of $\rho$ for larger particle activity in consistent with previously reported simulation results. The pictures for different values of $\tau_{p}=0.05$ and 0.167 are shown in Fig.6(b) and (c), respectively. The results are similar to those in (a), with the values of $\rho_{c}$ shifting to relatively larger values for larger $\tau_{p}$.

In Fig.5(b), it is shown that the relaxation time $\tau_{\alpha}$ increases when the system approaches the glass transition and it diverges at the glass transition point. Therefore, one may also study the glass transition by investigating the behavior of $\tau_{\alpha}$ as a function of $\rho$. The results are depicted in Fig.7 with the same parameter setting as in Fig.6. Obviously, $\tau_{\alpha}$ increases fastly with $\rho$ for fixed values of $v_{0}$ ($T_{\text{eff}}$) and $\tau_{p}$ and it diverges at some critical value $\rho_{c}$. For a very small $\tau_{p}=0.001$, it seems that changing $v_{0}$ does not affect very much the values of $\tau_{\alpha}$ as shown in Fig.7(a). The influence becomes more considerable when $\tau_{p}$ gets larger as demonstrated in 7(b) and (c), and the value of $\rho_{c}$ also shifts to larger values in consistent with Fig.6.

Surely, the value of $\rho_{c}$ should be the same either obtained by $f_{q}$ or $\tau_{\alpha}$ within reasonable fluctuations. In Fig.8(a) and (b), the dependence of $\rho_{c}$, obtained from both$f_{q}$ and $\tau_{\alpha}$, on $v_{0}$ and $T_{\text{eff }}$ are presented for different given values of $\tau_{p}$. Clearly, $\rho_{c}$ increases with both $v_{0}$ and $T_{\text{eff}}$ as expected. Interestingly, $\rho_{c}$ shows a nearly linear dependence on $v_{0}^{2}$, wherein the slope increases with $\tau_{p}$. This linear dependence was also observed in the simulation work of Ni. We also note that $\rho_{c}$ increases with $\tau_{p}$ for fixed $v_{0}$, whereas it decreases with $\tau_{p}$ for fixed $T_{\text{eff}}$ according to the data presented in Fig.8(a) and (b). This is shown more clearly in Fig.8(c) , where we have also plotted $\rho_{c}$ as a function of $\tau_{p}$ for different $T_{\text{eff}}$. For comparison, the dashed line gives the value of $\rho_{c}^{B}$ for the corresponding passive Brownian system with $T=T_{\text{eff}}$, which is obtained by setting $D_{t}=D_{0}$ and zero self-propulsion force $f_{i}=0$ in Eq.(1). Clearly, $\rho_{c}$ approaches $\rho_{c}^{B}$ in the limit $\tau_{p}\to 0$ as expected.

For equilibrium systems, glass transition is often studied in terms of the critical temperature $T_{c}$, below which the system enters the glassy state. In the present work, we may study the non-equilibrium glass transition in the same spirit by calculating the value of critical effective temperature $T_{\text{eff}}^{c}$ with fixed number density $\rho$. In Fig.9 (a) to (c), the results of $\tau_{\alpha}$ are presented as functions of $T_{\text{eff}}$ for different $\tau_{p}$ and $\rho$. For fixed value of $\tau_{p}$ and $\rho$, $\tau_{\alpha}$ decreases monotonically with $T_{\text{eff}}$. Below some critical values of $T_{\text{eff}}$ corresponding to $T_{\text{eff}}^{c}$, $\tau_{\alpha}$ diverges indicating the occurrence of glass transition. In Fig.9 (d), the dependence of $T_{\text{eff}}^{c}$ on $\tau_{p}$ for different $\rho$ is shown. One can see that $T_{\text{eff}}^{c}$ increases monotonically with $\tau_{p}$, and it approaches a constant value in the limit $\tau_{p}\to 0$. Such a constant value corresponds to the one for a passive Brownian system with $T_{c}=T_{\text{eff}}^{c}$. We also note that $T_{\text{eff}}^{c}$ increases with the number density $\rho$, indicating that a denser system enters glass transition at a higher critical temperature as expected.

IV Conclusions

In summary, we have developed a promising mode coupling theory to study the nonequilibrium glassy dynamics of a general model system of self-propelled particles. The self-propulsion force is given by a colored noise described by OU process, and thermal noises in the environment are also considered. Our work mainly contains two parts. By using Fox approximation method for Langevin systems with colored noise, an approximate Smoluchowski equation can be obtained, governing the time evolution of the distribution function of the particles’ positions. This effective SE is expected to be exactly valid for not large persistence time $\tau_{p}$ of the propulsion force, and it thus serves as a promising starting point to study the system’s relaxation or glassy dynamics. The SE involves a configuration dependent instantaneous diffusion function $D_{j}\left(\mathbf{r}^{N}\right)$ which is related to the gradient of force subjected to particle $j$. With this SE, we are able to derive a memory function equation for the time dependent behavior of the collective or self- intermediate scattering function $F_{q}\left(t\right)$ or $F_{q}^{s}\left(t\right)$ in the nonequilibrium steady state. Applying the basic assumption that macroscopic currents vanish in the steady state and using standard mode coupling approximation, we have obtained the expressions for the irreducible memory function as well as frequency terms. Particularly, we find that the dynamics are mainly determined by an effective diffusion coefficient $\bar{D}$, which is the ensemble average of $D_{j}\left(\mathbf{r}_{N}\right)$ in the nonequilibrium steady state, and a pseudo steady state structure factor $S_{2}\left(k\right)$, which involves the coupling between $D_{j}\left(\mathbf{r}^{N}\right)$ and density fluctuations. $\bar{D}$ enters the frequency term and thus governs the short time dynamics, whereas both enter the vortex for memory function and influence the long time dynamics. By direct simulations, we find that $\bar{D}$ increases with the single effective temperature $T_{\text{eff}}$ as well as the magnitude $v_{0}$ of propulsion force, while it decreases with $\tau_{p}$ for fixed $T_{\text{eff}}$ or $v_{0}$. The structure function $S_{2}\left(k\right)$ simply decouples into the product of $\bar{D}/D_{0}$ and $S\left(k\right)$, with $S\left(k\right)$ the nonequilibrium static structure factor and $D_{0}$ the single particle diffusion coefficient in the limit $\tau_{p}\to 0$, for relatively large values of $k$, whereas it shows apparent deviations from $\left(\bar{D}/D_{0}\right)\cdot S\left(k\right)$ for small $k$s. Our theory makes it feasible for us to investigate the glassy dynamics of the system, by investigating the time behavior of $F_{q}\left(t\right)$ or $F_{q}^{s}\left(t\right)$ in the long time limit, chosen the persistence time $\tau_{p}$, the effective temperature $T_{\text{eff}}$, as well as the number density $\rho$ as free parameters. We find that the critical density $\rho_{c}$ for glass transition shifts to larger values with increasing $T_{\text{eff}}$ or $v_{0}$ if $\tau_{p}$ is fixed, in good qualitative accordance with the simulation results of active Brownian particles and related systems. In addition, we have also investigated how the critical density $\rho_{c}$ changes with $\tau_{p}$ for a fixed $T_{\text{eff}}$, finding that $\rho_{c}$ decreases with $\tau_{p}$ monotonically and it approaches the value for the corresponding passive Brownian system in the limit $\tau_{p}\to 0$ as expected. We have also calculated the critical temperature $T_{\text{eff}}^{c}$ for glass transition at fixed density, finding that it increases monotonically with $\tau_{p}$ and also approaches the Brownian particle limit for $\tau_{p}\to 0$.

In future work, we would like to extend the present method to more complex systems such as mixtures of self-propelled particles but with different driving forces or to mixtures of active-passive particlesDing et al. (2015). As mentioned in the main text, the relaxation time $\tau_{\alpha}$ may show nontrivial dependence on the persistence time $\tau_{p}$, which would be also an interesting topic to address for system with both propulsion force and thermal noise. In addition, our results demonstrate that only in the limit $\tau_{p}\to 0$, the glass transition point $\rho_{c}$ or $T_{\text{eff}}^{c}$ approaches that of a Brownian system, indicating that the ’collective’ effective temperature with respect to the nonequilibrium glass transition is different from the single particle one Levis and Berthier (2015), which may deserve more detailed study. In a word, we believe that our work presents a useful theoretical framework to study the nonequilibrium dynamics of dense active particles system from the microscopic level which could find many applications in future works.

Acknowledgements.

This work is supported by National Basic Research Program of China(Grant No. 2013CB834606), by National Science Foundation of China (Grant Nos. 21673212, 21521001, 21473165, 21403204), by the Ministry of Science and Technology of China (Grant No.ￜ2016YFA0400904), and by the Fundamental Research Funds for the Central Universities (Grant Nos. WK2060030018, 2030020028,2340000074).

Appendix A Derivation of the Smoluchowski Equation

Generally, for a LE with colored noise, one can get the FPE within Fox approximationFox (1986). For illustration, consider a simple one dimensional over damped LE

[TABLE]

where $G(x)$ denotes the external or internal force and $\chi\left(t\right)$ is the stochastic noise with correlation

[TABLE]

Define a probability distribution function

[TABLE]

where $D\left[\chi\right]$ denotes integration over the noisy path of $\chi\left(t\right)$ and $P\left[\chi\right]$ is the distribution functional of $\chi$ which is assumed to be Gaussian. One can then obtain the FPE governing the evolution of $P\left(y,t\right)$ as follows Farage et al. (2015)

[TABLE]

Note that if $\chi\left(t\right)$ is white noise, $C\left(t-s\right)=D_{0}\delta\left(t-s\right)$, then the second term is just

[TABLE]

which recovers the standard FPE. For a colored noise with

[TABLE]

one can obtain the FPE approximately as

[TABLE]

by assuming

[TABLE]

For a general multi-variable case,

[TABLE]

where

[TABLE]

with $i,j=1,2,\cdots,N$, the FPE for distribution function

[TABLE]

reads

[TABLE]

Then, if $C_{ij}\left(t-s\right)=\delta_{ij}C\left(t-s\right)=\delta_{ij}\frac{D}{\tau_{0}}\exp\left(-\frac{\left|t-s\right|}{\tau_{0}}\right)$, using the assumption mentioned before we can get

[TABLE]

For our system described by Eq.(1), we have correspondingly $x\to\mathbf{r}^{N}$, $G\left(x\right)\to\gamma^{-1}\mathbf{F}\left(\mathbf{r}^{N}\right)=\beta D_{t}\mathbf{F}\left(\mathbf{r}^{N}\right)$, $\gamma^{-1}\mathbf{f}_{i}\left(t\right)\to\mathbf{\chi}_{i}\left(t\right)$. According to Eq.(5), the variable $D$ in Eq. (A.12) is $D_{f}\tau_{p}^{2}/\gamma^{2}$ and $\tau_{0}$ is $\tau_{p}$. Note that the white noise term $\mathbf{\eta}_{i}\left(t\right)$ in Eq.(1) contributes a normal diffusion term to the FPE. Thus we finally obtain

[TABLE]

Write

[TABLE]

and

[TABLE]

is the effective force. Finally, we have

[TABLE]

which corresponds exactly to Eqs. (11) to (14) in the main text.

Appendix B Derivation of the General Langevin Equation

B.1 Memory Function Equation

Here we present the derivation of the memory function equations, namely Eqs. (24) to (29) in the main text, for the scattering function F_{q}\left(t\right)=\frac{1}{N}\left\langle\text{\rho}_{\mathbf{q}}^{*}e^{\hat{\Omega}t}\rho_{\mathbf{q}}\right\rangle. This is most easily done in the Laplace domain, even for the complex Smoluchowski operator $\hat{\Omega}$ shown in Eq.(12) which contains the instantaneous effective diffusion constant $D_{j}\left(\mathbf{r}^{N}\right)$ given by Eq.(13). The main steps are similar to the derivation of MCT equations for passive colloidal systems, following Mori-Zwanzig projection operator procedures.

We start from Laplace transform of the scattering function

[TABLE]

where $\mathcal{LT}$ stands for Laplace transformation and $A_{\mathbf{q}}=\rho_{\mathbf{q}}/\sqrt{N}$. One can define a projection operator on the density

[TABLE]

which has the property

[TABLE]

Accordingly, we can define $\mathcal{Q}=\mathcal{I}-\mathcal{P}$, which satisfies $\mathcal{Q}A_{\mathbf{q}}=0$, $\mathcal{PQ}=0$, and $\mathcal{Q}^{n}=\mathcal{Q}$. Then for the operator $\left[z-\Omega\right]^{-1}$, one has the following identity

[TABLE]

which is known as Dyson decomposition and can be easily checked by right-multiplying both sides by $z-\hat{\Omega}$.

wherein the operator $\Omega$ acts on all the functions to its right side. In Laplace domain, this reads

[TABLE]

Using the definition of $\mathcal{P}$, the first term is

[TABLE]

While the second term is, using the identity(B.3),

[TABLE]

Note that $\mathcal{Q}A_{\mathbf{q}}=0$, hence \frac{1}{z-\hat{\mbox{\Omega}}\mathcal{Q}}A_{\mathbf{q}}=0 and

[TABLE]

where we have used the definition of $\mathcal{P}$ in the second equality and the fact $\mathcal{Q}\mathcal{Q}=\mathcal{Q}$ in the third equality. We may introduce

[TABLE]

which is Eq.(25) and define

[TABLE]

Then the time-evolution of $F\left(q,t\right)$ in Laplace domain reads

[TABLE]

Therefore

[TABLE]

For colloidal systems, there is a so-called irreducible issue Szamel (2016); Kawasaki (1995) . Following the procedure in Nägele (1996) one needs to introduce an irreducible memory function $\tilde{M}^{{\rm irr}}\left(q,z\right)$, which is related to$\tilde{M}\left(q,z\right)$ according to

[TABLE]

and

[TABLE]

Herein, $\hat{\Omega}^{\text{irr}}$ denotes an irreducible Smoluchowski operator of which the detailed form is not relevant within the MCT approximation below. Consequently, this leads to

[TABLE]

corresponding to in the time domain

[TABLE]

which is exactly Eq.(24) in the main text. Note that above derivations are quite general and do not depend on the explicit form of the operator $\hat{\Omega}$, whereas the expressions for $\omega_{q}$ and $M^{\text{irr}}$ should certainly depend on $\hat{\Omega}$.

B.2 Frequency $\omega_{q}$

We now substitute $A_{{\bf q}}=\rho_{\mathbf{q}}/\sqrt{N}=\sum_{j}e^{-i{\bf q}\cdot{\bf r}_{j}}/\sqrt{N}$ to calculate $\omega_{q}$. Note that

[TABLE]

where as mentioned before, the operator $\hat{\Omega}$ acts on all the functions to its right side including the steady-state distribution function $P_{s}\left(\mathbf{r}^{N}\right)$. In the steady state, the summation of all the currents $\mathbf{J}_{j}^{s}$, given by Eq.(17), is zero according to $\hat{\Omega}P_{s}\left(\mathbf{r}^{N}\right)=-\sum_{j}\mathbf{J}_{j}^{s}=0$. To proceed and as many authors have done, we assume more strongly that $\mathbf{J}_{j}^{s}=0$, i.e.,

[TABLE]

Therefore, one can obtain that

[TABLE]

which is the counterpart of Yvon theoremHansen and McDonald (1990) in this nonequilibrium system. Using this result, one has

[TABLE]

where the second equality results from partial integration and we have used $\nabla_{j}A_{\mathbf{q}}=-i\mathbf{q}\exp\left(-i\mathbf{q}\cdot\mathbf{r}_{j}\right)/\sqrt{N}$ in the third one. $\left\langle D_{j}\right\rangle=\int d\mathbf{r}^{N}D_{j}\left(\mathbf{r}^{N}\right)P_{s}\left(\mathbf{r}^{N}\right)$ denotes the averaged instantaneous diffusion function of particle $j$ and $\bar{D}=N^{-1}\sum_{j}\left\langle D_{j}\right\rangle$. Therefore, the effective frequency $\omega_{q}$ reads

[TABLE]

which is Eq.(25) in the main text.

B.3 Memory Function $M^{\text{irr}}\left(q,t\right)$

In the time domain, the irreducible memory function $M^{\text{irr}}\left(q,t\right)$ is given by

[TABLE]

Using the adjoint operator $\hat{\Omega}^{\dagger}$, the first term is

[TABLE]

where $\mathcal{QQ=Q}$ is used in the third equality and we have introduced

[TABLE]

which is a type of “random force”.

It is in this step one needs to introduce the mode-coupling approximation. The memory function is assumed to be dominated by the projection onto the coupling density modes Götze (2008). One can then define a second-order projection operator

[TABLE]

and make the approximation:

[TABLE]

Here we have to factories the static and dynamic four-point correlation functions into products of two-point functions

[TABLE]

and simultaneously replace the projected operator $\mathcal{Q}\hat{\Omega}^{{\rm irr}}\mathcal{Q}$ by the full Smoluchowski operator $\hat{\Omega}$ in the propagator governing the time evolution of the correlation function Hansen and McDonald (1990)

[TABLE]

We now need to calculate $\left\langle\rho_{\mathbf{p}}^{*}\rho_{\mathbf{k}}^{*}R_{\mathbf{q}}\right\rangle$, which is given by

[TABLE]

The first term in the bracket is

[TABLE]

where the second equality is simply a result of partial integration, and the fourth equality results from translational invariance. For short of notation, we have introduced a function $S_{2}\left(k\right)$ in the fourth equality defined as

[TABLE]

in accordance with Eq.(29) in the main text. If $D_{j}$ is a constant, such as$D_{j}=D_{0}$ in the $\tau_{p}\rightarrow 0$ limit, it can be drawn out of the bracket such that $S_{2}\left(k\right)=N^{-1}\left\langle\sum_{j,l}e^{-i\mathbf{k}\cdot\left(\mathbf{r}_{j}-\mathbf{r}_{l}\right)}\right\rangle=S\left(k\right)$ which is exactly the static structure factor. Nevertheless, in our present case, $D_{j}$ depends on the instantaneous configuration $\mathbf{r}^{N}$, such that $S_{2}\left(k\right)$ may show abundant features different from $S\left(k\right)$. It is interesting to note that for large $k$, $D_{j}$ seems to be decoupled from the Fourier components $\exp\left(-i\mathbf{k}\cdot\mathbf{r}_{l}\right)$, and $S_{2}\left(k\right)$ can be approximated by

[TABLE]

as shown in Fig.3 in the main text.

The second term in Eq.(B.13) can be calculated using the so-called convolution approximation, which is assumed to be still appropriate in nonequilibrium situationSzamel (2014),

[TABLE]

Therefore, we can get

[TABLE]

Substituting this into Eq.(B.12), we obtain

[TABLE]

where the vortex function $V_{\mathbf{q}}\left(\mathbf{k},\mathbf{p}\right)$ is defined as

[TABLE]

Now it is instructive for us to define a pseudo “direct correlation function” as

[TABLE]

Note that if Eq.(B.14) becomes a equality, namely, $S_{2}(k)=\bar{D}S(k)/D_{0}$, then

[TABLE]

which has the same form as the conventional direct correlation function. With this notation, the vortex can be written as

[TABLE]

and

[TABLE]

Consequently, we get the irreducible memory function

[TABLE]

Changing to integration by using $\sum_{\mathbf{k}}\to\left(2\pi\right)^{-3}V\int d^{3}\mathbf{k}$, one has

[TABLE]

where $\mathbf{\hat{q}}=\mathbf{q}/q$ is the unit vector in the direction of $\mathbf{q}$. This is Eq. (27) in the main text.

B.4 Tagged Particle Dynamics

We now consider tagged particle dynamics. The relevant variable is $A_{\mathbf{q}}=\rho_{\mathbf{q}}^{s}=e^{-i\mathbf{q}\cdot\mathbf{r}_{s}}$ ( the subscript or superscript ’s’ stands for the single particle) and the self-scattering function reads $F_{q}^{s}\left(t\right)=\left\langle\rho_{-\mathbf{q}}^{s}e^{\hat{\Omega}t}\rho_{\mathbf{q}}^{s}\right\rangle$ with $F_{q}^{s}\left(0\right)=1$. The derivation of the memory function equation for $F_{q}^{s}\left(t\right)$ are similar to that of $F_{q}\left(t\right)$, except that some relevant calculations are different. Briefly, we have

[TABLE]

where

[TABLE]

In the third equality of above equation for $R_{\mathbf{q}}^{s}$, we have used the result $\left\langle\rho_{-\mathbf{q}}^{s}\hat{\Omega}\rho_{\mathbf{q}}^{s}\right\rangle=-q^{2}\left\langle D_{s}\right\rangle$.

To calculate $\left\langle R_{\mathbf{q}}^{s*}e^{\mathcal{Q}\hat{\Omega}^{{\rm irr}}\mathcal{Q}t}R_{\mathbf{q}}^{s}\right\rangle$, one projects $R_{\mathbf{q}}^{s}$ onto the product of single and collective modes $\rho_{\mathbf{k}}\rho_{\mathbf{p}}^{s}$ with projection operator

[TABLE]

Then

[TABLE]

where

[TABLE]

Now we need to calculate $\left\langle\left(\rho_{\mathbf{k}}\rho_{\mathbf{p}}^{s}\right)^{*}R_{\mathbf{q}}^{s}\right\rangle$, which is given by

[TABLE]

The second term is

[TABLE]

where the $\delta$ symbol results from translational invariance and the result $\left\langle\rho_{\mathbf{k}}^{*}\rho_{\mathbf{k}}^{s}\right\rangle=\rho c\left(k\right)S\left(k\right)=S\left(k\right)-1$Götze (2008).

The first term is

[TABLE]

where we have used partial integration and Yvon theorem in the second equality.

[TABLE]

[TABLE]

Notice if $S_{2}(k)=\bar{D}S(k)/D_{0}$, $V_{q}\left(\mathbf{k},\mathbf{p}\right)=\delta_{\mathbf{q},\mathbf{k}+\mathbf{p}}\frac{\rho\bar{D}}{N}\left[\mathbf{k}\cdot\mathbf{q}c(k)\right]$, which is the equilibrium result. Next

[TABLE]

and

[TABLE]

Finally, the memory function equation for the self-scattering function $F_{q}^{s}\left(t\right)$ is given by

[TABLE]

where

[TABLE]