Mode Coupling Theory for Nonequilibrium Glassy Dynamics of Thermal Self-Propelled Particles
Mengkai Feng Zhonghuai Hou






Corresponding Author: hzhlj@ustc.edu.cn

Department of Chemical Physics & Hefei National Laboratory for Physical
Sciences at Microscales, iChEM, University of Science and Technology
of China, Hefei, Anhui 230026, China



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,
wherein the nonequilibrium character of the active system is manifested
through an averaged diffusion coefficient D̅ and a nontrivial
structural function S_2(q) with q the magnitude
of wave vector 𝐪. D̅ and S_2(q)
enter the frequency term and the vortex term for the memory function,
thus influence both the short time and the long time dynamics of the
system. 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 τ_α
as functions of the persistence time τ_p of self-propulsion,
the single particle effective temperature T_eff as well
as the number density ρ. Consequently, we find the critical
density ρ_c for given τ_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_eff^c for fixed ρ as well as its dependence
on τ_p. We find that T_eff^c increases with
τ_p and in the limit τ_p→0, 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. 

Mode Coupling Theory for Nonequilibrium Glassy Dynamics of Thermal
Self-Propelled Particles
    Mengkai Feng,  Zhonghuai Hou
    December 30, 2023
===========================================================================================





§ 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 <cit.>.
A wealth of new nonequilibrium phenomena have been reported, such
as active swarming, large scale vortex formation<cit.>,
phase separation<cit.>,
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 systems<cit.>.
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
increases<cit.>, amorphous solidification
process was found in collective motion of a cellular monolayer<cit.>,
and glassy behaviors could even be found for ant aggregates in large
scales<cit.>, 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. <cit.>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_eff can be introduced to quantify
the strength of the self-propulsion force, and a persistence time
τ_p controls the duration of the persistent self-propelled
motion. Berthier and co-workers<cit.>
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.<cit.> 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 values<cit.>
and particularly, hydrodynamic interactions can further enhance this
effect<cit.>. 

Besides experimental and simulation studies mentioned above, on the
theoretical side, important progresses have also been achieved in
recent years<cit.>. Starting from a generalized Langevin
equation with colored non-thermal noise, Berthier and Kurchan<cit.>
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
model<cit.>, 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 density<cit.>.
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
<cit.>. Very recently, Szamel et
al.<cit.>
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 particles<cit.>. This effective
SE allows us to derive a memory function equation for the nonequilibrium
steady state collective intermediate scattering function F_q(t),
as well as that for the self-intermediate scattering function F_q^s(t).
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 D̅ and a nontrivial steady state structure
function S_2(q), both depending on the effective temperature
T_eff, the persistence time τ_p as well as the
number density ρ. With D̅, S_2(q) and
the nonequilibrium structure factor S(q) as inputs,
we can calculate F_q(t) and F_q^s(t)
for different parameter settings and investigate the glass transition
behaviors. We calculate the critical density ρ_c for glass
transition as a function of the effective temperature T_eff,
the magnitude of propulsion force v_0 and the persistence time
τ_p, by investigating the Debye-Waller factor f_q as
well as the relaxation time τ_α. Consequently, ρ_c
shifts to larger values with increasing propulsion force v_0
or effective temperature T_eff, in good accordance with
previous simulation works. The theory also facilitates us to study
the critical temperature T_eff^c for fixed number density
ρ, as well as its dependence on τ_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. 



§ MODEL AND THEORY





 §.§ 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 γ 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


    𝐫̇_i(t)     =γ^-1[𝐅_i(t)+𝐟_i(t)]+ξ_i(t)

where 𝐫_i denotes the position vector of particle i,
the force 𝐅_i=-∑_j i∇_iu(r_ij)
originates from the interactions where is the pair-potential u(r_ij),
𝐟_i is the self-propulsion force, and ξ_i(t)
is the thermal noise with zero mean and variance 

    ⟨ξ_i(t)ξ_j(t')⟩ =2D_t1δ_ijδ(t-t'),

where D_t=k_BT/γ with k_B the Boltzmann constant
and T the ambient temperature, and 1 denotes the unit tensor.
The equations of motion for the self-propulsion force 𝐟_i
are given by


    𝐟̇_̇i̇(t)=-τ_p^-1𝐟_i(t)+η_i(t)

where τ_p is the persistence time of self-propulsion and η_i(t)
is a Gaussian white noise with zero mean and variance D_f


    ⟨η_i(t)η_j(t')⟩ =2D_f1δ_ijδ(t-t')

