Nonequilibrium dynamics of mixtures of active and passive colloidal particles

We develop a mesoscopic field theory for the collective nonequilibrium dynamics of multicomponent mixtures of interacting active (i.e., motile) and passive (i.e., nonmotile) colloidal particles with isometric shape in two spatial dimensions. By a stability analysis of the field theory, we obtain equations for the spinodal that describes the onset of a motility-induced instability leading to cluster formation in such mixtures. The prediction for the spinodal is found to be in good agreement with particle-resolved computer simulations. Furthermore, we show that in active-passive mixtures the spinodal instability can be of two different types. One type is associated with a stationary bifurcation and occurs also in one-component active systems, whereas the other type is associated with a Hopf bifurcation and can occur only in active-passive mixtures. Remarkably, the Hopf bifurcation leads to moving clusters. This explains recent results from simulations of active-passive particle mixtures, where moving clusters and interfaces that are not seen in the corresponding one-component systems have been observed.


Introduction
Mixtures of active (i.e., motile) and passive (i.e., nonmotile) colloidal particles exhibit an interesting set of collective behaviors that is quite different from the dynamics of the corresponding one-component systems. Many studies of active-passive mixtures have focused on the dilute regime, where long-range hydrodynamic flows generated by active particles affect the diffusion of passive tracer particles [1][2][3]. At higher densities, where short-range interactions such as excluded-volume interactions become important, other fascinating phenomena arise. Among them are effective depletion-like attractions between passive objects in an activeparticle suspension [4][5][6], crystallization and melting of hard-sphere glasses by doping with active particles [7,8], and mesoscale turbulence mediated by passive particles [9]. Furthermore, the intriguing phenomenon of motility-induced phase separation (MIPS) [10], whereby purely repulsive active particles spontaneously segregate into dense and dilute phases, has recently been numerically shown to occur also for mixtures of active and passive particles [11,12]. For purely active systems, where MIPS has already been extensively studied [13][14][15][16][17][18][19][20][21][22][23][24][25][26], this transition has been theoretically rationalized as a spinodal-like instability occurring for sufficiently large densities and propulsion speeds. A striking feature observed in active-passive mixtures undergoing MIPS is the emergence of persistently moving interfaces and clusters due to a spontaneous breaking of spatial symmetry [11,12], which is qualitatively different from the stationary, albeit fluctuating, clusters observed in purely active suspensions.
To gain a deeper theoretical understanding of MIPS in mixtures of active and passive particles, as well as mixtures of active particles with different properties, we here derive a mesoscopic field theory that describes the collective dynamics of multicomponent mixtures of such particles. In order to restrict the study to the basic features of their nonequilibrium dynamics, we consider unbounded systems in two spatial dimensions and particles with isometric shapes (i.e., disks or spheres in a plane). Furthermore, we neglect hydrodynamic interactions between the particles. The different species therefore differ only in their radii and motilities. By using appropriate approximations for the full (multidimensional) pair-distribution functions of the particles, we derive spinodal criteria based on coefficients that can be readily evaluated from computer simulations in the one-phase parameter regime. The predicted spinodals are in good accordance with the instability region observed in simulations of active-passive mixtures. We furthermore show how two different types of instabilities arise, depending on the properties of the mixtures: a stationary bifurcation, corresponding to MIPS in onecomponent active systems, and a Hopf bifurcation, which can only occur in mixtures and is associated with a steady state with moving clusters. This distinction is likely to be responsible for the intriguing collective dynamics observed in computer simulations of active-passive mixtures undergoing MIPS [11,12].
The article is organized as follows: in section 2 we derive dynamic equations for multicomponent mixtures of particles with different motilities as well as useful approximations of these equations. Based on these, in section 3 we address the stability of mixtures of active and passive particles and derive a condition for the onset of a dynamical instability in such mixtures. The good agreement between our analytic results and particle-resolved computer simulations of mixtures of active and passive particles is demonstrated in section 4. We finally conclude in section 5.

