Chinese Journal of Chemical Physics  2018, Vol. 31 Issue (4): 584-594

The article information

Meng-kai Feng, Zhong-huai Hou

Mode-Coupling Theory for Glass Transition of Active-Passive Binary Mixture

Chinese Journal of Chemical Physics, 2018, 31(4): 584-594

http://dx.doi.org/10.1063/1674-0068/31/cjcp1806148

Article history

Accepted on: July 22, 2018
Mode-Coupling Theory for Glass Transition of Active-Passive Binary Mixture
Meng-kai Feng, Zhong-huai Hou
Dated: Received on June 21, 2018; Accepted on July 22, 2018
Department of Chemical Physics & Hefei National Laboratory for Physical Sciences at the Microscales, iChEM, University of Science and Technology of China, Hefei 230026, China
*Author to whom correspondence should be addressed. Zhong-huai Hou, E-mail:hzhlj@ustc.edu.cn
Abstract: Collective behaviours of active particle systems have gained great research attentions in recent years. Here we present a mode-coupling theory (MCT) framework to study the glass transition of a mixture system of active and passive Brownian particles. The starting point is an effective Smoluchowski equation, which governs the dynamics of the probability distribution function in the position phase space. With the assumption of the existence of a nonequilibrium steady state, we are able to obtain dynamic equations for the intermediate scattering functions (ISFs), wherein an irreducible memory function is introduced which in turn can be written as functions of the ISFs based on standard mode-coupling approximations. The effect of particle activity is included through an effective diffusion coefficient which can be obtained via short time simulations. By calculating the long-time limit of the ISF, the Debye-Waller (DW) factor, one can determine the critical packing fraction ηc of glass transition. We find that for active-passive (AP) mixtures with the same particle sizes, ηc increases as the partial fraction of active particle xA increases, which is in agreement with previous simulation works. For system with different active/passive particle sizes, we find an interesting reentrance behaviour of glass transition, i.e., ηc shows a non-monotonic dependence on xA. In addition, such a reentrance behaviour would disappear if the particle activity is large enough. Our results thus provide a useful theoretical scheme to study glass transition behaviour of active-passive mixture systems in a promising way.
Key words: Glass transition    Mode-coupling theory    Active particles    Active-passive particle mixture
Ⅰ. INTRODUCTION

Understanding the collective behaviours of active particles has gained extensive research interest due to its importance for both fundamental physics and biological application [1-3]. In recent years, a number of active particle systems have been studied both experimentally and theoretically (including simulation), and a lot of interesting collective behaviours have been reported, such as large scale vertex formation, active swarming [4, 5], phase separation [3, 6-10], and glass transition [11-13] etc. Specifically, mixture system consisting of active and passive particles has become a hot field recently [14-18]. For instance, Cates et al. [18] studied the activity induced phase separation and self-assembly of active-passive mixtures, Wittmann et al. [16] studied the effective interactions between different active particles driven by coloured noises, Weber et al. [17] studied the demixing generated by different diffusivities in hot-cold mixture, and so on. Very recently, glassy dynamics and glass transition in dense active particle systems have also gained lots of attention [11, 19-22]. Experiments show that active particle systems exist many dynamic slowing down processes such as jamming and dynamic arrest. For example, amorphous solidification process has been found in collective motion of a cellular monolayer [23], and glassy behaviours have been found for ant aggregates in large scales [24]. Using extensive simulation, Ran et al. [11] found that the activity of particles pushes the glass transition critical point to a higher density until the random close packing. Besides, there have also been some purely theoretical studies about glassy dynamics of active particle systems. Voigtmann et al. [25] presented a mode-coupling theory for the slow dynamics of two-dimensional spherical active Brownian particles. Berthier and Kurchan [19] studied the nonequilibrium glass transitions in driven and active matter. Szamel et al. [12] studied the glassy dynamics of athermal self-propelled particles, both in computer simulations and in a nonequilibrium microscopic MCT-like theory. Nevertheless, theoretical study for glassy dynamics of mixture systems including active particles is still lacking, to the best of our knowledge.

In our recent work [26], we have developed a mode-coupling theory framework to study the glassy dynamics of a type of active particle systems, the so-called active particles driven by an Ornstein-Uhlenbeck (OU) process with thermal noise included (AOU-T). Our starting point was an effective Smoluchowski equation (SE) governing the dynamics of the probability distribution function in the position phase space. We were able to obtain a dynamic equation for the intermediate scattering function (ISF) of the system, from which one can calculate the Debye-Waller (DW) factor, the value of ISF in the long time limit. Glass transition takes place when the DW factor reaches a non-zero value upon change of the volume fraction or temperature. Our results there well reproduced the simulation results reported in the literature that activity induced a shift of the critical density for glass transition to higher values [11]. However, the theory was only for a monocomponent system and how it can be applied to mixture systems is still open.

Motivated by this, in the present work, we have extended our theory to a mixture binary system of active-passive particles. The active component is described by the AOU-T model as in our previous work. Following Fox's method of dealing with coloured noises, we have obtained an effective SE for the probability distribution function $\Psi\left(\boldsymbol{\textbf{r}}^{N}, t\right)$, the probability density that the system is at a specific configuration $\mathbf{\textbf{r}}^{N}$=$\left(\mathbf{r}_{1}, \mathbf{r}_{2}, \dots, \mathbf{r}_{N}\right)$ at time $t$. By using standard Zwanzig-Mori projection operator method, we can obtain the time evolution equations for the ISFs $F_{\alpha\beta}\left(q, t\right)$, where $\alpha, \beta$ stands for active or passive component and $q$ is a modulus of vector in the Fourier space. With these equations, we are able to calculate the DW factor and determine how the glass transition depends on the system parameters, such as the volume fraction, the particle activity, as well as the number fraction of the active component $x_{\textrm{A}}$. For system with equal active/passive particle sizes, we find that $\eta_{\textrm{c}}$ shifts to larger values with the increase of $x_{\textrm{A}}$ or particle activity, in accordance with previous simulation studies. For a system with different particle sizes, an interesting reentrance behaviour of glass transition is observed, i.e., $\eta_{\textrm{c}}$ shows non-monotonic dependence on $x_{\textrm{A}}$ if particle activity is not large. If the particle activity is large enough, however, such a reentrance behaviour is suppressed and $\eta_{\textrm{c}}$ increases monotonically with $x_{\textrm{A}}$ as in the case of equal sizes. These results indicate some interesting interplay between particle activity and mixing effect in such complex systems.

Ⅱ. MODEL AND THEORY A. Equations of motion

Consider a mixture system with total particle number $N$ and volume $V$, consisting of active particles with number fraction $x_{\textrm{A}}$=$N_{\textrm{A}}/N$ and passive particles with number fraction $x_{\textrm{P}}$=1-$x_{\textrm{A}}$=$N_{\textrm{P}}/N$. The interacting forces among all particles are pairwise and purely repulsive, namely, the Weeks-Chandler-Andersen (WCA) potential

 $\begin{eqnarray} u_{\alpha\beta}(r)\hspace{-0.1cm}=\hspace{-0.1cm}\begin{cases} 4\epsilon_{\alpha\beta}\left[\left(\displaystyle{\frac{\sigma_{\alpha\beta}}{r}}\right)^{12}\hspace{-0.1cm}-\hspace{-0.1cm}\left(\displaystyle{\frac{\sigma_{\alpha\beta}}{r}}\right)^{6}\right]\hspace{-0.05cm}+\hspace{-0.05cm}\epsilon_{\alpha\beta}, &\hspace{-0.15cm}r\hspace{-0.05cm}\le\hspace{-0.05cm}2^{1/6}\sigma_{\alpha\beta} \\ 0, &\hspace{-0.15cm} r\hspace{-0.05cm}>\hspace{-0.05cm}2^{1/6}\sigma_{\alpha\beta} \end{cases} \end{eqnarray}$

where $\alpha, \beta$=$\textrm{A}\text{ or }\textrm{P}$ (denoting active or passive particle, respectively), $\epsilon_{\alpha\beta}$ is the strength of potential, $\sigma_{\alpha\beta}$ is the interacting scale between $\alpha$ and $\beta$ particles.

In general forms, the dynamic equations of motion for the particles read

 $\begin{eqnarray} \dot{\boldsymbol{r}}_{i}^{\alpha}\left(t\right)=\gamma_{\alpha}^{-1}\left[\boldsymbol{f}_{i}^{\alpha}\left(t\right)+\boldsymbol{F}_{i}^{\alpha}\right]+\boldsymbol{\xi}_{i}^{\alpha}\left(t\right)\label{eq:Langevin} \end{eqnarray}$

where $\boldsymbol{r}_{i}^{\alpha}$ is the position vector of $i$-th particle of type-$\alpha$, $\boldsymbol{F}_{i}^{\alpha}$=$\displaystyle{-\sum_{\beta}\sum_{j=1, j\neq i}^{N_{\beta}}\nabla_{i}u_{\alpha\beta}\left(r_{ij}\right)}$ is the configurational force, $\boldsymbol{f}_{i}^{\alpha}$ is the self-propulsion force which only exerts on active particles (for simplicity we write $\boldsymbol{f}_{i}^{\textrm{A}}$=$\boldsymbol{f}_{i}$). We assume $\mathit{\boldsymbol{f}}_{i}$ is governed by an OU-process,

 $\begin{eqnarray} \dot{\boldsymbol{f}_{i}}\left(t\right)=-\tau_{\textrm{p}}^{-1}\boldsymbol{f}_{i}\left(t\right)+\mathit{\boldsymbol{\eta}}{}_{i}(t)\label{eq:OU_process} \end{eqnarray}$ (2)

with persistent time $\tau_{\textrm{p}}$ and a white noise term $\mathit{\boldsymbol{\eta}}_{i}$, obeying $\left\langle \mathit{\boldsymbol{\eta}}_{i}(t)\mathit{\boldsymbol{\eta}}_{j}(t')\right\rangle _{{\rm noise}}$=$2D_{f}\delta_{ij}\delta$($t$-$t'$)$\boldsymbol{1}$, where $\left\langle \cdots\right\rangle _{{\rm noise}}$ denotes averaging over noise distribution, and $\boldsymbol{1}$ denotes the unit tensor. Besides, $\mathit{\boldsymbol{\xi}}_{i}^{\alpha}\left(t\right)$ denotes the thermal noise with zero mean and white noise correlation

 $\begin{eqnarray} \left\langle \mathit{\boldsymbol{\xi}}_{i}^{\alpha}\left(t\right)\mathit{\boldsymbol{\xi}}_{j}^{\beta}\left(t'\right)\right\rangle _{{\rm noise}}=2D_{t}^{\alpha}\mathit{\boldsymbol{1}}\delta_{\alpha\beta}\delta_{ij}\delta\left(t-t'\right) \label{eq:thermal_noise} \end{eqnarray}$ (3)