Accordingly, the correlation function of the force 𝐟_i
reads

    ⟨𝐟_i(t)𝐟_j(t')⟩ =D_fτ_pe^-|t'-t|/τ_p1δ_ij


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

    ⟨δ r^2(t)⟩ =6D_fτ_p^2/γ^2[t+τ_p(e^-t/τ_p-1)]+6D_tt

In the short time limit t≪τ_p, the particle's motion is
diffusive with ⟨δ r^2(t)⟩ =6D_tt,
which is at variance with the AOU system, wherein the term 6D_tt
is absent and the motion is ballistic ⟨δ r^2(t)⟩ =3D_fτ_pγ^-2t^2
for t≪τ_p. For long time, the motion is also diffusive but
with ⟨δ r^2(t)⟩ =6(D_t+D_fτ_p^2/γ^2)t,
implying that the long diffusion coefficient is given by 

    D_0=D_t+D_fτ_p^2/γ^2

This allows us to introduce a single-particle effective temperature

    T_eff=T+k_BD_fτ_p^2/γ=(D_0/D_t)T.

In the limit of vanishing τ_p→0, ⟨𝐟_i(t)𝐟_j(t')⟩→2D_fτ_p^21δ(t-t')δ_ij
and the system becomes equivalent to a Brownian system at the effective
temperature T_eff with diffusion coefficient D_0.
For the AOU model where the thermal noise is absent, the effective
temperature is simply given by T_ eff=D_fτ_p^2/γ<cit.>
. 

We note here that the AOU-T model can be mapped onto the RAB model
at a coarse-grained level<cit.>. The
equation of motion for the RAB particles is ṙ_i=v_0 p_i+γ^-1 F_i+ξ_i,
where v_0 denotes the magnitude of the propulsion force and 𝐩_i
is the unit vector of the direction of particle i. 𝐩_i
changes randomly with time via rotational diffusion, ṗ_i=ζ_i× p_i,
where ζ_i is also a Gaussian white noise with zero mean
and correlation ⟨ζ_i(t)ζ_j(t')⟩ =2D_r1δ_ijδ(t-t')
where D_r denotes the rotational diffusion coefficient. Recently,
it was shown that 𝐩_i can be approximated by a colored
noise with persistence time τ_p=(2D_r)^-1 if
one average over the angular degree of freedom, i.e., 

    ⟨𝐩_i(t)𝐩_j(t')⟩≃1/3e^-2D_r|t-t'| 1δ_ij

Comparing Eq.(<ref>) with (<ref>), we see
that v_0𝐩_i has same correlation property as 𝐟_i
by the mapping τ_p→(2D_r)^-1 and D_fτ_p→ v_0^2/3.

In the studies of the AOU model, the authors usually used T_eff
as one of the independent parameters together with the persistence
time τ_p and number density ρ=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.<cit.>,
they have adopted Stokes-Einstein relation to set the rotational diffusion
coefficient D_r=3D_t/σ^2 where σ is the particle
diameter. In the dimensionless version by setting D_t=1, γ=1
and σ=1, this means that D_r is fixed and the persistence
time is τ_p=(2D_r)^-1=1/6. In more general
cases, the coupling of D_t and D_r may not hold and one
can thus set τ_p as a free independent parameter. In the present
work, we will set v_0=√(3D_fτ_p) or T_eff,
τ_p, and ρ 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

    u(r)=
    4ε[(σ/r)^12-(σ/r)^6]+ε    r≤2^1/6σ
    
    0     r>2^1/6σ

where ε is the strength of the potential. Here we set
ϵ=1k_BT, where k_BT=D_tγ is the unit of energy.
Moreover, σγ/k_BT is the unit of time. The number density
ρ 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. 




 §.§ Effective Smoluchowski Equation


The AOU-T model described in Eq.(<ref>) is non-Markovian
due to the colored noise term 𝐟_i. Consequently, it
is not possible to derive an exact Fokker-Planck equation (FPE) for
the time evolution of the probability distribution function Ψ(𝐫^N,t),
which gives the probability that the system is at a specific configuration
𝐫^N=(𝐫_1,𝐫_2,…,𝐫_N)
at time t. Nevertheless, one may obtain an approximate FPE for
such a colored noise system by applying the method introduced by Fox<cit.>,
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 Ψ(𝐫^N,t),
and analyzed the effective interaction among active particles in the
low density limit. Since the FPE only involves the distribution in
the configuration space 𝐫^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.<cit.>
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 


    ∂/∂ tΨ(𝐫^N,t)=Ω̂Ψ(𝐫^N,t)

where Ω̂ denotes the effective Smoluchowski operator
given by

    Ω̂=∑_j=1^N∇_j· D_j(𝐫^N)[∇_j-β𝐅_j^ eff(𝐫^N)].

Herein, D_j(𝐫^N) is a configuration-dependent
instantaneous diffusion coefficient of particle j which is given
by 

    D_j(𝐫^N)=D_t+D_fτ_p^2/γ^2/1-τ_pβ D_t∇_j·𝐅_j^.

with β=(k_BT)^-1. 𝐅_j^eff(𝐫^N)
defines an instantaneous effective force subject to particle j,
which also depends on the configuration, given by


    𝐅_j^(𝐫^N)=D_t/D_j(𝐫^N)[𝐅_j(𝐫^N)-1/β D_t∇_jD_j(𝐫^N)]


Note that for passive particles, one has D_fτ_p=0 such that
D_j(𝐫^N)=D_t and 𝐅_j^eff(𝐫^N)=𝐅_j(𝐫^N)
as expected. In the limit τ_p→0,corresponding to a white
noise 𝐟_i, we have D_j(𝐫^N)=D_t+D_fτ_p^2/γ^2=D_0
and 𝐅_j^eff=(D_t/D_0)𝐅_j.
In this latter case, the effective Smoluchowski operator is given
by<cit.>

    Ω̂_τ_p→0=D_0∑_j=1^N∇_j(∇_j-βD_t/D_0𝐅_j)=D_0∑_j=1^N∇_j(∇_j-β_eff𝐅_j)

where β_eff=(k_BT_eff)^-1,
and the system reduces to N interacting Brownian particles at the
effective temperature T_eff. 

We assume that the system will reach a nonequilibrium steady state
(NESS) P_s(𝐫^N) in the long time limit,
which satisfies 

    Ω̂P_s(𝐫^N)=-∑_i∇_i·𝐉_i^s=0

where the steady state current 𝐉_i^s is given by 

    𝐉_i^s=-D_j(𝐫^N)[∇_j-β𝐅_j^ eff(𝐫^N)]P_s(𝐫^N)

For a passive system, P_s(𝐫^N) will be
given by the canonical equilibrium distribution P_s^eq(𝐫^N)=exp(-β U(𝐫^N))/Z
where U(𝐫^N)=1/2∑_j iu(r_ij)
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(𝐫^N) is hard to obtain. Nevertheless,
in the case τ_p→0, P_s^τ_p→0(𝐫^N)∼exp(-β_effU(𝐫^N))
satisfies Ω̂_τ_p→0P_s^τ_p→0(𝐫^N)=0
indicating that the system can be described by an effective equilibrium
distribution at an effective temperature T_eff. 

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

    Ω̂^†=∑_j=1^N(∇_j+β𝐅_j^ eff)D_j(𝐫^N)·∇_j

which satisfies ∫ d𝐫^Nf^*(Ω̂g)=∫ d𝐫^N(Ω̂^†f)^*g
for any functions f(𝐫^N) and g(𝐫^N).
For the collective dynamic behaviors of the system, one can then
define the collective intermediate scattering function as<cit.>


    F_q(t)=1/N⟨_𝐪^*(e^Ω̂^†tρ_𝐪)⟩ =1/N⟨_-𝐪(e^Ω̂^†tρ_𝐪)⟩

where 

    ρ_ q=∑_j=1^Ne^-i q· r_j

is the Fourier transform with wave vector 𝐪 of the density
variable ρ(𝐫)=∑_j=1^Nδ(𝐫-𝐫_j)
and q=|𝐪|. In particular, one must emphasize
that the brackets ⟨⟩ in Eq.(<ref>)
denotes the ensemble average over the NESS distribution P_s(𝐫^N),
rather than over the equilibrium one P_s^eq(𝐫^N).
For t=0, F_q(t) is related to the non-equilibrium
static structure factor 

    F_q(0)=1/N⟨_-𝐪ρ_𝐪⟩ =S(q)

where again ⟨⋯⟩ denotes averaging
over the NESS. Nevertheless, for the non-equilibrium system studied
here, S(q) can not be calculated by analytical methods
like the Ornstein-Zernike (OZ) equations and must be obtained by direct
simulations. 

Note that F_q(t) can also be written as 

    F_q(t)=1/N⟨_-𝐪e^Ω̂tρ_𝐪⟩

wherein the operator Ω̂ acts on all the functions on
its right side including P_s(𝐫^N), while
in Eq.(<ref>) the adjoint operator Ω̂^†
only acts on ρ_𝐪. We also consider a closely related
function, F_q^s(t), called self-intermediate scattering
function

    F_q^s(t)     =    ⟨ρ_-q^sρ_q^s(t)⟩ =⟨ e^-i𝐪·(𝐫_s(t)-𝐫_s(0))⟩
        =    1/N∑_j=1^N⟨ e^-i𝐪·(𝐫_j(t)-𝐫_j(0))⟩
        =    1/N∑_j=1^N⟨_-𝐪^je^Ω̂tρ_𝐪^j⟩

where ρ_q^s is the Fourier transform of microscopic tagged
particle (tracer) density ρ_q^s=e^-i𝐪·𝐫_s. 



 §.§ Memory Function Equations


In this subsection, we derive a formal expression for the collective
(and self-) intermediate scattering functions Eqs. (<ref>)
and (<ref>) 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(t) is given by

    ∂/∂ tF_q(t)+ω_qF_q(t)+∫_0^tduM̃^ irr(q,t-u)∂/∂ uF_q(u)=0

where 

    ω_q=-⟨ρ_𝐪^*(^†ρ_𝐪)⟩⟨ρ_𝐪^*ρ_𝐪⟩ ^-1=q^2∑_j⟨ D_j(𝐫^N)⟩/NS(q)=q^2D̅/S(q)

is the frequency term where D_j(𝐫^N) is
given by Eq.(<ref>) and

    D̅=N^-1∑_j⟨ D_j(𝐫^N)⟩ .

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

The irreducible memory function M̃^ irr(q,t)
is given by 

    M̃^ irr(q,t)=ρD̅/16π^3∫ d𝐤[(q̂·𝐤)C_2(𝐪;𝐤)+(q̂·𝐩)C_2(𝐪;𝐤)]^2F_k(t)F_p(t).

with 𝐩=𝐪-𝐤, p=|𝐩|.
Herein, a pseudo-correlation function C_2(𝐪,𝐤)
is introduced which is defined as 

    C_2(𝐪;𝐤)=ρ^-1[1-D_0/D̅S_2(p)/S(p)S^-1(k)]

where

    S_2(k)=1/D_0N⟨∑_i,jD_j(𝐫^N)e^-i𝐤·𝐫_j+i𝐤·𝐫_i⟩

denotes a static structure function involving the coupling of the
instantaneous diffusion coefficient D_j(𝐫^N)
and density fluctuation e^-i𝐤·(𝐫_j-𝐫_i). 

Eqs. (<ref>) to (<ref>) contribute to the main
theoretical results of the present paper. The equation for F_q(t),
(<ref>), has the same form as that for an equilibrium colloid
system<cit.>. However, Eqs.(<ref>)
to (<ref>) contain important new features that are specific
to the AOU-T system. The frequency ω_q depends on the parameter
D̅, which denotes an averaged effective diffusion coefficient
of a particle in the NESS. The irreducible memory function, Eq.(<ref>),
has similar form as that for passive colloid systems, except that
a new pseudo-direct correlation function C_2(𝐪;𝐤)
is introduced in replace of the usual direct correlation function
c(k)=ρ^-1[1-S^-1(k)]. The definition
of C_2(𝐪;𝐤) now involves another
function S_2(k), defined by Eq.(<ref>), which
resembles the structure factor S(k) but with D_j(𝐫^N)
involved. Since D_j is a configuration-dependent function, it
cannot be drawn out of the summation ∑_i,j in Eq.(<ref>).
Interestingly, for a homogeneous passive system, D_j=D_t=D_0,
hence ω_q=q^2D_tS^-1(q), S_2(k)=N^-1⟨∑_i,je^-i𝐤·𝐫_j+i𝐤·𝐫_i⟩≡ S(k)
and C_2(𝐪;𝐤) simply reduces to c(k).
In this case, Eq. (<ref>) becomes 

    M̃^irr(q,t)=ρ D_t/16π^3∫ d𝐤[(q̂·𝐤)c(k)+(q̂·𝐩)c(k)]^2F_k(t)F_p(t)

which reduces exactly to that of a passive colloid system<cit.>.
Note that in the limit τ_p→0, we have D_j=D=D_0
such that S_2(k)=S(k) and C_2(𝐪;𝐤)=c(k)
also hold. In this case, ω_q=q^2D_0S^-1(q)
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(k)
in Eq.(<ref>). Interestingly, if we approximately replaces
D_j by its ensemble average value ⟨ D_j⟩
in the summation, we can obtain

    S_2(k)    ≃1/D_0N⟨∑_i,j⟨ D_j⟩ e^-i𝐤·𝐫_j+i𝐤·𝐫_i⟩
        =D̅/D_0N⟨∑_i,je^-i𝐤·𝐫_j+i𝐤·𝐫_i⟩ =D̅/D_0S(k)

in the second equality, we use the fact that ⟨ D_j⟩ =D̅
for a homogenous system. We find then 

    C_2(𝐪,𝐤)=ρ^-1[1-D_0/D̅S_2(p)/S(p)S^-1(k)]≃ c(k)

and the memory function Eq.(<ref>) reduces to that
for a passive system with effective diffusion coefficient D̅.
Note that this approximation holds if the coupling of D_j and
density fluctuation e^-i k·( r_j-𝐫_l)
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(k)≃(D̅/D_0)S(k)
is a very good approximation when k is large, whereas for small
k S_2(k)D_0/D̅S(k) 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, γ=D_t=k_BT=1
and we choose D_f, τ_p together with the number density
ρ as adjustable parameters. As already discussed in the model
description part, now the effective temperature is given by T_eff=1+D_fτ_p^2,
while the amplitude of active force is quantified by v_0=√(3D_fτ_p).
To compare our results with those simulation works of Ni and others,
we will choose v_0 and τ_p as independent variables together
with ρ. Nevertheless, we will also study the behavior of the
system by choosing T_eff and τ_p as independent
free parameters since it has been a regular choice in recent studies<cit.>.
To begin, we will run the system until it reaches the steady state
from which we can get the parameter D̅ and the function S_2(k),
with which the memory function Eq. (<ref>) can be numerically
calculated. We can then investigate the time dependence of F_q(t)
to address the glass transition issue. 

For the self-intermediate scattering function F_q^s(t),
the memory function equation reads(see the Appendix <ref>)

    ∂/∂ tF_q^s(t)+ω_q^sF_q^s(t)+∫_0^tduM_s^ irr(q,t-u)∂/∂ uF_q^s(u)=0

where ω_q^s=q^2D̅ and 

    M̃_s^ irr(q;t)=ρD̅/(2π)^3∫ d^3𝐤[(𝐤·𝐪̂)c(k)+(𝐩·𝐪̂)1/ρ(1-D_0S_2(k)/D̅S(k))]^2F_k(t)F_p^s(t)


If S_2(k)≃(D̅/D_0)S(k),
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 Brader<cit.> had tried to develop
a MCT for the RAB model in the limit τ_p→0. In this circumstance,
the effective Smoluchowski operator is actually given by Eq.(<ref>).
Starting from this effective operator, they obtained a memory function
for the collective scattering function F_q^eq(t),
but defined for the equilibrium distribution. In our work, the effective
Smoluchowski operator is now extended to finite (small) τ_p,
and importantly, the scattering function is now defined for the NESS
which is more relevant for the active system as pointed out by Szamel<cit.>.
The extension to finite τ_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 D̅,
S_2(k) in the present work must be obtained from simulations,
which is currently not avoidable. 

On the other hand, Szamel et al.<cit.>
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 𝐟_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(t). Importantly,
their theory involved a function ω_||(q) which
highlights the role of the velocity correlations. In particular, this
theory reproduced a nontrivial non-monotonic dependence of the relaxation
time τ_α with τ_p if T_eff is fixed
which was observed in their simulations for a standard LJ system,
although the theory apparently overestimated τ_α in the
τ_p→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.(<ref>) 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 D̅ and a
static structure function S_2(k) wherein both involves
the instantaneous diffusion coefficient D_j(𝐫^N).
Interestingly, although our method are quite different with that of
Szamel, we note that ω_||(q)τ_p in their
work plays the same role as D̅ in ours. We also note that
in a recent paper, Marconi et al.<cit.>
had studied the velocity correlations in the AOU-model, finding an
expression very similar to D_j(𝐫^N) under
mean-field approximation. 



§ NUMERICAL RESULTS AND DISCUSSIONS




 §.§ Static Properties 


As discussed above, to solve F_q(t), we must obtain
the parameter D̅ and the pseudo-structure factor S_2(k)
in the NESS via direct numerical simulations. In Fig.<ref>(a),
the dependence of D̅ on the effective temperature T_eff
is presented, for different fixed values of τ_p and ρ.
As can be seen, D̅ increases monotonically with T_eff,
which is reasonable since D̅ denotes a kind of averaged diffusion
coefficient that should be larger for a higher temperature. If T_eff
is fixed, D̅ decreases with increasing τ_p and the
variation of D̅ with T_eff becomes less sharp,
i.e., (∂D̅/∂ T_eff) decreases.
Such qualitative behaviors are robust with the change of number density
ρ, despite that the value of D̅ decreases slightly with
increasing ρ for given values of τ_p and T_eff.
In Fig.<ref>(b), we have also plotted D̅ as a function
of persistent time τ_p for different values of T_ eff
and ρ. In this case, we can see that D̅ decreases with
τ_p, tending to approach D_0=k_BT_ eff/γ
when τ_p→0 and close to D_t at a large τ_p
value. Besides, at the lower density D̅ has a slightly larger
value as shown in (a). 



In Fig.<ref>(a), the non-equilibrium static structure factors
S(k) obtained from direct simulations are drawn for
different particle activities v_0 for fixed ρ=1.12 and
τ_p=0.167. The value of ρ is chosen such that the system
is close to the glass transition and that of τ_p is consistent
with the setting in the simulation work of Ni<cit.>.
It seems that S(k) 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 τ_p, the effective temperature T_eff∼1+D_fτ_p^2
changes in the same tendency as v_0, such that Fig.<ref>(a)
also shows the change of S(k) with T_eff.
In Fig.<ref>(b), S(k) for different τ_p,
but with fixed T_eff have been presented. In this case,
one can see that the main peak increases apparently with increasing
τ_p and also shifts a little bit to smaller values of k.
Since T_eff∼1+v_0^2τ_p/3, increasing τ_p
with fixed T_eff corresponds to decreasing v_0, this
observation is consistent with Fig.<ref>(a). The second and
third peak also show observable differences with the variation of
τ_p in that the peak gets higher and moves to smaller values
of k with increasing τ_p. 



As discussed in the last section, an important new feature of our
theory is the introduction of the function S_2(k),
which couples the instantaneous diffusion coefficient D_j and
the density fluctuations. It is therefore instructive for us to investigate
how S_2(k) looks like. In Fig.<ref>, we have
plotted S_2(k) with the same parameter settings as
in Fig.<ref>. As shown in Fig.<ref>(a), the particle
activity (or effective temperature) drastically influences S_2(k),
with the main peak decreasing considerably with increasing v_0
or T_eff. Compared to Fig.<ref>(a), the value
of S_2(k) is much smaller than S(k),
which reflects the fact that D_j is generally less than D_0.
The behaviors of S_2(k) for fixed T_eff
but with variant τ_p are shown in Fig.<ref>(b). In
this latter case, we see that the main peak now slightly reduces with
increasing τ_p and it seems to saturate for large τ_p,
which are at variance with the observations in Fig.<ref>(b).
The apparent discrepancies between S_2(k) and S(k)
indicate that our theory may show interesting new features.



Another new feature of our theory is the pseudo-correlation function
C_2(𝐪,𝐤), which plays a similar
role to c(k) in the irreducible memory function M^irr(q,t).
As discussed above, C_2(𝐪,𝐤) reduces
to c(k) if S_2(p)D_0/D̅S(p)≃1.
In Fig.<ref>, the dependence of S_2(k)D_0/D̅S(k)
on k has been presented, for fixed ρ with varying T_ eff
and τ_p. Interestingly, we find that it is approximately one
if k is larger than 2π/σ which is approximately the peak
position for S(k). Nevertheless, for small values of
k, S_2(k)D_0/D̅S(k) shows some
structures. Specifically, S_2(k)D_0/D̅S(k)
becomes much less than one for small τ_p if T_eff
is fixed. Such a feature may lead to enhancement of the irreducible
memory function M_s^irr(q,t) shown in Eq.(<ref>)
with decreasing τ_p→0 if T_eff is fixed. This
would lead to the increment of τ_α,if other effects
are not accounted for. Note that, however, D̅ increases with
decreasing τ_p with constant T_eff as shown in
Fig.<ref>, such that τ_α would decrease with
decreasing τ_p in terms of this effect. Therefore, it might
be possible that the relaxation time τ_α shows some re-entrance
behaviors in the small τ_p region, similar to that reported
for the AOU model<cit.>. 






 §.§ Intermediate Scattering Function 


With the static properties obtained above, particularly D̅
and S_2(k), we are ready to investigate the behavior
of the intermediate scattering function F_q(t) by
numerically solving the memory functions Eqs. (<ref>) and
(<ref>). In Fig.<ref>(a), the normalized scattering
functions ϕ_q(t)=F_q(t)/S(q)
are shown for different values of v_0 (or T_eff)
and number density ρ, wherein we have chosen q=7.5 which
is around the first peak of S(q). The results for two
densities ρ=1.05 and 1.10 are plotted, and the value of τ_p
is fixed to be 0.167. For the higher density ρ=1.07, one can
see that F_q(t) finally reaches a plateau in the long
time limit for v_0=0 (or T_eff=1), indicating that
the system reaches the glassy state. For a nonzero value of v_0
as shown in the figure, F_q(t) 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 ρ=1.05, the system is
in the liquid state for v_0=0 and the relaxation of F_q(t)
also becomes faster with increasing v_0 or T_eff.
The behaviors of the self-scattering function F_q^s(t)
are similar as shown in Fig.<ref>(b). While for v_0=0
the tracer particle is trapped and F_q^s(t) reaches
a non-zero value for t→∞, it relaxes to zero for v_0=10
and 20 with the relaxation time τ_α decreases apparently
with increasing v_0. 



The limiting value f_q=lim_t→∞ϕ_q(t)
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 ρ, f_q may change from zero to an
apparent nonzero value, and the value ρ_c thus corresponds
to the glass transition point. One may also fix ρ but vary T_eff,
then f_q may become nonzero for T_eff less than some
critical value T_eff^c , which defines a critical temperature
for glass transition. According to the MCT Eq. (<ref>),
the Debye-Waller factor f_q follows

    f_q=m_q(∞)/1+m_q(∞)

where 

    m_q(∞)=ρD̅/16π^3q^2∫ d^3𝐤[(𝐪̂·𝐤)C_2(𝐪;𝐤)+(𝐪̂·𝐩)C_2(𝐪;𝐩)]^2S(k)S(q)S(p)f_kf_p

This equation can be solved numerically and self-consistently to get
f_q for given control parameters v_0 (or T_eff),
τ_p as well as ρ.

In Fig.<ref>(a), the dependence of f_q on the number
density ρ is presented for different v_0 (or T_eff)
with given τ_p=0.001. Clearly, f_q changes abruptly from
zero to a large nonzero value at a critical density ρ_c, indicated
for example by the vertical dashed line for v_0=0 at about ρ_c≃1.064.
For given τ_p, the curve shifts to larger values of ρ
with increasing v_0, indicating that that glass transition is
pushed to higher values of ρ for larger particle activity in
consistent with previously reported simulation results. The pictures
for different values of τ_p=0.05 and 0.167 are shown in Fig.<ref>(b)
and (c), respectively. The results are similar to those in (a), with
the values of ρ_c shifting to relatively larger values for
larger τ_p. 



In Fig.<ref>(b), it is shown that the relaxation time τ_α
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 τ_α as a
function of ρ. The results are depicted in Fig.<ref>
with the same parameter setting as in Fig.<ref>. Obviously,
τ_α increases fastly with ρ for fixed values of
v_0 (T_eff) and τ_p and it diverges at some
critical value ρ_c. For a very small τ_p=0.001, it
seems that changing v_0 does not affect very much the values
of τ_α as shown in Fig.<ref>(a). The influence
becomes more considerable when τ_p gets larger as demonstrated
in <ref>(b) and (c), and the value of ρ_c also
shifts to larger values in consistent with Fig.<ref>. 





Surely, the value of ρ_c should be the same either obtained
by f_q or τ_α within reasonable fluctuations. In
Fig.<ref>(a) and (b), the dependence of ρ_c, obtained
from bothf_q and τ_α, on v_0 and T_eff
are presented for different given values of τ_p. Clearly,
ρ_c increases with both v_0 and T_eff as
expected. Interestingly, ρ_c shows a nearly linear dependence
on v_0^2, wherein the slope increases with τ_p. This
linear dependence was also observed in the simulation work of Ni.
We also note that ρ_c increases with τ_p for fixed
v_0, whereas it decreases with τ_p for fixed T_eff
according to the data presented in Fig.<ref>(a) and (b).
This is shown more clearly in Fig.<ref>(c) , where we
have also plotted ρ_c as a function of τ_p for different
T_eff. For comparison, the dashed line gives the value
of ρ_c^B for the corresponding passive Brownian system with
T=T_eff, which is obtained by setting D_t=D_0 and
zero self-propulsion force f_i=0 in Eq.(<ref>). Clearly,
ρ_c approaches ρ_c^B in the limit τ_p→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_eff^c with fixed number density
ρ. In Fig.<ref> (a) to (c), the results of
τ_α are presented as functions of T_eff for
different τ_p and ρ. For fixed value of τ_p
and ρ, τ_α decreases monotonically with T_eff.
Below some critical values of T_eff corresponding to T_eff^c,
τ_α diverges indicating the occurrence of glass transition.
In Fig.<ref> (d), the dependence of T_eff^c
on τ_p for different ρ is shown. One can see that T_eff^c
increases monotonically with τ_p, and it approaches a constant
value in the limit τ_p→0. Such a constant value corresponds
to the one for a passive Brownian system with T_c=T_eff^c.
We also note that T_eff^c increases with the number
density ρ, indicating that a denser system enters glass transition
at a higher critical temperature as expected. 




§ 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 τ_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(𝐫^N) 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(t) or F_q^s(t) 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 D̅,
which is the ensemble average of D_j(𝐫_N)
in the nonequilibrium steady state, and a pseudo steady state structure
factor S_2(k), which involves the coupling between
D_j(𝐫^N) and density fluctuations. 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 D̅
increases with the single effective temperature T_eff
as well as the magnitude v_0 of propulsion force, while it decreases
with τ_p for fixed T_eff or v_0. The structure
function S_2(k) simply decouples into the product
of D̅/D_0 and S(k), with S(k)
the nonequilibrium static structure factor and D_0 the single
particle diffusion coefficient in the limit τ_p→0, for relatively
large values of k, whereas it shows apparent deviations from (D̅/D_0)· S(k)
for small ks. Our theory makes it feasible for us to investigate
the glassy dynamics of the system, by investigating the time behavior
of F_q(t) or F_q^s(t) in the long
time limit, chosen the persistence time τ_p, the effective
temperature T_eff, as well as the number density ρ
as free parameters. We find that the critical density ρ_c
for glass transition shifts to larger values with increasing T_eff
or v_0 if τ_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
ρ_c changes with τ_p for a fixed T_eff,
finding that ρ_c decreases with τ_p monotonically
and it approaches the value for the corresponding passive Brownian
system in the limit τ_p→0 as expected. We have also calculated
the critical temperature T_eff^c for glass transition
at fixed density, finding that it increases monotonically with τ_p
and also approaches the Brownian particle limit for τ_p→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 particles<cit.>.
As mentioned in the main text, the relaxation time τ_α
may show nontrivial dependence on the persistence time τ_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 τ_p→0, the glass transition
point ρ_c or T_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 <cit.>,
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. 

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). 


41
natexlab#1#1bibnamefont#1#1bibfnamefont#1#1citenamefont#1#1url<#>1urlprefixURL 


[Vicsek and
  Zafeiris(2012)]2012_PhysRep_Vicsek_CollectiveMotion
authorT. Vicsek and
  authorA. Zafeiris,
  journalPhys. Rep. volume517,
  pages71 (year2012).

[Marchetti and
  Joanny(2013)]2013_RMP_Marchetti_HydroSoftActi
authorM. C. Marchetti
  and authorJ. F.
  Joanny, journalRev. Mod. Phys.
  volume85, pages1147 (year2013).

[Bechinger et al.(2016)Bechinger,
  Di Leonardo, Löwen, Reichhardt, Volpe, and
  Volpe]2016_RMP_Bechinger_ActiParti_in_complex_enviro
authorC. Bechinger,
  authorR. Di Leonardo,
  authorH. Löwen,
  authorC. Reichhardt,
  authorG. Volpe, and
  authorG. Volpe,
  journalReviews of Modern Physics volume88,
  pages045006 (year2016).

[Schaller et al.(2010)Schaller, Weber,
  Semmrich, Frey, and Bausch]2010_Nature_Schaller_SoftFilament
authorV. Schaller,
  authorC. Weber,
  authorC. Semmrich,
  authorE. Frey, and
  authorA. R. Bausch,
  journalNature volume467, pages73
  (year2010).

[Sumino1 et al.(2012)Sumino1, Nagai,
  Shitaka, Tanaka, Yoshikawa, Chate, and
  Oiwa]2012_Nature_Suminol_HardFilament
authorY. Sumino1,
  authorK. H. Nagai,
  authorY. Shitaka,
  authorD. Tanaka,
  authorK. Yoshikawa,
  authorH. Chate, and
  authorK. Oiwa,
  journalNature volume483, pages48
  (year2012).

[Schwarz-Linek et al.(2012)Schwarz-Linek,
  Valeriani, Cacciuto, Cates, Marenduzzo, Morozov, and
  Poon]2012_PNAS_Linek_PhaseSepa
authorJ. Schwarz-Linek,
  authorC. Valeriani,
  authorA. Cacciuto,
  authorM. E. Cates,
  authorD. Marenduzzo,
  authorA. N. Morozov,
  and authorW. C. K.
  Poon, journalPNAS volume109,
  pages4052 (year2012).

[Fily and
  Marchetti1(2012)]2012_PRL_Fily_AthermalPhaseSepa
authorY. Fily and
  authorM. C. Marchetti1,
  journalPhys. Rev. Lett. volume108,
  pages235702 (year2012).

[Buttinoni et al.(2013)Buttinoni,
  Bialké, Kümmel, Löwen, Bechinger, and
  Speck]2013_PRL_Buttinoni_ClustPhasSepa
authorI. Buttinoni,
  authorJ. Bialké,
  authorF. Kümmel,
  authorH. Löwen,
  authorC. Bechinger,
  and authorT. Speck,
  journalPhysical review letters volume110,
  pages238301 (year2013).

[Wittkowski et al.(2014)Wittkowski,
  Tiribocchi, Stenhammar, Allen, Marenduzzo, and
  Cates]2014_NatComm_Cates_Chi4Theory
authorR. Wittkowski,
  authorA. Tiribocchi,
  authorJ. Stenhammar,
  authorR. J. Allen,
  authorD. Marenduzzo,
  and authorM. E. Cates,
  journalNature Communications volume5,
  pages4351 (year2014).

[Das et al.(2014)Das, Egorov, Trefz,
  Virnau, and Binder]2014_PRL_Das_AP_Mix
authorS. K. Das,
  authorS. A. Egorov,
  authorB. Trefz,
  authorP. Virnau, and
  authorK. Binder,
  journalPhys. Rev. Lett. volume112,
  pages198301 (year2014).

[Berthier and
  Kurchan(2013)]2013_NP_Berthier_NonEquiGlasTrans
authorL. Berthier and
  authorJ. Kurchan,
  journalNature Physics volume9,
  pages310 (year2013).

[Ni et al.(2013)Ni, Stuart, and
  Dijkstra]2013_NC_NiRan_PushGlasTrans
authorR. Ni,
  authorM. A. C. Stuart,
  and authorM. Dijkstra,
  journalNature Communications volume4,
  pages1 (year2013).

[Berthier(2013)]2013_NatPhys_Kurchan_NonEquiGlasTrans
authorL. Berthier,
  journalNature Physics volume9,
  pages310 (year2013).

[Mandal et al.(2014)Mandal, Lang, Gross,
  Oettel, Raabe, Franosch, and Varnik]2014_NatComm_Mandal_ReGT
authorS. Mandal,
  authorS. Lang,
  authorM. Gross,
  authorM. Oettel,
  authorD. Raabe,
  authorT. Franosch, and
  authorF. Varnik,
  journalNature Communications volume5,
  pages1 (year2014).

[Berthier(2014)]2014_PRL_Berthier_GTActiveHardDisk
authorL. Berthier,
  journalPhys. Rev. Lett. volume112,
  pages220602 (year2014).

[Sharma et al.(2017)Sharma, Wittmann, and
  Brader]2017_PRE_Brader
authorA. Sharma,
  authorR. Wittmann, and
  authorJ. M. Brader,
  journalPhys. Rev. E volume95,
  pages012115 (year2017).

[Angelini et al.(2011)Angelini, Hannezo,
  Trepat, Marquez, Fredberg, and Weitz]2011_PNAS_Angelini_Exp_Cells
authorT. E. Angelini,
  authorE. Hannezo,
  authorX. Trepat,
  authorM. Marquez,
  authorJ. J. Fredberg,
  and authorD. A. Weitz,
  journalPNAS volume108, pages4714
  (year2011).

[Krakoviack(2005)]2005_PRL_Krakoviack_porous
authorV. Krakoviack,
  journalPhys. Rev. Lett. volume94,
  pages065703 (year2005).

[Whitesides and
  Grzybowski(2002)]2002_Science_Whitesides_Self-assemble
authorG. M. Whitesides
  and
  authorB. Grzybowski,
  journalScience volume295,
  pages2418 (year2002).

[ten Hagen et al.(2011)ten Hagen, van
  Teeffelen, and Löwen]2011_JPCM_Hagen_ActBrown
authorB. ten Hagen,
  authorS. van Teeffelen,
  and
  authorH. Löwen,
  journalJournal of Physics: Condensed Matter
  volume23, pages194119
  (year2011).

[Ni et al.(2014)Ni, Stuart, Dijkstra, and
  Bolhuis]2014_SoftMatt_NiRan_CrystHardSpheGlas
authorR. Ni,
  authorM. A. C. Stuart,
  authorM. Dijkstra, and
  authorP. G. Bolhuis,
  journalSoft Matter volume10,
  pages6609 (year2014).

[Levis and
  Berthier(2014)]2014_PRE_Berthier_ClustKMCSelfPropHardSph
authorD. Levis and
  authorL. Berthier,
  journalPhys. Rev. E volume89,
  pages062301 (year2014).

[Flenner et al.(2016)Flenner, Szamel, and
  Berthier]2016_arXiv_Flenner_NonEquiGlasTrans
authorE. Flenner,
  authorG. Szamel, and
  authorL. Berthier,
  journalarXiv preprint arXiv:1606.00641  (year2016).

[Julian Bialke(2012)]2012_PRL_Bialke_CrysDenseSuspension
authora. H. L. Julian Bialke,
  Thomas Speck, journalPhys. Rev. Lett.
  volume108, pages168301
  (year2012).

[Li et al.(2015)Li, Jiang, and
  Hou]2015_SoftMatt_LiShuxian
authorS. Li,
  authorH. Jiang, and
  authorZ. Hou, journalSoft
  matter volume11, pages5712
  (year2015).

[Ding et al.(2015)Ding, Feng, Jiang, and
  Hou]2015_our_Arxiv
authorH. Ding,
  authorM. Feng,
  authorH. Jiang, and
  authorZ. Hou,
  journalarXiv preprint arXiv:1506.02754  (year2015).

[Farage and
  Brader(2014)]2014_arXiv_Farage_MCTActiveGlass
authorT. F. F. Farage
  and authorJ. M.
  Brader, journalarXiv pp. pages1–5
  (year2014).

[Farage et al.(2015)Farage, Krinninger,
  and Brader]2015_PRE_Farage_EffectIntera
authorT. F. Farage,
  authorP. Krinninger,
  and authorJ. M.
  Brader, journalPhysical Review E
  volume91, pages042310
  (year2015).

[Nandi(2016)]2016_arXiv_Nandi_Activity
authorS. K. Nandi,
  journalarXiv preprint arXiv:1605.06073  (year2016).

[Szamel et al.(2015)Szamel, Flenner, and
  Berthier]2015_PRE_Szamel_GlasDyna
authorG. Szamel,
  authorE. Flenner, and
  authorL. Berthier,
  journalPhysical Review E volume91,
  pages062304 (year2015).

[Szamel(2014)]2014_PRE_Szamel_EffecTemp
authorG. Szamel,
  journalPhysical Review E volume90,
  pages012111 (year2014).

[Fox(1986)]1986_PRA_Fox
authorR. F. Fox,
  journalPhysical Review A volume33,
  pages467 (year1986).

[Nägele(1996)]1996_PhysRep_Nagele_charged_suspensions
authorG. Nägele,
  journalPhys. Rep. volume272,
  pages215 (year1996).

[Cichocki and
  Hess(1987)]1987_PhysA_Cichocki_irreducible
authorB. Cichocki and
  authorW. Hess,
  journalPhys. A volume141,
  pages475 (year1987).

[Kawasaki(1994)]1994_PhysA_Kawasaki_irreducible
authorK. Kawasaki,
  journalPhys. A volume208,
  pages35 (year1994).

[Szamel(2016)]2016_PRE_Szamel_AthermalActi
authorG. Szamel,
  journalPhysical Review E volume93,
  pages012603 (year2016).

[Maggi et al.(2015)Maggi, Marconi, Gnan,
  and Di Leonardo]2015_ScitiRep_Maggi_Multidimen_Stationary
authorC. Maggi,
  authorU. M. B. Marconi,
  authorN. Gnan, and
  authorR. Di Leonardo,
  journalScientific reports volume5
  (year2015).

[Levis and
  Berthier(2015)]2015_EPL_Levis_Single-collective-effecTemp
authorD. Levis and
  authorL. Berthier,
  journalEPL (Europhysics Letters) volume111,
  pages60006 (year2015).

[Kawasaki(1995)]1995_Kawasaki_irreducible
authorK. Kawasaki,
  journalPhysica A: Statistical Mechanics and its Applications
  volume215, pages61 (year1995).

[Hansen and McDonald(1990)]book_Hansen_TSL
authorJ.-P. Hansen and
  authorI. R. McDonald,
  titleTheory of simple liquids
  (publisherElsevier, year1990).

[Götze(2008)]book_gotze2008complex
authorW. Götze,
  titleComplex dynamics of glass-forming liquids: A
  mode-coupling theory, vol. volume143 (publisherOUP
  Oxford, year2008).




equationsection 





§ DERIVATION OF THE SMOLUCHOWSKI EQUATION


Generally, for a LE with colored noise, one can get the FPE within
Fox approximation<cit.>. For illustration, consider
a simple one dimensional over damped LE 

    ẋ(t)=G(x)+χ(t)

where G(x) denotes the external or internal force and χ(t)
is the stochastic noise with correlation 

    ⟨χ(t)χ(s)⟩ =C(t-s)

Define a probability distribution function 

    P(y,t)=∫ D[χ]P[χ]δ(y-x(t))

where D[χ] denotes integration over the noisy path
of χ(t) and P[χ] is the distribution
functional of χ which is assumed to be Gaussian. One can then
obtain the FPE governing the evolution of P(y,t) as
follows <cit.>

    ∂/∂ tP(y,t)=-∂/∂ y[G(y)P(y,t)]+∂^2/∂ y^2{∫_0^tds'C(t-s')∫ D[χ]P[χ]e^∫_s'^tdsG'(x(s))δ(y-x(t))}