Dynamic equations for mixtures
We consider a multicomponent mixture of colloidal particles, each with different but constant magnitudes of self-propulsion. In addition to the deterministic self-propulsion, the positions and orientations of the particles (so-called 'active Brownian particles' [27]) undergo translational and rotational diffusion, respectively. In the following, n denotes the number of different species of active or passive colloidal particles and m N is the total number of particles of species is the orientation of the ith particle of species μ at time t. The interactions between a particle of species μ and a particle of species ν are described by the pairinteraction potentials - can be different 5 and   · is the Euclidean norm.

Exact dynamic equations
The microscopic active Brownian motion of an n-component mixture of these colloidal particles with different constant motilities is described by the Langevin equations 6 Here, an overdot denotes differentiation with respect to time and b = ( ) k T 1 B is the inverse thermal energy with the Boltzmann constant k B and the absolute temperature T of the implicit solvent. Furthermore, where ⊗ denotes the dyadic product and 1 is the identity matrix. The collective dynamics of the n-component mixture can equivalently be described by the Smoluchowski equation (see [30,31] for details) 5 This is the case especially for colloidal particles that catalyze chemical reactions on their surfaces and break the action-reaction symmetry stated by Newton's third law. As an interesting consequence of such asymmetric interactions, anisotropic clusters of passive particles can become self-propelled [29]. However, in this article we do not further study such systems. Our simulations described in section 4 use the same interaction potential for all particles. 6 In the limit of a one-component system, these Langevin equations reduce to the simpler ones that have been used in [16,19].
i i . Many-particle densities like the two-particle densities can be obtained analogously [33]. Applying an integration as in equation (5) . The dynamics of the one-particle densities r f m  ( ) r t , , depends on the two-particle densities r , which can be decomposed into the one-particle densities Note that this standard definition of where the term d mn m N is usually negligibly small. Since the one-particle densities in equation (7) depend on orientation and time, the commonly used orientation-and time-independent pair-distribution functions ¢ mn   ( ) g r r , cannot be simply obtained by an angular integration and a time average of f f ¢ ¢ mn   ( ) g r r t , , , , . Instead, they need to be calculated from an equation analogous to equation (7), where all functions are replaced by their orientation-integrated and time-averaged counterparts.
Inserting equation (7) into equation (6) results in the exact dynamic equations for the one-particle densities The pair-distribution functions f f ¢ ¢ mn   ( ) g r r t , , , , , however, are unknown and have to be approximated appropriately.

Approximate dynamic equations
In order to approximate the pair-distribution functions, we first assume global translational invariance of the system, leading to figure 1) this can be rewritten as f f f¢ mn ( ) g r t , , , , R . Assuming also invariance of the system with respect to global rigid rotations, the number of variables of the pair-distribution functions can be reduced further, yielding , , R . In addition, we assume that the pair-distribution functions are not explicitly timedependent: R . These common approximations become exact when the system is in a homogeneous steady state. In accordance with [16], we now further approximate , R with the symmetry property q q -= mn mn , . Notice that the pair-distribution functions f f mn ( ) g r, R still depend on the driving forces m F A of the active particles and on the constant average we assume short-range interactions between the particles. This allows us to simplify the final term of equation (9) through a gradient expansion [34][35][36].
With these approximations, equation (9) simplifies to the dynamic equations for n-component mixtures of active and passive particles with the effective speeds of the active particles T res that describe the slow-down of active particles due to collisions [15,16]. Here, the resistive forces   with the coefficients The four terms on the right-hand side of equation (10) describe entropic translational and rotational diffusion, drift caused by the particles' motility, and diffusive transport resulting from interactions of the particles, respectively. When the particles are passive, the third term vanishes and the function f ¡ m  ( ) r t , , reduces to the excess chemical potential corresponding to the particles of species μ. The dynamic equations (10) are applicable even for active systems far from thermodynamic equilibrium and constitute the first main result of this article.
So far, we approximated the pair-distribution functions by assuming that they are invariant with respect to translations and rotations of the coordinate system and that they have no explicit time-dependence, as well as by neglecting the dependence of , R on f f ¢ -. Furthermore, we assumed short-range particle interactions that allow fast convergence of the gradient expansions (12) and (13). Most of these approximations are common when dealing with passive particles. The only exception is the negligence of the f f ¢ --dependence of the pair-distribution functions, which depend only on r when the motilities of all particles vanish. With additional assumptions on the pair-distribution functions, one can reduce the number of  (12) and (13). These possible further approximations of the pairdistribution functions are described in appendix A, but not used in this article.
We proceed by expanding the orientation-dependent one-particle densities r f m  ( ) r t , , with respect to the orientation f, neglecting terms of third or higher order in gradients, and performing a quasi-stationary approximation (see appendix B for details). This results in the dynamic equations Here and in the following, ¶ i denotes the ith element of the del symbol , and Einstein's summation convention is used for Latin but not for Greek indices. Furthermore, the density-dependent diffusion tensor r m ({ }) D has the elements The division of mn Especially for one-component systems, this constitutes an alternative to determining the pair-distribution functions in simulations and calculating mn ( ) a 0 and mn ( ) a 1 using equations (14), (18), and (19).

