Mode-coupling theory for mixtures of athermal self-propelled particles

Dense or glassy active matter, as a result of its remarkable resemblance to passive glass-forming materials, is enjoying increasing scientific interest. To better grasp the subtle effect of active motion on the process of vitrification, a number of active mode-coupling theories (MCTs) have recently been developed. These have proven capable of qualitatively predicting important parts of the active glassy phenomenology. However, most efforts so far have only considered single-component materials, and their derivations are arguably more complex than the standard MCT case, which might hinder broader usage. Here we present a detailed derivation of a distinct active MCT for mixtures of athermal self-propelled particles that is more transparent than previously introduced versions. The key insight is that we can follow a similar strategy for our overdamped active system as is typically used for passive underdamped MCT. Interestingly, when only considering one particle species, our theory gives the exact same result as the one obtained in previous work which employed a highly different mode-coupling strategy. Moreover, we assess the quality of the theory and its novel extension to multi-component materials by using it to predict the dynamics of a Kob-Andersen mixture of athermal active Brownian quasi-hard spheres. We demonstrate that our theory is able to capture all qualitative features, most notably the location of the optimum of the dynamics when the persistence length and cage length coincide, for each combination of particle types.


Introduction
Motivated by its frequent occurrence in the biological realm, active matter has enjoyed much scientific interest during the past two decades [1][2][3] and is now firmly positioned at the forefront of soft matter and biophysical research.Initial active matter studies, e.g., to explain the swimming behavior of microorganisms, have been predominantly geared towards self-propelled particles in the now reasonably well-understood dilute limit [1].By contrast, the high-density regime has only recently seen a significant upsurge of interest [4,5].Dense active matter has already been shown to exhibit a myriad of exciting non-equilibrium phenomena, such as velocity ordering [6,7], motility induced phase separation (MIPS) [8][9][10][11][12][13], collective actuation [14], and activity-induced crystallization [15,16].So-called active glassy materials constitute an arguably even more extreme example of dense active matter.These are disordered assemblies of self-propelled particles that, irrespective of their intrinsic driving, have undergone kinetic arrest.Experimental examples can be found in the context of cell layers [17][18][19][20], individual cells [21][22][23], granular matter [24], and colloids [25,26], and they have also been witnessed in many theoretical and simulation studies [12,.As a result of the increasingly strong role of particleparticle interactions at high densities, active glassy materials show a remarkable resemblance to their passive * l.m.c.janssen@tue.nlcounterparts [4,50,51].However, this does not imply that activity is completely overshadowed by interactions as it can indeed influence glassy dynamics in surprising ways [29, 32, 33, 35-37, 47, 48].Grasping the subtle effect of active motion on the process of vitrification has therefore surfaced as an encouraging line of research and one that unites the comprehensive fields of glassy physics and active matter [4].
For passive materials, one of the few first-principlesbased theories that is capable of making reasonable predictions about their glassy dynamics is mode-coupling theory (MCT) [52][53][54][55][56] (and its hierarchical extension GMCT [57][58][59][60][61][62][63][64][65][66][67]).Requiring only the static structure factor as input, it can predict the full relaxation dynamics via the intermediate scattering function with qualitative, and in few instances even semi-quantitative, agreement.In an attempt to better comprehend active glassy dynamics, and inspired by the previous successes of MCT, several works have recently set out to extend the theory by also including self-propelled particle motion [27, 28, 30, 31, 39-41, 68, 69].This has resulted in a number of so-called active MCTs that qualitatively capture part of the active glassy phenomenology.However, most efforts have so far only considered monodisperse systems, and their mathematical derivations are rather involved, which might discourage a broader usage of the theories.
Here, to add to our theoretical understanding of active glassy dynamics, we present an active MCT for mixtures of athermal self-propelled particles that is distinct from, and conceptually more transparent than, previously introduced versions.We consider an overdamped active system for which we can follow a similar strategy as typically used for the well-established passive underdamped MCT [55].This helps make the derivation more insightful and, since we do not have to specify the time evolution of our active force, allows it to be more general.Interestingly, we show that for a monodisperse system our theory yields the exact same equation for the intermediate scattering function as has been previously derived in another (seemingly more involved) mode-coupling theory [28,30].Finally, we test our theory's predictive capabilities for multicomponent systems by comparing it to simulation results obtained for a Kob-Andersen mixture of athermal active Brownian quasi-hard spheres.We find that for all particle types our theory gives qualitatively consistent predictions and thus captures a non-trivial maximum of the dynamics when the persistence length equals the socalled cage length.

Theory
As our model system we consider an athermal N -particle active fluid of volume V which consists of m different species α with component number densities ρ α = N α /V .In the overdamped limit, the motion in time t of each individual particle is described via the following equation [31,33,70] Here, r α i denotes the position of the ith particle of type α, the over-dot the derivative with respect to time, ζ α the (species-dependent) friction constant, and F α i and f α i the interaction and self-propulsion force acting on particle i respectively.Normally, one proceeds by also introducing a time-evolution equation for the self-propulsion force.This will turn out not to be necessary for our derivation and we will therefore refrain from doing so to improve the generality of our derivation.Examples of popular model systems include so-called active Ornstein Uhlenbeck (AOUPs) [29-31, 34, 36, 37] and active Brownian particles (ABPs) [27,28].
Let us now introduce the velocity of each particle, i.e., p α i = ṙα i , and use it to complement the particle positions as our degrees of freedom (essentially replacing the active forces).Note that for thermal particles this is not feasible since p α i would become discontinuous due to the thermal noise term.The joint N -particle probability distribution of positions and velocities P N (Γ; t) then evolves in time via where Γ = ({r α i }, {p α i }) represents the degrees of freedom and Ω is the Smoluchowski operator which can be inferred from the equations of motion for the positions and velocities (see [71] for an example with AOUPs).Now we assume that our system can reach a steady-state characterized by a probability distribution P ss N (Γ) that obeys [28,30] ΩP ss N (Γ) = 0. (3) Using the steady-state distribution, one can define the time-correlation function C(t) of any dynamical vector A(t) whose elements A i (t) = [A(t)] i are functions of the degrees of freedom as Here averages . . .are taken with respect to P ss N (Γ), A = A(0), the asterisk denotes complex conjugation, and the adjoint (or backward) Smoluchowski operator Ω † is working on everything to its right except the probability distribution.Moreover, by taking the derivative with respect to time and setting it to zero, we also obtain the useful property The strategy of mode-coupling theory is then to choose for the elements of the dynamical vector slow or quasiconserved quantities of the system.We will focus on the conventional ones usually considered in the derivation of underdamped passive MCT, i.e., the density modes and their respective time-derivatives or current modes [53,55,62].Note that we thus assume that current modes retain their slow character, which might not necessarily be true in an active matter setting.In particular, we have where √ N α and j α k = −i ρα k are the density and current modes respectively, N α the number of particles of type α, and k a wave vector which probes the (inverse) length scale of interest.
Having specified our vector we can then employ the Mori-Zwanzig formalism [72,73] to develop an equation of motion for its time-correlation function.We define a projection operator onto the subspace spanned by A as , where we have introduced G = A * A .The superscript −1 denotes the inverse matrix of the respective quantity, i.e.
ij ensures the idempotence of P. Following standard procedure one may find (6) where Given the division of our dynamical vector into density and current modes, it is now convenient to segment C(t) into four separate (m × m) matrices as Similar to previous work in active MCT we assume that current densities vanish after integrating out the velocities (or active forces) [28,30].In other words, we assume that Consequently, we have (9) Here, we have introduced the partial static structure factor S αβ (k) and the static velocity-structure correlation function (both of which are symmetric in their species label), Interestingly, the latter function is identical to the one introduced in a different active MCT approach [30] and has subsequently been studied in several numerical works [34,36,37].It serves to quantify correlations between the velocities of individual particles and represents a distinct nonequilibrium feature.That is, it is a constant in the passive limit and develops significant oscillations upon increasing the persistence of particles.
To find a suitable expression for the frequency matrix we will also assume that the time-correlated version of ω(k) decays exponentially over a timescale equal to the persistence time τ p of the active force (assumed to be the same for all species).This corresponds to the fact that the velocities decorrelate over a similar characteristic time as the active forces and implies that we have and thus Note that this approximation likely becomes progressively worse for larger persistence as velocities are probably decorrelating on shorter time scales due to collisions with other particles.Finally, we realize that Combining the above derived results and focusing on the lower-left term of C(t), whose elements are proportional to the time derivative of the intermediate scattering function F αβ (k, t) = ρ α * k ρ β k (t) , allows us to write down the following dynamical equation for F αβ (k, t): The most involved term in this equation is the memory kernel, which is formally written as with |R r as defined in eq. ( 13) and To proceed and make analytical progress we need to approximate this term.Therefore, we apply standard techniques from conventional MCT and replace its projected with full dynamics, while also projecting on density doublets [62,67,74].In other words, we have where we have used the projection operator, (17) Our aim is now to calculate explicit expressions for the so-called vertices, i.e., R α l ρ µ1 q1 ρ µ2 q2 and ρ ν3 * q3 ρ ν4 * q4 R β r .For this we first note that due to our assumption of vanishing currents [see eq. ( 8)], terms of the form j α * k ρ β q1 ρ γ q2 are equal to zero.Furthermore, we will employ the convolution approximation [74] with δ ij a Kronecker delta, and make use of the fact that due to time-translational invariance we may rewrite jα * The last expression we require is for terms of the form ρα * k ρβ q1 ρ γ q2 .In particular, we propose a multi-component extension of a previously introduced convolution approximation for correlation functions involving active particle velocities (see Ref. [30]).This yields ρα * where we have introduced We emphasize that for a single component system (m = 1) eq. ( 19) reduces to the one presented in Ref. [30].
Using these results one can then show that the memory kernel simplifies to with the vertices given by which in turn are described by a modified direct correlation function, To make our equation self-consistent (and thus solvable) we factorize the four-point density correlation function, i.e., ρ µ * q ρ ν * k−q e Ω † t ρ µ q ρ ν k−q ≈ F µµ (q, t)F νν (|k − q| , t) δ q,q + F µν (q, t)F νµ (|k − q| , t) δ k−q,q .
(24) so that we have and, taking the thermodynamic limit, we finally obtain Interestingly, when only considering one particle type (m = 1), the equations of motion for the intermediate scattering function and specifically the derived memory kernel are identical to the ones for AOUPs and ABPs (neglecting thermal noise) obtained in previous work, which employed a highly different (and seemingly more involved) mode-coupling strategy [28,30].Moreover, in the passive limit where τ p ω αβ (k) = D t δ αβ (with D t the translational diffusion coefficient) our equation reduces to the well-known MCT equation derived for mixtures of Brownian particles [67,74].

MCT Numerics
To self-consistently solve the derived active MCT equations, one needs to complement the theory with a numerical scheme.For this we invoke the rotational symmetry of our system to rewrite the three-dimensional integral over q in eq. ( 26) in terms of the bipolar coordinates q = |q| and p = |k − q|.The single wavenumber integrals are in turn approximated by a Riemann sum on an equidistant grid kσ = [0.2,0.6, . . ., 39.8] where σ is the unit of length [75].The integration over time in eq. ( 14) is evaluated by means of Fuchs' algorithm [76].In particular, we evaluate the first N t /2 = 16 points in time via a Taylor expansion with a step size ∆t = 10 −6 , numerically integrate eq. ( 26) for the next N t /2 points in time, duplicate the timestep, and repeat this integration process until the long-time limit is reached.

Simulation Details
Given the extension of our theory to multi-component materials, we test its predictions for a model binary active glassformer as a proof of principle.In particular, we choose to study the dynamics of a Kob-Andersen (KA) mixture [77] consisting of N A = 800 and N B = 200 quasi-hard active spheres.The evolution of each particle i is described by eq. ( 1), with the interaction force F i derived from a steep repulsive power-law potential 36 [78,79].The corresponding interaction parameters are given by AA = 1, AB = 1.5, BB = 0.5, σ AA = 1, σ AB = 0.8, σ BB = 0.88, which, combined with setting ζ αβ = 1, allow for glassy behavior by suppressing crystallization [77,80].For the self-propulsion force f α i we use the active Brownian particle (ABP) model [27,28].That is, the absolute value of the force f remains constant over time, f α i = f e α i , while the orientation e α i undergoes random diffusion on the unit sphere, i.e., Here, χ α i represents a Gaussian noise process with zero mean and variance χ α i (t)χ β j (t ) noise = 2D r Iδ ij δ αβ δ(t− t ), with I the unit matrix and D r the rotational diffusion coefficient (taken to be the same for each particle type).
In the dilute limit when particle-particle interactions are assumed to be absent, each of our particles performs a persistent random walk (PRW).As a result, their mean square displacement (MSD) is given by [31] δr 2 (t) = 6T eff τ p (e −t/τp − 1) + t , (28) where we have introduced the persistence time τ p = (2D r ) −1 and an effective temperature T eff = f 2 τ p /3.A closer look at eq. ( 28) reveals that for t τ p , the motion is ballistic δr 2 (t) ≈ 3T eff t 2 /τ p , while diffusive motion, δr 2 (t) ≈ 6T eff t, is obtained for long times (t τ p ).We can therefore conclude that in the limit τ p → 0 (with T eff ∼ constant), our active system reduces to a Brownian one at a temperature T = T eff .It is thus convenient to introduce T eff as our control parameter, which we complement with the persistence length l p = f τ p as a measure for how far we are from the Brownian limit [47].
Each individual simulation consists of solving the overdamped equation [eq.( 1)] in time with a forward Euler scheme using LAMMPS [81].We set the cutoff radius of the repulsive potential at r c = 2.5σ αβ and fix the size of the cubic periodic simulation box to L = 9.41, such that the number density is ρ = 1.2.We run the system sufficiently long (typically between 200 and 1000 time units) to prevent aging, and afterwards track the particles over time for at least twice the initialization time.All results are presented in reduced units where σ AA , AA , AA /k B , and ζσ 2 AA / AA represent the units of length, energy, temperature, and time respectively [82].To correct for the influence of diffusive center-of-mass motion, all particle positions are retrieved relative to the momentary center of mass [82].

Results & Discussion
In previous work involving the same model glassformer, it has been shown that for a fixed effective temperature T eff , the dynamics exhibits a nonmonotonic dependence on the persistence length l p [47].As an initial assessment of the quality of our theory, it is interesting to see whether it is capable of predicting this nontrivial behavior.Before we test our theoretical prediction, however, we first want to verify the nonmonotonic dynamics.For this we have extracted the long-time diffusion coefficient D = lim t→∞ δr 2 (t) /6t as a function of the persistence length l p at a fixed value of T eff = 4.0.The resulting values are shown in fig. 1 and clearly illustrate nonmonotonic dynamics for both particle species (A and B).Moreover, we find that, consistent with literature [47], the maximum of the dynamics corresponds to the point where the persistence length is approximately equal to the cage length, i.e., l p ∼ 0.1.Having benchmarked our simulation results we now proceed to the active MCT predictions.Based on the retrieved particle trajectories we have calculated the static structure factors S αβ (k) and velocity-structure correlation functions ω αβ (k), which in turn have been rewritten in terms on an equidistant grid via cubic spline.Using these as input for our active MCT, we have calculated the predicted intermediate scattering function F αβ (k, t).The AA contribution (normalized by the static structure factor) is plotted as a function of time for a subset of persistence lengths in fig. 2. Interestingly, this scattering function decays to zero fastest at an intermediate persistence length, and hence the theory seems able to capture the nonmonotonic dependence of the relaxation dynamics.
To quantify the observed behavior in more detail, we have also retrieved the MCT-predicted alpha-relaxation time via F αβ (k, τ α )/S αβ (k) = e −1 where the wavenumber k corresponds to the main peak of S αβ (k).The results are presented in fig. 3 and show clear nonmonotonic and almost identical behavior for all combinations of particle types (either AA, BB, or AB).Note that the fastest relaxation dynamics corresponds to the smallest value of τ α (in contrast to the largest value of D).Moreover, the optimum of the dynamics again coincides with the point where the persistence length is approximately equal to the cage length, i.e., l p ∼ 0.1.Our theory is thus capable of accurately depicting the qualitative behavior of the relaxation dynamics upon increasing the persistence of the constituent particles.
We finalize our discussion by highlighting two noticeable quantitative features.First we find that in all cases our theory predicts faster relaxation dynamics for the smaller type B particles, which is consistent with our simulation results (see fig. 1) and is intuitively to be expected [82].More strikingly, we also observe that our active MCT predicts a dramatic speedup of the dynamics (orders of magnitude decrease of the relaxation time) compared to that obtained from standard MCT for an analogous passive Brownian system at the same effective temperature (l p = 0.0, T = T eff ).At first glance this might seem surprising, and while it is probably influenced by the assumptions made in the theory, we argue that this behavior is at least partly to be anticipated.
To illustrate this we have calculated, based on simulation data, the self-intermediate scattering function, i.e., F s α (k, t) = e −ik•r α j (0) e ik•r α j (t) , for the majority type A particles and the corresponding alpha relaxation time τ s α defined via F s α (k, τ s α ) = e −1 .The results for different effective temperatures at a fixed persistence length (l p = 0.1, i.e., on the order of the cage length corresponding to the optimum of the dynamics) are shown in fig. 4. For comparison, we have also added the values obtained for a passive Brownian system (l p = 0.0, T = T eff ).We can see that as we lower the (effective) temperature, the relaxation time in both cases increases, but the relative difference between the passive and optimum active dynamics is simultaneously being amplified and can reach differences of several orders of magnitude.Thus, as we approach dynamical arrest by lowering T eff , an optimal active system becomes relatively much more dynamic (i.e., more liquid-like) than its Brownian counterpart.Given that our passive MCT predicts an extremely large relaxation time, it is thus to some degree expected that the optimal active dynamics yields a relaxation time that is orders of magnitude smaller.The alpha relaxation time τ s α as a function of the inverse effective temperature for a passive Brownian system (lp = 0.0, T = T eff ) and an active system with a persistence length on the order of the cage length (lp = 0.1).The results are extracted from the self-intermediate scattering of the majority A species F s A (k, t) which has been directly calculated from the simulation data.

Conclusion
To conclude, we have presented a fully time-dependent microscopic mode-coupling theory for mixtures of athermal self-propelled particles.The crucial insight for our derivation is that, since we neglect thermal diffusion, the total velocity of each particle is well-behaved; Therefore we can introduce these velocities (instead of the active forces) as our degrees of freedom complementing the particle positions.This then allows us to follow a similar strategy for our overdamped active system as is typically used for passive underdamped MCT [55].Moreover, it also enables us to leave the time-evolution of the active force unspecified, thereby adding to the generality of the theory.The main result consists of an equation of motion for the (partial) intermediate scattering function, which can be self-consistently solved using the static structure factor and a distinctly non-equilibrium static velocitystructure correlation function as input.Remarkably, for a monodisperse system this equation turns out to be exactly equal to one that has been derived in a previous (and possibly more convoluted) active mode-coupling theory [28,30].
As an initial assessment of the quality of the theory and especially to test its novel extension to multi-component materials, we have used it to predict the dynamics of a Kob-Andersen mixture of athermal active Brownian particles.Such particles exhibit nonmonotonic behavior for increasing particle persistence and thus form a stringent test for the theory.Our theory is indeed able to capture all qualitative features, most notably the location of the optimum of the dynamics when the persistence length and cage length coincide, for each combination of particle types.On a quantitative level active MCT predicts (upon approaching dynamical arrest) a dramatic enhancement of the dynamics (multiple orders of magnitude) compared to that obtained from standard MCT for an analogous passive Brownian system.Though surprising, we show that this effect can in fact be anticipated from simulations.Given the success of our theoretical framework to give qualitatively consistent results, it would be interesting to see whether the analogy between overdamped athermal active systems and underdamped passive ones can be further exploited to better understand the intriguing phenomenology of active glassy matter.

FIG. 1 .
FIG.1.The long-time diffusion coefficient D of both the Aand B-type particles as a function of the persistence length lp at a constant effective temperature T eff = 4.0.The results are directly obtained from the simulation data and the value retrieved from passive Brownian dynamics simulations (lp = 0) is added as a reference (dashed lines).

FIG. 2 .FIG. 3 .
FIG. 2. The intermediate scattering function of the majorityA species FAA(k, t) [normalized by the static structure factor SAA(k)] as a function of time t for different persistence lengths lp at a constant effective temperature T eff = 4.0.The results are obtained using active MCT for a wavenumber k corresponding to the main peak of the static structure factor.
FIG. 4. The alpha relaxation time τ sα as a function of the inverse effective temperature for a passive Brownian system (lp = 0.0, T = T eff ) and an active system with a persistence length on the order of the cage length (lp = 0.1).The results are extracted from the self-intermediate scattering of the majority A species F s A (k, t) which has been directly calculated from the simulation data.