where $D_{t}^{\alpha}$ is the translational diffusion coefficient satisfying fluctuation-dissipation theorem $D_{t}^{\alpha}\gamma_{\alpha}$=$k_{\textrm{B}}T$, with $\gamma_{\alpha}$ being the friction constant for ${{\alpha}}$-type particles. For long time, the motion of isolated active particle is diffusive, with a long time diffusion coefficient $D_{0}$=$D_{t}^{\textrm{A}}+D_{f}\tau_{\textrm{p}}^{2}/\gamma_{\textrm{A}}^{2}$. Since larger activity leads to larger diffusivity, and $D_{t}^{\textrm{A}}$ and $\gamma_{\textrm{A}}$ are constants here, $D_{f}\tau_{\textrm{p}}^{2}$ is a reasonable parameter to describe the activity.

In this work, we use dimensionless units by setting $k_{\textrm{B}}T$=1, $\sigma_{\textrm{AA}}$=1 and $D_{t}^{\textrm{A}}$=1. Time unit is then $\sigma_{\textrm{AA}}^{2}/D_{t}^{A}$. For potential parameters, we set $\epsilon_{\textrm{AA}}$=$\epsilon_{\textrm{AP}}$=$\epsilon_{\textrm{PP}}$=1, and assume $2\sigma_{\textrm{AP}}$=$\sigma_{\textrm{AA}}$+$\sigma_{\textrm{PP}}$. For active particles, we use $\tau_{\textrm{p}}$=0.02, and change $D_{f}$ to adjust the activity. Simulations are performed in a periodic boundary three dimensional box with 4096 particles and volume fraction is changed via change of the box size. Initial conditions are prepared with random particle positions and random self-propulsion forces satisfying the steady-state distribution of OU process, and the simulation method is standard Brownian dynamics algorithm.

B. Effective Smoluchowski equation

One notes that the Langevin equation (Eq.(1)) for active particles is non-Markovian due to the coloured noise term $\boldsymbol{f}_{i}$, such that it is not possible to derive an exact SE (which is actually Fokker-Planck equation) for the time evolution of the probability distribution function $\Psi\left(\mathbf{r}^{N}, t\right)$. Nevertheless, it is possible to obtain an approximate SE for $\Psi\left(\mathbf{r}^{N}, t\right)$, as already shown in many literatures. For instance, Fox had introduced a method using functional calculus which allows one to obtain a nearly accurate SE for short correlation times [27, 28], even being valid when the correlation time is long [29]. Using this method, Brader et al. [29] had presented a theory to study the effective interactions among active Brownian particles which facilitated them to understand activity-induced phase separation in a novel way. In our recent work [26], we have also used the same method to obtain the SE for a single-component active particle system described by the AOU-T model.

In the present work, we have applied similar approach to the mixture system described by Eq.(1). Details are given in Appendix A. Consequently, the effective SE reads,

 $\begin{eqnarray} \partial_{t}\Psi\left(\boldsymbol{\textbf{r}}^{N}, t\right)=\Omega\Psi=\left(\Omega_{\textrm{A}}+\Omega_{\textrm{P}}\right)\Psi\left(\boldsymbol{\textbf{r}}^{N}, t\right) \end{eqnarray}$ (4)

wherein the Smoluchowski operator $\Omega$ can be divided into two additive parts, an active part $\Omega_{\textrm{A}}$ and a passive one $\Omega_{\textrm{P}}$, given by

 $\begin{eqnarray} \Omega_{\textrm{P}}\hspace{-0.15cm}&=&\hspace{-0.15cm}\sum\limits_{i=1}^{N_{\textrm{P}}}\nabla_{i}\cdot D_{t}^{\textrm{P}}\left[\nabla_{i}-\beta\boldsymbol{F}_{i}^{\textrm{P}}\left(\boldsymbol{\textbf{r}}^{N}\right)\right]\label{eq:PassSO} \end{eqnarray}$ (5)
 $\begin{eqnarray} \Omega_{\textrm{A}} \hspace{-0.15cm}& = &\hspace{-0.15cm} \sum\limits_{j=1}^{N_{\textrm{A}}}\nabla_{j}\cdot D_{j}^{\textrm{A}}\left[\nabla_{j}-\beta\boldsymbol{F}_{j}^{{\rm eff}}\left(\boldsymbol{\textbf{r}}^{N}\right)\right]\label{eq:ActSO} \end{eqnarray}$ (6)