Dynamical instability
For certain combinations of the free swim speeds m v 0 and the average particle number densities r m the mixture becomes unstable, particle clusters form, and the system phase separates [11]. In the following, we use the dynamic equations (15) to derive the spinodal for the onset of this activity-induced dynamical instability by means of a linear stability analysis.

Multicomponent mixtures
To derive the spinodal, we consider a system with slightly perturbed concentration fields r  (23) can be written as In the next subsection, we apply the stability criterion (25) to binary mixtures.

Binary mixtures
We now focus on binary mixtures with two species and the spinodal is the outer surface that encloses the two nonoverlapping 8  As can be seen from the eigenvalues (27) and (28) and accompanying eigenvectors  (29), where a complex conjugate pair of eigenvalues ofD passes through the imaginary axis, is associated with a Hopf bifurcation [37][38][39]. This type of bifurcation is possible for a binary mixture, but not for a one-component system of active or passive colloidal particles. In contrast, under the condition (30), l 1 and l 2 are real and one of these eigenvalues is zero. The surface described by equation (30), where one of the two real eigenvalues ofD passes through zero, is associated with a stationary bifurcation [37][38][39]. This type of bifurcation is possible also for one-component systems. Hence, there is a fundamental difference between the already widely studied MIPS of one-component systems of active particles, which we related to a stationary bifurcation, and the nonequilibrium dynamics of active-passive mixtures, which exhibits also a Hopf bifurcation. 8 Combining equations (29) and (30)  To study the features of these bifurcations, we now insert the plane-wave ansatz 1 2 , imaginary unit i, wave vectors  Î   k k , 1 2 2 , and angular frequencies  w w Î , 1 2 into equation (22) From this equation it is obvious that-in accordance with the stability criterion (25)-the amplitudes of small spatially harmonic perturbations are exponentially growing with time when at least one of the eigenvalues ofD has a negative real part. Furthermore, we see that the perturbations are traveling plane waves for l ¹ ( ) Im 0 1 , i.e., for the Hopf bifurcation associated with equation (29), whereas the perturbations are static for l = ( ) Im 0 1 , as holds for the stationary bifurcation associated with equation (30). This is an important difference between the two types of bifurcations and helps to distinguish the two underlying instabilities in simulations and experiments. The speed of a traveling spatially harmonic perturbation with wave number k is l ( )k Im 2 . Note that this expression holds only for small k, since the continuum description by the field theory (15) breaks down for where R A and R B are the radii of particles of species A and B, respectively. In the parameter regime where traveling perturbations occur, an additional condition for k, which can restrict the applicability of the field theory (15) to even smaller k, results from the quasi-stationary approximation (see section 2.2). Since the applicability of the quasi-stationary approximation requires that the time scale D 1 R , on which the orientational order-parameter fields described in appendix B relax, is much smaller than the period p w 2 of the traveling waves, where ω is their angular frequency, the condition p w  with the effective swim speeds .