Note that if χ(t) is white noise, C(t-s)=D_0δ(t-s),
then the second term is just 

    D_0∂^2/∂ y^2[∫ D[χ]P[χ]δ(y-x(t))]≡ D_0∂^2/∂ y^2P(y,t)

which recovers the standard FPE. For a colored noise with 

    C(t-s)=D/τ_0exp(-|t-s|/τ_0)

one can obtain the FPE approximately as

    ∂/∂ tP(y,t)=-∂/∂ y[G(y)P(y,t)]+D∂^2/∂ y^2[1/1-τ_0G'(y)P(y,t)]

by assuming 

    ∫_s'^tdsG'(x(s))≈ G(x(t))(t-s')


For a general multi-variable case, 

    dx_i(t)/dt=G_i({x_i})+χ_i(t)

where

    ⟨χ_i(t)χ_j(s)⟩ =C_ij(t-s)

with i,j=1,2,⋯,N, the FPE for distribution function 

    P(𝐲,t)=∫ D[χ]P[χ]δ(𝐲-𝐱(t))

reads

    ∂/∂ tP(𝐲,t)     =     -∑_i∂_i[G_i(𝐲)P(𝐲,t)]
           +∑_ij∂_i{∑_l∫_0^tds'C_il(t-s')∂_j∫ D[χ]P[χ]exp[∫_s'^tds∂/∂ x_lG_j(𝐱(s))δ_jl]δ(𝐲-𝐱(t))}
        =     -∑_i∂_i[G_i(𝐲)P(𝐲,t)]
           +∑_ij∂_i{∫_0^tds'C_ij(t-s')∂_j∫ D[χ]P[χ]exp[∫_s'^tds∂_jG_j(x(s))]δ(𝐲-𝐱(t))}


Then, if C_ij(t-s)=δ_ijC(t-s)=δ_ijD/τ_0exp(-|t-s|/τ_0),
using the assumption mentioned before we can get

    ∂/∂ tP(𝐲,t)=-∑_i∂_i[G_i(𝐲)P(𝐲,t)]+∑_iD∂_i^2{[1/1-τ_0∂_iG_i(𝐲)]P(𝐲,t)}


For our system described by Eq.(<ref>), we have correspondingly
x→𝐫^N, G(x)→γ^-1𝐅(𝐫^N)=β D_t𝐅(𝐫^N),
γ^-1𝐟_i(t)→χ_i(t).
According to Eq.(<ref>), the variable D in Eq. (<ref>)
is D_fτ_p^2/γ^2 and τ_0 is τ_p.
Note that the white noise term η_i(t)
in Eq.(<ref>) contributes a normal diffusion term to the
FPE. Thus we finally obtain

    ∂/∂ tΨ(𝐫^N,t)     =     -∑_i∇_i·[β D_t𝐅_i(𝐫^N)-D_t∇_i]Ψ(𝐫^N,t)
           +∑_i∇_i^2[D_fτ_p^2/γ^2/1-τ_p·β D_t∇_i𝐅_i(𝐫^N,t)]Ψ(𝐫^N,t)
        =     +∑_i∇_i·{ D_t+[D_fτ_p^2/γ^2/1-τ_p·β D_t∇_i𝐅_i(𝐫^N)]}·∇_iΨ(𝐫^N,t)
           -∑_i∇_i·{β D_t𝐅_i(𝐫^N)-∇_i[D_fτ_p^2/γ^2/1-τ_p·β D_t∇_i𝐅_i(𝐫^N)]}Ψ(𝐫^N,t)

Write

    D_i(𝐫^N)=D_t+[D_fτ_p^2/γ^2/1-τ_p·β D_t∇_i𝐅_i(𝐫^N)]

and 

    𝐅_i^(𝐫^N)     =    D_t/D_i(𝐫^N){𝐅_i(𝐫^N)-1/β D_t∇_i[D_fτ_p^2/γ^2/1-τ_p·β D_t∇_i𝐅_i(𝐫^N)]}
        =    D_t/D_i(𝐫^N)[𝐅_i(𝐫^N)-1/β D_t∇_iD_i(𝐫^N)]

is the effective force. Finally, we have 

    ∂/∂ tΨ(𝐫^N,t)=-∑_i∇_i· D_i(𝐫^N)·[∇_i-β𝐅_i^(𝐫^N)]Ψ(𝐫^N,t)

which corresponds exactly to Eqs. (<ref>) to (<ref>)
in the main text. 



§ DERIVATION OF THE GENERAL LANGEVIN EQUATION




 §.§ Memory Function Equation 


Here we present the derivation of the memory function equations, namely
Eqs. (<ref>) to (<ref>) in the main text, for
the scattering function F_q(t)=1/N⟨_𝐪^*e^Ω̂tρ_𝐪⟩.
This is most easily done in the Laplace domain, even for the complex
Smoluchowski operator Ω̂ shown in Eq.(<ref>)
which contains the instantaneous effective diffusion constant D_j(𝐫^N)
given by Eq.(<ref>). 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 

    F̃(q,z)=ℒ𝒯[F_q(t)]=⟨ A_-𝐪1/z-Ω̂A_𝐪⟩

where ℒ𝒯 stands for Laplace transformation and A_𝐪=ρ_𝐪/√(N).
One can define a projection operator on the density 

    𝒫(⋯)=A_𝐪⟩⟨ A_𝐪A_-𝐪⟩ ^-1⟨ A_-𝐪(⋯)⟩

which has the property

    𝒫A_𝐪=A_𝐪𝒫𝒫=𝒫,𝒫^n=𝒫.

Accordingly, we can define 𝒬=ℐ-𝒫,
which satisfies 𝒬A_𝐪=0, 𝒫𝒬=0,
and 𝒬^n=𝒬. Then for the operator [z-Ω]^-1,
one has the following identity

    1/z-Ω̂=1/z-𝒬+1/z-𝒬𝒫1/z-Ω̂

which is known as Dyson decomposition and can be easily checked by
right-multiplying both sides by z-Ω̂. 

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


    ℒ𝒯[∂_tF(q,t)](z)     =     zF̃(q,z)-F(q,t=0)=⟨ A_-𝐪Ω̂1/z-Ω̂A_𝐪⟩
        =    ⟨ A_-𝐪Ω̂𝒫1/z-Ω̂A_𝐪⟩ +⟨ A_-𝐪Ω̂𝒬1/z-Ω̂A_𝐪⟩

Using the definition of 𝒫, the first term is 

    ⟨ A_-𝐪Ω̂𝒫1/z-Ω̂A_𝐪⟩    =    ⟨ A_-𝐪Ω̂A_𝐪⟩⟨ A_-𝐪A_𝐪⟩ ^-1⟨ A_-𝐪1/z-Ω̂A_𝐪⟩
        =    ⟨ A_-𝐪Ω̂A_𝐪⟩⟨ A_-𝐪A_𝐪⟩ ^-1F̃(q,z)

While the second term is, using the identity(<ref>), 

    ⟨ A_-𝐪Ω̂𝒬1/z-Ω̂A_𝐪⟩ =⟨ A_-𝐪Ω̂𝒬[1/z-𝒬+1/z-𝒬𝒫1/z-Ω̂]A_𝐪⟩

Note that 𝒬A_𝐪=0, hence 1/z-𝒬A_𝐪=0
and 

    ⟨ A_-𝐪Ω̂𝒬1/z-Ω̂A_𝐪⟩    =    ⟨ A_-𝐪Ω̂𝒬1/z-𝒬𝒫1/z-Ω̂A_𝐪⟩
        =    ⟨ A_-𝐪Ω̂𝒬1/z-𝒬A_𝐪⟩⟨ A_-𝐪A_𝐪⟩ ^-1⟨ A_-𝐪1/z-Ω̂A_𝐪⟩
        =    ⟨ A_-𝐪Ω̂𝒬1/z-𝒬𝒬𝒬A_𝐪⟩⟨ A_-𝐪A_𝐪⟩ ^-1F̃(q,z)

where we have used the definition of 𝒫 in the second
equality and the fact 𝒬𝒬=𝒬 in the
third equality. We may introduce 

    ω_q=-⟨ A_-𝐪Ω̂A_𝐪⟩⟨ A_-𝐪A_𝐪⟩ ^-1

which is Eq.(<ref>) and define 

    M̃(q,z)=-⟨ A_-𝐪Ω̂𝒬1/z-𝒬𝒬𝒬A_𝐪⟩⟨ A_-𝐪Ω̂A_𝐪⟩ ^-1

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

    ℒ𝒯[∂_tF(q,t)](z)=zF̃(q,z)-F(q,t=0)=-ω_q[1-M̃(q,z)]F̃(q,z)

Therefore 

    F̃(q,z)=F(q,t=0)/z+ω_q[1-M̃(q,z)]


For colloidal systems, there is a so-called irreducible issue <cit.>
. Following the procedure in <cit.>
 one needs to introduce an irreducible memory function M̃^ irr(q,z),
which is related toM̃(q,z) according to 

    M̃(q,z)=M̃^ irr(q,z)[1+M̃^ irr(q,z)]^-1

and 

    M̃^ irr(q,z)=-⟨ A_-𝐪Ω̂𝒬1/z-𝒬^ irr𝒬𝒬A_𝐪⟩⟨ A_-𝐪Ω̂A_𝐪⟩ ^-1

Herein, Ω̂^irr denotes an irreducible Smoluchowski
operator of which the detailed form is not relevant within the MCT
approximation below. Consequently, this leads to

    F̃(q,z)=F(q,t=0)/z+ω_q/1+M̃^ irr(q,z)

corresponding to in the time domain 

    ∂/∂ tF_q(t)+ω_qF_q(t)+∫_0^tduM^ irr(q,t-u)∂/∂ uF_q(u)=0,

which is exactly Eq.(<ref>) in the main text. Note that
above derivations are quite general and do not depend on the explicit
form of the operator Ω̂, whereas the expressions for
ω_q and M^irr should certainly depend on Ω̂. 



 §.§ Frequency ω_q


We now substitute A_ q=ρ_𝐪/√(N)=∑_je^-i q· r_j/√(N)
to calculate ω_q. Note that 

    ⟨ A_-𝐪Ω̂A_𝐪⟩    =∫ d𝐫^NA_-𝐪∑_j∇_j· D_j[∇_j-β𝐅_j^eff]A_𝐪P_s(𝐫^N)
        =∫ d𝐫^NA_-𝐪∑_j∇_j· D_j{∇_j[A_𝐪P_s(𝐫^N)]-β𝐅_j^effA_𝐪P_s(𝐫^N)}

where as mentioned before, the operator Ω̂ acts on all
the functions to its right side including the steady-state distribution
function P_s(𝐫^N). In the steady state,
the summation of all the currents 𝐉_j^s, given by
Eq.(<ref>), is zero according to Ω̂P_s(𝐫^N)=-∑_j𝐉_j^s=0.
To proceed and as many authors have done, we assume more strongly
that 𝐉_j^s=0, i.e., 

    𝐉_j^s=-D_j(𝐫^N)[∇_j-β𝐅_j^ eff(𝐫^N)]P_s(𝐫^N)=0.

Therefore, one can obtain that 

    ∇_jP_s(𝐫^N)=β𝐅_j^ eff(𝐫^N)P_s(𝐫^N)

which is the counterpart of Yvon theorem<cit.> in
this nonequilibrium system.  Using this result, one has 

    ⟨ A_-𝐪Ω̂A_𝐪⟩    =∫ d𝐫^NA_-𝐪∑_j∇_j·[D_j(∇_jA_𝐪)P_s(𝐫^N)]
        =-∑_j∫ d𝐫^N[(∇_jA_-q)·(∇_jA_𝐪)]D_jP_s(𝐫^N)
        =-N^-1∑_j∫ d𝐫^Nq^2D_jP_s(𝐫^N)
        =-q^2∑_j⟨ D_j⟩ /N=-q^2D̅

where the second equality results from partial integration and we
have used ∇_jA_𝐪=-i𝐪exp(-i𝐪·𝐫_j)/√(N)
in the third one. ⟨ D_j⟩ =∫ d𝐫^ND_j(𝐫^N)P_s(𝐫^N)
denotes the averaged instantaneous diffusion function of particle
j and D̅=N^-1∑_j⟨ D_j⟩.
Therefore, the effective frequency ω_q reads 

    ω_q=-⟨ A_-𝐪Ω̂A_𝐪⟩/⟨ A_-𝐪A_𝐪⟩=q^2D̅/S(q)

which is Eq.(<ref>) in the main text.



 §.§ Memory Function M^irr(q,t)


In the time domain, the irreducible memory function M^irr(q,t)
is given by 

    M̃^ irr(q,t)=-⟨ A_-𝐪Ω̂𝒬e^𝒬Ω̂^ irr𝒬t𝒬A_𝐪⟩⟨ A_-𝐪Ω̂A_𝐪⟩ ^-1

Using the adjoint operator Ω̂^†, the first term
is 

    ⟨ A_-𝐪Ω̂𝒬e^𝒬Ω̂^ irr𝒬t𝒬Ω̂A_𝐪⟩    =    ⟨ A_-𝐪Ω̂𝒬e^𝒬Ω̂^ irr𝒬t𝒬(Ω̂^†A_𝐪)⟩
        =    ⟨(Ω̂^†A_-𝐪)𝒬e^𝒬Ω̂^ irr𝒬t𝒬(Ω̂^†A_𝐪)⟩
        =    ⟨(𝒬Ω̂^†A_-𝐪)𝒬e^𝒬Ω̂^ irr𝒬t𝒬(𝒬Ω̂^†A_𝐪)⟩
        =    ⟨ R_𝐪^*𝒬e^𝒬Ω̂^ irr𝒬t𝒬R_𝐪⟩

where 𝒬𝒬=𝒬 is used in the third equality and we have
introduced

    R_𝐪    =    𝒬(Ω̂^†A_𝐪)=(Ω̂^†A_𝐪)-𝒫(Ω^†A_𝐪)
        =     (Ω̂^†A_𝐪)-⟨ A_-𝐪(Ω̂^†A_𝐪)⟩/⟨ A_-𝐪A_𝐪⟩A_𝐪
        =     =(Ω̂^†A_𝐪)-⟨ A_-𝐪Ω̂A_𝐪⟩/⟨ A_-𝐪A_𝐪⟩A_𝐪
        =     (Ω̂^†A_𝐪)+ω_qA_𝐪

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 <cit.>. One can
then define a second-order projection operator

    𝒫_2≡1/2∑_ k, p|A_𝐩A_𝐤⟩⟨ A_𝐩^*A_𝐤^*A_𝐩A_𝐤⟩ ^-1⟨ A_𝐩A_𝐤|

and make the approximation:

    ⟨ R_𝐪^*e^𝒬Ω̂^ irr𝒬tR_𝐪⟩   ≈   ⟨ R_𝐪^*𝒫_2e^𝒬Ω̂^ irr𝒬t𝒫_2R_𝐪⟩
        =    1/4∑_ k, p∑_ k', p'⟨ R_𝐪^*A_𝐩A_𝐤⟩⟨ A_𝐩^*A_𝐤^*A_𝐩A_𝐤⟩ ^-1
          ×⟨ A_𝐩'^*A_𝐤'^*R_𝐪⟩⟨ A_𝐩'^*A_𝐤'^*A_𝐩'A_𝐤'⟩ ^-1
          ×⟨ A_𝐩A_𝐤e^𝒬Ω̂^ irr𝒬tA_𝐩'A_𝐤'⟩
       ≈   1/4∑_ k, p∑_ k', p'⟨ρ_𝐩'^*ρ_𝐤'^*R_𝐪⟩/NS(k')S(p')⟨ R_𝐪^*ρ_𝐩ρ_𝐤⟩/NS(k)S(p)
          ×1/N^2[δ_ pp'δ_ kk'⟨ρ_𝐩^*e^Ω̂tρ_𝐩'⟩⟨ρ_𝐤^*e^Ω̂tρ_𝐤'⟩.
           +.δ_ pk'δ_ kp'⟨ρ_𝐩^*e^Ω̂tρ_𝐤'⟩⟨ρ_𝐤^*e^Ω̂tρ_𝐩'⟩]
        =    1/2∑_ k, p|⟨ρ_𝐩^*ρ_𝐤^*R_𝐪⟩|^2/[N^2S(k)S(p)]^2⟨ρ_𝐩^*e^Ω̂tρ_𝐩⟩⟨ρ_𝐤^*e^Ω̂tρ_𝐤⟩

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

    ⟨ρ_𝐩^*ρ_𝐤^*ρ_𝐩'ρ_𝐤'⟩   ≈   ⟨ρ_𝐩^*ρ_𝐤'⟩⟨ρ_𝐤^*ρ_𝐩'⟩ +⟨ρ_𝐩^*ρ_𝐩'⟩⟨ρ_𝐤^*ρ_𝐤'⟩
        =    δ_ p, k'δ_ k, p'N^2S(p)S(k)+δ_ p, p'δ_ k, k'N^2S(p)S(k)

and simultaneously replace the projected operator 𝒬Ω̂^ irr𝒬
by the full Smoluchowski operator Ω̂ in the propagator
governing the time evolution of the correlation function <cit.>

    ⟨ A_𝐩A_𝐤e^𝒬Ω̂^ irr𝒬tA_𝐩'A_𝐤'⟩≈δ_ pp'δ_ kk'⟨ρ_𝐩^*e^Ω̂tρ_𝐩'⟩⟨ρ_𝐤^*e^Ω̂tρ_𝐤'⟩ +δ_ pk'δ_ kp'⟨ρ_𝐩^*e^Ω̂tρ_𝐤'⟩⟨ρ_𝐤^*e^Ω̂tρ_𝐩'⟩


We now need to calculate ⟨ρ_𝐩^*ρ_𝐤^*R_𝐪⟩,
which is given by 

    ⟨ρ_𝐩^*ρ_𝐤^*R_𝐪⟩ =1/√(N)[⟨(ρ_𝐩ρ_𝐤)^*(Ω̂^†ρ_𝐪)⟩ +ω_q⟨(ρ_𝐩ρ_𝐤)^*ρ_𝐪⟩]

The first term in the bracket is 

    ⟨(ρ_𝐩ρ_𝐤)^*(Ω̂^†ρ_𝐪)⟩    =    ⟨(ρ_𝐩ρ_𝐤)^*Ω̂ρ_𝐪⟩
        =    ∑_j⟨(-i𝐩e^i𝐩·𝐫_jρ_𝐤-i𝐤e^i𝐤·𝐫_jρ_𝐩)· D_j(-i𝐪e^-i𝐪·𝐫_j)⟩
        =     -𝐪·𝐩⟨∑_j,lD_je^-i(𝐪-𝐩)·𝐫_j+i𝐤·𝐫_l⟩ -𝐪·𝐤⟨∑_j,lD_je^-i(𝐪-𝐤)·𝐫_j+i𝐩·𝐫_l⟩
        =     -𝐪·𝐩δ_𝐪,𝐤+𝐩[⟨∑_j,lD_je^-i𝐤·𝐫_j+i𝐤·𝐫_l⟩ +⟨∑_j,lD_je^-i𝐩·𝐫_j+i𝐩·𝐫_l⟩]
        =     -ND_0δ_𝐪,𝐤+𝐩[(𝐪·𝐩)S_2(k)+(𝐪·𝐤)S_2(p)]

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(k)
in the fourth equality defined as 

    S_2(k)=1/ND_0⟨∑_j,lD_je^-i𝐤·(𝐫_j-𝐫_l)⟩

in accordance with Eq.(<ref>) in the main text. If D_j
is a constant, such asD_j=D_0 in the τ_p→0
limit, it can be drawn out of the bracket such that S_2(k)=N^-1⟨∑_j,le^-i𝐤·(𝐫_j-𝐫_l)⟩ =S(k)
which is exactly the static structure factor. Nevertheless, in our
present case, D_j depends on the instantaneous configuration
𝐫^N, such that S_2(k) may show abundant
features different from S(k). It is interesting to note
that for large k, D_j seems to be decoupled from the Fourier
components exp(-i𝐤·𝐫_l), and
S_2(k) can be approximated by 

    S_2(k)≃D̅/ND_0⟨∑_j,le^-i𝐤·(𝐫_j-𝐫_l)⟩ =D̅/D_0S(k)

as shown in Fig.<ref> in the main text.

The second term in Eq.(<ref>) can be calculated using
the so-called convolution approximation, which is assumed to be still
appropriate in nonequilibrium situation<cit.>,

    ⟨ρ_𝐩^*ρ_𝐤^*ρ_𝐪⟩≈δ_𝐤+𝐩,𝐪NS(q)S(p)S(k)

Therefore, we can get

    ⟨ρ_𝐩^*ρ_𝐤^*R_𝐪⟩    =     -√(N)D_0δ_𝐪,𝐤+𝐩[𝐪·𝐩S_2(k)+𝐪·𝐤S_2(p)-q^2D̅/D_0S(p)S(k)]

Substituting this into Eq.(<ref>), we obtain

    ⟨ R_𝐪^*e^𝒬Ω̂^ irr𝒬tR_𝐪⟩   ≈1/2∑_ k, p|⟨ρ_𝐩^*ρ_𝐤^*R_𝐪⟩|^2/[N^2S(k)S(p)]^2⟨ρ_𝐩^*e^Ω̂tρ_𝐩⟩⟨ρ_𝐤^*e^Ω̂tρ_𝐤⟩
        =1/2N∑_ k, p|V_𝐪(𝐤,𝐩)|^2⟨ρ_𝐩^*e^Ω̂tρ_𝐩⟩⟨ρ_𝐤^*e^Ω̂tρ_𝐤⟩

where the vortex function V_𝐪(𝐤,𝐩)
is defined as


    V_𝐪(𝐤,𝐩)     =    √(N)⟨ρ_𝐩^*ρ_𝐤^*R_𝐪⟩[N^2S(k)S(p)]^-1
        =     -δ_𝐤+𝐩,𝐪D_0/N{(𝐪·𝐤)S_2(p)/S(p)S(k)+(𝐪·𝐩)S_2(k)/S(p)S(k)-q^2D̅/D_0}
        =     -δ_𝐤+𝐩,𝐪D̅/N{(𝐪·𝐤)[D_0S_2(p)/D̅S(p)S(k)-1]+(𝐪·𝐩)[D_0S_2(k)/D̅S(p)S(k)-1]}

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

    C_2(𝐪;𝐤)=δ_𝐤+𝐩,𝐪/ρ[1-D_0S_2(p)/D̅S(p)S(k)]≡1/ρ[1-D_0S_2(|𝐪-𝐤|)/D̅S(|𝐪-𝐤|)S(k)]

Note that if Eq.(<ref>) becomes a equality, namely,
S_2(k)=D̅S(k)/D_0, then

    C_2( q; k)=1/ρ[1-1/S(k)]=c(k)

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

    V_𝐪(𝐤,𝐩)=ρD̅/N[(𝐪·𝐤)C_2(𝐪;𝐤)+(𝐪·𝐩)C_2(𝐪;𝐤)]

and 

    ⟨ R_𝐪^*e^𝒬Ω̂^ irr𝒬tR_𝐪⟩    =1/2N∑_ k, p|V_𝐪(𝐤,𝐩)|^2⟨ρ_𝐩^*e^Ω̂tρ_𝐩⟩⟨ρ_𝐤^*e^Ω̂tρ_𝐤⟩
       ≈1/2∑_𝐤ρ^2D̅^2/N[(𝐪·𝐤)C_2( q; k)+(𝐪·𝐩)C_2( q; p)]^2F_k(t)F_p(t)

Consequently, we get the irreducible memory function 

    M̃^ irr(q;t)    ≈    -⟨ R_𝐪^*e^𝒬Ω̂^ irr𝒬tR_𝐪⟩⟨ A_-𝐪Ω̂A_𝐪⟩ ^-1
        =    ρ^2D̅/2q^2N∑_𝐤[(𝐪·𝐤)C_2( q; k)+(𝐪·𝐩)C_2( q; p)]^2F_k(t)F_p(t)

Changing to integration by using ∑_𝐤→(2π)^-3V∫ d^3𝐤,
one has

    M̃^ irr(q;t)=ρD̅/16π^3∫ d^3𝐤[(𝐪̂·𝐤)C_2( q; k)+(𝐪̂·𝐩)C_2( q; p)]^2F_k(t)F_p(t)

where 𝐪̂=𝐪/q is the unit vector in the direction
of 𝐪. This is Eq. (<ref>) in the main text.



 §.§ Tagged Particle Dynamics


We now consider tagged particle dynamics. The relevant variable is
A_𝐪=ρ_𝐪^s=e^-i𝐪·𝐫_s
( the subscript or superscript 's' stands for the single particle)
and the self-scattering function reads F_q^s(t)=⟨ρ_-𝐪^se^Ω̂tρ_𝐪^s⟩
with F_q^s(0)=1. The derivation of the memory function
equation for F_q^s(t) are similar to that of F_q(t),
except that some relevant calculations are different. Briefly, we
have 

    M_s^irr(q,t)=-⟨ R_𝐪^s*e^𝒬Ω̂^ irr𝒬tR_𝐪^s⟩⟨ρ_-𝐪^sΩ̂ρ_𝐪^s⟩ ^-1

where 

    R_𝐪^s    =    𝒬(Ω^†ρ_𝐪^s)=(Ω^†ρ_𝐪^s)-𝒫(Ω^†ρ_𝐪^s)
        =     (Ω^†ρ_𝐪^s)-⟨ρ_-𝐪^s(Ω̂^†ρ_𝐪^s)⟩/⟨ρ_-𝐪^sρ_𝐪^s⟩ρ_𝐪^s
        =     (Ω^†ρ_𝐪^s)+q^2D̅ρ_𝐪^s=(Ω^†ρ_𝐪^s)

In the third equality of above equation for R_𝐪^s,
we have used the result ⟨ρ_-𝐪^sΩ̂ρ_𝐪^s⟩ =-q^2⟨ D_s⟩. 

To calculate ⟨ R_𝐪^s*e^𝒬Ω̂^ irr𝒬tR_𝐪^s⟩,
one projects R_𝐪^s onto the product of single and
collective modes ρ_𝐤ρ_𝐩^s with projection
operator

    𝒫_2^s=∑_𝐤,𝐩.ρ_𝐤ρ_𝐩^s⟩⟨(ρ_𝐤ρ_𝐩^s)^*(ρ_𝐤ρ_𝐩^s)⟩ ^-1⟨(ρ_𝐤ρ_𝐩^s)^*.

Then

    ⟨ R_-𝐪^se^𝒬Ω̂^ irr𝒬tR_𝐪^s⟩   ≈   ⟨ R_-𝐪^s𝒫_2^se^𝒬Ω̂^ irr𝒬t𝒫_2^sR_𝐪^s⟩
        =    ∑_𝐤,𝐩∑_𝐤',𝐩'⟨ R_-𝐪^s(ρ_𝐤ρ_𝐩^s)⟩⟨(ρ_𝐤ρ_𝐩^s)^*(ρ_𝐤ρ_𝐩^s)⟩ ^-1
          ×⟨(ρ_𝐤'ρ_𝐩'^s)^*R_𝐪^s⟩⟨(ρ_𝐤'ρ_𝐩'^s)^*(ρ_𝐤'ρ_𝐩'^s)⟩ ^-1⟨(ρ_𝐤ρ_𝐩^s)^*e^𝒬Ω̂^ irr𝒬t(ρ_𝐤'ρ_𝐩'^s)⟩
       ≃   ∑_𝐤,𝐩∑_𝐤',𝐩'⟨ R_-𝐪^s(ρ_𝐤ρ_𝐩^s)⟩/NS(k)⟨(ρ_𝐤'ρ_𝐩'^s)^*R_𝐪^s⟩/NS(k')δ_𝐤𝐤'δ_𝐩𝐩'⟨ρ_𝐤^*e^Ω̂tρ_𝐤'⟩⟨ρ_𝐩^s*e^Ω̂tρ_𝐩'^s⟩
        =    ∑_𝐤,𝐩|V_q^s(𝐤,𝐩)|^2⟨ρ_-𝐤e^Ω̂tρ_𝐤⟩⟨ρ_-𝐩^se^Ω̂tρ_𝐩^s⟩

where


    V_𝐪^s(𝐤,𝐩)=⟨(ρ_𝐤ρ_𝐩^s)^*R_𝐪^s⟩/NS(k)

Now we need to calculate ⟨(ρ_𝐤ρ_𝐩^s)^*R_𝐪^s⟩,
which is given by 

    ⟨(ρ_𝐤ρ_𝐩^s)^*R_𝐪^s⟩ =⟨(ρ_𝐤ρ_𝐩^s)^*(Ω^†ρ_𝐪^s)⟩ +q^2D̅⟨(ρ_𝐤ρ_𝐩^s)^*ρ_𝐪^s⟩


The second term is 

    ⟨(ρ_𝐤ρ_𝐩^s)^*ρ_𝐪^s⟩    =⟨ρ_-𝐤ρ_𝐪-𝐩^s⟩ =δ_𝐪,𝐤+𝐩ρ c(k)S(k)=δ_𝐪,𝐤+𝐩[S(k)-1]

where the δ symbol results from translational invariance and
the result ⟨ρ_𝐤^*ρ_𝐤^s⟩ =ρ c(k)S(k)=S(k)-1<cit.>.

The first term is

    ⟨(ρ_𝐤ρ_𝐩^s)^*(Ω^†ρ_𝐪^s)⟩    =    ⟨(ρ_𝐤ρ_𝐩^s)^*Ω̂ρ_𝐪^s⟩
        =    ∑_j=1^N⟨(-i𝐤(∑_l≠ sδ_jle^i𝐤·𝐫_l)e^i𝐩·𝐫_s-δ_jsi𝐩ρ_-𝐤e^i𝐩·𝐫_s)D_j·(-δ_jsi𝐪e^-i𝐪·𝐫_s)⟩
        =     -⟨(𝐩·𝐪)D_sρ_-𝐤e^i(𝐩-𝐪)·𝐫_s⟩
        =     -δ_𝐪,𝐤+𝐩(𝐩·𝐪)(D_0S_2(k)-D̅)

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


    ⟨(ρ_kρ_p^s)^*R_𝐪^s⟩    =     -δ_𝐪,𝐤+𝐩[𝐩·𝐪(D_0S_2(k)-D̅)]+q^2D̅δ_𝐪,𝐤+𝐩(S(k)-1)
        =    δ_𝐪,𝐤+𝐩D̅[𝐤·𝐪(S(k)-1)+𝐩·𝐪(S(k)-D_0S_2(k)/D̅)]


    V_q^s(𝐤,𝐩)     =    ⟨(ρ_𝐤ρ_𝐩^s)^*R_𝐪^s⟩/NS(k)=δ_𝐪,𝐤+𝐩D̅/N[𝐤·𝐪(1-1/S(k))+𝐩·𝐪(1-D_0S_2(k)/D̅S(k))]
        =    δ_𝐪,𝐤+𝐩ρD̅/N[𝐤·𝐪c(k)+𝐩·𝐪1/ρ(1-D_0S_2(k)/D̅S(k))]

Notice if S_2(k)=D̅S(k)/D_0, V_q(𝐤,𝐩)=δ_𝐪,𝐤+𝐩ρD̅/N[𝐤·𝐪c(k)],
which is the equilibrium result. Next

    ⟨ R_𝐪^s*e^𝒬Ω̂^ irr𝒬tR_𝐪^s⟩    =∑_𝐤,𝐩|V_q^s(𝐤,𝐩)|^2⟨ρ_-𝐤e^Ω̂tρ_𝐤⟩⟨ρ_-𝐩^se^Ω̂tρ_𝐩^s⟩
       ≈∑_𝐤ρ^2D̅^2/N[𝐤·𝐪c(k)+𝐩·𝐪1/ρ(1-D_0S_2(k)/D̅S(k))]^2F_k(t)F_p^s(t)

and

    M̃_s^ irr(q;t)     =     -⟨ R_𝐪^s*e^𝒬Ω̂^ irr𝒬tR_𝐪^s⟩⟨ρ_-𝐪^sΩ̂ρ_𝐪^s⟩ ^-1
       ≈   ρ^2D̅/q^2N∑_𝐤[𝐤·𝐪c(k)+𝐩·𝐪1/ρ(1-D_0S_2(k)/D̅S(k))]^2F_k(t)F_p^s(t)
        =    ρD̅/(2π)^3∫ d^3𝐤[𝐤·𝐪̂c(k)+𝐩·𝐪̂1/ρ(1-D_0S_2(k)/D̅S(k))]^2F_k(t)F_p^s(t)

Finally, the memory function equation for the self-scattering function
F_q^s(t) is given by

    ∂/∂ tF_q^s(t)+ω_q^sF_q^s(t)+∫_0^tduM_s^ irr(q,t-u)∂/∂ uF_q^s(u)=0

where 

    ω_q^s=-⟨ρ_-𝐪^s(Ω̂^†ρ_𝐪^s)⟩/⟨ρ_-𝐪^sρ_𝐪^s⟩=q^2D̅.