where

 $\begin{eqnarray} D_{j}^{\textrm{A}}=D_{t}^{\textrm{A}}+\left[\frac{D_{f}\tau_{\textrm{p}}^{2}/\gamma^{2}}{1-\tau_{\textrm{p}}\beta D_{t}^{\textrm{A}}\nabla_{j}\cdot\boldsymbol{F}'_{j}\left(\mathbf{r}^{N}\right)}\right]\label{eq:InsDiffCoe} \end{eqnarray}$ (7)

is a configuration-dependent instantaneous diffusion coefficient of active particles, and $\boldsymbol{F}_{j}^{{\rm eff}}$ denotes an instantaneous effective force satisfying

 $\begin{eqnarray} \beta D_{j}^{\textrm{A}}\boldsymbol{F}_{j}^{{\rm eff}}=\beta D_{t}^{\textrm{A}}\boldsymbol{F}_{j}^{\textrm{A}}-\nabla_{i}\frac{D_{f}\tau_{\textrm{p}}^{2}/\gamma^{2}}{1-\tau_{\textrm{p}}\beta D_{t}^{\textrm{A}}\nabla_{j}\cdot\boldsymbol{F}'_{j}\left(\mathbf{r}^{N}\right)} \end{eqnarray}$ (8)

Herein, $\boldsymbol{F}'_{j}$ denotes the force exerted on particle $j$ by all other active particles, i.e., $\boldsymbol{F}'_{j}$=$\displaystyle{\sum_{i=1, j\neq i}^{N_{\textrm{A}}}} \boldsymbol{F}_{ij}$= $-\nabla_{j}\displaystyle{\sum_{i=1, j\neq i}^{N_{\textrm{A}}}}u_{\textrm{AA}}(r_{ij})$. In the long time limit, we assume that the system will reach a nonequilibrium steady-state $P_{\textrm{s}}\left(\mathbf{r}^{N}\right)$, satisfying

 $\begin{eqnarray} \Omega P_{\textrm{s}}\left(\mathbf{r}^{N}\right)=-\sum\limits_{\alpha}\sum\limits_{i=1}^{N_{\alpha}}\nabla_{i, \alpha}\cdot\mathbf{J}_{i, \alpha}^{s}=0\label{eq:Omega_Pss} \end{eqnarray}$ (9)

where $\mathbf{J}_{i, \alpha}^{\textrm{s}}$ denotes the probability density current given by

 $\begin{eqnarray} \mathbf{J}_{i, \textrm{A}}^{\textrm{s}}\hspace{-0.15cm}&=&\hspace{-0.15cm}-D_{i}^{\textrm{A}}\left(\mathbf{r}^{N}\right)\left[\nabla_{i}-\beta\boldsymbol{F}_{i}^{{\rm eff}}\left(\boldsymbol{\textbf{r}}^{N}\right)\right]P_{\textrm{s}}\left(\boldsymbol{\textbf{r}}^{N}\right) \end{eqnarray}$ (10)
 $\begin{eqnarray} \mathbf{J}_{i, \textrm{P}}^{\textrm{s}}\hspace{-0.15cm}&=&\hspace{-0.15cm}-D_{i}^{\textrm{P}}\left(\mathbf{r}^{N}\right)\left[\nabla_{i}-\beta\boldsymbol{F}_{i}^{P}\left(\boldsymbol{\textbf{r}}^{N}\right)\right]P_{\textrm{s}}\left(\boldsymbol{\textbf{r}}^{N}\right) \end{eqnarray}$ (11)

To proceed, here we assume that in the steady state all the currents vanish [27], i.e., $\mathbf{J}_{i, \alpha}^{\textrm{s}}$=0. As argued in Ref.[30], there was no clear evidence of non-vanishing system current in dense active particle systems close to glass transition. With this assumption, we then have $\nabla_{j}\ln P_{\textrm{s}}\left(\mathbf{\textbf{r}}^{N}\right)$=$\beta\boldsymbol{F}_{j}^{{\rm eff}}$ for active particles and $\nabla_{i}\ln P_{\textrm{s}}\left(\mathbf{\textbf{r}}^{N}\right)$=$\beta\boldsymbol{F}_{i}^{\textrm{P}}$ for passive particles.

To study the collective dynamic behaviours of the system, we define a matrix of collective intermediate scattering functions (ISFs) $\mathbb{S}(q, t)$ with each element

 $\begin{eqnarray} S_{\alpha\beta}\left(q, t\right)=\left\langle \rho_{-\boldsymbol{q}}^{\alpha}\textrm{e}^{\Omega t}\rho_{\boldsymbol{q}}^{\beta}\right\rangle \label{eq:ISF} \end{eqnarray}$ (12)

where $\left\langle \cdots\right\rangle =\int \textrm{d}\mathbf{r}^{N}\cdots P_{\textrm{s}}\left(\mathbf{r}^{N}\right)$ means the steady-state average, and

 $\begin{eqnarray} \rho_{\boldsymbol{q}}^{\alpha}=\frac{1}{\sqrt{N_{\alpha}}}\sum\limits_{j=1}^{N_{\alpha}}\textrm{e}^{-i\boldsymbol{q}\cdot\boldsymbol{r}_{j}^{\alpha}}\label{eq:DenFlu} \end{eqnarray}$ (13)

is the Fourier transform of density fluctuation $\rho^{\alpha}\left(\mathbf{r}\right)$=$\displaystyle{\sum_{j=1}^{N_{\alpha}}}\delta\left(\mathbf{r}-\mathbf{r}_{j}^{\alpha}\right)$. At $t$=0, $S_{\alpha\beta}\left(q, t\right)$ is the non-equilibrium static structure factor $\mathbb{S}(q)$ with elements

 $\begin{eqnarray} S_{\alpha\beta}(q)=\left\langle \rho_{-\boldsymbol{q}}^{\alpha}\rho_{\boldsymbol{q}}^{\beta}\right\rangle \label{eq:NE_Sq} \end{eqnarray}$ (14)

Clearly, one has $S_{\alpha\beta}\left(q\right)$=$S_{\beta\alpha}\left(q\right)$. Note that for the non-equilibrium system studied here, it is not feasible to calculate $\mathbb{S}\left(q\right)$ through theoretical methods like the Ornstein-Zernike (OZ) equations, while direct simulations must be used to calculate them in the steady state.

C. Memory function equation

To study glassy dynamics, one needs to study the long time relaxation behaviour of $S_{\alpha\beta}\left(q, t\right)$. Generally, $S_{\alpha\beta}\left(q, t\right)$ will relax to zero in the long time limit for a liquid system, while it will reach an apparent non-zero value for a glass state. Following standard projection operator procedure as shown in detail in Appendix B, we can obtain the dynamic evolution equation for the matrix $\mathbb{S}(q, t)$, given by

 $\begin{eqnarray} \hspace{-0.3cm}&&\partial_{t}\mathbb{S}\left(q, t\right)\hspace{-0.1cm}+\hspace{-0.1cm} \mathbb{W}\left(q\right)\mathbb{S}^{-1}\left(q\right)\mathbb{S}\left(q, t\right) \hspace{-0.1cm} \\ &+&\hspace{-0.1cm}\int_{0}^{t}\textrm{d}u\mathbb{M}^{{\rm irr}}\left(q, u\right)\partial_{u}\mathbb{S}\left(q, u\right)=0\label{eq:GLE} \end{eqnarray}$ (15)

In this equation, $\mathbb{W}\left(q\right)$ denotes a frequency term given by

 $\begin{eqnarray} W_{\alpha\beta}\left(q\right)\hspace{-0.1cm}&=&\hspace{-0.1cm}-\left\langle \rho_{-\boldsymbol{q}}^{\alpha}\Omega\rho_{\boldsymbol{q}}^{\beta}\right\rangle \\ &=&\hspace{-0.1cm}\frac{\delta_{\alpha\beta}q^{2}}{N_{\alpha}}\sum\limits_{j=1}^{N_{\alpha}}\left\langle D_{j}^{\alpha}\right\rangle \\ &\equiv&\hspace{-0.1cm}\delta_{\alpha\beta}q^{2}\bar{D}^{\alpha}\label{eq:Freq} \end{eqnarray}$ (16)

where $D_{j}^{\alpha}$ is given by Eq.(7) and hence

 $\begin{eqnarray} \bar{D}^{\textrm{A}}=N_{\textrm{A}}^{-1}\sum\limits_{j=1}^{N_{\textrm{A}}}\left\langle D_{j}^{\textrm{A}}\right\rangle \label{eq:D-bar} \end{eqnarray}$ (17)

denotes an averaged single-particle diffusion coefficient of active particles, and $\bar{D}^{\textrm{P}}$=$D_{t}^{\textrm{P}}$.

$\mathbb{M}^{{\rm irr}}\left(q, u\right)$ is the so-called irreducible memory function, given by

 $\begin{eqnarray} M_{\alpha\beta}^{{\rm irr}}\left(q, t\right) \frac{1}{2}\sum\limits_{\boldsymbol{k}}\sum\limits_{\nu}\frac{\delta_{\nu\beta}}{q^{2}\bar{D}^{\nu}}\sum\limits_{\alpha'\beta'\gamma\delta}V_{\alpha}^{\gamma\delta}\left(\boldsymbol{q};\boldsymbol{k}\right)\cdot \\ S_{\gamma\alpha'}\left(k, t\right)S_{\delta\beta'}\left(\left|\boldsymbol{q}-\boldsymbol{k}\right|, t\right)V_{\nu}^{\alpha'\beta'}\left(\boldsymbol{q};\boldsymbol{k}\right)\label{eq:M_irr} \end{eqnarray}$ (18)

where $\boldsymbol{p}$=$\boldsymbol{q}$-$\boldsymbol{k}$, $q$=$\left|\boldsymbol{q}\right|$, and the vertex term reads

 $\begin{array}{l} V_\alpha ^{\gamma \delta }\left( {\mathit{\boldsymbol{q;k}}} \right)\frac{{{{\bar D}^\alpha }}}{{\sqrt {{N_\alpha }} }}\left[ {\mathit{\boldsymbol{q}} \cdot \mathit{\boldsymbol{k}}{\delta _{\alpha \delta }}{C_{\alpha \gamma }}\left( k \right) + } \right.\\ \mathit{\boldsymbol{q}} \cdot \left( {\mathit{\boldsymbol{q - k}}} \right){\delta _{\alpha \gamma }}{C_{\alpha \delta }}\left( {\left| {\mathit{\boldsymbol{q - k}}} \right|} \right) \end{array}$ (19)

with $\mathbb{S}^{-1}(k)$=$1-\mathbb{C}(k)$ wherein $C_{\alpha\beta}\left(k\right)$=$\sqrt{\rho_{\alpha}\rho_{\beta}}c_{\alpha\beta}\left(k\right)$ is the direct correlation function, $\rho_\alpha$=$N_\alpha$/$V$ denotes the number density of $\alpha$-species.

Since the goal of the present work is to determine whether the system is at liquid or glass state, one can just calculate the infinite-time limit of $\mathbb{S}\left(q, t\right)$, the so-called DW factor

 $\begin{eqnarray} \boldsymbol{f}\left(q\right)=\lim\limits_{t\rightarrow\infty}\mathbb{S}\left(q, t\right)\label{eq:DW_def} \end{eqnarray}$ (20)

since a non-zero value of $\boldsymbol{f}\left(q\right)$ indicates that the system is in glass state. According to Eq.(15), we can have

 $\begin{eqnarray} \hspace{-0.1cm}\mathbb{W}\left(q\right)\mathbb{S}^{-1}\left(q\right)\boldsymbol{f}\left(q\right)\hspace{-0.1cm}+\hspace{-0.1cm}\mathbb{M}^{\textrm{irr}}\left(q, \infty\right)\left[\boldsymbol{f}\left(q\right)-\mathbb{S}\left(q\right)\right]\hspace{-0.1cm}=\hspace{-0.1cm}0\label{eq:DW_equ} \end{eqnarray}$ (21)

where

 $\begin{eqnarray} M_{\alpha\beta}^{{\rm irr}}\left(q, \infty\right) \hspace{-0.15cm}& \approx&\hspace{-0.15cm} \frac{1}{2}\sum\limits_{\boldsymbol{k}}\sum\limits_{\nu}\frac{\delta_{\nu\beta}}{q^{2}\bar{D}^{\nu}}\sum\limits_{\alpha'\beta'\gamma\delta}V_{\alpha}^{\gamma\delta}\left(\boldsymbol{q};\boldsymbol{k}\right)\cdot \\ &&S_{\gamma\alpha'}\left(k, t\right)S_{\delta\beta'}\left(\left|\boldsymbol{q}-\boldsymbol{k}\right|, t\right)V_{\nu}^{\alpha'\beta'}\left(\boldsymbol{q};\boldsymbol{k}\right) \\ \label{eq:M_inf} \end{eqnarray}$ (22)

These two equations can be solved by iteration method using $\mathbb{S}\left(q\right)$ as the initial value of $\boldsymbol{f}(q)$.

The above equations (Eqs.(15)-(19)) allow us to study the glassy dynamics of the system. Usually, one can perform numerical calculations to obtain the intermediate scattering functions $\mathbb{S}(q, t)$, given the nonequilibrium steady-state values of $\bar{D}^{\alpha}$ and the nonequilirium structure factors $\mathbb{S}\left(q\right)$ as inputs, as we have done in our previous work about MCT for purely active systems [27]. For a single component passive system, the ISF close to glass transition often shows a two-step relaxation. However, for a single component active system, simulations showed weak or no two-step relaxation [11], which was also reproduced in our previous MCT work [27]. For the mixture system considered in the present work, surely it would be an interesting and important topic to investigate the relaxation dynamics of $S_{\alpha\beta}\left(q, t\right)$. However, as a first step, the main purpose of the present work is to study the glass transition point. This can be achieved by calculating the long-time value of ISFs, i.e., the DW factors, in a much easier way by solving Eq.(21) iteratively, rather than numerically solving Eq.(15) which is technically much harder.

Ⅲ. RESULTS AND DISCUSSION A. Structure information

As stated above, one needs to calculate $\boldsymbol{f}(q)$ to study the glass transition behaviour of our mixture system. According to Eq.(21) and Eq.(22), one must use $\bar{D}^{\textrm{A}}$ and $\mathbb{S}(q)$ as inputs. Since the system is in nonequilibrium, the structure factor $\mathbb{S}(q)$ cannot be obtained through OZ equations, and theoretical calculations of $\bar{D}^{A}$ is not available either. Therefore, one must use direct simulations to obtain these parameters in the steady state. Combining these relatively cheaper simulations with our MCT equations above, one can feasibly study the very long time glassy dynamics of the system, which is usually very expensive to study directly by simulations.

We are interested in how the system's dynamics change with the volume fraction $\eta$ of the system. For a single-component hard sphere system, it is simply defined as $\eta$=$\displaystyle{\frac{\pi}{6}\rho d^{3}}$ where $\rho$=$N/V$ is the number density and $d$ is the diameter of the particle. In our present work, however, the particle is not hard sphere and the definition of packing fraction for such a soft potential system is not straightforward. In addition, we are considering a mixture system wherein active or passive particles could have different parameters $\sigma_{\textrm{AA}}$ or $\sigma_{\textrm{PP}}$. Generally, one may write $\eta$=$\left(N_{\textrm{A}}V_{\textrm{A}}+N_{\textrm{P}}V_{\textrm{P}}\right)/V$=$NV_{\textrm{A}}\left(x_{\textrm{A}}+x_{\textrm{P}}\theta\right)/V$, wherein $V_{\textrm{A, P}}$ denotes the effective volume of an active (passive) particle and $\theta$=$V_{\textrm{P}}/V_{\textrm{A}}$ is just the volume ratio. In our present work, we set $V_{\textrm{A}}$=$\displaystyle{\frac{\pi}{6}}\left(d_{\text{eff}}^{\textrm{A}}\right)^{3}$ where $d_{{\rm eff}}^{\textrm{A}}$ denotes an effective diameter given by $\displaystyle{\int_{0}^{\infty}}\left(1-\textrm{e}^{-\beta u_{\textrm{AA}}(r)}\right)\textrm{d}r$=$d_{{\rm eff}}^{\textrm{A}}$ (note that $d_{{\rm eff}}^{\textrm{A}}$ recovers the particle diameter for a hard sphere potential [30]). To determine $\theta$, however, we require that a binary {passive} mixture system would have same critical packing fraction $\eta_{\textrm{c}}$ of glass transition for either $x_{\textrm{A}}$=1 (however with zero activity) or $x_{\textrm{P}}$=1 if the particles are of different sizes. For instance, for parameters $\sigma_{\textrm{PP}}/\sigma_{\textrm{AA}}$=0.8, our data shows that $\theta$$\approx0.6733. Surely for \sigma_{\textrm{PP}}/\sigma_{\textrm{AA}}=1.0, \theta is just one. In FIG. 1(a), we show the partial structure factor \mathbb{S}(q) for different volume fraction \eta as examples, with constant particle activity 3.0 and fraction of active component x_{\textrm{A}}=0.2. The diagonal elements S_{\textrm{AA}}(q) and S_{\textrm{PP}}(q) are shown in black and blue lines, respectively. The peak heights are positively correlated with the fraction x_{\textrm{A}} or x_{\textrm{P}} of corresponding component. Since x_{\textrm{A}}=0.2, the peaks for S_{\textrm{PP}}(q) (blue lines) are much higher than those for S_{\textrm{AA}}(q) (black lines). In addition, we have also shown the non-diagonal elements S_{\textrm{AP}}\left(q\right) (red lines), which are approximately zero at large q-values and can even take negative values. The insets show the zoom-in of the peaks, wherein the peak heights for all three elements increase with the increment of packing fraction \eta.  FIG. 1 The dependence of partial structure factor S_{\alpha\beta}(q) and averaged diffusion coefficient \bar{D}^{\textrm{A}} on the effective packing fraction \eta, and activity D_{f}\tau_{\textrm{p}}^{2}. All of the information is gained from simulations, with x_{\textrm{A}}=0.2, \tau_{\textrm{p}}=0.02, and the same particle size. (a) shows the maximum {peak} of partial structure factor rises with the density of particles, with the constant activity; the insets show the zoom-out of the peaks. (b) \bar{D}^{\textrm{A}} monotonically decreases with effective packing fraction \eta. FIG. 1(b) shows examples for another input parameter \bar{D}^{\textrm{A}} as functions of the packing fraction \eta, for a few different values of particle activity D_{f}\tau_{\textrm{p}}^{2}=1.0, 3.0 and 10.0, respectively, with fixed x_{\textrm{A}}=0.2. Generally, larger activity leads to larger \bar{D}^{\textrm{A}}. Besides, with constant activity, \bar{D}^{\textrm{A}} slightly decreases with increasing packing fraction due to the shrinking of average distance between particles, similar to those found in active particles systems [26]. B. Debye-Waller factor Using the iteration method mentioned above, we get the DW factors \mathit{\boldsymbol{f}}\left(q\right) as shown in FIG. 2(a), with the same parameters as in FIG. 1(a). Blue and black lines are diagonal elements f_{PP}(q) and f_{\textrm{AA}}(q), respectively. At large qs, all f_{\alpha\beta}\left(q\right) are close to zero, indicating that small scale fluctuations relax approximately to zero. For a relatively low packing fraction, say \eta=0.6, all elements f_{\alpha\beta}\left(q\right) are zero for all qs, corresponding to a liquid state. With slightly increasing \eta to 0.601, however, clear structures show up for all elements of \mathit{\boldsymbol{f}}\left(q\right), with peaks at certain values of q. The apparent non-zero values of these DW factors clearly demonstrate that the system is in glassy state. With further increase of \eta to a larger value 0.614, the values of \mathit{\boldsymbol{f}}$$\left(q\right)$ also increase. Obviously, the system undergoes a sharp glass transition when the packing fraction $\eta$ changes from 0.6 to 0.601.

 FIG. 2 The dependence of Debye-Waller factor on packing fractions and activities. (a) shows every element of $\boldsymbol{f}(q)$, all of them increases with the packing fraction. Select the maximum value of $f_{\textrm{AA}}\left(q\right)$ and draw figure (b), glass phase transitions correspond to the step of $f_{\textrm{AA}}\left(q\right)_{\textrm{max}}$ from 0 to a finite value. Parameters set as $x_{\textrm{A}}$=0.2, $\sigma_{\textrm{PP}}/\sigma_{\textrm{AA}}$=1.0.

In FIG. 2(b), we depict the first peak value of $f_{\textrm{AA}}\left(q\right)$ as functions of $\eta$ for three different values of particle activity, $D_{f}\tau_{\textrm{p}}^{2}$=1.0, 3.0 and 10.0, respectively. All three curves show abrupt jumps at a certain critical value $\eta_{\textrm{c}}$, corresponding to the glass transition point. With increasing particle activity, one can see that the critical value $\eta_{\textrm{c}}$ shifts to higher values, indicating that particle activity hinders glass transition. This is in qualitatively accordance with the simulation results reported in Ran's work [11].

C. Glass transition of active-passive mixture

In the above subsection, we show the method of determining the critical packing fraction of glass transition $\eta_{\textrm{c}}$ with constant activity and constant $x_{\textrm{A}}$. In FIG. 3(a), we show the dependence of critical packing fraction $\eta_{\textrm{c}}$ on number fraction of active component $x_{\textrm{A}}$, for different values of the activity $D_{f}\tau_{\textrm{p}}^{2}$=0, 1.0, 3.0 and 10.0. Generally, larger activity leads to larger shift of critical packing fraction $\Delta\eta_{\textrm{c}}$. For small activity such as $D_{f}\tau_{\textrm{p}}^{2}$=1.0, the increases of critical packing fraction $\Delta\eta_{\textrm{c}}$ are almost proportional to fraction of active component $x_{\textrm{A}}$. However for large activity, since the packing fractions are obviously bounded, the increases of critical packing fractions are strongly nonlinear. Interestingly, for hard sphere systems [11], one had found that there is a definite upper bond of critical packing fraction of glass transition, i.e. the random closest packing, approximately \eta_{\textrm{max}}\approx0.64. The activity of particles can only push the critical packing fraction close to this value. However for WCA potential particles, despite the critical packing fraction is obviously bounded in physical region, it is quite difficult to find out the exact upper bound. Nevertheless, we have observed the suppression of this upper bound, analogous to the phenomenon in hard sphere systems.  FIG. 3 Phase diagram of active-passive mixture systems, phase points above the lines are glass state, otherwise are liquids. In (a), each line denotes an activity, with the equal size of active and passive particles. While for (b), the size ratio changes to \sigma_{\textrm{PP}}/\sigma_{\textrm{AA}}=0.8. For particles with soft core interactions, as we know, the critical density \rho_{\textrm{c}} of glass transition increases with the temperature in equilibrium systems. For nonequilibrium active systems, it is suggested that one can define an effective temperature T_{\text{eff}}. In our previous work for single-component active system described by the AOU-T model [26], we found that the critical density of glass transition increases with an effective temperature T_{\text{eff}}=k_{\rm{B}}^{ - 1}\gamma \left( {{D_t}} \right.+\left. {{D_f}\tau _{\rm{p}}^2/{\gamma ^2}} \right)). For the AP mixture system considered in the present work, it is yet not clear to define an effective temperature in a consistent way. Nevertheless, such an effective temperature, if can be well-defined, should increase with the increase of particle activity D_{f}\tau_{\textrm{p}}^{2} or the fraction of active component x_{\textrm{A}}. Therefore, the critical volume fraction \eta_{\textrm{c}} should also increase with D_{f}\tau_{\textrm{p}}^{2} or x_{\textrm{A}}, as demonstrated in FIG. 3(a). In FIG. 3(b), we also calculate the phase diagram with different size ratio of active to passive particles \sigma_{\textrm{PP}}/\sigma_{\textrm{AA}}=0.8, but with the same other parameters as FIG. 3(a). The lowest (black) line in this figure denotes a passive binary mixture, where the critical packing fraction \eta_{\textrm{c}} shows a non-monotonic dependence on x_{\textrm{A}}. Such a reentrance behaviour comes from the mixing effect [31]. When activity is present, the critical packing fraction \eta_{\textrm{c}} shifts to a higher value. {The shift value \Delta\eta_\textrm{c} increases with the fraction of active component x_{\textrm{A}}, as well as the activity D_{f}\tau_{\textrm{p}}^{2}. }When the activity is large enough such as D_{f}\tau_{\textrm{p}}^{2}=10.0, the reentrance behaviour disappears. Ⅳ. CONCLUSION In summary, we have presented a mode-coupling theory for active-passive binary mixture, and studied the glass transition behaviour of this system. The active forces propelling the active particles are described by the OU process. Using the Fox's method to handle this coloured noise, we obtain an approximate Smoluchowski equation. With this SE, we derive the mode-coupling equation for the time evolution of intermediate scattering functions, and calculate Debye-Waller factor since we focus on phase behaviour in present work. To obtain the DW factor, we only need effective diffusion coefficient of active particles \bar{D}^{\textrm{A}} and partial structure factor \mathbb{S}\left(q\right) as input information, both of them are obtained through direct simulations. According to the behaviour of DW factors, one can determine the critical point of glass transition. We find that activity pushes the critical effective packing fraction of glass transition to a higher value. For the systems consisting of active and passive particles with the same size, the critical packing fraction \eta_{\textrm{c}} increases monotonically with the increasing of number fraction of active component x_{\textrm{A}}. For AP mixture with different particle sizes, we find that \eta_{\textrm{c}} bypasses a maximum value with the increasing of x_{\textrm{A}}, demonstrating an interesting type of reentrance behavior. In addition, such reentrance can only be observed for small particle activity. In next step, we will calculate ISFs to study the relaxation behaviour near glass transition. We will also study the situation with smaller active particle and larger passive particle, to seek out the behaviour in so-called active bath systems [32]. In summary, we believe that our MCT framework for AP mixture is applicable to study glass transition or other behaviours in dense systems from the theoretical point of view, and could find many developments and applications in the future. Ⅴ. ACKNOWLEDGEMENT This work is supported by the Ministry of Science and Technology of China (No.2016YFA0400904 and No.2013CB834606), the National Natural Science Foundation of China (No.21673212, No.21521001, No.21473165, No.21403204), and by the Fundamental Research Funds for the Central Universities (No.WK2030020028 and No.2340000074). APPENDIX A: Effective Smoluchowski equation Previously, the method of achieving the approximate Smoluchowski equation from a Langevin equation with coloured noise has been presented by Fox [28, 29]. Later, Brader et al. [30] and our group [27] have used this method to calculate the approximate Smoluchowski equation of active particle systems. Here, we concentrate on handling the mixture system and more details can be found in these papers. For a general multi-variable OU noise case,  \begin{eqnarray} \frac{\textrm{d}x_{i}^{\alpha}\left(t\right)}{\textrm{d}t}=G_{i}^{\alpha}\left(\{x_{i}\}\right)+\chi_{i}^{\alpha}\left(t\right)+\xi_{i}^{\alpha}\left(t\right) \end{eqnarray} (A1) where \left\langle \xi_{i}^{\alpha}\left(t\right)\xi_{j}^{\beta}\left(t'\right)\right\rangle =2D_{t}^{\alpha}\delta_{\alpha\beta}\delta_{ij}\delta\left(t-t'\right), \chi_{i}^{\textrm{P}}=0, G_{i}^{\alpha}=\gamma_{\alpha}^{-1}F_{i}=-\gamma_{\alpha}^{-1}\nabla_{i}\displaystyle{\sum_{j}}u_{ij}, and  \begin{eqnarray} \left\langle \chi_{i}^{\textrm{A}}\left(t\right)\chi_{j}^{\textrm{A}}\left(s\right)\right\rangle \hspace{-0.15cm}&=&\hspace{-0.15cm}C_{il}\left(t-s\right) \\ \hspace{-0.15cm}&=&\hspace{-0.15cm}D_{f}\tau_{\textrm{p}}\textrm{e}^{-\left|s-t\right|/\tau_{\textrm{p}}}\delta_{ij} \end{eqnarray} (A2) the Smoluchowski equation writes  \begin{array}{l} \frac{\partial }{{\partial t}}P\left( {\mathit{\boldsymbol{y}},t} \right) - \sum\limits_\alpha {\sum\limits_{i = 1}^{{N_\alpha }} {{\partial _i}} } \left[ {G_i^\alpha \left( \mathit{\boldsymbol{y}} \right)P\left( {\mathit{\boldsymbol{y}},t} \right)} \right]\\ + \sum\limits_\alpha {\sum\limits_{i = 1}^{{N_\alpha }} {D_t^\alpha } } \partial _i^2P\left( {\mathit{\boldsymbol{y}},t} \right)\\ + \sum\limits_{i,j = 1}^{{N_A}} {{\partial _i}} \sum\limits_{l = 1}^{{N_A}} {\int_0^t {\rm{d}} } s'{C_{il}}\left( {t - s'} \right){\partial _j}\\ \int D \left[ \chi \right]P\left[ \chi \right] \times \exp \int_{s'}^t {\rm{d}} s\frac{\partial }{{\partial {x_l}}}G_j^A\\ \left( {\mathit{\boldsymbol{x}}\left( s \right)} \right){\delta _{jl}}]\delta \left( {\mathit{\boldsymbol{y - x}}\left( t \right)} \right)\} \end{array} where P\left(\boldsymbol{y}, t\right) is the probability distribution function. Using the assumption presented in Ref.[28], we get  \begin{array}{l} \frac{\partial }{{\partial t}}P\left( {\mathit{\boldsymbol{y}},t} \right) = - \sum\limits_\alpha {} \sum\limits_{i = 1}^{{N_\alpha }} {\{ {\partial _i}\left[ {\gamma _\alpha ^{ - 1}{F_i}\left( \mathit{\boldsymbol{y}} \right)P\left( {\mathit{\boldsymbol{y}},t} \right)} \right] - } \\ D_t^\alpha \partial _i^2P\left( {\mathit{\boldsymbol{y}},t} \right)\} + \sum\limits_{i = 1}^{{N_{\rm{A}}}} {{D_f}} \partial _i^2\\ \left\{ {\left[ {\frac{1}{{1 - {\tau _{\rm{p}}}\gamma _{\rm{A}}^{ - 1}{\partial _i}{{F'}_i}\left( \mathit{\boldsymbol{y}} \right)}}} \right]P\left( {\mathit{\boldsymbol{y}},t} \right)} \right\} \end{array} (A3) where \boldsymbol{F}'_{j}=\displaystyle{\sum_{i=1, j\neq i}^{N_{\textrm{A}}}}, \boldsymbol{F}_{ij}=-\nabla_{j}\displaystyle{\sum_{i=1, j\neq i}^{N_{\textrm{A}}}}u_{\textrm{AA}}(r_{ij}), notice that the summation is only for active particles. APPENDIX B: The mode-coupling equation It is convenient to achieve the memory function equation (Eq.(15)) in Laplace domain. At the beginning, we introduce the Laplace transform of ISF  \begin{eqnarray} \tilde{S}_{\alpha\beta}\left(q, z\right)=\mathcal{LT}\left[S_{\alpha\beta}\left(q, t\right)\right]=\left\langle \rho_{-\boldsymbol{q}}^{\alpha}\frac{1}{z-\Omega}\rho_{\boldsymbol{q}}^{\beta}\right\rangle \end{eqnarray} (B1) where the Laplace transform and corresponding inversed transform are defined as  \mathcal{LT}\left[\cdots\right]=\int_{0}^{\infty}\cdots \textrm{e}^{-zt}\textrm{d}t  \mathcal{LT}^{-1}\left[\cdots\right]=\frac{1}{2\pi i}\int_{\sigma-i\infty}^{\sigma+i\infty}\cdots \textrm{e}^{zt}\textrm{d}z Using the Zwanzig-Mori projection operator method, define a projection operator on the density fluctuation \rho_{\boldsymbol{q}}^{\alpha}  {\cal P} \equiv \sum\limits_{\alpha ,\beta = {\rm{A}},{\rm{P}}} {} \rho _\mathit{\boldsymbol{q}}^\alpha \rangle {\left[ {\mathbb{S}{^{ - 1}}(q)} \right]_{\alpha \beta }}\langle \rho _{ - \mathit{\boldsymbol{q}}}^\beta (B2) Clearly we have \mathcal{PP}=\mathcal{P} and \mathcal{P}\left.\rho_{\boldsymbol{q}}^{\gamma}\right\rangle =\left.\rho_{\boldsymbol{q}}^{\gamma}\right\rangle . Introducing a counterpart operator \mathcal{Q}=1-\mathcal{P}, we also have \mathcal{QQ}=\mathcal{Q}, \mathcal{PQ}=\mathcal{QP}=0, \mathcal{Q}\rho_{\boldsymbol{q}}^{\gamma}\rangle =0. Using Dyson decomposition  \begin{eqnarray} \frac{1}{z-\Omega}=\frac{1}{z-\Omega\mathcal{Q}}+\frac{1}{z-\Omega\mathcal{Q}}\Omega\mathcal{P}\frac{1}{z-\Omega} \end{eqnarray} (B3) (notice: the operator function can be understood as series \displaystyle{\frac{1}{z-\Omega}}=\sum_{i=0}^{\infty}\frac{\Omega^{i}}{z^{i+1}}), we have  \begin{eqnarray*} \mathcal{LT}\left[\partial_{t}S_{\mu\nu}\left(q, t\right)\right]&=&z\tilde{S}_{\mu\nu}\left(q, z\right)-S_{\mu\nu}\left(q, t=0\right) \\ &=&\left\langle \rho_{-\boldsymbol{q}}^{\mu}\Omega\frac{1}{z-\Omega}\rho_{\boldsymbol{q}}^{\nu}\right\rangle \\ &=&\left\langle \rho_{-\boldsymbol{q}}^{\mu}\Omega\mathcal{P}\frac{1}{z-\Omega}\rho_{\boldsymbol{q}}^{\nu}\right\rangle \\ &&+\left\langle \rho_{-\boldsymbol{q}}^{\mu}\Omega\mathcal{Q}\frac{1}{z-\Omega}\rho_{\boldsymbol{q}}^{\nu}\right\rangle \end{eqnarray*} where  \begin{eqnarray*} \left\langle \rho_{-\boldsymbol{q}}^{\mu}\Omega\mathcal{P}\frac{1}{z-\Omega}\rho_{\boldsymbol{q}}^{\nu}\right\rangle &=&\sum\limits_{\alpha\beta}\left\langle \rho_{-\boldsymbol{q}}^{\mu}\Omega\rho_{\boldsymbol{q}}^{\alpha}\right\rangle \left[\mathbb{S}^{-1}(q)\right]_{\alpha\beta}\cdot \\ &&\left\langle \rho_{-\boldsymbol{q}}^{\beta}\frac{1}{z-\Omega}\rho_{\boldsymbol{q}}^{\nu}\right\rangle \\ &=&\sum\limits_{\alpha\beta}\left\langle \rho_{-\boldsymbol{q}}^{\mu}\Omega\rho_{\boldsymbol{q}}^{\alpha}\right\rangle \left[\mathbb{S}^{-1}(q)\right]_{\alpha\beta}\tilde{S}_{\beta\nu}\left(q, z\right) \end{eqnarray*} and  \begin{align} &\left\langle \rho _{-\mathit{\boldsymbol{q}}}^{\mu }\Omega \mathcal{Q}\frac{1}{z-\Omega }\rho _{\mathit{\boldsymbol{q}}}^{\nu } \right\rangle =\left\langle \rho _{-\mathit{\boldsymbol{q}}}^{\mu }\Omega \mathcal{Q}\left[ \frac{1}{z-\Omega \mathcal{Q}}+ \right. \right. \\ &\left. \frac{1}{z-\Omega \mathcal{Q}}\Omega \mathcal{P}\frac{1}{z-\Omega }\left. {} \right]\rho _{\mathit{\boldsymbol{q}}}^{\nu } \right\rangle \\ \end{align} Since \mathcal{Q}\left.\rho_{\boldsymbol{q}}^{\gamma}\right\rangle=0 we have  \begin{eqnarray} \left\langle \rho_{-\boldsymbol{q}}^{\mu}\Omega\mathcal{Q}\frac{1}{z-\Omega\mathcal{Q}}\rho_{\boldsymbol{q}}^{\nu}\right\rangle =0 \end{eqnarray} and the rest part  \begin{eqnarray} & &\hspace{-0.15cm} \left\langle \rho_{-\boldsymbol{q}}^{\mu}\Omega\mathcal{Q}\frac{1}{z-\Omega}\rho_{\boldsymbol{q}}^{\nu}\right\rangle \\ &= &\hspace{-0.15cm} \left\langle \rho_{-\boldsymbol{q}}^{\mu}\Omega\mathcal{Q}\frac{1}{z-\Omega\mathcal{Q}}\Omega\mathcal{P}\frac{1}{z-\Omega}\rho_{\boldsymbol{q}}^{\nu}\right\rangle \\ &= &\hspace{-0.15cm} \sum\limits_{\alpha\beta}\left\langle \rho_{-\boldsymbol{q}}^{\mu}\Omega\mathcal{Q}\frac{1}{z-\Omega\mathcal{Q}}\Omega\rho_{\boldsymbol{q}}^{\alpha}\right\rangle \left[\mathbb{S}^{-1}(q)\right]_{\alpha\beta}\cdot \\ &&\left\langle \rho_{-\boldsymbol{q}}^{\beta}\frac{1}{z-\Omega}\rho_{\boldsymbol{q}}^{\nu}\right\rangle \\ &= &\hspace{-0.15cm} \sum\limits_{\alpha\beta}\left\langle \rho_{-\boldsymbol{q}}^{\mu}\Omega\mathcal{Q}\frac{1}{z-\Omega\mathcal{Q}}\Omega\rho_{\boldsymbol{q}}^{\alpha}\right\rangle\cdot \\ &&\left[\mathbb{S}^{-1}(q)\right]_{\alpha\beta}\tilde{S}_{\beta\nu}\left(q, z\right) \end{eqnarray} and  \begin{eqnarray} &&\left\langle \rho_{-\boldsymbol{q}}^{\alpha}\Omega\mathcal{Q}\frac{1}{z-\Omega\mathcal{Q}}\Omega\rho_{\boldsymbol{q}}^{\beta}\right\rangle \\ && = \left\langle \rho_{-\boldsymbol{q}}^{\alpha}\Omega\mathcal{Q}\frac{1}{z-\mathcal{Q}\Omega\mathcal{Q}}\mathcal{Q}\Omega\rho_{\boldsymbol{q}}^{\beta}\right\rangle \\ &&= \left\langle \rho_{-\boldsymbol{q}}^{\alpha}\Omega\mathcal{Q}\frac{1}{z-\mathcal{Q}\Omega\mathcal{Q}}\left(\mathcal{Q}\Omega^{\dagger}\rho_{\boldsymbol{q}}^{\beta}\right)\right\rangle \\ &&=\left\langle \left(\mathcal{Q}\Omega^{\dagger}\rho_{-\boldsymbol{q}}^{\alpha}\right)\frac{1}{z-\mathcal{Q}\Omega\mathcal{Q}}\left(\mathcal{Q}\Omega^{\dagger}\rho_{\boldsymbol{q}}^{\beta}\right)\right\rangle \\ &&=\left\langle R_{-\boldsymbol{q}}^{\alpha}\frac{1}{z-\mathcal{Q}\Omega\mathcal{Q}}R_{\boldsymbol{q}}^{\beta}\right\rangle \end{eqnarray} where R_{\boldsymbol{q}}^{\beta}=\mathcal{Q}\Omega^{\dagger}\rho_{\boldsymbol{q}}^{\beta}=(1-\mathcal{P})\Omega^{\dagger}\rho_{\boldsymbol{q}}^{\beta} is a kind of random force. Here, the operator \Omega^{\dagger} means the adjoint of Smoluchowski operator \Omega, defined as \displaystyle{\int} f\Omega^{\dagger}g\textrm{d}\Gamma=\displaystyle{\int} g\Omega f\textrm{d}\Gamma, and therefore  \begin{eqnarray} \Omega^{\dagger}=\sum\limits_{\alpha=\textrm{A, P}}\sum\limits_{j=1}^{N_{\alpha}}\left(\nabla_{j}^{\alpha}+\beta\boldsymbol{F}_{j}^{\alpha}\right)D_{j}^{\alpha}\left(\boldsymbol{r}^{N}\right)\cdot\nabla_{j}^{\alpha} \end{eqnarray} (B4) Then we define the frequency term  \begin{eqnarray} W_{\alpha\beta}\left(q\right)=-\left\langle \rho_{-\boldsymbol{q}}^{\alpha}\Omega\rho_{\boldsymbol{q}}^{\beta}\right\rangle \end{eqnarray} (B5) and  \begin{eqnarray*} \left\langle \rho_{-\boldsymbol{q}}^{\alpha}\Omega\rho_{\boldsymbol{q}}^{\beta}\right\rangle \hspace{-0.15cm}&= &\hspace{-0.15cm} \sum\limits_{\gamma}\sum\limits_{j=1}^{N_{\gamma}}\left\langle \rho_{-\boldsymbol{q}}^{\alpha}\nabla_{j}^{\gamma}\cdot D_{j}^{\gamma}\left(\nabla_{j}^{\gamma}-\beta\boldsymbol{F}_{j}^{\gamma}\right)\rho_{\boldsymbol{q}}^{\beta}\right\rangle \\ \hspace{-0.15cm}&= &\hspace{-0.15cm} \sum\limits_{\gamma}\sum\limits_{j=1}^{N_{\gamma}}\left\langle -\left(\nabla_{j}^{\gamma}\rho_{-\boldsymbol{q}}^{\alpha}\right)\cdot D_{j}^{\gamma}\left(\nabla_{j}^{\gamma}-\beta\boldsymbol{F}_{j}^{\gamma}\right)\rho_{\boldsymbol{q}}^{\beta}\right\rangle \\ \hspace{-0.15cm} &= &\hspace{-0.15cm} -\frac{1}{\sqrt{N_{\alpha}N_{\beta}}}\sum\limits_{\gamma}\hspace{-0.1cm} \sum\limits_{j=1}^{N_{\gamma}}q^{2}\hspace{-0.1cm}\left\langle \delta_{\gamma\alpha}\textrm{e}^{i\boldsymbol{q}\cdot\boldsymbol{r}_{j}^{\gamma}}D_{j}^{\gamma}\delta_{\gamma\beta}\textrm{e}^{-i\boldsymbol{q}\cdot\boldsymbol{r}_{j}^{\gamma}}\right\rangle \\ &=&\hspace{-0.15cm} -\frac{\delta_{\alpha\beta}q^{2}}{N_{\alpha}}\sum\limits_{j=1}^{N_{\alpha}}\left\langle D_{j}^{\alpha}\right\rangle \\ \hspace{-0.15cm} &\equiv&\hspace{-0.15cm} -\delta_{\alpha\beta}q^{2}\bar{D}^{\alpha} \end{eqnarray*} therefore  \begin{eqnarray} \mathbb{W}(q)=q^{2}\left[\begin{array}{cc} \bar{D}^{\textrm{A}}&0\\ 0&D_{t}^{\textrm{P}} \end{array}\right] \end{eqnarray} (B6) Considering the inverse matrix of structure factor  \begin{align} & {{\left[ \begin{matrix} {{S}_{\text{AA}}}\left( q \right)\quad {{S}_{\text{AP}}}\left( q \right) \\ {{S}_{\text{PA}}}\left( q \right)\quad {{S}_{\text{PP}}}\left( q \right) \\ \end{matrix} \right]}^{-1}} \\ & =\frac{1}{{{S}_{\text{AA}}}\left( q \right){{S}_{\text{PP}}}\left( q \right)-{{S}_{\text{AP}}}{{\left( q \right)}^{2}}}\cdot \\ & \left[ \begin{matrix} {{S}_{\text{PP}}}\left( q \right)\quad -{{S}_{\text{AP}}}\left( q \right) \\ -{{S}_{\text{AP}}}\left( q \right)\quad {{S}_{\text{AA}}}\left( q \right) \\ \end{matrix} \right] \\ \end{align} we calculate the random force term  \begin{eqnarray*} \mathcal{P}\Omega^{\dagger}\rho_{\boldsymbol{q}}^{\beta} \hspace{-0.15cm}&= &\hspace{-0.15cm} \sum\limits_{\mu\nu}\left.\rho_{\boldsymbol{q}}^{\mu}\right\rangle \left[\mathbb{S}^{-1}(q)\right]_{\mu\nu}\left\langle \rho_{-\boldsymbol{q}}^{\nu}\left(\Omega^{\dagger}\rho_{\boldsymbol{q}}^{\beta}\right)\right\rangle \\ \hspace{-0.15cm}&= &\hspace{-0.15cm} -\sum\limits_{\mu\nu}\left.\rho_{\boldsymbol{q}}^{\mu}\right\rangle \left[\mathbb{S}^{-1}(q)\right]_{\mu\nu}W_{\nu\beta}(q)\\ R_{\boldsymbol{q}}^{\beta}&=&\left(1-\mathcal{P}\right)\Omega^{\dagger}\rho_{\boldsymbol{q}}^{\beta}\\ \hspace{-0.15cm} &= &\hspace{-0.15cm} \Omega^{\dagger}\rho_{\boldsymbol{q}}^{\beta}+\sum\limits_{\mu\nu}\left.\rho_{\boldsymbol{q}}^{\mu}\right\rangle \left[\mathbb{S}^{-1}(q)\right]_{\mu\nu}W_{\nu\beta}(q)\\ \hspace{-0.15cm}&= &\hspace{-0.15cm} \Omega^{\dagger}\rho_{\boldsymbol{q}}^{\beta}+\sum\limits_{\mu}\left.\rho_{\boldsymbol{q}}^{\mu}\right\rangle \left[\mathbb{S}^{-1}(q)\right]_{\mu\beta}W_{\beta\beta}(q) \end{eqnarray*} Define the memory function  \begin{eqnarray} \tilde{M}_{\alpha\beta}\left(q, z\right) \hspace{-0.15cm}&= &\hspace{-0.15cm} -\sum\limits_{\gamma}\left\langle \rho_{-\boldsymbol{q}}^{\alpha}\Omega\mathcal{Q}\frac{1}{z-\Omega\mathcal{Q}}\Omega\rho_{\boldsymbol{q}}^{\gamma}\right\rangle\cdot \\ &&\hspace{-0.15cm}\left[\left\langle \rho_{-\boldsymbol{q}}^{\mu}\Omega\rho_{\boldsymbol{q}}^{\nu}\right\rangle \right]_{\gamma\beta}^{-1} \\ &= &\hspace{-0.15cm} \sum\limits_{\gamma}\left\langle \rho_{-\boldsymbol{q}}^{\alpha}\Omega\mathcal{Q}\frac{1}{z-\Omega\mathcal{Q}}\Omega\rho_{\boldsymbol{q}}^{\gamma}\right\rangle \left[\mathbb{W}\left(q\right)\right]_{\gamma\beta}^{-1} \\ \end{eqnarray} (B7) and therefore, we have an evolution function of scattering function in Laplace space  \begin{align} & z{{{\tilde{S}}}_{\mu \nu }}\left( q,z \right)-{{S}_{\mu \nu }}\left( q,0 \right) \\ & =-\sum\limits_{\alpha \beta \gamma }{\left[ {{\delta }_{\mu \alpha }}-{{{\tilde{M}}}_{\mu \alpha }}(q,z \right)}\left. {} \right]\cdot \\ & {{W}_{\alpha \beta }}\left( q \right){{\left[ {{\mathbb{S}}^{-1}}\left( q \right) \right]}_{\beta \gamma }}{{{\tilde{S}}}_{\gamma \nu }}\left( q,z \right) \\ \end{align} or in matrix form  \begin{eqnarray} \tilde{\mathbb{S}}\left(q, z\right)=\left\{ z+\left(\boldsymbol{1}-\tilde{\mathbb{M}}\left(q, z\right)\right)\mathbb{W}\left(q\right)\mathbb{S}^{-1}\left(q\right)\right\} ^{-1}\mathbb{S}\left(q, 0\right) \end{eqnarray} Now introduce the so-called irreducible memory function \mathbb{M}^{{\rm irr}}\left(q, t\right), satisfying  \tilde{\mathbb{M}}\left(q, z\right)=\tilde{\mathbb{M}}^{{\rm irr}}\left(q, z\right)\left[1+\tilde{\mathbb{M}}^{{\rm irr}}\left(q, z\right)\right]^{-1} where  \tilde{M}_{\alpha\beta}^{{\rm irr}}\left(q, z\right)=-\sum\limits_{\gamma}\left\langle R_{-\boldsymbol{q}}^{\alpha}\frac{1}{z-\mathcal{Q}\Omega^{{\rm irr}}\mathcal{Q}}R_{\boldsymbol{q}}^{\gamma}\right\rangle \left[\left\langle \rho_{-\boldsymbol{q}}^{\mu}\Omega\rho_{\boldsymbol{q}}^{\nu}\right\rangle \right]_{\gamma\beta}^{-1}  \Omega^{{\rm irr}}=\sum\limits_{\alpha}\sum\limits_{j=1}^{N}\nabla_{j}^{\alpha}\cdot D_{j}^{\alpha}Q_{j}^{\alpha}\left(\nabla_{j}^{\alpha}-\beta\boldsymbol{F}_{j}^{\alpha}\right) Here, the operator \Omega^{{\rm irr}} is the irreducible Smoluchowski operator, where Q_{j}^{\alpha}=1-P_{j}^{\alpha} and P_{j}^{\alpha}\equiv\left.\textrm{e}^{-i\boldsymbol{q}\cdot\boldsymbol{r}_{j}^{\alpha}}\right\rangle \left\langle \textrm{e}^{i\boldsymbol{q}\cdot\boldsymbol{r}_{j}^{\alpha}}\right.. This step is a standard operation in MCT. Then the evolution of ISF writes  \begin{align} & \widetilde{\mathbb{S}}\left( q,z \right)=\left\{ {} \right.z+{{\left[ \mathbf{1}+{{\widetilde{\mathbb{M}}}^{\text{irr}}}\left( q,z \right) \right]}^{-1}} \\ & \mathbb{W}\left( q \right){{\mathbb{S}}^{-1}}\left( q \right){{\}}^{-1}}\mathbb{S}\left( q,0 \right) \\ \end{align} and in real time space  \begin{eqnarray} &&\partial_{t}\mathbb{S}\left(q, t\right)+\mathbb{W}\left(q\right)\mathbb{S}^{-1}\left(q\right)\mathbb{S}\left(q, t\right)+ \\ &&\int_{0}^{t}\textrm{d}u\mathbb{M}^{{\rm irr}}\left(q, u\right)\partial_{u}\mathbb{S}\left(q, u\right)=0 \end{eqnarray} Next, we need to introduce the mode-coupling approximation to calculate the irreducible memory function. Firstly define a second-order projection operator  \begin{array}{l} {{\cal P}_2} \equiv \frac{1}{2}\sum\limits_\mathit{\boldsymbol{k}} {} \sum\limits_{\alpha \beta \gamma \delta } {} \rho _\mathit{\boldsymbol{k}}^\alpha \rho _{\mathit{\boldsymbol{q - k}}}^\beta \left. {} \right\rangle \left[ {\mathbb{S}\left( k \right)} \right]_{\alpha \gamma }^{ - 1}\left[ {\mathbb{S}\left( {\left| {\mathit{\boldsymbol{q - k}}} \right|} \right)} \right]_{\beta \delta }^{ - 1}\\ \left\langle {} \right.\rho _{ - \mathit{\boldsymbol{k}}}^\gamma \rho _{\mathit{\boldsymbol{k - q}}}^\delta \end{array} (B8) and insert the operator into the irreducible memory function  \begin{eqnarray*} &&\left\langle R_{-\boldsymbol{q}}^{\mu}\textrm{e}^{\mathcal{Q}\Omega^{\textrm{irr}}\mathcal{Q}t}R_{\boldsymbol{q}}^{\nu}\right\rangle \\ & \approx&\left\langle R_{-\boldsymbol{q}}^{\mu}\mathcal{P}_{2}\textrm{e}^{\mathcal{Q}\Omega^{\textrm{irr}}\mathcal{Q}t}\mathcal{P}_{2}R_{\boldsymbol{q}}^{\nu}\right\rangle \\ &=&\frac{1}{4}\sum\limits_{\boldsymbol{k}, \boldsymbol{k}'}\sum\limits_{\alpha\beta\gamma\delta}\sum\limits_{\alpha'\beta'\gamma'\delta'}\left\langle R_{-\boldsymbol{q}}^{\mu}\rho_{\boldsymbol{k}}^{\alpha}\rho_{\boldsymbol{q}-\boldsymbol{k}}^{\beta}\right\rangle \left[\mathbb{S}\left(k\right)\right]_{\alpha\gamma}^{-1} \times\\ &&\left[\mathbb{S}\left(\left|\boldsymbol{q}-\boldsymbol{k}\right|\right)\right]_{\beta\delta}^{-1}\left\langle \rho_{-\boldsymbol{k}}^{\gamma}\rho_{\boldsymbol{k}-\boldsymbol{q}}^{\delta}\textrm{e}^{\mathcal{Q}\Omega^{irr}\mathcal{Q}t}\rho_{\boldsymbol{k}'}^{\alpha'}\rho_{\boldsymbol{q}-\boldsymbol{k}'}^{\beta'}\right\rangle \times\\ &&\left[\mathbb{S}\left(k'\right)\right]_{\alpha'\gamma'}^{-1}\left[\mathbb{S}\left(\left|\boldsymbol{q}-\boldsymbol{k}'\right|\right)\right]_{\beta'\delta'}^{-1}\left\langle \rho_{-\boldsymbol{k}'}^{\gamma'}\rho_{\boldsymbol{k}'-\boldsymbol{q}}^{\delta'}R_{\boldsymbol{q}}^{\nu}\right\rangle \end{eqnarray*} The mode-coupling approximation means that approximate the four-point correlation function to the products of two-point correlation functions.  \begin{eqnarray} &&\left\langle \rho_{-\boldsymbol{k}}^{\gamma}\rho_{\boldsymbol{k}-\boldsymbol{q}}^{\delta}\textrm{e}^{\mathcal{Q}\Omega^{irr}\mathcal{Q}t}\rho_{\boldsymbol{k}'}^{\alpha'}\rho_{\boldsymbol{q}-\boldsymbol{k}'}^{\beta'}\right\rangle \\ \hspace{-0.15cm}& \approx &\hspace{-0.15cm} \delta_{\boldsymbol{kk}'}\left\langle \rho_{-k}^{\gamma}\textrm{e}^{\Omega t}\rho_{\boldsymbol{k}'}^{\alpha'}\right\rangle \left\langle \rho_{\boldsymbol{k}-\boldsymbol{q}}^{\delta}\textrm{e}^{\Omega t}\rho_{\boldsymbol{q}-\boldsymbol{k}'}^{\beta'}\right\rangle \\ & &\hspace{-0.15cm} +\delta_{\boldsymbol{k}, \boldsymbol{q}-\boldsymbol{k}'}\left\langle \rho_{-k}^{\gamma}\textrm{e}^{\Omega t}\rho_{\boldsymbol{q}-\boldsymbol{k}'}^{\beta'}\right\rangle \left\langle \rho_{\boldsymbol{k}-\boldsymbol{q}}^{\delta}\textrm{e}^{\Omega t}\rho_{\boldsymbol{k}'}^{\alpha'}\right\rangle \\ \hspace{-0.15cm}&=&\hspace{-0.15cm} \delta_{\boldsymbol{kk}'}S_{\gamma\alpha'}\left(k, t\right)S_{\delta\beta'}\left(\left|\boldsymbol{q}-\boldsymbol{k}\right|, t\right) \\ &&\hspace{-0.15cm} +\delta_{\boldsymbol{k}, \boldsymbol{q}-\boldsymbol{k}'}S_{\gamma\beta'}\left(k, t\right)S_{\delta\alpha'}\left(\left|\boldsymbol{q}-\boldsymbol{k}\right|, t\right) \end{eqnarray} (B9) Calculate several terms in irreducible memory function  \begin{eqnarray*} &&\hspace{-0.15cm}\left\langle R_{-\boldsymbol{q}}^{\mu}\rho_{\boldsymbol{k}}^{\alpha}\rho_{\boldsymbol{q}-\boldsymbol{k}}^{\beta}\right\rangle \\ &= &\hspace{-0.15cm} \left\langle \rho_{-\boldsymbol{q}}^{\mu}\Omega\mathcal{Q}\rho_{\boldsymbol{k}}^{\alpha}\rho_{\boldsymbol{q}-\boldsymbol{k}}^{\beta}\right\rangle \\ &= &\hspace{-0.15cm} \left\langle \rho_{-\boldsymbol{q}}^{\mu}\Omega\rho_{\boldsymbol{k}}^{\alpha}\rho_{\boldsymbol{q}-\boldsymbol{k}}^{\beta}\right\rangle - \\ &&\sum\limits_{\gamma\delta}\left\langle \rho_{-\boldsymbol{q}}^{\mu}\Omega\rho_{\boldsymbol{q}}^{\gamma}\right\rangle \left[\mathbb{S}^{-1}\left(q\right)\right]_{\gamma\delta}\left\langle \rho_{-\boldsymbol{q}}^{\delta}\rho_{\boldsymbol{k}}^{\alpha}\rho_{\boldsymbol{q}-\boldsymbol{k}}^{\beta}\right\rangle \end{eqnarray*}  \begin{eqnarray*} &&\left\langle \rho_{-\boldsymbol{q}}^{\mu}\Omega\rho_{\boldsymbol{k}}^{\alpha}\rho_{\boldsymbol{q}-\boldsymbol{k}}^{\beta}\right\rangle \\ & = &\sum\limits_{\gamma}\sum\limits_{j=1}^{N_{\gamma}}\left\langle \rho_{-\boldsymbol{q}}^{\mu}\nabla_{j}^{\gamma}\cdot D_{j}^{\gamma}\left(\nabla_{j}^{\gamma}-\beta\boldsymbol{F}_{j}^{\gamma}\right)\rho_{\boldsymbol{k}}^{\alpha}\rho_{\boldsymbol{q}-\boldsymbol{k}}^{\beta}\right\rangle \\ &=&-\delta_{\mu\alpha}\frac{\boldsymbol{q}\cdot\boldsymbol{k}}{\sqrt{N_{\mu}}}D_{t}^{\alpha}S_{2}^{\alpha\beta}\left(\left|\boldsymbol{q}-\boldsymbol{k}\right|\right)- \\ &&\delta_{\mu\beta}\frac{\boldsymbol{q}\cdot\left(\boldsymbol{q}-\boldsymbol{k}\right)}{\sqrt{N_{\mu}}}D_{t}^{\beta}S_{2}^{\beta\alpha}\left(k\right) \end{eqnarray*} where  \begin{eqnarray} S_{2}^{\alpha\beta}\left(q\right)=\frac{1}{D_{t}^{\alpha}}\sum\limits_{j=1}^{N_{\alpha}}\left\langle D_{j}^{\alpha}\rho_{-\boldsymbol{q}}^{\alpha}\rho_{\boldsymbol{q}}^{\beta}\right\rangle \end{eqnarray} (B10) is a pseudo "structure factor", which is discussed in Ref.[27]. In practice, we find that the correlations between D_{j}^{\alpha} and \rho_{-\boldsymbol{q}}^{\alpha}\rho_{\boldsymbol{q}}^{\beta} are weak and can be neglected in mode-coupling equation. Therefore, we use the approximation S_{2}^{\alpha\beta}\left(q\right)\approx$\displaystyle{\frac{\bar{D}^{\alpha}}{D_{t}^{\alpha}}}S_{\alpha\beta}\left(q\right)$ in the following derivation. Then we need the convolution approximation

 $\begin{eqnarray} &&\hspace{-0.3cm}\left\langle \rho_{-\boldsymbol{q}}^{\delta}\rho_{\boldsymbol{k}}^{\alpha}\rho_{\boldsymbol{q}-\boldsymbol{k}}^{\beta}\right\rangle \\ \hspace{-1cm}&\approx&\hspace{-0.3cm}\sum\limits_{\lambda=\textrm{A, P}}\frac{1}{\sqrt{N_{\lambda}}}S_{\delta\lambda}\left(q\right)S_{\alpha\lambda}\left(k\right)S_{\beta\lambda}\left(\left|\boldsymbol{q}-\boldsymbol{k}\right|\right) \end{eqnarray}$ (B11)

and therefore

 $\begin{eqnarray} &&\sum\limits_{\alpha\beta}\left\langle R_{-\boldsymbol{q}}^{\mu}\rho_{\boldsymbol{k}}^{\alpha}\rho_{\boldsymbol{q}-\boldsymbol{k}}^{\beta}\right\rangle \left[\mathbb{S}\left(k\right)\right]_{\alpha\gamma}^{-1}\left[\mathbb{S}\left(\left|\boldsymbol{q}-\boldsymbol{k}\right|\right)\right]_{\beta\delta}^{-1} \\ & \approx&\frac{\bar{D}^{\mu}}{\sqrt{N_{\mu}}}\left[\boldsymbol{q}\cdot\boldsymbol{k}\delta_{\mu\delta}C_{\mu\gamma}\left(k\right)\right. \\ &&+\left.\boldsymbol{q}\cdot\left(\boldsymbol{q}-\boldsymbol{k}\right)\delta_{\mu\gamma}C_{\mu\delta}\left(\left|\boldsymbol{q}-\boldsymbol{k}\right|\right)\right] \\ &=&V_{\mu}^{\gamma\delta}\left(\boldsymbol{q};\boldsymbol{k}\right) \end{eqnarray}$ (B12)

where $\mathbb{S}^{-1}(k)$=1-$\mathbb{C}(k)$, $C_{\alpha\beta}\left(k\right)$=$\sqrt{\rho_{\alpha}\rho_{\beta}}c_{\alpha\beta}\left(k\right)$ and proportion to direct correlation function. Finally, we can calculate the memory term

 $\begin{eqnarray} &\left\langle R_{-\boldsymbol{q}}^{\mu}\mathcal{P}_{2}\textrm{e}^{\mathcal{Q}\Omega^{\textrm{irr}}\mathcal{Q}t}\mathcal{P}_{2}R_{\boldsymbol{q}}^{\nu}\right\rangle \approx\displaystyle{\frac{1}{2}}\sum\limits_{\boldsymbol{k}}\sum\limits_{\alpha'\beta'\gamma\delta}V_{\mu}^{\gamma\delta}\left(\boldsymbol{q};\boldsymbol{k}\right)\times \\ &S_{\gamma\alpha'}\left(k, t\right)S_{\delta\beta'}\left(\left|\boldsymbol{q}-\boldsymbol{k}\right|, t\right)V_{\nu}^{\alpha'\beta'}\left(\boldsymbol{q};\boldsymbol{k}\right) \end{eqnarray}$ (B13)

and irreducible memory function writes

 $\begin{eqnarray} M_{\alpha\beta}^{{\rm irr}}\left(q, t\right) \hspace{-0.15cm}& \approx &\hspace{-0.15cm} \frac{1}{2}\sum\limits_{\boldsymbol{k}}\sum\limits_{\nu}\frac{\delta_{\nu\beta}}{q^{2}\bar{D}^{\nu}}\sum\limits_{\alpha'\beta'\gamma\delta}V_{\alpha}^{\gamma\delta}\left(\boldsymbol{q};\boldsymbol{k}\right)\times \\ &&\hspace{-1.5cm}S_{\gamma\alpha'}\left(k, t\right)S_{\delta\beta'}\left(\left|\boldsymbol{q}-\boldsymbol{k}\right|, t\right)V_{\nu}^{\alpha'\beta'}\left(\boldsymbol{q};\boldsymbol{k}\right) \end{eqnarray}$ (B14)
Reference
 [1] T. Vicsek, and A. Zafeiris, Phys. Rep. 517 , 71 (2012). DOI:10.1016/j.physrep.2012.03.004 [2] M. C. Marchetti, . J. F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. D. Rao, and R. A. Simha, Rev. Mod. Phys. 85 , 1147 (2013). [3] C. Bechinger, R. Di Leonardo, H. Löwen, C. Reichhardt, G. Volpe, and G. Volpe, Rev. Mod. Phys. 88 , 045006 (2016). [4] V. Schaller, C. Weber, C. Semmrich, E. Frey, and A. R. Bausch, Nature 467 , 73 (2010). [5] Y. Sumino, K. H. Nagai, Y. Shitaka, D. Tanaka, K. Yoshikawa, H. Chat, and K. Oiwa, Nature 483 , 48 (2012). [6] J. Schwarz-Linek, C. Valeriani, A. Cacciuto, M. E. Cates, D. Marenduzzo, A. N. Morozov, and W. C. K. Poon, Proc. Nat. Acad. Sci. USA 109 , 4052 (2012). DOI:10.1073/pnas.1116334109 [7] Y. Fily, and M. C. Marchetti, Phys. Rev. Lett. 108 , 235702 (2012). DOI:10.1103/PhysRevLett.108.235702 [8] I. Buttinoni, J. Bialk, F. Kmmel, H. Löwen, C. Bechinger, and T. Speck, Phys. Rev. Lett. 110 , 238301 (2013). DOI:10.1103/PhysRevLett.110.238301 [9] R. Wittkowski, A. Tiribocchi, J. Stenhammar, R. J. Allen, D. Marenduzzo, and M. E. Cates, Nat. Commun. 5 , 4351 (2014). [10] S. K. Das, S. A. Egorov, B. Trefz, P. Virnau, and K. Binder, Phys. Rev. Lett. 112 , 198301 (2014). DOI:10.1103/PhysRevLett.112.198301 [11] R. Ni, M. A. Cohen Stuart, and M. Dijkstra, Nat. Commun. 4 , 2704 (2013). DOI:10.1038/ncomms3704 [12] G. Szamel, E. Flenner, and L. Berthier, Phys. Rev. E 91 , 062304 (2015). [13] E. Flenner, G. Szamel, and L. Berthier, The Nonequilibrium Glassy Dynamics of Self-Propelled Particles, arXiv preprint arXiv:1606. , 00641 (2016). [14] Z. Ma, Q. L Lei, and R. Ni, Mobility Induced Demixing in the Mixtures of Passive and Eccentric Active Brownian Particles, arXiv preprint arXiv:1706. , 00202 (2017). [15] M. E. Cates, and E. Tjhung, J. Fluid Mech. 836, P , 1 (2018). [16] R. Wittmann, J. M. Brader, A. Sharma, and U. M. B. Marconi, Phys. Rev. E 97 , 012601 (2018). DOI:10.1103/PhysRevE.97.012601 [17] S. N. Weber, C. A. Weber, and E. Frey, Phys. Rev. Lett. 116 , 058301 (2016). DOI:10.1103/PhysRevLett.116.058301 [18] J. Stenhammar, R. Wittkowski, D. Marenduzzo, and M. E. Cates, Phys. Rev. Lett. 114 , 018301 (2015). DOI:10.1103/PhysRevLett.114.018301 [19] L. Berthier, and J. Kurchan, Nat. Phys. 9 , 310 (2013). DOI:10.1038/nphys2592 [20] S. Mandal, S. Lang, M. Gross, M. Oettel, D. Raabe, T. Franosch, and F. Varnik, Nat. Commun. 5 , 4435 (2014). DOI:10.1038/ncomms5435 [21] L. Berthier, Phys. Rev. Lett. 112 , 220602 (2014). DOI:10.1103/PhysRevLett.112.220602 [22] A. Sharma, R. Wittmann, and J. M. Brader, Phys. Rev. E 95 , 012115 (2017). [23] V. Krakoviack, Phys. Rev. Lett. 94 , 065703 (2005). DOI:10.1103/PhysRevLett.94.065703 [24] G. M. Whitesides, and B. Grzybowski, Science 295 , 2418 (2002). [25] A. Liluashvili, J. Ónody, and T. Voigtmann, Phys. Rev. E 96 , 062608 (2017). DOI:10.1103/PhysRevE.96.062608 [26] M. K. Feng, and Z. H. Hou, Soft Matter 13 , 4464 (2017). DOI:10.1039/C7SM00852J [27] R. F. Fox, Phys. Rev. A 34 , 4525 (1986). DOI:10.1103/PhysRevA.34.4525 [28] R. F. Fox, Phys. Rev. A 33 , 467 (1986). DOI:10.1103/PhysRevA.33.467 [29] T. F. F. Farage, P. Krinninger, and J. M. Brader, Phys. Rev. E 91 , 042310 (2015). [30] J. P. Hansen and I. R. McDonald, Theory of Simple Liquids, Amsterdam: Elsevier, (1990). [31] T. Voigtmann, Europhys. Lett. 96 , 36006 (2011). DOI:10.1209/0295-5075/96/36006 [32] X. L. Wu, and A. Libchaber, Phys. Rev. Lett. 84 , 3017 (2000).