Large v 0
We now consider large values  ¥ v 0 of the self-propulsion speed of a free particle v 0 , but not too large particle densities r A and r P so that  ¥ v . Assuming that the coefficients ( ) ( ) ( ) ( ) ( ) a a a a a , ,  which is a simple condition for the low-density branch of the spinodal, as found previously in [11].

One-component systems
Another limiting case of equations (47) and (48) is that of a one-component system of active particles where r = 0 P . Also in this case, equation (47) has no solution so that the Hopf bifurcation with its moving perturbations cannot occur. The stationary bifurcation, in contrast, is still possible. For this bifurcation we obtain from equation (48) the spinodal condition with the effective speed of the active particles , equation (51) becomes equivalent to the simpler spinodal condition that has been proposed in [16]. This approximation is, however, in general not applicable. The term proportional to ( ) a 1 AA in the diffusion coefficientD AA takes the density-dependence of the particles' collective diffusion into account that originates from their interactions and is present also if the particles are passive ( = v 0 0 ). Further below we show that neglecting this term has a strong influence on the predicted spinodal.
If all particles interact purely repulsively, the onset of the dynamical instability requires a sufficiently large value of the low-density self-propulsion speed v 0 . The threshold value of v 0 depends on r A : for small values of r A , the threshold value of v 0 is infinite and the instability cannot occur at all; for larger r A , the threshold value is finite, and for a certain value of r A , which we denote as r A,min , the threshold value attains its minimal possible value v 0,min . Considering ( ) a 0 AA and ( ) a 1 AA as nonzero finite constants, which is-as we will see in section 4.2only an approximation, v 0,min and the corresponding r A,min are given by r r r = + -+ -  and has been shown to lead to too small values for v 0,min compared to particle-resolved simulations [18,20]. Note that, within the approximations mentioned above, r (¯) v , A,min 0,min is the critical point of the one-component system of active particles.
Next, we consider the spinodal condition (51) in the limit  ¥ v and assume that ( ) a 0 AA and ( ) a 1 AA are finite in this limit. This results in the well-known condition for the low-density spinodal branch of a suspension of active Brownian particles r = ( )

Simulations
To confirm our analytical results from section 3, we carried out Brownian dynamics computer simulations of binary mixtures of isometric active and passive colloidal particles in two spatial dimensions. These particleresolved simulations were based on solving the Langevin equations (1) and (2) numerically and made use of the LAMMPS molecular dynamics package [40]. are the average area fractions of the active and passive particles, respectively. The Péclet number describes the ratio between motility and thermal diffusion of a particle. To ensure that the effective particle radius does not depend on the propulsion speed of a free active particle v 0 , we varied Pe by changing D T for all particles, while we kept v 0 constant (see [20,41]

Pair-distribution functions
Since the pair-distribution functions 0, 500 and F Î [ ] 0, 0.8 . However, we did not include parameter combinations where the isotropic state of the simulated system was not stable. The resulting pair-distribution functions for = Pe 50and F = 0.6 are shown in figure 2.
For fixed f f -R , the shapes of these functions of r are qualitatively similar to the typical shape of a radial distribution function for repulsive disks or spheres [43]. However, the maxima and minima of  figure 2(c)). This is a consequence of the different type of the reference particle, which is active for ( ) g r, 0 AA and ( ) g r, 0 AP but passive for ( ) g r, 0 PA and ( ) g r, 0 PP . When the reference particle is active, it pushes other particles ahead of itself so that in front of it (i.e., at f f = R ) on average over time the concentration of particles is larger and the particle separations are smaller than around a passive reference particle.

Coefficients
Using equation ( . While for small Pe we observed a strong dependence of these coefficients on Pe, we found that their dependence on Pe is negligible for > Pe 50. Under the condition > Pe 50, the calculated coefficients are approximately given by the expressions Interestingly, the coefficient ( ) a 0 AA is independent of the total particle area fraction Φ and ( ) a 0 AP has only a moderate dependence on Φ, whereas the other coefficients depend strongly on Φ.
With equations (56) and (57) the effective swim speed of an active particle (46) becomes

Dynamical state diagram
Based on our simulations, we now consider the dynamical state diagram of a binary mixture of isometric active and passive colloidal particles with the same size and the same purely repulsive interactions in two spatial dimensions. For small motilities or low concentrations of active particles, the steady state of such a mixture is homogeneous. In contrast, when the Péclet number Pe and the average area fraction of the active particles F A are sufficiently large, the homogeneous state becomes unstable and particle clusters form [11].  In the stable region (blue) the system has a homogeneous steady state, whereas in the unstable region (red) particle clusters can be observed. The predicted spinodals for the onset of a dynamical instability in the system correspond to equation (47)  region where the total particle area fraction Φ is near the close-packing area fraction p F = » ( ) 2 3 0.907 cp or even larger so that effective motion of the active particles is not possible.
In the case of an active-passive mixture with c = 0.5 A , clustering of the particles requires > Pe 50and F > 0.4 (see figure 3(b)). Evidently, the predicted spinodal for the onset of a dynamical instability in the mixture, which is in general given by equations (47) and (48) with the parameters (56)-(61), is in accordance with the border between the stable and unstable regions that we observed in the simulations. Here, only equation (47) describes the spinodal curve, whereas equation (48) has no solution in the parameter range of the state diagram. The appearance of the predicted spinodal clearly inside the actual one can partly be explained by the fact that-due to the high noise level in the simulations-in large parts of the simulated instability region clustering occurs through nucleation and growth instead of spinodal decomposition.
The state diagram of one-component systems of active particles (c = 1 A ) is qualitatively similar to that of active-passive mixtures with c = 0.5 A (see figure 3(c)). Also in the case of solely active particles, the occurrence of a motility-induced instability requires > Pe 50, but in this case an instability can already be observed for Φ down to F » 0.4 [20]. For c = 1 A , equation (47) has no solution and the spinodal curve is entirely described by equation (48) or, equivalently, equation (51). The predicted spinodal is in rather good agreement with that directly obtained from the simulations. Especially the minimal values of Pe required for the instability to occur, that are associated with the critical points of the spinodals, are very close together. When 6.9and the agreement with the simulation results is significantly worse. This shows that the density-dependence of the particles' collective diffusion has a strong influence on the predicted spinodal and should not be neglected. Setting in the case of a mixture with c = 0.5 A leads to an even stronger deterioration of the predictions. Then the predicted spinodal not only includes regions of too small Pe, but also corresponds to the wrong type of instability.
Remarkably, through equations (47) and (48), which describe the spinodals for the active-passive mixtures with c = 0.5 A and for the purely active systems with c = 1 A , respectively, also these spinodals are associated with different types of instabilities, including moving perturbations in the first case and static perturbations in the second case. This is in line with the collective dynamics of mixtures of active and passive particles being more violent, with traveling interfaces and strongly dynamic clusters (see, e.g., the supplemental movies associated with [11]), than that of purely active systems, where clusters only move through particle diffusion [11,12]. Given the different nature of the instabilities for c = 0.5 A and c = 1 A , and the fact that the mixtures exhibit significantly larger fluctuations, also the ease with which nucleation occurs in the metastable region is likely to be different for one-and two-component systems. In the case c = 0.5 A , one has to take into account that in the parameter regime where moving perturbations occur the validity of the quasi-stationary approximation, which is involved in the derivation of equations (47) and (48), is ensured only for perturbations with small wave numbers  k k max , where s » k 0.1 max near the predicted spinodal shown in figure 3(b). Both a facilitated nucleation and a reduced accuracy of the quasi-stationary approximation for c = 0.5 A could, at least partly, explain, why the agreement of the predicted spinodals with the simulation results is better for c = 1 A than for c = 0.5 A . For c = 0.5 A the agreement of the analytical prediction and the simulation results for the spinodal could therefore be improved by avoiding the quasi-stationary approximation through carrying out a stability analysis of equations (B8)-(B10) instead of equation (15). A general improvement of the field theory given by equations (B8)-(B10) and thus of equation (15) is possible by truncating the gradient expansions in equations (B8)-(B10) at higher than second order. This would lead to spinodal conditions that are-for all values of c A -more accurate than those presented in the current work.

Conclusions
We studied multicomponent mixtures of interacting colloidal particles with different but constant motilities and isometric shapes that are suspended in an atomic or molecular solvent. For this purpose, we derived a mesoscopic field theory that describes the collective nonequilibrium dynamics of such particles in two spatial dimensions. By carrying out a stability analysis for this field theory, we obtained equations that describe the onset of a motility-induced dynamical instability in such mixtures. This instability occurs for sufficiently strong selfpropulsion and a sufficiently large concentration of the active particles and leads to cluster formation.
After a general treatment of multicomponent mixtures we focused on the important special case of a binary mixture of interacting active and passive colloidal particles. We found that the instability leading to cluster formation in such an active-passive mixture can be fundamentally different from that occurring in a onecomponent system of active particles. While in the latter type of systems we found only an instability that is associated with a stationary bifurcation and a steady state of static clusters, in an active-passive mixture an instability that is associated with a Hopf bifurcation and moving clusters can also occur. This finding is in excellent agreement with the different dynamics of clusters in both systems that have already been reported before and explains why these systems can behave in very different ways. To complement our analytical calculations, we carried out Brownian dynamics computer simulations of mixtures of active and passive particles. The results of these simulations are well in line with our analytical predictions for the spinodal curves.
While our mesoscopic field theory accurately predicts the qualitatively different dynamics of clusters in active-passive mixtures and one-component systems of active particles, respectively, it does not give direct insights into the microscopic mechanism that is responsible for the different cluster dynamics. It is likely that this mechanism is analogous to the one that causes moving interfaces in active-passive mixtures reported in [12]. Different from the situation in one-component active systems, where clusters have a rather homogeneous composition of active particles with their local polarization pointing inwards the cluster [17], in active-passive mixtures the clusters are destabilized by an inhomogeneous distribution of passive particles [11]. We assume that this inhomogeneous composition of the clusters locally reduces the stabilizing inwards-pointing polarization, which causes a force and particle flux imbalance at the interfaces that leads to dynamic evaporation and deposition of particles at the interfaces and sets the clusters into motion. For reasons discussed in [12], such a flux imbalance can persist for considerable time.
An interesting feature of our theory is that it is applicable for both attractive and repulsive particle interactions and even if the interactions are different for particles of different species. Our theory therefore constitutes a general framework that can be applied to a particular system by choosing the number of particle species, interaction potentials, motilities, etc appropriately. In the future it would be useful to extend this theory towards mixtures of particles with anisometric shapes and thus to make its scope of application even wider.  . They can still depend on the speeds m v 0 and the average particle number densities r m , but are not functions of r or f. Furthermore, due to their close relationship to the pair-distribution functions f mn ( ) g r, it is possible to estimate at least lower and upper bounds for the values of + . It should also be possible to make good estimates for their dependence on m v 0 and r m .
indices. 12 Notice that the concentration fields r m  ( ) r t , are conserved quantities, whereas the local polarizations  m   ( ) r t , and the nematic tensors  m  ( ) r t , are not conserved. 13 Setting n=1 and = =